Faulhaber's triangle: Difference between revisions

Content added Content deleted

Inline

Latest revision as of 11:46, 4 December 2023

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:

\sum _{k=1}^{n}k^{p}={1 \over p+1}\sum _{j=0}^{p}{p+1 \choose j}B_{j}n^{p+1-j}

where $B_{n}$ 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:

$\sum _{k=1}^{n}k^{2}={1 \over 6}n+{1 \over 2}n^{2}+{1 \over 3}n^{3}$

Task

show the first 10 rows of Faulhaber's triangle.
using the 18th row of Faulhaber's triangle, compute the sum: $\sum _{k=1}^{1000}k^{17}$ (extra credit).

See also

ALGOL 68

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32

Using code from the Algol 68 samples for the Arithmetic/Rational and Bernoulli numbers tasks and the Algol W sample for the Evaluate binomial coefficients task.
Note that in the Bernoulli numbers task, the Algol 68 sample returns -1/2 for B(1) - this is modified here so B(1) is 1/2.
Assumes LONG LONG INT is long enough to calculate the 17th power sum, the default precision of LONG LONG INT in ALGOL 68G is large enough.

BEGIN # show some rows of Faulhaber's triangle                               #

    # utility operators                                                      #
    OP   LENGTH = ( STRING a )INT: ( UPB a - LWB a ) + 1;
    PRIO PAD    = 9;
    OP   PAD    = ( INT width, STRING v )STRING: # left blank pad v to width #
         IF LENGTH v >= width THEN v ELSE ( " " * ( width - LENGTH v ) ) + v FI;

    MODE INTEGER   = LONG LONG INT; # mode for FRAC numberator & denominator #
    OP   TOINTEGER = ( INT n )INTEGER: n;      # force widening n to INTEGER #

    # Code from the Arithmetic/Rational task                                 #
  
    MODE FRAC = STRUCT( INTEGER num #erator#,  den #ominator#);

    PROC gcd = (INTEGER a, b) INTEGER: # greatest common divisor #
       (a = 0 | b |: b = 0 | a |: ABS a > ABS b  | gcd(b, a MOD b) | gcd(a, b MOD a));
 
    PROC lcm = (INTEGER a, b)INTEGER: # least common multiple #
       a OVER gcd(a, b) * b;
 
    PRIO // = 9; # higher then the ** operator #
    OP // = (INTEGER num, den)FRAC: ( # initialise and normalise #
       INTEGER common = gcd(num, den);
       IF den < 0 THEN
         ( -num OVER common, -den OVER common)
       ELSE
         ( num OVER common, den OVER common)
       FI
     );

    OP + = (FRAC a, b)FRAC: (
       INTEGER common = lcm(den OF a, den OF b);
       FRAC result := ( common OVER den OF a * num OF a + common OVER den OF b * num OF b, common );
       num OF result//den OF result
    );
 
    OP - = (FRAC a, b)FRAC: a + -b,
       * = (FRAC a, b)FRAC: (
           INTEGER num = num OF a * num OF b,
           den = den OF a * den OF b;
           INTEGER common = gcd(num, den);
           (num OVER common) // (den OVER common)
         );
 
    OP - = (FRAC frac)FRAC: (-num OF frac, den OF frac);

    # end code from the Arithmetic/Rational task                             #

    # alternative // operator for standard size INT values                   #
    OP // = (INT num, den)FRAC: TOINTEGER num // TOINTEGER den;
    # returns a * b                                                          #
    OP *   = ( INT     a, FRAC b )FRAC: ( num OF b * a ) // den OF b;
    OP *   = ( INTEGER a, FRAC b )FRAC: ( num OF b * a ) // den OF b;
    # sets a to a + b and returns a                                          #
    OP +:= = ( REF FRAC a, FRAC b )FRAC: a := a + b;
    # sets a to - a and returns a                                            #
    OP -=: = ( REF FRAC a )FRAC: BEGIN num OF a := - num OF a; a END;

    # returns the nth Bernoulli number, n must be >= 0                       #
    OP   BERNOULLI = ( INT n )FRAC:
         IF n < 0
         THEN # n is out of range # 0 // 1
         ELSE # n is valid        #
            [ 0 : n ]FRAC a;
            FOR m FROM 0 TO n DO 
                a[ m ] := 1 // ( m + 1 );
                FOR j FROM m BY -1 TO 1 DO
                    a[ j - 1 ] := j * ( a[ j - 1 ] - a[ j ] )
                OD
            OD;
            IF n = 1 THEN - a[ 0 ] ELSE a[ 0 ] FI
         FI # BERNOULLI # ;

    # returns n! / k!                                                        #
    PROC factorial over factorial = ( INT n, k )INTEGER:
        IF     k > n THEN 0
        ELIF   k = n THEN 1
        ELSE # k < n #
            INTEGER f := 1;
            FOR i FROM k + 1 TO n DO f *:= i OD;
            f
        FI # factorial over Factorial # ;

    # returns n!                                                             #
    PROC factorial = ( INT n )INTEGER:
         BEGIN
            INTEGER f := 1;
            FOR i FROM 2 TO n DO f *:= i OD;
            f
         END # factorial # ;

    # returns the binomial coefficient of (n k)                              #
    PROC binomial coefficient = ( INT n, k )INTEGER:
         IF n - k > k
         THEN factorial over factorial( n, n - k ) OVER factorial(   k   )
         ELSE factorial over factorial( n,   k   ) OVER factorial( n - k )
         FI # binomial coefficient # ;

    # returns a string representation of a                                   #
    OP   TOSTRING = ( FRAC a )STRING:
         whole( num OF a, 0 ) + IF den OF a = 1 THEN "" ELSE "/" + whole( den OF a, 0 ) FI;

    # returns the pth row of Faulhaber's triangle                            #
    OP   FAULHABER = ( INT p )[]FRAC:
         BEGIN
            FRAC q := -1 // ( p + 1 );
            [ 0 : p ]FRAC coeffs;
            FOR j FROM 0 TO p DO
                coeffs[ p - j ] := binomial coefficient( p + 1, j ) * BERNOULLI j * -=: q
            OD;
            coeffs
         END # faulhaber # ;

    FOR i FROM 0 TO 9 DO               # show the triabngle's first 10 rows #
        []FRAC frow =  FAULHABER i;
        FOR j FROM LWB frow TO UPB frow DO
            print( ( " ", 6 PAD TOSTRING frow[ j ] ) )
        OD;
        print( ( newline ) )
    OD;
    BEGIN # compute the sum of k^17 for k = 1 to 1000 using triangle row 18 #
        []FRAC  frow = FAULHABER 17;
        FRAC    sum := 0 // 1;
        INTEGER kn  := 1;
        FOR j FROM LWB frow TO UPB frow DO
            VOID( sum +:= ( kn *:= 1000 ) * frow[ j ] )
        OD;
        print( ( TOSTRING sum, newline ) )
    END
