Faulhaber's triangle
You are encouraged to solve this task according to the task description, using any language you may know.
Named after Johann Faulhaber, the rows of Faulhaber's triangle are the coefficients of polynomials that represent sums of integer powers, which are extracted from Faulhaber's formula:
where is the nth-Bernoulli number.
The first 5 rows of Faulhaber's triangle, are:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
Using the third row of the triangle, we have:
- Task
- show the first 10 rows of Faulhaber's triangle.
- using the 18th row of Faulhaber's triangle, compute the sum: (extra credit).
- See also
- Bernoulli numbers
- Evaluate binomial coefficients
- Faulhaber's formula (Wikipedia)
- Faulhaber's triangle (PDF)
C
<lang c>#include <stdbool.h>
- include <stdio.h>
- include <stdlib.h>
- include <string.h>
int binomial(int n, int k) {
int num, denom, i;
if (n < 0 || k < 0 || n < k) return -1; if (n == 0 || k == 0) return 1;
num = 1;
for (i = k + 1; i <= n; ++i) {
num = num * i;
}
denom = 1;
for (i = 2; i <= n - k; ++i) {
denom *= i;
}
return num / denom;
}
int gcd(int a, int b) {
int temp;
while (b != 0) {
temp = a % b;
a = b;
b = temp;
}
return a;
}
typedef struct tFrac {
int num, denom;
} Frac;
Frac makeFrac(int n, int d) {
Frac result; int g;
if (d == 0) {
result.num = 0;
result.denom = 0;
return result;
}
if (n == 0) {
d = 1;
} else if (d < 0) {
n = -n;
d = -d;
}
g = abs(gcd(n, d));
if (g > 1) {
n = n / g;
d = d / g;
}
result.num = n; result.denom = d; return result;
}
Frac negateFrac(Frac f) {
return makeFrac(-f.num, f.denom);
}
Frac subFrac(Frac lhs, Frac rhs) {
return makeFrac(lhs.num * rhs.denom - lhs.denom * rhs.num, rhs.denom * lhs.denom);
}
Frac multFrac(Frac lhs, Frac rhs) {
return makeFrac(lhs.num * rhs.num, lhs.denom * rhs.denom);
}
bool equalFrac(Frac lhs, Frac rhs) {
return (lhs.num == rhs.num) && (lhs.denom == rhs.denom);
}
bool lessFrac(Frac lhs, Frac rhs) {
return (lhs.num * rhs.denom) < (rhs.num * lhs.denom);
}
void printFrac(Frac f) {
char buffer[7]; int len;
if (f.denom != 1) {
snprintf(buffer, 7, "%d/%d", f.num, f.denom);
} else {
snprintf(buffer, 7, "%d", f.num);
}
len = 7 - strlen(buffer);
while (len-- > 0) {
putc(' ', stdout);
}
printf(buffer);
}
Frac bernoulli(int n) {
Frac a[16]; int j, m;
if (n < 0) {
a[0].num = 0;
a[0].denom = 0;
return a[0];
}
for (m = 0; m <= n; ++m) {
a[m] = makeFrac(1, m + 1);
for (j = m; j >= 1; --j) {
a[j - 1] = multFrac(subFrac(a[j - 1], a[j]), makeFrac(j, 1));
}
}
if (n != 1) {
return a[0];
}
return negateFrac(a[0]);
}
void faulhaber(int p) {
Frac q, *coeffs; int j, sign;
coeffs = malloc(sizeof(Frac)*(p + 1));
q = makeFrac(1, p + 1);
sign = -1;
for (j = 0; j <= p; ++j) {
sign = -1 * sign;
coeffs[p - j] = multFrac(multFrac(multFrac(q, makeFrac(sign, 1)), makeFrac(binomial(p + 1, j), 1)), bernoulli(j));
}
for (j = 0; j <= p; ++j) {
printFrac(coeffs[j]);
}
printf("\n");
free(coeffs);
}
int main() {
int i;
for (i = 0; i < 10; ++i) {
faulhaber(i);
}
return 0;
}</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
C#
<lang csharp>using System;
namespace FaulhabersTriangle {
internal class Frac {
private long num;
private long denom;
public static readonly Frac ZERO = new Frac(0, 1);
public static readonly Frac ONE = new Frac(1, 1);
public Frac(long n, long d) {
if (d == 0) {
throw new ArgumentException("d must not be zero");
}
long nn = n;
long dd = d;
if (nn == 0) {
dd = 1;
}
else if (dd < 0) {
nn = -nn;
dd = -dd;
}
long g = Math.Abs(Gcd(nn, dd));
if (g > 1) {
nn /= g;
dd /= g;
}
num = nn;
denom = dd;
}
private static long Gcd(long a, long b) {
if (b == 0) {
return a;
}
return Gcd(b, a % b);
}
public static Frac operator -(Frac self) {
return new Frac(-self.num, self.denom);
}
public static Frac operator +(Frac lhs, Frac rhs) {
return new Frac(lhs.num * rhs.denom + lhs.denom * rhs.num, rhs.denom * lhs.denom);
}
public static Frac operator -(Frac lhs, Frac rhs) {
return lhs + -rhs;
}
public static Frac operator *(Frac lhs, Frac rhs) {
return new Frac(lhs.num * rhs.num, lhs.denom * rhs.denom);
}
public static bool operator <(Frac lhs, Frac rhs) {
double x = (double)lhs.num / lhs.denom;
double y = (double)rhs.num / rhs.denom;
return x < y;
}
public static bool operator >(Frac lhs, Frac rhs) {
double x = (double)lhs.num / lhs.denom;
double y = (double)rhs.num / rhs.denom;
return x > y;
}
public static bool operator ==(Frac lhs, Frac rhs) {
return lhs.num == rhs.num && lhs.denom == rhs.denom;
}
public static bool operator !=(Frac lhs, Frac rhs) {
return lhs.num != rhs.num || lhs.denom != rhs.denom;
}
public override string ToString() {
if (denom == 1) {
return num.ToString();
}
return string.Format("{0}/{1}", num, denom);
}
public override bool Equals(object obj) {
var frac = obj as Frac;
return frac != null &&
num == frac.num &&
denom == frac.denom;
}
public override int GetHashCode() {
var hashCode = 1317992671;
hashCode = hashCode * -1521134295 + num.GetHashCode();
hashCode = hashCode * -1521134295 + denom.GetHashCode();
return hashCode;
}
}
class Program {
static Frac Bernoulli(int n) {
if (n < 0) {
throw new ArgumentException("n may not be negative or zero");
}
Frac[] a = new Frac[n + 1];
for (int m = 0; m <= n; m++) {
a[m] = new Frac(1, m + 1);
for (int j = m; j >= 1; j--) {
a[j - 1] = (a[j - 1] - a[j]) * new Frac(j, 1);
}
}
// returns 'first' Bernoulli number
if (n != 1) return a[0];
return -a[0];
}
static int Binomial(int n, int k) {
if (n < 0 || k < 0 || n < k) {
throw new ArgumentException();
}
if (n == 0 || k == 0) return 1;
int num = 1;
for (int i = k + 1; i <= n; i++) {
num = num * i;
}
int denom = 1;
for (int i = 2; i <= n - k; i++) {
denom = denom * i;
}
return num / denom;
}
static Frac[] FaulhaberTriangle(int p) {
Frac[] coeffs = new Frac[p + 1];
for (int i = 0; i < p + 1; i++) {
coeffs[i] = Frac.ZERO;
}
Frac q = new Frac(1, p + 1);
int sign = -1;
for (int j = 0; j <= p; j++) {
sign *= -1;
coeffs[p - j] = q * new Frac(sign, 1) * new Frac(Binomial(p + 1, j), 1) * Bernoulli(j);
}
return coeffs;
}
static void Main(string[] args) {
for (int i = 0; i < 10; i++) {
Frac[] coeffs = FaulhaberTriangle(i);
foreach (Frac coeff in coeffs) {
Console.Write("{0,5} ", coeff);
}
Console.WriteLine();
}
}
}
}</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
C++
Uses C++ 17 <lang cpp>#include <exception>
- include <iomanip>
- include <iostream>
- include <numeric>
- include <sstream>
- include <vector>
class Frac { public:
Frac() : num(0), denom(1) {}
Frac(int n, int d) { if (d == 0) { throw std::runtime_error("d must not be zero"); }
int sign_of_d = d < 0 ? -1 : 1; int g = std::gcd(n, d);
num = sign_of_d * n / g;
denom = sign_of_d * d / g;
}
Frac operator-() const { return Frac(-num, denom); }
Frac operator+(const Frac& rhs) const { return Frac(num*rhs.denom + denom * rhs.num, rhs.denom*denom); }
Frac operator-(const Frac& rhs) const { return Frac(num*rhs.denom - denom * rhs.num, rhs.denom*denom); }
Frac operator*(const Frac& rhs) const { return Frac(num*rhs.num, denom*rhs.denom); }
Frac operator*(int rhs) const { return Frac(num * rhs, denom); }
friend std::ostream& operator<<(std::ostream&, const Frac&);
private: int num; int denom; };
std::ostream & operator<<(std::ostream & os, const Frac &f) { if (f.num == 0 || f.denom == 1) { return os << f.num; }
std::stringstream ss; ss << f.num << "/" << f.denom; return os << ss.str(); }
Frac bernoulli(int n) { if (n < 0) { throw std::runtime_error("n may not be negative or zero"); }
std::vector<Frac> a; for (int m = 0; m <= n; m++) { a.push_back(Frac(1, m + 1)); for (int j = m; j >= 1; j--) { a[j - 1] = (a[j - 1] - a[j]) * j; } }
// returns 'first' Bernoulli number if (n != 1) return a[0]; return -a[0]; }
int binomial(int n, int k) { if (n < 0 || k < 0 || n < k) { throw std::runtime_error("parameters are invalid"); } if (n == 0 || k == 0) return 1;
int num = 1; for (int i = k + 1; i <= n; i++) { num *= i; }
int denom = 1; for (int i = 2; i <= n - k; i++) { denom *= i; }
return num / denom; }
std::vector<Frac> faulhaberTraingle(int p) { std::vector<Frac> coeffs(p + 1);
Frac q{ 1, p + 1 }; int sign = -1; for (int j = 0; j <= p; j++) { sign *= -1; coeffs[p - j] = q * sign * binomial(p + 1, j) * bernoulli(j); }
return coeffs; }
int main() {
for (int i = 0; i < 10; i++) { std::vector<Frac> coeffs = faulhaberTraingle(i); for (auto frac : coeffs) { std::cout << std::right << std::setw(5) << frac << " "; } std::cout << std::endl; }
return 0; }</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
D
<lang D>import std.algorithm : fold; import std.conv : to; import std.exception : enforce; import std.format : formattedWrite; import std.numeric : cmp, gcd; import std.range : iota; import std.stdio; import std.traits;
auto abs(T)(T val) if (isNumeric!T) {
if (val < 0) {
return -val;
}
return val;
}
struct Frac {
long num; long denom;
enum ZERO = Frac(0, 1); enum ONE = Frac(1, 1);
this(long n, long d) in {
enforce(d != 0, "Parameter d may not be zero.");
} body {
auto nn = n;
auto dd = d;
if (nn == 0) {
dd = 1;
} else if (dd < 0) {
nn = -nn;
dd = -dd;
}
auto g = gcd(abs(nn), abs(dd));
if (g > 1) {
nn /= g;
dd /= g;
}
num = nn;
denom = dd;
}
auto opBinary(string op)(Frac rhs) const {
static if (op == "+" || op == "-") {
return mixin("Frac(num*rhs.denom"~op~"denom*rhs.num, rhs.denom*denom)");
} else if (op == "*") {
return Frac(num*rhs.num, denom*rhs.denom);
}
}
auto opUnary(string op : "-")() const {
return Frac(-num, denom);
}
int opCmp(Frac rhs) const {
return cmp(cast(real) this, cast(real) rhs);
}
bool opEquals(Frac rhs) const {
return num == rhs.num && denom == rhs.denom;
}
void toString(scope void delegate(const(char)[]) sink) const {
if (denom == 1) {
formattedWrite(sink, "%d", num);
} else {
formattedWrite(sink, "%d/%s", num, denom);
}
}
T opCast(T)() const if (isFloatingPoint!T) {
return cast(T) num / denom;
}
}
auto abs(Frac f) {
if (f.num >= 0) {
return f;
}
return -f;
}
auto bernoulli(int n) in {
enforce(n >= 0, "Parameter n must not be negative.");
} body {
Frac[] a;
a.length = n+1;
a[0] = Frac.ZERO;
foreach (m; 0..n+1) {
a[m] = Frac(1, m+1);
foreach_reverse (j; 1..m+1) {
a[j-1] = (a[j-1] - a[j]) * Frac(j, 1);
}
}
if (n != 1) {
return a[0];
}
return -a[0];
}
auto binomial(int n, int k) in {
enforce(n>=0 && k>=0 && n>=k);
} body {
if (n==0 || k==0) return 1; auto num = iota(k+1, n+1).fold!"a*b"(1); auto den = iota(2, n-k+1).fold!"a*b"(1); return num / den;
}
Frac[] faulhaberTriangle(int p) {
Frac[] coeffs;
coeffs.length = p+1;
coeffs[0] = Frac.ZERO;
auto q = Frac(1, p+1);
auto sign = -1;
foreach (j; 0..p+1) {
sign *= -1;
coeffs[p - j] = q * Frac(sign, 1) * Frac(binomial(p+1, j), 1) * bernoulli(j);
}
return coeffs;
}
void main() {
foreach (i; 0..10) {
auto coeffs = faulhaberTriangle(i);
foreach (coeff; coeffs) {
writef("%5s ", coeff.to!string);
}
writeln;
}
writeln;
}</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
F#
The Function
<lang fsharp> // Generate Faulhaber's Triangle. Nigel Galloway: May 8th., 2018 let Faulhaber=let fN n = (1N - List.sum n)::n
let rec Faul a b=seq{let t = fN (List.mapi(fun n g->b*g/BigRational.FromInt(n+2)) a)
yield t
yield! Faul t (b+1N)}
Faul [] 0N
</lang>
The Task
<lang fsharp> Faulhaber |> Seq.take 10 |> Seq.iter (printfn "%A") </lang>
- Output:
[1N] [1/2N; 1/2N] [1/6N; 1/2N; 1/3N] [0N; 1/4N; 1/2N; 1/4N] [-1/30N; 0N; 1/3N; 1/2N; 1/5N] [0N; -1/12N; 0N; 5/12N; 1/2N; 1/6N] [1/42N; 0N; -1/6N; 0N; 1/2N; 1/2N; 1/7N] [0N; 1/12N; 0N; -7/24N; 0N; 7/12N; 1/2N; 1/8N] [-1/30N; 0N; 2/9N; 0N; -7/15N; 0N; 2/3N; 1/2N; 1/9N] [0N; -3/20N; 0N; 1/2N; 0N; -7/10N; 0N; 3/4N; 1/2N; 1/10N]
Factor
<lang factor>USING: kernel math math.combinatorics math.extras math.functions math.ranges prettyprint sequences ;
- faulhaber ( p -- seq )
1 + dup recip swap dup 0 (a,b] [ [ nCk ] [ -1 swap ^ ] [ bernoulli ] tri * * * ] 2with map ;
10 [ faulhaber . ] each-integer</lang>
- Output:
{ 1 }
{ 1/2 1/2 }
{ 1/6 1/2 1/3 }
{ 0 1/4 1/2 1/4 }
{ -1/30 0 1/3 1/2 1/5 }
{ 0 -1/12 0 5/12 1/2 1/6 }
{ 1/42 0 -1/6 0 1/2 1/2 1/7 }
{ 0 1/12 0 -7/24 0 7/12 1/2 1/8 }
{ -1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9 }
{ 0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10 }
FreeBASIC
<lang freebasic>' version 12-08-2017 ' compile with: fbc -s console ' uses GMP
- Include Once "gmp.bi"
- Define i_max 17
Dim As UInteger i, j, x Dim As String s Dim As ZString Ptr gmp_str : gmp_str = Allocate(100)
Dim As Mpq_ptr n, tmp1, tmp2, sum, one, zero n = Allocate(Len(__mpq_struct)) : Mpq_init(n) tmp1 = Allocate(Len(__mpq_struct)) : Mpq_init(tmp1) tmp2 = Allocate(Len(__mpq_struct)) : Mpq_init(tmp2) sum = Allocate(Len(__mpq_struct)) : Mpq_init(sum) zero = Allocate(Len(__mpq_struct)) : Mpq_init(zero) one = Allocate(Len(__mpq_struct)) : Mpq_init(one) Mpq_set_ui(zero, 0, 0) ' 0/0 = 0 Mpq_set_ui(one , 1, 1) ' 1/1 = 1
Dim As Mpq_ptr Faulhaber_triangle(0 To i_max, 1 To i_max +1) ' only initialize the variables we need For i = 0 To i_max
For j = 1 To i +1
Faulhaber_triangle(i, j) = Allocate(Len(__Mpq_struct))
Mpq_init(Faulhaber_triangle(i, j))
Next
Next
Mpq_set(Faulhaber_triangle(0, 1), one)
' we calculate the first 18 rows For i = 1 To i_max
Mpq_set(sum, zero)
For j = i +1 To 2 Step -1
Mpq_set_ui(tmp1, i, j) ' i / j
Mpq_set(tmp2, Faulhaber_triangle(i -1, j -1))
Mpq_mul(Faulhaber_triangle(i, j), tmp2, tmp1)
Mpq_canonicalize(Faulhaber_triangle(i, j))
Mpq_add(sum, sum, Faulhaber_triangle(i, j))
Next
Mpq_sub(Faulhaber_triangle(i, 1), one, sum)
Next
Print "The first 10 rows" For i = 0 To 9
For j = 1 To i +1
Mpq_get_str(gmp_str, 10, Faulhaber_triangle(i, j))
s = Space(6) + *gmp_str + Space(6)
x = InStr(s,"/")
If x = 0 Then x = 7 ' in case of 0 or 1
Print Mid(s, x -3, 7);
Next
Print
Next print
' using the 17'the row Mpq_set(sum, zero) Mpq_set_ui(n, 1000, 1) ' 1000/1 = 1000 Mpq_set(tmp2, n) For j = 1 To 18
Mpq_mul(tmp1, n, Faulhaber_triangle(17, j)) Mpq_add(sum, sum, tmp1) Mpq_mul(n, n, tmp2)
Next
Mpq_get_str(gmp_str, 10, sum) Print *gmp_str
' free memory DeAllocate(gmp_str) Mpq_clear(tmp1) : Mpq_clear(tmp2) : Mpq_clear(n) Mpq_clear(zero) : Mpq_clear(one) : Mpq_clear(sum)
For i = 0 To i_max
For j = 1 To i +1
Mpq_clear(Faulhaber_triangle(i, j))
Next
Next
' empty keyboard buffer While Inkey <> "" : Wend Print : Print "hit any key to end program" Sleep End</lang>
- Output:
The first 10 rows 1 1/2 1/2 1/6 1/2 1/3 0 1/4 1/2 1/4 -1/30 0 1/3 1/2 1/5 0 -1/12 0 5/12 1/2 1/6 1/42 0 -1/6 0 1/2 1/2 1/7 0 1/12 0 -7/24 0 7/12 1/2 1/8 -1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9 0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10 56056972216555580111030077961944183400198333273050000
Fōrmulæ
Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text. Moreover, there can be multiple visual representations of the same program. Even though it is possible to have textual representation —i.e. XML, JSON— they are intended for storage and transfer purposes more than visualization and edition.
Programs in Fōrmulæ are created/edited online in its website, However they run on execution servers. By default remote servers are used, but they are limited in memory and processing power, since they are intended for demonstration and casual use. A local server can be downloaded and installed, it has no limitations (it runs in your own computer). Because of that, example programs can be fully visualized and edited, but some of them will not run if they require a moderate or heavy computation/memory resources, and no local server is being used.
In this page you can see the program(s) related to this task and their results.
Go
Except that there is no need to roll our own Frac type when we can use the big.Rat type from the Go standard library. <lang go>package main
import (
"fmt" "math/big"
)
func bernoulli(n uint) *big.Rat {
a := make([]big.Rat, n+1)
z := new(big.Rat)
for m := range a {
a[m].SetFrac64(1, int64(m+1))
for j := m; j >= 1; j-- {
d := &a[j-1]
d.Mul(z.SetInt64(int64(j)), d.Sub(d, &a[j]))
}
}
// return the 'first' Bernoulli number
if n != 1 {
return &a[0]
}
a[0].Neg(&a[0])
return &a[0]
}
func binomial(n, k int) int64 {
if n <= 0 || k <= 0 || n < k {
return 1
}
var num, den int64 = 1, 1
for i := k + 1; i <= n; i++ {
num *= int64(i)
}
for i := 2; i <= n-k; i++ {
den *= int64(i)
}
return num / den
}
func faulhaberTriangle(p int) []big.Rat {
coeffs := make([]big.Rat, p+1)
q := big.NewRat(1, int64(p)+1)
t := new(big.Rat)
u := new(big.Rat)
sign := -1
for j := range coeffs {
sign *= -1
d := &coeffs[p-j]
t.SetInt64(int64(sign))
u.SetInt64(binomial(p+1, j))
d.Mul(q, t)
d.Mul(d, u)
d.Mul(d, bernoulli(uint(j)))
}
return coeffs
}
func main() {
for i := 0; i < 10; i++ {
coeffs := faulhaberTriangle(i)
for _, coeff := range coeffs {
fmt.Printf("%5s ", coeff.RatString())
}
fmt.Println()
}
fmt.Println()
// get coeffs for (k + 1)th row
k := 17
cc := faulhaberTriangle(k)
n := int64(1000)
nn := big.NewRat(n, 1)
np := big.NewRat(1, 1)
sum := new(big.Rat)
tmp := new(big.Rat)
for _, c := range cc {
np.Mul(np, nn)
tmp.Set(np)
tmp.Mul(tmp, &c)
sum.Add(sum, tmp)
}
fmt.Println(sum.RatString())
}</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
56056972216555580111030077961944183400198333273050000
Groovy
<lang groovy>import java.math.MathContext import java.util.stream.LongStream
class FaulhabersTriangle {
private static final MathContext MC = new MathContext(256)
private static long gcd(long a, long b) {
if (b == 0) {
return a
}
return gcd(b, a % b)
}
private static class Frac implements Comparable<Frac> {
private long num
private long denom
public static final Frac ZERO = new Frac(0, 1)
Frac(long n, long d) {
if (d == 0) throw new IllegalArgumentException("d must not be zero")
long nn = n
long dd = d
if (nn == 0) {
dd = 1
} else if (dd < 0) {
nn = -nn
dd = -dd
}
long g = Math.abs(gcd(nn, dd))
if (g > 1) {
nn /= g
dd /= g
}
num = nn
denom = dd
}
Frac plus(Frac rhs) {
return new Frac(num * rhs.denom + denom * rhs.num, rhs.denom * denom)
}
Frac negative() {
return new Frac(-num, denom)
}
Frac minus(Frac rhs) {
return this + -rhs
}
Frac multiply(Frac rhs) {
return new Frac(this.num * rhs.num, this.denom * rhs.denom)
}
@Override
int compareTo(Frac o) {
double diff = toDouble() - o.toDouble()
return Double.compare(diff, 0.0)
}
@Override
boolean equals(Object obj) {
return null != obj && obj instanceof Frac && this == (Frac) obj
}
@Override
String toString() {
if (denom == 1) {
return Long.toString(num)
}
return String.format("%d/%d", num, denom)
}
double toDouble() {
return (double) num / denom
}
BigDecimal toBigDecimal() {
return BigDecimal.valueOf(num).divide(BigDecimal.valueOf(denom), MC)
}
}
private static Frac bernoulli(int n) {
if (n < 0) throw new IllegalArgumentException("n may not be negative or zero")
Frac[] a = new Frac[n + 1]
Arrays.fill(a, Frac.ZERO)
for (int m = 0; m <= n; ++m) {
a[m] = new Frac(1, m + 1)
for (int j = m; j >= 1; --j) {
a[j - 1] = (a[j - 1] - a[j]) * new Frac(j, 1)
}
}
// returns 'first' Bernoulli number
if (n != 1) return a[0]
return -a[0]
}
private static long binomial(int n, int k) {
if (n < 0 || k < 0 || n < k) throw new IllegalArgumentException()
if (n == 0 || k == 0) return 1
long num = LongStream.rangeClosed(k + 1, n).reduce(1, { a, b -> a * b })
long den = LongStream.rangeClosed(2, n - k).reduce(1, { acc, i -> acc * i })
return num / den
}
private static Frac[] faulhaberTriangle(int p) {
Frac[] coeffs = new Frac[p + 1]
Arrays.fill(coeffs, Frac.ZERO)
Frac q = new Frac(1, p + 1)
int sign = -1
for (int j = 0; j <= p; ++j) {
sign *= -1
coeffs[p - j] = q * new Frac(sign, 1) * new Frac(binomial(p + 1, j), 1) * bernoulli(j)
}
return coeffs
}
static void main(String[] args) {
for (int i = 0; i <= 9; ++i) {
Frac[] coeffs = faulhaberTriangle(i)
for (Frac coeff : coeffs) {
printf("%5s ", coeff)
}
println()
}
println()
// get coeffs for (k + 1)th row
int k = 17
Frac[] cc = faulhaberTriangle(k)
int n = 1000
BigDecimal nn = BigDecimal.valueOf(n)
BigDecimal np = BigDecimal.ONE
BigDecimal sum = BigDecimal.ZERO
for (Frac c : cc) {
np = np * nn
sum = sum.add(np * c.toBigDecimal())
}
println(sum.toBigInteger())
}
}</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
Haskell
<lang haskell>import Data.Ratio (Ratio, denominator, numerator, (%))
FAULHABER -----------------------
faulhaber :: Int -> Rational -> Rational faulhaber p n =
sum $ zipWith ((*) . (n ^)) [1 ..] (faulhaberTriangle !! p)
faulhaberTriangle :: Rational
faulhaberTriangle =
tail $
scanl
( \rs n ->
let xs = zipWith ((*) . (n %)) [2 ..] rs
in 1 - sum xs : xs
)
[]
[0 ..]
TEST -------------------------
main :: IO () main = do
let triangle = take 10 faulhaberTriangle
widths = maxWidths triangle
mapM_
putStrLn
[ unlines
( (justifyRatio widths 8 ' ' =<<)
<$> triangle
),
(show . numerator) (faulhaber 17 1000)
]
DISPLAY ------------------------
justifyRatio ::
(Int, Int) -> Int -> Char -> Rational -> String
justifyRatio (wn, wd) n c nd =
go $
[numerator, denominator] <*> [nd]
where
-- Minimum column width, or more if specified.
w = max n (wn + wd + 2)
go [num, den]
| 1 == den = center w c (show num)
| otherwise =
let (q, r) = quotRem (w - 1) 2
in concat
[ justifyRight q c (show num),
"/",
justifyLeft (q + r) c (show den)
]
justifyLeft :: Int -> a -> [a] -> [a] justifyLeft n c s = take n (s <> replicate n c)
justifyRight :: Int -> a -> [a] -> [a] justifyRight n c = (drop . length) <*> (replicate n c <>)
center :: Int -> a -> [a] -> [a] center n c s =
let (q, r) = quotRem (n - length s) 2
pad = replicate q c
in concat [pad, s, pad, replicate r c]
maxWidths :: Rational -> (Int, Int) maxWidths xss =
let widest f xs = maximum $ fmap (length . show . f) xs
in ((,) . widest numerator <*> widest denominator) $
concat xss</lang>
- Output:
1 1/2 1/2 1/6 1/2 1/3 0 1/4 1/2 1/4 -1/30 0 1/3 1/2 1/5 0 -1/12 0 5/12 1/2 1/6 1/42 0 -1/6 0 1/2 1/2 1/7 0 1/12 0 -7/24 0 7/12 1/2 1/8 -1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9 0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10 56056972216555580111030077961944183400198333273050000
J
faulhaberTriangle=: ([: %. [: x: (1 _2 (p.) 2 | +/~) * >:/~ * (!~/~ >:))@:i.
faulhaberTriangle 10
1 0 0 0 0 0 0 0 0 0
1r2 1r2 0 0 0 0 0 0 0 0
1r6 1r2 1r3 0 0 0 0 0 0 0
0 1r4 1r2 1r4 0 0 0 0 0 0
_1r30 0 1r3 1r2 1r5 0 0 0 0 0
0 _1r12 0 5r12 1r2 1r6 0 0 0 0
1r42 0 _1r6 0 1r2 1r2 1r7 0 0 0
0 1r12 0 _7r24 0 7r12 1r2 1r8 0 0
_1r30 0 2r9 0 _7r15 0 2r3 1r2 1r9 0
0 _3r20 0 1r2 0 _7r10 0 3r4 1r2 1r10
NB.--------------------------------------------------------------
NB. matrix inverse (%.) of the extended precision (x:) product of
(!~/~ >:)@:i. 5 NB. Pascal's triangle
1 1 0 0 0
1 2 1 0 0
1 3 3 1 0
1 4 6 4 1
1 5 10 10 5
NB. second task in progress
(>:/~)@: i. 5 NB. lower triangle
1 0 0 0 0
1 1 0 0 0
1 1 1 0 0
1 1 1 1 0
1 1 1 1 1
(1 _2 p. (2 | +/~))@:i. 5 NB. 1 + (-2) * (parity table)
1 _1 1 _1 1
_1 1 _1 1 _1
1 _1 1 _1 1
_1 1 _1 1 _1
1 _1 1 _1 1
Java
<lang Java>import java.math.BigDecimal; import java.math.MathContext; import java.util.Arrays; import java.util.stream.LongStream;
public class FaulhabersTriangle {
private static final MathContext MC = new MathContext(256);
private static long gcd(long a, long b) {
if (b == 0) {
return a;
}
return gcd(b, a % b);
}
private static class Frac implements Comparable<Frac> {
private long num;
private long denom;
public static final Frac ZERO = new Frac(0, 1);
public Frac(long n, long d) {
if (d == 0) throw new IllegalArgumentException("d must not be zero");
long nn = n;
long dd = d;
if (nn == 0) {
dd = 1;
} else if (dd < 0) {
nn = -nn;
dd = -dd;
}
long g = Math.abs(gcd(nn, dd));
if (g > 1) {
nn /= g;
dd /= g;
}
num = nn;
denom = dd;
}
public Frac plus(Frac rhs) {
return new Frac(num * rhs.denom + denom * rhs.num, rhs.denom * denom);
}
public Frac unaryMinus() {
return new Frac(-num, denom);
}
public Frac minus(Frac rhs) {
return this.plus(rhs.unaryMinus());
}
public Frac times(Frac rhs) {
return new Frac(this.num * rhs.num, this.denom * rhs.denom);
}
@Override
public int compareTo(Frac o) {
double diff = toDouble() - o.toDouble();
return Double.compare(diff, 0.0);
}
@Override
public boolean equals(Object obj) {
return null != obj && obj instanceof Frac && this.compareTo((Frac) obj) == 0;
}
@Override
public String toString() {
if (denom == 1) {
return Long.toString(num);
}
return String.format("%d/%d", num, denom);
}
public double toDouble() {
return (double) num / denom;
}
public BigDecimal toBigDecimal() {
return BigDecimal.valueOf(num).divide(BigDecimal.valueOf(denom), MC);
}
}
private static Frac bernoulli(int n) {
if (n < 0) throw new IllegalArgumentException("n may not be negative or zero");
Frac[] a = new Frac[n + 1];
Arrays.fill(a, Frac.ZERO);
for (int m = 0; m <= n; ++m) {
a[m] = new Frac(1, m + 1);
for (int j = m; j >= 1; --j) {
a[j - 1] = a[j - 1].minus(a[j]).times(new Frac(j, 1));
}
}
// returns 'first' Bernoulli number
if (n != 1) return a[0];
return a[0].unaryMinus();
}
private static long binomial(int n, int k) {
if (n < 0 || k < 0 || n < k) throw new IllegalArgumentException();
if (n == 0 || k == 0) return 1;
long num = LongStream.rangeClosed(k + 1, n).reduce(1, (a, b) -> a * b);
long den = LongStream.rangeClosed(2, n - k).reduce(1, (acc, i) -> acc * i);
return num / den;
}
private static Frac[] faulhaberTriangle(int p) {
Frac[] coeffs = new Frac[p + 1];
Arrays.fill(coeffs, Frac.ZERO);
Frac q = new Frac(1, p + 1);
int sign = -1;
for (int j = 0; j <= p; ++j) {
sign *= -1;
coeffs[p - j] = q.times(new Frac(sign, 1)).times(new Frac(binomial(p + 1, j), 1)).times(bernoulli(j));
}
return coeffs;
}
public static void main(String[] args) {
for (int i = 0; i <= 9; ++i) {
Frac[] coeffs = faulhaberTriangle(i);
for (Frac coeff : coeffs) {
System.out.printf("%5s ", coeff);
}
System.out.println();
}
System.out.println();
// get coeffs for (k + 1)th row
int k = 17;
Frac[] cc = faulhaberTriangle(k);
int n = 1000;
BigDecimal nn = BigDecimal.valueOf(n);
BigDecimal np = BigDecimal.ONE;
BigDecimal sum = BigDecimal.ZERO;
for (Frac c : cc) {
np = np.multiply(nn);
sum = sum.add(np.multiply(c.toBigDecimal()));
}
System.out.println(sum.toBigInteger());
}
}</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
JavaScript
ES6
JavaScript is probably not the right instrument to choose for this task, which requires both a ratio number type and arbitrary precision integers. JavaScript has neither – its only numeric datatype is the IEEE 754 double-precision floating-point format number, into which integers and all else must fit. (See the built-in JS name Number.MAX_SAFE_INTEGER)
This means that we can print Faulhaber's triangle (hand-coding some rudimentary ratio-arithmetic functions), but our only reward for evaluating faulhaber(17, 1000) is an integer overflow. With JS integers out of the box, we can get about as far as faulhaber(17, 8), or faulhaber(4, 1000).
(Further progress would entail implementing some hand-crafted representation of arbitrary precision integers – perhaps a bit beyond the intended scope of this task, and good enough motivation to use a different language) <lang JavaScript>(() => {
// Order of Faulhaber's triangle -> rows of Faulhaber's triangle // faulHaberTriangle :: Int -> Ratio Int const faulhaberTriangle = n => map(x => tail( scanl((a, x) => { const ys = map((nd, i) => ratioMult(nd, Ratio(x, i + 2)), a); return cons(ratioMinus(Ratio(1, 1), ratioSum(ys)), ys); }, [], enumFromTo(0, x)) ), enumFromTo(0, n));
// p -> n -> Sum of the p-th powers of the first n positive integers
// faulhaber :: Int -> Ratio Int -> Ratio Int
const faulhaber = (p, n) =>
ratioSum(map(
(nd, i) => ratioMult(nd, Ratio(raise(n, i + 1), 1)),
last(faulhaberTriangle(p))
));
// RATIOS -----------------------------------------------------------------
// (Max numr + denr widths) -> Column width -> Filler -> Ratio -> String
// justifyRatio :: (Int, Int) -> Int -> Char -> Ratio Integer -> String
const justifyRatio = (ws, n, c, nd) => {
const
w = max(n, ws.nMax + ws.dMax + 2),
[num, den] = [nd.num, nd.den];
return all(Number.isSafeInteger, [num, den]) ? (
den === 1 ? center(w, c, show(num)) : (() => {
const [q, r] = quotRem(w - 1, 2);
return concat([
justifyRight(q, c, show(num)),
'/',
justifyLeft(q + r, c, (show(den)))
]);
})()
) : "JS integer overflow ... ";
};
// Ratio :: Int -> Int -> Ratio
const Ratio = (n, d) => ({
num: n,
den: d
});
// ratioMinus :: Ratio -> Ratio -> Ratio
const ratioMinus = (nd, nd1) => {
const
d = lcm(nd.den, nd1.den);
return simpleRatio({
num: (nd.num * (d / nd.den)) - (nd1.num * (d / nd1.den)),
den: d
});
};
// ratioMult :: Ratio -> Ratio -> Ratio
const ratioMult = (nd, nd1) => simpleRatio({
num: nd.num * nd1.num,
den: nd.den * nd1.den
});
// ratioPlus :: Ratio -> Ratio -> Ratio
const ratioPlus = (nd, nd1) => {
const
d = lcm(nd.den, nd1.den);
return simpleRatio({
num: (nd.num * (d / nd.den)) + (nd1.num * (d / nd1.den)),
den: d
});
};
// ratioSum :: [Ratio] -> Ratio
const ratioSum = xs =>
simpleRatio(foldl((a, x) => ratioPlus(a, x), {
num: 0,
den: 1
}, xs));
// ratioWidths :: Ratio -> {nMax::Int, dMax::Int} const ratioWidths = xss => { return foldl((a, x) => { const [nw, dw] = ap( [compose(length, show)], [x.num, x.den] ), [an, ad] = ap( [curry(flip(lookup))(a)], ['nMax', 'dMax'] ); return { nMax: nw > an ? nw : an, dMax: dw > ad ? dw : ad }; }, { nMax: 0, dMax: 0 }, concat(xss)); };
// simpleRatio :: Ratio -> Ratio
const simpleRatio = nd => {
const g = gcd(nd.num, nd.den);
return {
num: nd.num / g,
den: nd.den / g
};
};
// GENERIC FUNCTIONS ------------------------------------------------------
// all :: (a -> Bool) -> [a] -> Bool const all = (f, xs) => xs.every(f);
// A list of functions applied to a list of arguments
// <*> :: [(a -> b)] -> [a] -> [b]
const ap = (fs, xs) => //
[].concat.apply([], fs.map(f => //
[].concat.apply([], xs.map(x => [f(x)]))));
// Size of space -> filler Char -> Text -> Centered Text
// center :: Int -> Char -> Text -> Text
const center = (n, c, s) => {
const [q, r] = quotRem(n - s.length, 2);
return concat(concat([replicate(q, c), s, replicate(q + r, c)]));
};
// compose :: (b -> c) -> (a -> b) -> (a -> c) const compose = (f, g) => x => f(g(x));
// concat :: a -> [a] | [String] -> String const concat = xs => xs.length > 0 ? (() => { const unit = typeof xs[0] === 'string' ? : []; return unit.concat.apply(unit, xs); })() : [];
// cons :: a -> [a] -> [a] const cons = (x, xs) => [x].concat(xs);
// 2 or more arguments
// curry :: Function -> Function
const curry = (f, ...args) => {
const go = xs => xs.length >= f.length ? (f.apply(null, xs)) :
function () {
return go(xs.concat(Array.from(arguments)));
};
return go([].slice.call(args, 1));
};
// enumFromTo :: Int -> Int -> [Int]
const enumFromTo = (m, n) =>
Array.from({
length: Math.floor(n - m) + 1
}, (_, i) => m + i);
// flip :: (a -> b -> c) -> b -> a -> c const flip = f => (a, b) => f.apply(null, [b, a]);
// foldl :: (b -> a -> b) -> b -> [a] -> b const foldl = (f, a, xs) => xs.reduce(f, a);
// gcd :: Integral a => a -> a -> a
const gcd = (x, y) => {
const _gcd = (a, b) => (b === 0 ? a : _gcd(b, a % b)),
abs = Math.abs;
return _gcd(abs(x), abs(y));
};
// head :: [a] -> a const head = xs => xs.length ? xs[0] : undefined;
// intercalate :: String -> [a] -> String const intercalate = (s, xs) => xs.join(s);
// justifyLeft :: Int -> Char -> Text -> Text
const justifyLeft = (n, cFiller, strText) =>
n > strText.length ? (
(strText + cFiller.repeat(n))
.substr(0, n)
) : strText;
// justifyRight :: Int -> Char -> Text -> Text
const justifyRight = (n, cFiller, strText) =>
n > strText.length ? (
(cFiller.repeat(n) + strText)
.slice(-n)
) : strText;
// last :: [a] -> a const last = xs => xs.length ? xs.slice(-1)[0] : undefined;
// length :: [a] -> Int const length = xs => xs.length;
// lcm :: Integral a => a -> a -> a
const lcm = (x, y) =>
(x === 0 || y === 0) ? 0 : Math.abs(Math.floor(x / gcd(x, y)) * y);
// lookup :: Eq a => a -> [(a, b)] -> Maybe b
const lookup = (k, pairs) => {
if (Array.isArray(pairs)) {
let m = pairs.find(x => x[0] === k);
return m ? m[1] : undefined;
} else {
return typeof pairs === 'object' ? (
pairs[k]
) : undefined;
}
};
// map :: (a -> b) -> [a] -> [b] const map = (f, xs) => xs.map(f);
// max :: Ord a => a -> a -> a const max = (a, b) => b > a ? b : a;
// min :: Ord a => a -> a -> a const min = (a, b) => b < a ? b : a;
// quotRem :: Integral a => a -> a -> (a, a) const quotRem = (m, n) => [Math.floor(m / n), m % n];
// raise :: Num -> Int -> Num const raise = (n, e) => Math.pow(n, e);
// replicate :: Int -> a -> [a]
const replicate = (n, x) =>
Array.from({
length: n
}, () => x);
// scanl :: (b -> a -> b) -> b -> [a] -> [b]
const scanl = (f, startValue, xs) =>
xs.reduce((a, x) => {
const v = f(a.acc, x);
return {
acc: v,
scan: cons(a.scan, v)
};
}, {
acc: startValue,
scan: [startValue]
})
.scan;
// show :: a -> String
const show = (...x) =>
JSON.stringify.apply(
null, x.length > 1 ? [x[0], null, x[1]] : x
);
// tail :: [a] -> [a] const tail = xs => xs.length ? xs.slice(1) : undefined;
// unlines :: [String] -> String
const unlines = xs => xs.join('\n');
// TEST -------------------------------------------------------------------
const
triangle = faulhaberTriangle(9),
widths = ratioWidths(triangle);
return unlines(
map(row =>
concat(map(cell =>
justifyRatio(widths, 8, ' ', cell), row)), triangle)
) +
'\n\n' + unlines(
[
'faulhaber(17, 1000)',
justifyRatio(widths, 0, ' ', faulhaber(17, 1000)),
'\nfaulhaber(17, 8)',
justifyRatio(widths, 0, ' ', faulhaber(17, 8)),
'\nfaulhaber(4, 1000)',
justifyRatio(widths, 0, ' ', faulhaber(4, 1000)),
]
);
})();</lang>
- Output:
1 1/2 1/2 1/6 1/2 1/3 0 1/4 1/2 1/4 -1/30 0 1/3 1/2 1/5 0 -1/12 0 5/12 1/2 1/6 1/42 0 -1/6 0 1/2 1/2 1/7 0 1/12 0 -7/24 0 7/12 1/2 1/8 -1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9 0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10 faulhaber(17, 1000) JS integer overflow ... faulhaber(17, 8) 2502137235710736 faulhaber(4, 1000) 200500333333300
jq
Works with jq (*)
Works with gojq, the Go implementation of jq
The solution presented here requires a "rational arithmetic" package, such as the "Rational" module at Arithmetic/Rational#jq. As a reminder, the jq directive for including such a package appears as the first line of the program below.
Note also that the function `bernoulli` is defined here in a way that ensures B(1) is `1 // 2`.
(*) The C implementation of jq does not have sufficient numeric precision for the "extra credit" task. <lang jq>include "Rational";
- Preliminaries
def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;
- for gojq
def idivide($j):
. as $i | ($i % $j) as $mod | ($i - $mod) / $j ;
- use idivide for precision
def binomial(n; k):
if k > n / 2 then binomial(n; n-k) else reduce range(1; k+1) as $i (1; . * (n - $i + 1) | idivide($i)) end;
- Here we conform to the modern view that B(1) is 1 // 2
def bernoulli:
if type != "number" or . < 0 then "bernoulli must be given a non-negative number vs \(.)" | error
else . as $n
| reduce range(0; $n+1) as $i ([];
.[$i] = r(1; $i + 1)
| reduce range($i; 0; -1) as $j (.;
.[$j-1] = rmult($j; rminus(.[$j-1]; .[$j])) ) )
| .[0] # the modern view
end;
- Input: a non-negative integer, $p
- Output: an array of Rationals corresponding to the
- Faulhaber coefficients for row ($p + 1) (counting the first row as row 1).
def faulhabercoeffs:
def altBernoulli: # adjust B(1) for this task
bernoulli as $b
| if . == 1 then rmult(-1; $b) else $b end;
. as $p
| r(1; $p + 1) as $q
| { coeffs: [], sign: -1 }
| reduce range(0; $p+1) as $j (.;
.sign *= -1
| binomial($p + 1; $j) as $b
| .coeffs[$p - $j] = ([ .sign, $q, $b, ($j|altBernoulli) ] | rmult))
| .coeffs
- Calculate the sum for ($k|faulhabercoeffs)
def faulhabersum($n; $k):
($k|faulhabercoeffs) as $coe
| reduce range(0;$k+1) as $i ({sum: 0, power: 1};
.power *= $n
| .sum = radd(.sum; rmult(.power; $coe[$i]))
)
| .sum;
- pretty print a Rational assumed to have the {n,d} form
def rpp:
if .n == 0 then "0" elif .d == 1 then .n | tostring else "\(.n)/\(.d)" end;
def testfaulhaber:
(range(0; 10) as $i
| ($i | faulhabercoeffs | map(rpp | lpad(6)) | join(" "))),
"\nfaulhabersum(1000; 17):",
(faulhabersum(1000; 17) | rpp) ;
testfaulhaber</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
faulhabersum(1000; 17):
56056972216555580111030077961944183400198333273050000
Julia
<lang julia>function bernoulli(n)
A = Vector{Rational{BigInt}}(undef, n + 1)
for i in 0:n
A[i + 1] = 1 // (i + 1)
for j = i:-1:1
A[j] = j * (A[j] - A[j + 1])
end
end
return n == 1 ? -A[1] : A[1]
end
function faulhabercoeffs(p)
coeffs = Vector{Rational{BigInt}}(undef, p + 1)
q = Rational{BigInt}(1, p + 1)
sign = -1
for j in 0:p
sign *= -1
coeffs[p - j + 1] = bernoulli(j) * (q * sign) * Rational{BigInt}(binomial(p + 1, j), 1)
end
coeffs
end
faulhabersum(n, k) = begin coe = faulhabercoeffs(k); mapreduce(i -> BigInt(n)^i * coe[i], +, 1:k+1) end
prettyfrac(x) = (x.num == 0 ? "0" : x.den == 1 ? string(x.num) : replace(string(x), "//" => "/"))
function testfaulhaber()
for i in 0:9
for c in faulhabercoeffs(i)
print(prettyfrac(c), "\t")
end
println()
end
println("\n", prettyfrac(faulhabersum(1000, 17)))
end
testfaulhaber()
</lang>
- Output:
1 1/2 1/2 1/6 1/2 1/3 0 1/4 1/2 1/4 -1/30 0 1/3 1/2 1/5 0 -1/12 0 5/12 1/2 1/6 1/42 0 -1/6 0 1/2 1/2 1/7 0 1/12 0 -7/24 0 7/12 1/2 1/8 -1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9 0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10 56056972216555580111030077961944183400198333273050000
Kotlin
Uses appropriately modified code from the Faulhaber's Formula task: <lang scala>// version 1.1.2
import java.math.BigDecimal import java.math.MathContext
val mc = MathContext(256)
fun gcd(a: Long, b: Long): Long = if (b == 0L) a else gcd(b, a % b)
class Frac : Comparable<Frac> {
val num: Long val denom: Long
companion object {
val ZERO = Frac(0, 1)
val ONE = Frac(1, 1)
}
constructor(n: Long, d: Long) {
require(d != 0L)
var nn = n
var dd = d
if (nn == 0L) {
dd = 1
}
else if (dd < 0) {
nn = -nn
dd = -dd
}
val g = Math.abs(gcd(nn, dd))
if (g > 1) {
nn /= g
dd /= g
}
num = nn
denom = dd
}
constructor(n: Int, d: Int) : this(n.toLong(), d.toLong())
operator fun plus(other: Frac) =
Frac(num * other.denom + denom * other.num, other.denom * denom)
operator fun unaryMinus() = Frac(-num, denom)
operator fun minus(other: Frac) = this + (-other)
operator fun times(other: Frac) = Frac(this.num * other.num, this.denom * other.denom)
fun abs() = if (num >= 0) this else -this
override fun compareTo(other: Frac): Int {
val diff = this.toDouble() - other.toDouble()
return when {
diff < 0.0 -> -1
diff > 0.0 -> +1
else -> 0
}
}
override fun equals(other: Any?): Boolean {
if (other == null || other !is Frac) return false
return this.compareTo(other) == 0
}
override fun toString() = if (denom == 1L) "$num" else "$num/$denom"
fun toDouble() = num.toDouble() / denom
fun toBigDecimal() = BigDecimal(num).divide(BigDecimal(denom), mc)
}
fun bernoulli(n: Int): Frac {
require(n >= 0)
val a = Array(n + 1) { Frac.ZERO }
for (m in 0..n) {
a[m] = Frac(1, m + 1)
for (j in m downTo 1) a[j - 1] = (a[j - 1] - a[j]) * Frac(j, 1)
}
return if (n != 1) a[0] else -a[0] // returns 'first' Bernoulli number
}
fun binomial(n: Int, k: Int): Long {
require(n >= 0 && k >= 0 && n >= k)
if (n == 0 || k == 0) return 1
val num = (k + 1..n).fold(1L) { acc, i -> acc * i }
val den = (2..n - k).fold(1L) { acc, i -> acc * i }
return num / den
}
fun faulhaberTriangle(p: Int): Array<Frac> {
val coeffs = Array(p + 1) { Frac.ZERO }
val q = Frac(1, p + 1)
var sign = -1
for (j in 0..p) {
sign *= -1
coeffs[p - j] = q * Frac(sign, 1) * Frac(binomial(p + 1, j), 1) * bernoulli(j)
}
return coeffs
}
fun main(args: Array<String>) {
for (i in 0..9){
val coeffs = faulhaberTriangle(i)
for (coeff in coeffs) print("${coeff.toString().padStart(5)} ")
println()
}
println()
// get coeffs for (k + 1)th row
val k = 17
val cc = faulhaberTriangle(k)
val n = 1000
val nn = BigDecimal(n)
var np = BigDecimal.ONE
var sum = BigDecimal.ZERO
for (c in cc) {
np *= nn
sum += np * c.toBigDecimal()
}
println(sum.toBigInteger())
}</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
56056972216555580111030077961944183400198333273050000
Lua
<lang lua>function binomial(n,k)
if n<0 or k<0 or n<k then return -1 end if n==0 or k==0 then return 1 end
local num = 1
for i=k+1,n do
num = num * i
end
local denom = 1
for i=2,n-k do
denom = denom * i
end
return num / denom
end
function gcd(a,b)
while b ~= 0 do
local temp = a % b
a = b
b = temp
end
return a
end
function makeFrac(n,d)
local result = {}
if d==0 then
result.num = 0
result.denom = 0
return result
end
if n==0 then
d = 1
elseif d < 0 then
n = -n
d = -d
end
local g = math.abs(gcd(n, d))
if g>1 then
n = n / g
d = d / g
end
result.num = n result.denom = d return result
end
function negateFrac(f)
return makeFrac(-f.num, f.denom)
end
function subFrac(lhs, rhs)
return makeFrac(lhs.num * rhs.denom - lhs.denom * rhs.num, rhs.denom * lhs.denom)
end
function multFrac(lhs, rhs)
return makeFrac(lhs.num * rhs.num, lhs.denom * rhs.denom)
end
function equalFrac(lhs, rhs)
return (lhs.num == rhs.num) and (lhs.denom == rhs.denom)
end
function lessFrac(lhs, rhs)
return (lhs.num * rhs.denom) < (rhs.num * lhs.denom)
end
function printFrac(f)
local str = tostring(f.num)
if f.denom ~= 1 then
str = str.."/"..f.denom
end
for i=1, 7 - string.len(str) do
io.write(" ")
end
io.write(str)
return nil
end
function bernoulli(n)
if n<0 then
return {num=0, denom=0}
end
local a = {}
for m=0,n do
a[m] = makeFrac(1, m+1)
for j=m,1,-1 do
a[j-1] = multFrac(subFrac(a[j-1], a[j]), makeFrac(j, 1))
end
end
if n~=1 then
return a[0]
end
return negateFrac(a[0])
end
function faulhaber(p)
local q = makeFrac(1, p+1)
local sign = -1
local coeffs = {}
for j=0,p do
sign = -1 * sign
coeffs[p-j] = multFrac(multFrac(multFrac(q, makeFrac(sign, 1)), makeFrac(binomial(p + 1, j), 1)), bernoulli(j))
end
for j=0,p do
printFrac(coeffs[j])
end
print()
return nil
end
-- main for i=0,9 do
faulhaber(i)
end</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
-0 1/4 1/2 1/4
-1/30 -0 1/3 1/2 1/5
-0 -1/12 -0 5/12 1/2 1/6
1/42 -0 -1/6 -0 1/2 1/2 1/7
-0 1/12 -0 -7/24 -0 7/12 1/2 1/8
-1/30 -0 2/9 -0 -7/15 -0 2/3 1/2 1/9
-0 -3/20 -0 1/2 -0 -7/10 -0 3/4 1/2 1/10
Mathematica / Wolfram Language
<lang Mathematica>ClearAll[Faulhaber] bernoulliB[1] := 1/2 bernoulliB[n_] := BernoulliB[n] Faulhaber[n_, p_] := 1/(p + 1) Sum[Binomial[p + 1, j] bernoulliB[j] n^(p + 1 - j), {j, 0, p}] Table[Rest@CoefficientList[Faulhaber[n, t], n], {t, 0, 9}] // Grid Faulhaber[1000, 17]</lang>
- Output:
1 1/2 1/2 1/6 1/2 1/3 0 1/4 1/2 1/4 -(1/30) 0 1/3 1/2 1/5 0 -(1/12) 0 5/12 1/2 1/6 1/42 0 -(1/6) 0 1/2 1/2 1/7 0 1/12 0 -(7/24) 0 7/12 1/2 1/8 -(1/30) 0 2/9 0 -(7/15) 0 2/3 1/2 1/9 0 -(3/20) 0 1/2 0 -(7/10) 0 3/4 1/2 1/10 56056972216555580111030077961944183400198333273050000
Nim
For the task, we could use the standard module “rationals” but for the extra task we need big numbers (and big rationals). We use third party module “bignum” for this purpose. <lang Nim>import algorithm, math, strutils import bignum
type FaulhaberSequence = seq[Rat]
- ---------------------------------------------------------------------------------------------------
func bernoulli(n: Natural): Rat =
## Return nth Bernoulli coefficient.
var a = newSeq[Rat](n + 1)
for m in 0..n:
a[m] = newRat(1, m + 1)
for k in countdown(m, 1):
a[k - 1] = (a[k - 1] - a[k]) * k
result = if n != 1: a[0] else: -a[0]
- ---------------------------------------------------------------------------------------------------
func faulhaber(n: Natural): FaulhaberSequence =
## Return nth Faulhaber sequence (high degree first).
var a = newRat(1, n + 1) var sign = -1 for k in 0..n: sign = -sign result.add(a * sign * binom(n + 1, k) * bernoulli(k))
- ---------------------------------------------------------------------------------------------------
proc display(fs: FaulhaberSequence) =
## Return the string representing a Faulhaber sequence.
var str = ""
for i, coeff in reversed(fs):
str.addSep(" ", 0)
str.add(($coeff).align(6))
echo str
- ---------------------------------------------------------------------------------------------------
func evaluate(fs: FaulhaberSequence; n: int): Rat =
## Evaluate the polynomial associated to a sequence for value "n".
result = newRat(0) for coeff in fs: result = result * n + coeff result *= n
- ———————————————————————————————————————————————————————————————————————————————————————————————————
for n in 0..9:
display(faulhaber(n))
echo "" let fs18 = faulhaber(17) # 18th row. echo fs18.evaluate(1000)</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
56056972216555580111030077961944183400198333273050000
Pascal
A console application in Free Pascal, created with the Lazarus IDE.
Row numbering below is 0-based, so row r has r+1 elements. Rather than use a rational number type, the program scales up row r by (r+1)!, which means that all the entries are integers. <lang pascal> program FaulhaberTriangle; uses uIntX, uEnums, // units in the library IntXLib4Pascal
SysUtils;
// Convert a rational num/den to a string, right-justified in the given width. // Before converting, remove any common factor of num and den. // For this application we can assume den > 0. function RationalToString( num, den : TIntX;
minWidth : integer) : string;
var
num1, den1, divisor : TIntX; w : integer;
begin
divisor := TIntX.GCD( num, den);
// TIntx.Divide requires the caller to specifiy the division mode
num1 := TIntx.Divide( num, divisor, uEnums.dmClassic);
den1 := TIntx.Divide( den, divisor, uEnums.dmClassic);
result := num1.ToString;
if not den1.IsOne then result := result + '/' + den1.ToString;
w := minWidth - Length( result);
if (w > 0) then result := StringOfChar(' ', w) + result;
end;
// Main routine const
r_MAX = 17;
var
g : array [1..r_MAX + 1] of TIntX; r, s, k : integer; r_1_fac, sum, k_intx : TIntX;
begin
// Calculate rows 0..17 of Faulhaner's triangle, and show rows 0..9.
// For a given r, the subarray g[1..(r+1)] contains (r + 1)! times row r.
r_1_fac := 1; // (r + 1)!
g[1] := 1;
for r := 0 to r_MAX do begin
r_1_fac := r_1_fac * (r+1);
sum := 0;
for s := r downto 1 do begin
g[s + 1] := r*(r+1)*g[s] div (s+1);
sum := sum + g[s + 1];
end;
g[1] := r_1_fac - sum; // the scaled row must sum to (r + 1)!
if (r <= 9) then begin
for s := 1 to r + 1 do Write( RationalToString( g[s], r_1_fac, 7));
WriteLn;
end;
end;
// Use row 17 to sum 17th powers from 1 to 1000 sum := 0; for s := r_MAX + 1 downto 1 do sum := (sum + g[s]) * 1000; sum := TIntx.Divide( sum, r_1_fac, uEnums.dmClassic); WriteLn; WriteLn( 'Sum by Faulhaber = ' + sum.ToString);
// Check by direct calculation sum := 0; for k := 1 to 1000 do begin k_intx := k; sum := sum + TIntX.Pow( k_intx, r_MAX); end; WriteLn( 'by direct calc. = ' + sum.ToString);
end. </lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
Sum by Faulhaber = 56056972216555580111030077961944183400198333273050000
by direct calc. = 56056972216555580111030077961944183400198333273050000
Perl
<lang perl>use 5.010; use List::Util qw(sum); use Math::BigRat try => 'GMP'; use ntheory qw(binomial bernfrac);
sub faulhaber_triangle {
my ($p) = @_;
map {
Math::BigRat->new(bernfrac($_))
* binomial($p, $_)
/ $p
} reverse(0 .. $p-1);
}
- First 10 rows of Faulhaber's triangle
foreach my $p (1 .. 10) {
say map { sprintf("%6s", $_) } faulhaber_triangle($p);
}
- Extra credit
my $p = 17; my $n = Math::BigInt->new(1000); my @r = faulhaber_triangle($p+1); say "\n", sum(map { $r[$_] * $n**($_ + 1) } 0 .. $#r);</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
56056972216555580111030077961944183400198333273050000
Phix
with javascript_semantics include builtins\pfrac.e -- (0.8.0+) function bernoulli(integer n) sequence a = {} for m=0 to n do a = append(a,{1,m+1}) for j=m to 1 by -1 do a[j] = frac_mul({j,1},frac_sub(a[j+1],a[j])) end for end for if n!=1 then return a[1] end if return frac_uminus(a[1]) end function function binomial(integer n, k) if n<0 or k<0 or n<k then ?9/0 end if if n=0 or k=0 then return 1 end if atom num = 1, denom = 1 for i=k+1 to n do num *= i end for for i=2 to n-k do denom *= i end for return num / denom end function function faulhaber_triangle(integer p, bool asString=true) sequence coeffs = repeat(frac_zero,p+1) for j=0 to p do frac coeff = frac_mul({binomial(p+1,j),p+1},bernoulli(j)) coeffs[p-j+1] = iff(asString?sprintf("%5s",{frac_sprint(coeff)}):coeff) end for return coeffs end function for i=0 to 9 do printf(1,"%s\n",{join(faulhaber_triangle(i)," ")}) end for puts(1,"\n") if platform()!=JS then sequence row18 = faulhaber_triangle(17,false) frac res = frac_zero atom t1 = time()+1 integer lim = 1000 for k=1 to lim do bigatom nn = BA_ONE for i=1 to length(row18) do res = frac_add(res,frac_mul(row18[i],{nn,1})) nn = ba_mul(nn,lim) end for if time()>t1 then printf(1,"calculating, k=%d...\r",k) t1 = time()+1 end if end for printf(1,"%s \n",{frac_sprint(res)}) end if
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
56056972216555580111030077961944183400198333273050000
Note the extra credit takes about 90s, so I disabled it under pwa/p2js.
Prolog
<lang Prolog> ft_rows(Lz) :-
lazy_list(ft_row, [], Lz).
ft_row([], R1, R1) :- R1 = [1]. ft_row(R0, R2, R2) :-
length(R0, P), Jmax is 1 + P, numlist(2, Jmax, Qs), maplist(term(P), Qs, R0, R1), sum_list(R1, S), Bk is 1 - S, % Bk is Bernoulli number R2 = [Bk | R1].
term(P, Q, R, S) :- S is R * (P rdiv Q).
show(N) :-
ft_rows(Rs),
length(Rows, N), prefix(Rows, Rs),
forall(
member(R, Rows),
(format(string(S), "~w", [R]),
re_replace(" rdiv "/g, "/", S, T),
re_replace(","/g, ", ", T, U),
write(U), nl)).
sum(N, K, S) :- % sum I=1,N (I ** K)
ft_rows(Rows), drop(K, Rows, [Coefs|_]), reverse([0|Coefs], Poly), foldl(horner(N), Poly, 0, S).
horner(N, A, S0, S1) :-
S1 is N*S0 + A.
drop(N, Lz1, Lz2) :-
append(Pfx, Lz2, Lz1), length(Pfx, N), !.
</lang>
- Output:
?- show(10). [1] [1/2, 1/2] [1/6, 1/2, 1/3] [0, 1/4, 1/2, 1/4] [-1/30, 0, 1/3, 1/2, 1/5] [0, -1/12, 0, 5/12, 1/2, 1/6] [1/42, 0, -1/6, 0, 1/2, 1/2, 1/7] [0, 1/12, 0, -7/24, 0, 7/12, 1/2, 1/8] [-1/30, 0, 2/9, 0, -7/15, 0, 2/3, 1/2, 1/9] [0, -3/20, 0, 1/2, 0, -7/10, 0, 3/4, 1/2, 1/10] true. ?- sum(1000, 17, S). S = 56056972216555580111030077961944183400198333273050000.
Python
<lang python>Faulhaber's triangle
from itertools import accumulate, chain, count, islice from fractions import Fraction
- faulhaberTriangle :: Int -> Fraction
def faulhaberTriangle(m):
List of rows of Faulhaber fractions.
def go(rs, n):
def f(x, y):
return Fraction(n, x) * y
xs = list(map(f, islice(count(2), m), rs))
return [Fraction(1 - sum(xs), 1)] + xs
return list(accumulate(
[[]] + list(islice(count(0), 1 + m)),
go
))[1:]
- faulhaberSum :: Integer -> Integer -> Integer
def faulhaberSum(p, n):
Sum of the p-th powers of the first n
positive integers.
def go(x, y):
return y * (n ** x)
return sum(
map(go, count(1), faulhaberTriangle(p)[-1])
)
- ------------------------- TEST -------------------------
def main():
Tests
fs = faulhaberTriangle(9)
print(
fTable(__doc__ + ':\n')(str)(
compose(concat)(
fmap(showRatio(3)(3))
)
)(
index(fs)
)(range(0, len(fs)))
)
print()
print(
faulhaberSum(17, 1000)
)
- ----------------------- DISPLAY ------------------------
- fTable :: String -> (a -> String) ->
- (b -> String) -> (a -> b) -> [a] -> String
def fTable(s):
Heading -> x display function ->
fx display function -> f -> xs -> tabular string.
def gox(xShow):
def gofx(fxShow):
def gof(f):
def goxs(xs):
ys = [xShow(x) for x in xs]
w = max(map(len, ys))
def arrowed(x, y):
return y.rjust(w, ' ') + ' -> ' + (
fxShow(f(x))
)
return s + '\n' + '\n'.join(
map(arrowed, xs, ys)
)
return goxs
return gof
return gofx
return gox
- ----------------------- GENERIC ------------------------
- compose (<<<) :: (b -> c) -> (a -> b) -> a -> c
def compose(g):
Right to left function composition. return lambda f: lambda x: g(f(x))
- concat :: a -> [a]
- concat :: [String] -> String
def concat(xs):
The concatenation of all the elements
in a list or iterable.
def f(ys):
zs = list(chain(*ys))
return .join(zs) if isinstance(ys[0], str) else zs
return (
f(xs) if isinstance(xs, list) else (
chain.from_iterable(xs)
)
) if xs else []
- fmap :: (a -> b) -> [a] -> [b]
def fmap(f):
fmap over a list.
f lifted to a function over a list.
def go(xs):
return list(map(f, xs))
return go
- index (!!) :: [a] -> Int -> a
def index(xs):
Item at given (zero-based) index.
return lambda n: None if 0 > n else (
xs[n] if (
hasattr(xs, "__getitem__")
) else next(islice(xs, n, None))
)
- showRatio :: Int -> Int -> Ratio -> String
def showRatio(m):
Left and right aligned string
representation of the ratio r.
def go(n):
def f(r):
d = r.denominator
return str(r.numerator).rjust(m, ' ') + (
('/' + str(d).ljust(n, ' ')) if 1 != d else (
' ' * (1 + n)
)
)
return f
return go
- MAIN ---
if __name__ == '__main__':
main()</lang>
- Output:
Faulhaber's triangle: 0 -> 1 1 -> 1/2 1/2 2 -> 1/6 1/2 1/3 3 -> 0 1/4 1/2 1/4 4 -> -1/30 0 1/3 1/2 1/5 5 -> 0 -1/12 0 5/12 1/2 1/6 6 -> 1/42 0 -1/6 0 1/2 1/2 1/7 7 -> 0 1/12 0 -7/24 0 7/12 1/2 1/8 8 -> -1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9 9 -> 0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10 56056972216555580111030077961944183400198333273050000
Racket
<lang racket>#lang racket (require math/number-theory)
(define (second-bernoulli-number n)
(if (= n 1) 1/2 (bernoulli-number n)))
(define (faulhaber-row:formulaic p)
(let ((p+1 (+ p 1)))
(reverse
(for/list ((j (in-range p+1)))
(* (/ p+1) (second-bernoulli-number j) (binomial p+1 j))))))
(define (sum-k^p:formulaic p n)
(for/sum ((f (faulhaber-row:formulaic p)) (i (in-naturals 1))) (* f (expt n i))))
(module+ main
(map faulhaber-row:formulaic (range 10)) (sum-k^p:formulaic 17 1000))
(module+ test
(require rackunit)
(check-equal? (sum-k^p:formulaic 17 1000)
(for/sum ((k (in-range 1 (add1 1000)))) (expt k 17))))</lang>
- Output:
'((1) (1/2 1/2) (1/6 1/2 1/3) (0 1/4 1/2 1/4) (-1/30 0 1/3 1/2 1/5) (0 -1/12 0 5/12 1/2 1/6) (1/42 0 -1/6 0 1/2 1/2 1/7) (0 1/12 0 -7/24 0 7/12 1/2 1/8) (-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9) (0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10)) 56056972216555580111030077961944183400198333273050000
Raku
(formerly Perl 6)
<lang perl6># Helper subs
sub infix:<reduce> (\prev, \this) { this.key => this.key * (this.value - prev.value) }
sub next-bernoulli ( (:key($pm), :value(@pa)) ) {
$pm + 1 => [ map *.value, [\reduce] ($pm + 2 ... 1) Z=> 1 / ($pm + 2), |@pa ]
}
constant bernoulli = (0 => [1.FatRat], &next-bernoulli ... *).map: { .value[*-1] };
sub binomial (Int $n, Int $p) { combinations($n, $p).elems }
sub asRat (FatRat $r) { $r ?? $r.denominator == 1 ?? $r.numerator !! $r.nude.join('/') !! 0 }
- The task
sub faulhaber_triangle ($p) { map { binomial($p + 1, $_) * bernoulli[$_] / ($p + 1) }, ($p ... 0) }
- First 10 rows of Faulhaber's triangle:
say faulhaber_triangle($_)».&asRat.fmt('%5s') for ^10; say ;
- Extra credit:
my $p = 17; my $n = 1000; say sum faulhaber_triangle($p).kv.map: { $^value * $n**($^key + 1) }</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
56056972216555580111030077961944183400198333273050000
REXX
<lang rexx>Numeric Digits 100 Do r=0 To 20
ra=r-1
If r=0 Then
f.r.1=1
Else Do
rsum=0
Do c=2 To r+1
ca=c-1
f.r.c=fdivide(fmultiply(f.ra.ca,r),c)
rsum=fsum(rsum,f.r.c)
End
f.r.1=fsubtract(1,rsum)
End
End
Do r=0 To 9
ol= Do c=1 To r+1 ol=ol right(f.r.c,5) End Say ol End
Say x=0 Do c=1 To 18
x=fsum(x,fmultiply(f.17.c,(1000**c))) End
Say k(x) s=0 Do k=1 To 1000
s=s+k**17 End
Say s Exit
fmultiply: Procedure Parse Arg a,b Parse Var a ad '/' an Parse Var b bd '/' bn If an= Then an=1 If bn= Then bn=1 res=(abs(ad)*abs(bd))'/'||(an*bn) Return s(ad,bd)k(res)
fdivide: Procedure Parse Arg a,b Parse Var a ad '/' an Parse Var b bd '/' bn If an= Then an=1 If bn= Then bn=1 res=s(ad,bd)(abs(ad)*bn)'/'||(an*abs(bd)) Return k(res)
fsum: Procedure Parse Arg a,b Parse Var a ad '/' an Parse Var b bd '/' bn If an= Then an=1 If bn= Then bn=1 n=an*bn d=ad*bn+bd*an res=d'/'n Return k(res)
fsubtract: Procedure Parse Arg a,b Parse Var a ad '/' an Parse Var b bd '/' bn If an= Then an=1 If bn= Then bn=1 n=an*bn d=ad*bn-bd*an res=d'/'n Return k(res)
s: Procedure Parse Arg ad,bd s=sign(ad)*sign(bd) If s<0 Then Return '-'
Else Return
k: Procedure Parse Arg a Parse Var a ad '/' an Select
When ad=0 Then Return 0 When an=1 Then Return ad Otherwise Do g=gcd(ad,an) ad=ad/g an=an/g Return ad'/'an End End
gcd: procedure Parse Arg a,b if b = 0 then return abs(a) return gcd(b,a//b)</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
56056972216555580111030077961944183400198333273050000
56056972216555580111030077961944183400198333273050000
Ruby
<lang ruby>class Frac
attr_accessor:num attr_accessor:denom
def initialize(n,d)
if d == 0 then
raise ArgumentError.new('d cannot be zero')
end
nn = n
dd = d
if nn == 0 then
dd = 1
elsif dd < 0 then
nn = -nn
dd = -dd
end
g = nn.abs.gcd(dd.abs)
if g > 1 then
nn = nn / g
dd = dd / g
end
@num = nn
@denom = dd
end
def to_s
if self.denom == 1 then
return self.num.to_s
else
return "%d/%d" % [self.num, self.denom]
end
end
def -@
return Frac.new(-self.num, self.denom)
end
def +(rhs)
return Frac.new(self.num * rhs.denom + self.denom * rhs.num, rhs.denom * self.denom)
end
def -(rhs)
return Frac.new(self.num * rhs.denom - self.denom * rhs.num, rhs.denom * self.denom)
end
def *(rhs)
return Frac.new(self.num * rhs.num, rhs.denom * self.denom)
end
end
FRAC_ZERO = Frac.new(0, 1) FRAC_ONE = Frac.new(1, 1)
def bernoulli(n)
if n < 0 then
raise ArgumentError.new('n cannot be negative')
end
a = Array.new(n + 1) a[0] = FRAC_ZERO
for m in 0 .. n do
a[m] = Frac.new(1, m + 1)
m.downto(1) do |j|
a[j - 1] = (a[j - 1] - a[j]) * Frac.new(j, 1)
end
end
if n != 1 then
return a[0]
end
return -a[0]
end
def binomial(n, k)
if n < 0 then
raise ArgumentError.new('n cannot be negative')
end
if k < 0 then
raise ArgumentError.new('k cannot be negative')
end
if n < k then
raise ArgumentError.new('n cannot be less than k')
end
if n == 0 or k == 0 then
return 1
end
num = 1
for i in k + 1 .. n do
num = num * i
end
den = 1
for i in 2 .. n - k do
den = den * i
end
return num / den
end
def faulhaberTriangle(p)
coeffs = Array.new(p + 1)
coeffs[0] = FRAC_ZERO
q = Frac.new(1, p + 1)
sign = -1
for j in 0 .. p do
sign = -sign
coeffs[p - j] = q * Frac.new(sign, 1) * Frac.new(binomial(p + 1, j), 1) * bernoulli(j)
end
return coeffs
end
def main
for i in 0 .. 9 do
coeffs = faulhaberTriangle(i)
coeffs.each do |coeff|
print "%5s " % [coeff]
end
puts
end
end
main()</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
Scala
<lang scala>import java.math.MathContext
import scala.collection.mutable
abstract class Frac extends Comparable[Frac] {
val num: BigInt val denom: BigInt
def unary_-(): Frac = {
Frac(-num, denom)
}
def +(rhs: Frac): Frac = {
Frac(
num * rhs.denom + rhs.num * denom,
denom * rhs.denom
)
}
def -(rhs: Frac): Frac = {
Frac(
num * rhs.denom - rhs.num * denom,
denom * rhs.denom
)
}
def *(rhs: Frac): Frac = {
Frac(num * rhs.num, denom * rhs.denom)
}
override def compareTo(rhs: Frac): Int = {
val ln = num * rhs.denom
val rn = rhs.num * denom
ln.compare(rn)
}
def canEqual(other: Any): Boolean = other.isInstanceOf[Frac]
override def equals(other: Any): Boolean = other match {
case that: Frac =>
(that canEqual this) &&
num == that.num &&
denom == that.denom
case _ => false
}
override def hashCode(): Int = {
val state = Seq(num, denom)
state.map(_.hashCode()).foldLeft(0)((a, b) => 31 * a + b)
}
override def toString: String = {
if (denom == 1) {
return s"$num"
}
s"$num/$denom"
}
}
object Frac {
val ZERO: Frac = Frac(0) val ONE: Frac = Frac(1)
def apply(n: BigInt): Frac = new Frac {
val num: BigInt = n
val denom: BigInt = 1
}
def apply(n: BigInt, d: BigInt): Frac = {
if (d == 0) {
throw new IllegalArgumentException("Parameter d may not be zero.")
}
var nn = n var dd = d
if (nn == 0) {
dd = 1
} else if (dd < 0) {
nn = -nn
dd = -dd
}
val g = nn.gcd(dd)
if (g > 0) {
nn /= g
dd /= g
}
new Frac {
val num: BigInt = nn
val denom: BigInt = dd
}
}
}
object Faulhaber {
def bernoulli(n: Int): Frac = {
if (n < 0) {
throw new IllegalArgumentException("n may not be negative or zero")
}
val a = Array.fill(n + 1)(Frac.ZERO)
for (m <- 0 to n) {
a(m) = Frac(1, m + 1)
for (j <- m to 1 by -1) {
a(j - 1) = (a(j - 1) - a(j)) * Frac(j)
}
}
// returns 'first' Bernoulli number
if (n != 1) {
return a(0)
}
-a(0)
}
def binomial(n: Int, k: Int): Int = {
if (n < 0 || k < 0 || n < k) {
throw new IllegalArgumentException()
}
if (n == 0 || k == 0) {
return 1
}
val num = (k + 1 to n).product
val den = (2 to n - k).product
num / den
}
def faulhaberTriangle(p: Int): List[Frac] = {
val coeffs = mutable.MutableList.fill(p + 1)(Frac.ZERO)
val q = Frac(1, p + 1)
var sign = -1
for (j <- 0 to p) {
sign *= -1
coeffs(p - j) = q * Frac(sign) * Frac(binomial(p + 1, j)) * bernoulli(j)
}
coeffs.toList
}
def main(args: Array[String]): Unit = {
for (i <- 0 to 9) {
val coeffs = faulhaberTriangle(i)
for (coeff <- coeffs) {
print("%5s ".format(coeff))
}
println()
}
println()
}
}</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
Scheme
<lang scheme>; Return the first row-count rows of Faulhaber's Triangle as a vector of vectors. (define faulhabers-triangle
(lambda (row-count)
; Calculate and store the value of the first column of a row.
; The value is one minus the sum of all the rest of the columns.
(define calc-store-first!
(lambda (row)
(vector-set! row 0
(do ((col-inx 1 (1+ col-inx))
(col-sum 0 (+ col-sum (vector-ref row col-inx))))
((>= col-inx (vector-length row)) (- 1 col-sum))))))
; Generate the Faulhaber's Triangle one row at a time.
; The element at row i >= 0, column j >= 1 (both 0-based) is the product
; of the element at i - 1, j - 1 and the fraction ( i / ( j + 1 ) ).
; The element at column 0 is one minus the sum of all the rest of the columns.
(let ((tri (make-vector row-count)))
(do ((row-inx 0 (1+ row-inx)))
((>= row-inx row-count) tri)
(let ((row (make-vector (1+ row-inx))))
(vector-set! tri row-inx row)
(do ((col-inx 1 (1+ col-inx)))
((>= col-inx (vector-length row)))
(vector-set! row col-inx
(* (vector-ref (vector-ref tri (1- row-inx)) (1- col-inx))
(/ row-inx (1+ col-inx)))))
(calc-store-first! row))))))
- Convert elements of a vector to a string for display.
(define vector->string
(lambda (vec)
(do ((inx 0 (1+ inx))
(str "" (string-append str (format "~7@a" (vector-ref vec inx)))))
((>= inx (vector-length vec)) str))))
- Display a Faulhaber's Triangle.
(define faulhabers-triangle-display
(lambda (tri)
(do ((inx 0 (1+ inx)))
((>= inx (vector-length tri)))
(printf "~a~%" (vector->string (vector-ref tri inx))))))
- Task..
(let ((row-count 10))
(printf "The first ~a rows of Faulhaber's Triangle..~%" row-count) (faulhabers-triangle-display (faulhabers-triangle row-count)))
(newline) (let ((power 17)
(sum-to 1000))
(printf "Sum over k=1..~a of k^~a using Faulhaber's Triangle..~%" sum-to power)
(let* ((tri (faulhabers-triangle (1+ power)))
(coefs (vector-ref tri power)))
(printf "~a~%" (do ((inx 0 (1+ inx))
(term-expt sum-to (* term-expt sum-to))
(sum 0 (+ sum (* (vector-ref coefs inx) term-expt))))
((>= inx (vector-length coefs)) sum)))))</lang>
- Output:
The first 10 rows of Faulhaber's Triangle..
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
Sum over k=1..1000 of k^17 using Faulhaber's Triangle..
56056972216555580111030077961944183400198333273050000
Sidef
<lang ruby>func faulhaber_triangle(p) {
{ binomial(p, _) * bernoulli(_) / p }.map(p ^.. 0)
}
{ |p|
say faulhaber_triangle(p).map{ '%6s' % .as_rat }.join
} << 1..10
const p = 17 const n = 1000
say say faulhaber_triangle(p+1).map_kv {|k,v| v * n**(k+1) }.sum</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
56056972216555580111030077961944183400198333273050000
Alternative solution: <lang ruby>func find_poly_degree(a) {
var c = 0
loop {
++c
a = a.map_cons(2, {|n,k| n-k })
return 0 if a.is_empty
return c if a.all { .is_zero }
}
}
func faulhaber_triangle(n) {
var a = (0..(n+2) -> accumulate { _**n })
var c = find_poly_degree(a)
var A = c.of {|n|
c.of {|k| n**k }
}
A.msolve(a).slice(1)
}
10.times { say faulhaber_triangle(_) }</lang>
- Output:
[1] [1/2, 1/2] [1/6, 1/2, 1/3] [0, 1/4, 1/2, 1/4] [-1/30, 0, 1/3, 1/2, 1/5] [0, -1/12, 0, 5/12, 1/2, 1/6] [1/42, 0, -1/6, 0, 1/2, 1/2, 1/7] [0, 1/12, 0, -7/24, 0, 7/12, 1/2, 1/8] [-1/30, 0, 2/9, 0, -7/15, 0, 2/3, 1/2, 1/9] [0, -3/20, 0, 1/2, 0, -7/10, 0, 3/4, 1/2, 1/10]
Visual Basic .NET
<lang vbnet>Module Module1
Class Frac
Private ReadOnly num As Long
Private ReadOnly denom As Long
Public Shared ReadOnly ZERO = New Frac(0, 1)
Public Shared ReadOnly ONE = New Frac(1, 1)
Public Sub New(n As Long, d As Long)
If d = 0 Then
Throw New ArgumentException("d must not be zero")
End If
Dim nn = n
Dim dd = d
If nn = 0 Then
dd = 1
ElseIf dd < 0 Then
nn = -nn
dd = -dd
End If
Dim g = Math.Abs(Gcd(nn, dd))
If g > 1 Then
nn /= g
dd /= g
End If
num = nn
denom = dd
End Sub
Private Shared Function Gcd(a As Long, b As Long) As Long
If b = 0 Then
Return a
Else
Return Gcd(b, a Mod b)
End If
End Function
Public Shared Operator -(self As Frac) As Frac
Return New Frac(-self.num, self.denom)
End Operator
Public Shared Operator +(lhs As Frac, rhs As Frac) As Frac
Return New Frac(lhs.num * rhs.denom + lhs.denom * rhs.num, rhs.denom * lhs.denom)
End Operator
Public Shared Operator -(lhs As Frac, rhs As Frac) As Frac
Return lhs + -rhs
End Operator
Public Shared Operator *(lhs As Frac, rhs As Frac) As Frac
Return New Frac(lhs.num * rhs.num, lhs.denom * rhs.denom)
End Operator
Public Shared Operator <(lhs As Frac, rhs As Frac) As Boolean
Dim x = lhs.num / lhs.denom
Dim y = rhs.num / rhs.denom
Return x < y
End Operator
Public Shared Operator >(lhs As Frac, rhs As Frac) As Boolean
Dim x = lhs.num / lhs.denom
Dim y = rhs.num / rhs.denom
Return x > y
End Operator
Public Shared Operator =(lhs As Frac, rhs As Frac) As Boolean
Return lhs.num = rhs.num AndAlso lhs.denom = rhs.denom
End Operator
Public Shared Operator <>(lhs As Frac, rhs As Frac) As Boolean
Return lhs.num <> rhs.num OrElse lhs.denom <> rhs.denom
End Operator
Public Overrides Function ToString() As String
If denom = 1 Then
Return num.ToString
Else
Return String.Format("{0}/{1}", num, denom)
End If
End Function
Public Overrides Function Equals(obj As Object) As Boolean
Dim frac = CType(obj, Frac)
Return Not IsNothing(frac) AndAlso num = frac.num AndAlso denom = frac.denom
End Function
End Class
Function Bernoulli(n As Integer) As Frac
If n < 0 Then
Throw New ArgumentException("n may not be negative or zero")
End If
Dim a(n + 1) As Frac
For m = 0 To n
a(m) = New Frac(1, m + 1)
For j = m To 1 Step -1
a(j - 1) = (a(j - 1) - a(j)) * New Frac(j, 1)
Next
Next
' returns 'first' Bernoulli number
If n <> 1 Then
Return a(0)
Else
Return -a(0)
End If
End Function
Function Binomial(n As Integer, k As Integer) As Integer
If n < 0 OrElse k < 0 OrElse n < k Then
Throw New ArgumentException()
End If
If n = 0 OrElse k = 0 Then
Return 1
End If
Dim num = 1
For i = k + 1 To n
num *= i
Next
Dim denom = 1
For i = 2 To n - k
denom *= i
Next
Return num \ denom
End Function
Function FaulhaberTriangle(p As Integer) As Frac()
Dim coeffs(p + 1) As Frac
For i = 1 To p + 1
coeffs(i - 1) = Frac.ZERO
Next
Dim q As New Frac(1, p + 1)
Dim sign = -1
For j = 0 To p
sign *= -1
coeffs(p - j) = q * New Frac(sign, 1) * New Frac(Binomial(p + 1, j), 1) * Bernoulli(j)
Next
Return coeffs
End Function
Sub Main()
For i = 1 To 10
Dim coeffs = FaulhaberTriangle(i - 1)
For Each coeff In coeffs
Console.Write("{0,5} ", coeff)
Next
Console.WriteLine()
Next
End Sub
End Module</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
Vlang
<lang vlang>import math.fractions import math.big
fn bernoulli(n int) fractions.Fraction {
mut a := []fractions.Fraction{len: n+1}
for m,_ in a {
a[m] = fractions.fraction(1, i64(m+1))
for j := m; j >= 1; j-- {
mut d := a[j-1]
d = fractions.fraction(i64(j),i64(1)) * (d-a[j])
a[j-1]=d
}
}
// return the 'first' Bernoulli number
if n != 1 {
return a[0]
}
a[0] = a[0].negate()
return a[0]
}
fn binomial(n int, k int) i64 {
if n <= 0 || k <= 0 || n < k {
return 1
}
mut num, mut den := i64(1), i64(1)
for i := k + 1; i <= n; i++ {
num *= i64(i)
}
for i := 2; i <= n-k; i++ {
den *= i64(i)
}
return num / den
}
fn faulhaber_triangle(p int) []fractions.Fraction {
mut coeffs := []fractions.Fraction{len: p+1}
q := fractions.fraction(1, i64(p)+1)
mut t := fractions.fraction(1,1)
mut u := fractions.fraction(1,1)
mut sign := -1
for j,_ in coeffs {
sign *= -1
mut d := coeffs[p-j]
t=fractions.fraction(i64(sign),1)
u = fractions.fraction(binomial(p+1, j),1)
d=q*t
d*=u
d*=bernoulli(j)
coeffs[p-j]=d
}
return coeffs
}
fn main() {
for i in 0..10 {
coeffs := faulhaber_triangle(i)
for coeff in coeffs {
print("${coeff:5} ")
}
println()
}
println()
}</lang>
- Output:
1/1 1/2 1/2 1/6 1/2 1/3 0/1 1/4 1/2 1/4 -1/30 0/1 1/3 1/2 1/5 0/1 -1/12 0/1 5/12 1/2 1/6 1/42 0/1 -1/6 0/1 1/2 1/2 1/7 0/1 1/12 0/1 -7/24 0/1 7/12 1/2 1/8 -1/30 0/1 2/9 0/1 -7/15 0/1 2/3 1/2 1/9 0/1 -3/20 0/1 1/2 0/1 -7/10 0/1 3/4 1/2 1/10
Wren
<lang ecmascript>import "/fmt" for Fmt import "/math" for Int import "/big" for BigRat
var bernoulli = Fn.new { |n|
if (n < 0) Fiber.abort("Argument must be non-negative")
var a = List.filled(n+1, null)
for (m in 0..n) {
a[m] = BigRat.new(1, m+1)
var j = m
while (j >= 1) {
a[j-1] = (a[j-1] - a[j]) * BigRat.new(j, 1)
j = j - 1
}
}
return (n != 1) ? a[0] : -a[0] // 'first' Bernoulli number
}
var binomial = Fn.new { |n, k|
if (n < 0 || k < 0) Fiber.abort("Arguments must be non-negative integers")
if (n < k) Fiber.abort("The second argument cannot be more than the first.")
if (n == k) return 1
var prod = 1
var i = n - k + 1
while (i <= n) {
prod = prod * i
i = i + 1
}
return prod / Int.factorial(k)
}
var faulhaberTriangle = Fn.new { |p|
var coeffs = List.filled(p+1, null)
var q = BigRat.new(1, p+1)
var sign = -1
for (j in 0..p) {
sign = sign * -1
var b = BigRat.new(binomial.call(p+1, j), 1)
coeffs[p-j] = q * BigRat.new(sign, 1) * b * bernoulli.call(j)
}
return coeffs
}
BigRat.showAsInt = true for (i in 0..9) {
var coeffs = faulhaberTriangle.call(i)
for (coeff in coeffs) Fmt.write("$5s ", coeff)
System.print()
} System.print()
// get coeffs for (k + 1)th row
var k = 17 var cc = faulhaberTriangle.call(k) var n = BigRat.new(1000, 1) var np = BigRat.one var sum = BigRat.zero for (c in cc) {
np = np * n sum = sum + np*c
} System.print(sum)</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
56056972216555580111030077961944183400198333273050000
zkl
GNU Multiple Precision Arithmetic Library
Uses the code from Faulhaber's formula#zkl <lang zkl>foreach p in (10){
faulhaberFormula(p).apply("%7s".fmt).concat().println();
}
// each term of faulhaberFormula is BigInt/BigInt [1..].zipWith(fcn(n,rat){ rat*BN(1000).pow(n) }, faulhaberFormula(17)) .walk() // -->(0, -3617/60 * 1000^2, 0, 595/3 * 1000^4 ...) .reduce('+) // rat + rat + ... .println();</lang>
- Output:
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
56056972216555580111030077961944183400198333273050000