END

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

C

Translation of: 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;
}

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#

Translation of: Java

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();
            }
        }
    }
}

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++

Translation of: C#

Uses C++ 17

#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;
}

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

Translation of: Kotlin

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;
}

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

// 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

The Task

Faulhaber |> Seq.take 10 |> Seq.iter (printfn "%A")

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

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

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

Library: GMP

' 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

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.

In this page you can see and run the program(s) related to this task and their results. You can also change either the programs or the parameters they are called with, for experimentation, but remember that these programs were created with the main purpose of showing a clear solution of the task, and they generally lack any kind of validation.

Solution

The following function creates the Faulhaber's coefficients up to a given number of rows, according to the paper of of Mohammad Torabi Dashti:

(This is exactly the same as the task Faulhaber's formula)

Excecise 1. To show the first 11 rows (the first is the 0 row) of Faulhaber's triangle:

In order to show the previous result as a triangle:

The following function creates the sum of the p-th powers of the first n positive integers as a (p + 1)th-degree polynomial function of n:

Notes. The -1 index means the last element (-2 is the penultimate element, and so on). So it retrieves the last row of the triangle. |x| is the cardinality (number of elements) of x.

(This is exactly the same as the task Faulhaber's formula)

This function can be used for both symbolic or numeric computation of the polynomial:

Excecise 2. Using the 18th row of Faulhaber's triangle, compute the sum $\sum _{k=1}^{1000}k^{17}$

Verification:

Go

Translation of: Kotlin

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.

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())
}

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

Translation of: Java

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())
    }
}

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

Works with: GHC version 8.6.4

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

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

Translation of: Kotlin

Works with: Java version 8

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());
    }
}

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

Translation of: Haskell

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)

(() => {

    // 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)),
            ]
        );
})();

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.

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

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

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()

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:

// 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())
}

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

Translation of: C

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

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

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]

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

Maxima

faulhaber_fraction(n, k) :=
	    if n = 0 and k = 1 then 1
	    else if k >= 2 and k <= n + 1 then (n/k) * faulhaber_fraction(n-1, k-1)
	    else if k = 1 then 1 - sum(faulhaber_fraction(n, i), i, 2, n+1)
	    else 0$

faulhaber_row(n):=makelist(faulhaber_fraction(n,k),k,1,n+1)$
/* Example */
triangle_faulhaber_first_ten_rows:block(makelist(faulhaber_row(i),i,0,9),table_form(%%));

File:Faulhaber.png

Nim

Library: bignum

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.

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)

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

Library: IntXLib4Pascal

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.

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.

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

Library: ntheory

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);

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

Translation of: C#

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

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), !.

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

Translation of: Haskell

Works with: Python version 3.7

'''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()

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
(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))))

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)

Works with: Rakudo version 2017.05

Translation of: Sidef

# 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) }

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

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)

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

Translation of: D

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()

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

Translation of: Java

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()
  }
}

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

Works with: Chez 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)))))

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

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

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:

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(_) }

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

Translation of: C#

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

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

V (Vlang)

Translation of: Go

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('')
}

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

Translation of: Kotlin

Library: Wren-fmt

Library: Wren-math

Library: Wren-big

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)

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

Library: GMP

GNU Multiple Precision Arithmetic Library

Uses the code from Faulhaber's formula#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();

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