Faulhaber's formula: Difference between revisions

Content added Content deleted

Inline

Latest revision as of 13:06, 17 March 2024

In mathematics, Faulhaber's formula, named after Johann Faulhaber, expresses the sum of the p-th powers of the first n positive integers as a (p + 1)th-degree polynomial function of n, the coefficients involving Bernoulli numbers.

Task

Generate the first 10 closed-form expressions, starting with p = 0.

Related tasks

See also

The Wikipedia entry: Faulhaber's formula.
The Wikipedia entry: Bernoulli numbers.
The Wikipedia entry: binomial coefficients.

C

Translation of: Modula-2

#include <stdbool.h>
#include <stdio.h>
#include <stdlib.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) {
    printf("%d", f.num);
    if (f.denom != 1) {
        printf("/%d", f.denom);
    }
}

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 coeff, q;
    int j, pwr, sign;

    printf("%d : ", p);
    q = makeFrac(1, p + 1);
    sign = -1;
    for (j = 0; j <= p; ++j) {
        sign = -1 * sign;
        coeff = multFrac(multFrac(multFrac(q, makeFrac(sign, 1)), makeFrac(binomial(p + 1, j), 1)), bernoulli(j));
        if (equalFrac(coeff, makeFrac(0, 1))) {
            continue;
        }
        if (j == 0) {
            if (!equalFrac(coeff, makeFrac(1, 1))) {
                if (equalFrac(coeff, makeFrac(-1, 1))) {
                    printf("-");
                } else {
                    printFrac(coeff);
                }
            }
        } else {
            if (equalFrac(coeff, makeFrac(1, 1))) {
                printf(" + ");
            } else if (equalFrac(coeff, makeFrac(-1, 1))) {
                printf(" - ");
            } else if (lessFrac(makeFrac(0, 1), coeff)) {
                printf(" + ");
                printFrac(coeff);
            } else {
                printf(" - ");
                printFrac(negateFrac(coeff));
            }
        }
        pwr = p + 1 - j;
        if (pwr > 1) {
            printf("n^%d", pwr);
        } else {
            printf("n");
        }
    }
    printf("\n");
}

int main() {
    int i;

    for (i = 0; i < 10; ++i) {
        faulhaber(i);
    }

    return 0;
}

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

C#

Translation of: Java

using System;

namespace FaulhabersFormula {
    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 void Faulhaber(int p) {
            Console.Write("{0} : ", p);
            Frac q = new Frac(1, p + 1);
            int sign = -1;
            for (int j = 0; j <= p; j++) {
                sign *= -1;
                Frac coeff = q * new Frac(sign, 1) * new Frac(Binomial(p + 1, j), 1) * Bernoulli(j);
                if (Frac.ZERO == coeff) continue;
                if (j == 0) {
                    if (Frac.ONE != coeff) {
                        if (-Frac.ONE == coeff) {
                            Console.Write("-");
                        }
                        else {
                            Console.Write(coeff);
                        }
                    }
                }
                else {
                    if (Frac.ONE == coeff) {
                        Console.Write(" + ");
                    }
                    else if (-Frac.ONE == coeff) {
                        Console.Write(" - ");
                    }
                    else if (Frac.ZERO < coeff) {
                        Console.Write(" + {0}", coeff);
                    }
                    else {
                        Console.Write(" - {0}", -coeff);
                    }
                }
                int pwr = p + 1 - j;
                if (pwr > 1) {
                    Console.Write("n^{0}", pwr);
                }
                else {
                    Console.Write("n");
                }
            }
            Console.WriteLine();
        }

        static void Main(string[] args) {
            for (int i = 0; i < 10; i++) {
                Faulhaber(i);
            }
        }
    }
}

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

C++

Translation of: D

Uses C++17

#include <iostream>
#include <numeric>
#include <sstream>
#include <vector>

class Frac {
public:
	Frac(long n, long d) {
		if (d == 0) {
			throw new std::runtime_error("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 = abs(std::gcd(nn, dd));
		if (g > 1) {
			nn /= g;
			dd /= g;
		}

		num = nn;
		denom = dd;
	}

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

	bool operator==(const Frac& rhs) const {
		return num == rhs.num && denom == rhs.denom;
	}

	bool operator!=(const Frac& rhs) const {
		return num != rhs.num || denom != rhs.denom;
	}

	bool operator<(const Frac& rhs) const {
		if (denom == rhs.denom) {
			return num < rhs.num;
		}
		return num * rhs.denom < rhs.num * denom;
	}

	friend std::ostream& operator<<(std::ostream&, const Frac&);

	static Frac ZERO() {
		return Frac(0, 1);
	}

	static Frac ONE() {
		return Frac(1, 1);
	}

private:
	long num;
	long 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 new 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]) * Frac(j, 1);
		}
	}

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

void faulhaber(int p) {
	using namespace std;
	cout << p << " : ";

	auto q = Frac(1, p + 1);
	int sign = -1;
	for (int j = 0; j < p + 1; j++) {
		sign *= -1;
		auto coeff = q * Frac(sign, 1) * Frac(binomial(p + 1, j), 1) * bernoulli(j);
		if (coeff == Frac::ZERO()) {
			continue;
		}
		if (j == 0) {
			if (coeff == -Frac::ONE()) {
				cout << "-";
			} else if (coeff != Frac::ONE()) {
				cout << coeff;
			}
		} else {
			if (coeff == Frac::ONE()) {
				cout << " + ";
			} else if (coeff == -Frac::ONE()) {
				cout << " - ";
			} else if (coeff < Frac::ZERO()) {
				cout << " - " << -coeff;
			} else {
				cout << " + " << coeff;
			}
		}
		int pwr = p + 1 - j;
		if (pwr > 1) {
			cout << "n^" << pwr;
		} else {
			cout << "n";
		}
	}
	cout << endl;
}

int main() {
	for (int i = 0; i < 10; i++) {
		faulhaber(i);
	}

	return 0;
}

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

D

Translation of: Kotlin

import std.algorithm : fold;
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;
}

auto faulhaber(int p) {
    write(p, " : ");
    auto q = Frac(1, p+1);
    auto sign = -1;
    foreach (j; 0..p+1) {
        sign *= -1;
        auto coeff = q * Frac(sign, 1) * Frac(binomial(p+1, j), 1) * bernoulli(j);
        if (coeff == Frac.ZERO) continue;
        if (j == 0) {
            if (coeff == -Frac.ONE) {
                write("-");
            } else if (coeff != Frac.ONE) {
                write(coeff);
            }
        } else {
            if (coeff == Frac.ONE) {
                write(" + ");
            } else if (coeff == -Frac.ONE) {
                write(" - ");
            } else if (coeff > Frac.ZERO) {
                write(" + ", coeff);
            } else {
                write(" - ", -coeff);
            }
        }
        auto pwr = p + 1 - j;
        if (pwr > 1) {
            write("n^", pwr);
        } else {
            write("n");
        }
    }
    writeln;
}

void main() {
    foreach (i; 0..10) {
        faulhaber(i);
    }
}

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

EchoLisp

(lib 'math) ;; for bernoulli numbers
(string-delimiter "")

;; returns list of polynomial coefficients
(define (Faulhaber p)
	(cons 0
	(for/list ([k (in-range p -1 -1)])
		(* (Cnp (1+ p) k) (bernoulli k)))))

;; prints formal polynomial
(define (task (pmax 10))
    (for ((p pmax)) 
    (writeln p '→  (/ 1 (1+ p)) '* (poly->string 'n (Faulhaber p)))))
    
;; extra credit - compute sums
(define (Faulcomp n p)
	(printf "Σ(1..%d) n^%d = %d" n p (/  (poly n (Faulhaber p)) (1+ p) )))

Output:

(task)
0     →     1     *     n     
1     →     1/2     *     n^2 + n     
2     →     1/3     *     n^3 + 3/2 n^2 + 1/2 n     
3     →     1/4     *     n^4 + 2 n^3 + n^2     
4     →     1/5     *     n^5 + 5/2 n^4 + 5/3 n^3 -1/6 n     
5     →     1/6     *     n^6 + 3 n^5 + 5/2 n^4 -1/2 n^2     
6     →     1/7     *     n^7 + 7/2 n^6 + 7/2 n^5 -7/6 n^3 + 1/6 n     
7     →     1/8     *     n^8 + 4 n^7 + 14/3 n^6 -7/3 n^4 + 2/3 n^2     
8     →     1/9     *     n^9 + 9/2 n^8 + 6 n^7 -21/5 n^5 + 2 n^3 -3/10 n     
9     →     1/10     *     n^10 + 5 n^9 + 15/2 n^8 -7 n^6 + 5 n^4 -3/2 n^2     

(Faulcomp 100 2)
    Σ(1..100) n^2 = 338350
(Faulcomp 100 1)
    Σ(1..100) n^1 = 5050

(lib 'bigint)
(Faulcomp 100 9)
    Σ(1..100) n^9 = 10507499300049998000

;; check it ...
(for/sum ((n 101)) (expt n 9))
    → 10507499300049998500

Factor

USING: formatting kernel math math.combinatorics math.extras
math.functions regexp sequences ;

: faulhaber ( p -- seq )
    1 + dup recip swap dup <iota>
    [ [ nCk ] [ -1 swap ^ ] [ bernoulli ] tri * * * ] 2with map ;

: (poly>str) ( seq -- str )
    reverse [ 1 + "%un^%d" sprintf ] map-index reverse " + " join ;

: clean-up ( str -- str' )
    R/ n\^1\z/ "n" re-replace            ! Change n^1 to n.
    R/ 1n/ "n" re-replace                ! Change 1n to n.
    R/ \+ -/ "- " re-replace             ! Change + - to - .
    R/ [+-] 0n(\^\d+ )?/ "" re-replace ; ! Remove terms of zero.

: poly>str ( seq -- str ) (poly>str) clean-up ;

10 [ dup faulhaber poly>str "%d: %s\n" printf ] each-integer

Output:

0: n
1: 1/2n^2 + 1/2n
2: 1/3n^3 + 1/2n^2 + 1/6n
3: 1/4n^4 + 1/2n^3 + 1/4n^2 
4: 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5: 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2 
6: 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7: 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2 
8: 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9: 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

FreeBASIC

Translation of: C

Type Fraction
    As Integer num
    As Integer den
End Type

Function Binomial(n As Integer, k As Integer) As Integer
    If n < 0 Or k < 0 Then Print "Arguments must be non-negative integers": Exit Function
    If n < k Then Print "The second argument cannot be more than the first.": Exit Function
    If n = 0 Or k = 0 Then Return 1
    
    Dim As Integer i, num, den
    num = 1
    For i = k + 1 To n
        num *= i
    Next i
    
    den = 1
    For i = 2 To n - k
        den *= i
    Next i
    
    Return num / den
End Function

Function GCD(n As Integer, d As Integer) As Integer
    Return Iif(d = 0, n, GCD(d, n Mod d))
End Function

Function makeFrac(n As Integer, d As Integer) As Fraction
    Dim As Fraction result
    Dim As Integer g
    
    If d = 0 Then
        result.num = 0
        result.den = 0
        Return result
    End If
    
    If n = 0 Then
        d = 1
    Elseif d < 0 Then
        n = -n
        d = -d
    End If
    
    g = Abs(gcd(n, d))
    If g > 1 Then
        n /= g
        d /= g
    End If
    
    result.num = n
    result.den = d
    Return result
End Function

Function negateFrac(f As Fraction) As Fraction
    Return makeFrac(-f.num, f.den)
End Function

Function subFrac(lhs As Fraction, rhs As Fraction) As Fraction
    Return makeFrac(lhs.num * rhs.den - lhs.den * rhs.num, rhs.den * lhs.den)
End Function

Function multFrac(lhs As Fraction, rhs As Fraction) As Fraction
    Return makeFrac(lhs.num * rhs.num, lhs.den * rhs.den)
End Function

Function equalFrac(lhs As Fraction, rhs As Fraction) As Integer
    Return (lhs.num = rhs.num) And (lhs.den = rhs.den)
End Function

Function lessFrac(lhs As Fraction, rhs As Fraction) As Integer
    Return (lhs.num * rhs.den) < (rhs.num * lhs.den)
End Function

Sub printFrac(f As Fraction)
    Print Str(f.num);
    If f.den <> 1 Then Print "/" & f.den;
End Sub

Function Bernoulli(n As Integer) As Fraction
    If n < 0 Then Print "Argument must be non-negative": Exit Function
    Dim As Fraction a(16)
    Dim As Integer j, m
    
    If (n < 0) Then
        a(0).num = 0
        a(0).den = 0
        Return a(0)
    End If
    
    For m = 0 To n
        a(m) = makeFrac(1, m + 1)
        For j = m To 1 Step -1
            a(j - 1) = multFrac(subFrac(a(j - 1), a(j)), makeFrac(j, 1))
        Next j
    Next m
    
    If n <> 1 Then Return a(0)
    
    Return negateFrac(a(0))
End Function

Sub Faulhaber(p As Integer)
    Dim As Fraction coeff, q
    Dim As Integer j, pwr, sign
    
    Print p & " : ";
    q = makeFrac(1, p + 1)
    sign = -1
    For j = 0 To p
        sign = -1 * sign
        coeff = multFrac(multFrac(multFrac(q, makeFrac(sign, 1)), makeFrac(Binomial(p + 1, j), 1)), Bernoulli(j))
        If (equalFrac(coeff, makeFrac(0, 1))) Then Continue For
        If j = 0 Then
            If Not equalFrac(coeff, makeFrac(1, 1)) Then
                If equalFrac(coeff, makeFrac(-1, 1)) Then
                    Print "-";
                Else
                    printFrac(coeff)
                End If
            End If
        Else
            If equalFrac(coeff, makeFrac(1, 1)) Then
                Print " + ";
            Elseif equalFrac(coeff, makeFrac(-1, 1)) Then
                Print " - ";
            Elseif lessFrac(makeFrac(0, 1), coeff) Then
                Print " + ";
                printFrac(coeff)
            Else
                Print " - ";
                printFrac(negateFrac(coeff))
            End If
        End If
        pwr = p + 1 - j
        Print Iif(pwr > 1, "n^" & pwr, "n");
    Next j
    Print
End Sub

For i As Integer = 0 To 9
    Faulhaber(i)
Next i

Sleep

Output:

Same as C entry.

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 as than the task Faulhaber's 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 as than the task Faulhaber's triangle)

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

To generate the first 10 closed-form expressions, starting with p = 0:

GAP

Straightforward implementation using GAP polynomials, and two different formulas: one based on Stirling numbers of the second kind (sum1, see Python implementation below in this page), and the usual Faulhaber formula (sum2). No optimization is made (one could compute Stirling numbers row by row, or the product in sum1 may be kept from one call to the other). Notice the Bernoulli term in the first formula is here only to correct the value of sum1(0), which is off by one because sum1 computes sums from 0 to n.

n := X(Rationals, "n");
sum1 := p -> Sum([0 .. p], k -> Stirling2(p, k) * Product([0 .. k], j -> n + 1 - j) / (k + 1)) + 2 * Bernoulli(2 * p + 1);
sum2 := p -> Sum([0 .. p], j -> (-1)^j * Binomial(p + 1, j) * Bernoulli(j) * n^(p + 1 - j)) / (p + 1);
ForAll([0 .. 20], k -> sum1(k) = sum2(k));

for p in [0 .. 9] do
    Print(sum1(p), "\n");
od;

n
1/2*n^2+1/2*n
1/3*n^3+1/2*n^2+1/6*n
1/4*n^4+1/2*n^3+1/4*n^2
1/5*n^5+1/2*n^4+1/3*n^3-1/30*n
1/6*n^6+1/2*n^5+5/12*n^4-1/12*n^2
1/7*n^7+1/2*n^6+1/2*n^5-1/6*n^3+1/42*n
1/8*n^8+1/2*n^7+7/12*n^6-7/24*n^4+1/12*n^2
1/9*n^9+1/2*n^8+2/3*n^7-7/15*n^5+2/9*n^3-1/30*n
1/10*n^10+1/2*n^9+3/4*n^8-7/10*n^6+1/2*n^4-3/20*n^2

Go

package main

import (
	"fmt"
	"math/big"
)

func bernoulli(z *big.Rat, n int64) *big.Rat {
	if z == nil {
		z = new(big.Rat)
	}
	a := make([]big.Rat, n+1)
	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 z.Set(&a[0])
}

func main() {
	// allocate needed big.Rat's once
	q := new(big.Rat)
	c := new(big.Rat)      // coefficients
	be := new(big.Rat)     // for Bernoulli numbers
	bi := big.NewRat(1, 1) // for binomials

	for p := int64(0); p < 10; p++ {
		fmt.Print(p, " : ")
		q.SetFrac64(1, p+1)
		neg := true
		for j := int64(0); j <= p; j++ {
			neg = !neg
			if neg {
				c.Neg(q)
			} else {
				c.Set(q)
			}
			bi.Num().Binomial(p+1, j)
			bernoulli(be, j)
			c.Mul(c, bi)
			c.Mul(c, be)
			if c.Num().BitLen() == 0 {
				continue
			}
			if j == 0 {
				fmt.Printf(" %4s", c.RatString())
			} else {
				fmt.Printf(" %+2d/%-2d", c.Num(), c.Denom())
			}
			fmt.Print("×n")
			if exp := p + 1 - j; exp > 1 {
				fmt.Printf("^%-2d", exp)
			}
		}
		fmt.Println()
	}
}

Output:

0 :     1×n
1 :   1/2×n^2  -1/2 ×n
2 :   1/3×n^3  -1/2 ×n^2  +1/6 ×n
3 :   1/4×n^4  -1/2 ×n^3  +1/4 ×n^2 
4 :   1/5×n^5  -1/2 ×n^4  +1/3 ×n^3  -1/30×n
5 :   1/6×n^6  -1/2 ×n^5  +5/12×n^4  -1/12×n^2 
6 :   1/7×n^7  -1/2 ×n^6  +1/2 ×n^5  -1/6 ×n^3  +1/42×n
7 :   1/8×n^8  -1/2 ×n^7  +7/12×n^6  -7/24×n^4  +1/12×n^2 
8 :   1/9×n^9  -1/2 ×n^8  +2/3 ×n^7  -7/15×n^5  +2/9 ×n^3  -1/30×n
9 :  1/10×n^10 -1/2 ×n^9  +3/4 ×n^8  -7/10×n^6  +1/2 ×n^4  -3/20×n^2

Groovy

Translation of: Java

import java.util.stream.IntStream

class FaulhabersFormula {
    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 static final Frac ONE = new Frac(1, 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)
        }

        private double toDouble() {
            return (double) num / denom
        }
    }

    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 int binomial(int n, int k) {
        if (n < 0 || k < 0 || n < k) throw new IllegalArgumentException()
        if (n == 0 || k == 0) return 1
        int num = IntStream.rangeClosed(k + 1, n).reduce(1, { a, b -> a * b })
        int den = IntStream.rangeClosed(2, n - k).reduce(1, { acc, i -> acc * i })
        return num / den
    }

    private static void faulhaber(int p) {
        print("$p : ")
        Frac q = new Frac(1, p + 1)
        int sign = -1
        for (int j = 0; j <= p; ++j) {
            sign *= -1
            Frac coeff = q * new Frac(sign, 1) * new Frac(binomial(p + 1, j), 1) * bernoulli(j)
            if (Frac.ZERO == coeff) continue
            if (j == 0) {
                if (Frac.ONE != coeff) {
                    if (-Frac.ONE == coeff) {
                        print("-")
                    } else {
                        print(coeff)
                    }
                }
            } else {
                if (Frac.ONE == coeff) {
                    print(" + ")
                } else if (-Frac.ONE == coeff) {
                    print(" - ")
                } else if (coeff > Frac.ZERO) {
                    print(" + $coeff")
                } else {
                    print(" - ${-coeff}")
                }
            }
            int pwr = p + 1 - j
            if (pwr > 1) {
                print("n^$pwr")
            } else {
                print("n")
            }
        }
        println()
    }

    static void main(String[] args) {
        for (int i = 0; i <= 9; ++i) {
            faulhaber(i)
        }
    }
}

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

Haskell

Bernouilli polynomials

import Data.Ratio ((%), numerator, denominator)
import Data.List (intercalate, transpose)
import Data.Bifunctor (bimap)
import Data.Char (isSpace)
import Data.Monoid ((<>))
import Data.Bool (bool)
 
------------------------- FAULHABER ------------------------
faulhaber :: [[Rational]]
faulhaber =
  tail $
  scanl
    (\rs n ->
        let xs = zipWith ((*) . (n %)) [2 ..] rs
        in 1 - sum xs : xs)
    []
    [0 ..]
    
polynomials :: [[(String, String)]]
polynomials = fmap ((ratioPower =<<) . reverse . flip zip [1 ..]) faulhaber
    

---------------------------- TEST --------------------------
main :: IO ()
main = (putStrLn . unlines . expressionTable . take 10) polynomials

 
--------------------- EXPRESSION STRINGS -------------------

-- Rows of (Power string, Ratio string) tuples -> Printable lines
expressionTable :: [[(String, String)]] -> [String]
expressionTable ps =
  let cols = transpose (fullTable ps)
  in expressionRow <$>
     zip
       [0 ..]
       (transpose $
        zipWith
          (\(lw, rw) col ->
              fmap (bimap (justifyLeft lw ' ') (justifyLeft rw ' ')) col)
          (colWidths cols)
          cols)
          
-- Value pair -> String pair (lifted into list for use with >>=)
ratioPower :: (Rational, Integer) -> [(String, String)]
ratioPower (nd, j) =
  let (num, den) = ((,) . numerator <*> denominator) nd
      sn
        | num == 0 = []
        | (j /= 1) = ("n^" <> show j)
        | otherwise = "n"
      sr
        | num == 0 = []
        | den == 1 && num == 1 = []
        | den == 1 = show num <> "n"
        | otherwise = intercalate "/" [show num, show den]
      s = sr <> sn
  in bool [(sn, sr)] [] (null s)
 
-- Rows of uneven length -> All rows padded to length of longest
fullTable :: [[(String, String)]] -> [[(String, String)]]
fullTable xs =
  let lng = maximum $ length <$> xs
  in (<>) <*> (flip replicate ([], []) . (-) lng . length) <$> xs
 
justifyLeft :: Int -> Char -> String -> String
justifyLeft n c s = take n (s <> replicate n c)
 
-- (Row index, Expression pairs) -> String joined by conjunctions
expressionRow :: (Int, [(String, String)]) -> String
expressionRow (i, row) =
  concat
    [ show i
    , " ->  "
    , foldr
        (\s a -> concat [s, bool " + " " " (blank a || head a == '-'), a])
        []
        (polyTerm <$> row)
    ]
 
-- (Power string, Ratio String) -> Combined string with possible '*'
polyTerm :: (String, String) -> String
polyTerm (l, r)
  | blank l || blank r = l <> r
  | head r == '-' = concat ["- ", l, " * ", tail r]
  | otherwise = intercalate " * " [l, r]
 
blank :: String -> Bool
blank = all isSpace
 
-- Maximum widths of power and ratio elements in each column
colWidths :: [[(String, String)]] -> [(Int, Int)]
colWidths =
  fmap
    (foldr
       (\(ls, rs) (lMax, rMax) -> (max (length ls) lMax, max (length rs) rMax))
       (0, 0))
 
-- Length of string excluding any leading '-'
unsignedLength :: String -> Int
unsignedLength xs =
  let l = length xs
  in bool (bool l (l - 1) ('-' == head xs)) 0 (0 == l)

Output:

0 ->  n                                                 
1 ->  n^2  * 1/2  + n   * 1/2                                   
2 ->  n^3  * 1/3  + n^2 * 1/2 + n   * 1/6                            
3 ->  n^4  * 1/4  + n^3 * 1/2 + n^2 * 1/4                            
4 ->  n^5  * 1/5  + n^4 * 1/2 + n^3 * 1/3  - n   * 1/30                  
5 ->  n^6  * 1/6  + n^5 * 1/2 + n^4 * 5/12 - n^2 * 1/12                  
6 ->  n^7  * 1/7  + n^6 * 1/2 + n^5 * 1/2  - n^3 * 1/6  + n   * 1/42          
7 ->  n^8  * 1/8  + n^7 * 1/2 + n^6 * 7/12 - n^4 * 7/24 + n^2 * 1/12          
8 ->  n^9  * 1/9  + n^8 * 1/2 + n^7 * 2/3  - n^5 * 7/15 + n^3 * 2/9  - n   * 1/30 
9 ->  n^10 * 1/10 + n^9 * 1/2 + n^8 * 3/4  - n^6 * 7/10 + n^4 * 1/2  - n^2 * 3/20

J

Implementation:

Bsecond=:verb define"0
  +/,(<:*(_1^[)*!*(y^~1+[)%1+])"0/~i.1x+y
)

Bfirst=: Bsecond - 1&=

Faul=:adverb define
  (0,|.(%m+1x) * (_1x&^ * !&(m+1) * Bfirst) i.1+m)&p.
)

Task example:

   0 Faul
0 1x&p.
   1 Faul
0 1r2 1r2&p.
   2 Faul
0 1r6 1r2 1r3&p.
   3 Faul
0 0 1r4 1r2 1r4&p.
   4 Faul
0 _1r30 0 1r3 1r2 1r5&p.
   5 Faul
0 0 _1r12 0 5r12 1r2 1r6&p.
   6 Faul
0 1r42 0 _1r6 0 1r2 1r2 1r7&p.
   7 Faul
0 0 1r12 0 _7r24 0 7r12 1r2 1r8&p.
   8 Faul
0 _1r30 0 2r9 0 _7r15 0 2r3 1r2 1r9&p.
   9 Faul
0 0 _3r20 0 1r2 0 _7r10 0 3r4 1r2 1r10&p.

Double checking our work:

   Fcheck=: dyad def'+/(1+i.y)^x'"0
   9 Faul i.5
0 1 513 20196 282340
   9 Fcheck i.5
0 1 513 20196 282340
   2 Faul i.5
0 1 5 14 30
   2 Fcheck i.5
0 1 5 14 30

Java

Translation of: Kotlin

Works with: Java version 8

import java.util.Arrays;
import java.util.stream.IntStream;

public class FaulhabersFormula {
    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 static final Frac ONE = new Frac(1, 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);
        }

        private double toDouble() {
            return (double) num / denom;
        }
    }

    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 int binomial(int n, int k) {
        if (n < 0 || k < 0 || n < k) throw new IllegalArgumentException();
        if (n == 0 || k == 0) return 1;
        int num = IntStream.rangeClosed(k + 1, n).reduce(1, (a, b) -> a * b);
        int den = IntStream.rangeClosed(2, n - k).reduce(1, (acc, i) -> acc * i);
        return num / den;
    }

    private static void faulhaber(int p) {
        System.out.printf("%d : ", p);
        Frac q = new Frac(1, p + 1);
        int sign = -1;
        for (int j = 0; j <= p; ++j) {
            sign *= -1;
            Frac coeff = q.times(new Frac(sign, 1)).times(new Frac(binomial(p + 1, j), 1)).times(bernoulli(j));
            if (Frac.ZERO.equals(coeff)) continue;
            if (j == 0) {
                if (!Frac.ONE.equals(coeff)) {
                    if (Frac.ONE.unaryMinus().equals(coeff)) {
                        System.out.print("-");
                    } else {
                        System.out.print(coeff);
                    }
                }
            } else {
                if (Frac.ONE.equals(coeff)) {
                    System.out.print(" + ");
                } else if (Frac.ONE.unaryMinus().equals(coeff)) {
                    System.out.print(" - ");
                } else if (coeff.compareTo(Frac.ZERO) > 0) {
                    System.out.printf(" + %s", coeff);
                } else {
                    System.out.printf(" - %s", coeff.unaryMinus());
                }
            }
            int pwr = p + 1 - j;
            if (pwr > 1)
                System.out.printf("n^%d", pwr);
            else
                System.out.print("n");
        }
        System.out.println();
    }

    public static void main(String[] args) {
        for (int i = 0; i <= 9; ++i) {
            faulhaber(i);
        }
    }
}

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

jq

Works with: jq

Works with gojq, the Go implementation of jq, and with fq

The following assumes the following rational arithmetic functions as provided by the jq "rational" module at Arithmetic/Rational#jq:

r/2 (constructor)
requal/2 (equality)
rlessthan/1
rminus/1
rminus/2 (subtraction)
rmult/0 and rmult/2 (multiplication)

In addition, the pretty-printing function defined here (rpp) assumes rationals are JSON objects of the form {n,d}

include "rational" {search: "."}; # see [[Arithmetic/Rational#jq]]:

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

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

# The following definition reflects 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;


# The task
def faulhaber($p):
  # The traditional view:
  def bernouilli($n):
     $n | bernouilli | if $n==1 then rminus else . end;
   
  r(1; $p+1) as $q
  | { sign: -1 }
  | reduce range(0; 1+$p) as $j (.;
      .sign *= -1
      | r(binomial($p+1; $j); 1) as $b
      | ([$q, .sign, $b, ($j|bernoulli)] | rmult) as $coeff
      | if requal($coeff; r(0;1))|not
        then .emit +=
           (if $j == 0
            then (if   requal($coeff; r( 1;1)) then ""
                  elif requal($coeff; r(-1;1)) then "-"
                  else "\($coeff|rpp) " end)
            else (if requal($coeff; r(1;1)) then " + "
                  elif requal($coeff; r(-1;1)) then " - " 
                  elif r(0;1)|rlessthan($coeff) then " + \($coeff|rpp) "
                  else " - \($coeff|rminus|rpp) "
                  end)
            end)
        | ($p + 1 - $j) as $pwr
        | .emit += (if 1 < $pwr then "n^\($pwr)" else "n" end)
        else .
        end )
 | .emit ;

range(0;10) |  "\(.) : \(faulhaber(.))"

Output:

0 : n
1 : 1/2 n^2 - 1/2 n
2 : 1/3 n^3 - 1/2 n^2 + 1/6 n
3 : 1/4 n^4 - 1/2 n^3 + 1/4 n^2
4 : 1/5 n^5 - 1/2 n^4 + 1/3 n^3 - 1/30 n
5 : 1/6 n^6 - 1/2 n^5 + 5/12 n^4 - 1/12 n^2
6 : 1/7 n^7 - 1/2 n^6 + 1/2 n^5 - 1/6 n^3 + 1/42 n
7 : 1/8 n^8 - 1/2 n^7 + 7/12 n^6 - 7/24 n^4 + 1/12 n^2
8 : 1/9 n^9 - 1/2 n^8 + 2/3 n^7 - 7/15 n^5 + 2/9 n^3 - 1/30 n
9 : 1/10 n^10 - 1/2 n^9 + 3/4 n^8 - 7/10 n^6 + 1/2 n^4 - 3/20 n^2

Julia

Translation of: Kotlin

Module:

module Faulhaber

function bernoulli(n::Integer)
    n ≥ 0 || throw(DomainError(n, "n must be a positive-or-0 number"))
    a = fill(0 // 1, n + 1)
    for m in 1:n
        a[m] = 1 // (m + 1)
        for j in m:-1:2
            a[j - 1] = (a[j - 1] - a[j]) * j
        end
    end
    return ifelse(n != 1, a[1], -a[1])
end

const _exponents = collect(Char, "⁰¹²³⁴⁵⁶⁷⁸⁹")
toexponent(n) = join(_exponents[reverse(digits(n)) .+ 1])

function formula(p::Integer)
    print(p, ": ")
    q = 1 // (p + 1)
    s = -1
    for j in 0:p
        s *= -1
        coeff = q * s * binomial(p + 1, j) * bernoulli(j)
        iszero(coeff) && continue
        if iszero(j)
            print(coeff == 1 ? "" : coeff == -1 ? "-" : "$coeff")
        else
            print(coeff == 1 ? " + " : coeff == -1 ? " - " : coeff > 0 ? " + $coeff " : " - $(-coeff) ")
        end
        pwr = p + 1 - j
        if pwr > 1
            print("n", toexponent(pwr))
        else
            print("n")
        end
    end
    println()
end

end  # module Faulhaber

Main:

Faulhaber.formula.(1:10)

Output:

1:  + 1//2 n
2:  + 1//2 n² + 1//3 n
3:  + 1//2 n³ + 1//2 n² - 1//6 n
4:  + 1//2 n⁴ + 2//3 n³ - 1//3 n² + 1//30 n
5:  + 1//2 n⁵ + 5//6 n⁴ - 5//9 n³ + 1//12 n² + 1//30 n
6:  + 1//2 n⁶ + n⁵ - 5//6 n⁴ + 1//6 n³ + 1//10 n² - 1//42 n
7:  + 1//2 n⁷ + 7//6 n⁶ - 7//6 n⁵ + 7//24 n⁴ + 7//30 n³ - 1//12 n² - 1//42 n
8:  + 1//2 n⁸ + 4//3 n⁷ - 14//9 n⁶ + 7//15 n⁵ + 7//15 n⁴ - 2//9 n³ - 2//21 n² + 1//30 n
9:  + 1//2 n⁹ + 3//2 n⁸ - 2//1 n⁷ + 7//10 n⁶ + 21//25 n⁵ - 1//2 n⁴ - 2//7 n³ + 3//20 n² + 1//30 n
10:  + 1//2 n¹⁰ + 5//3 n⁹ - 5//2 n⁸ + n⁷ + 7//5 n⁶ - n⁵ - 5//7 n⁴ + 1//2 n³ + 1//6 n² - 5//66 n

Kotlin

As Kotlin doesn't have support for rational numbers built in, a cut-down version of the Frac class in the Arithmetic/Rational task has been used in order to express the polynomial coefficients as fractions.

// version 1.1.2

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 bernoulli(n: Int): Frac {
    require(n >= 0)
    val a = Array<Frac>(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): Int {
    require(n >= 0 && k >= 0 && n >= k) 
    if (n == 0 || k == 0) return 1
    val num = (k + 1..n).fold(1) { acc, i -> acc * i }
    val den = (2..n - k).fold(1) { acc, i -> acc * i }
    return num / den
}

fun faulhaber(p: Int) {
    print("$p : ")
    val q = Frac(1, p + 1)
    var sign = -1
    for (j in 0..p) {        
        sign *= -1
        val coeff = q * Frac(sign, 1) * Frac(binomial(p + 1, j), 1) * bernoulli(j)
        if (coeff == Frac.ZERO) continue
        if (j == 0) {
            print(when {
                coeff == Frac.ONE  -> ""
                coeff == -Frac.ONE -> "-"
                else               -> "$coeff"
            }) 
        }
        else { 
            print(when {
                coeff == Frac.ONE  -> " + "
                coeff == -Frac.ONE -> " - "
                coeff >  Frac.ZERO -> " + $coeff"
                else               -> " - ${-coeff}"
            })
        } 
        val pwr = p + 1 - j
        if (pwr > 1)
            print("n^${p + 1 - j}")
        else
            print("n")
    }
    println()
}


fun main(args: Array<String>) {    
    for (i in 0..9) faulhaber(i)
}

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

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)
    io.write(f.num)
    if f.denom ~= 1 then
        io.write("/"..f.denom)
    end
    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)
    io.write(p.." : ")
    local q = makeFrac(1, p+1)
    local sign = -1
    for j=0,p do
        sign = -1 * sign
        local coeff = multFrac(multFrac(multFrac(q, makeFrac(sign, 1)), makeFrac(binomial(p + 1, j), 1)), bernoulli(j))
        if not equalFrac(coeff, makeFrac(0, 1)) then
            if j==0 then
                if not equalFrac(coeff, makeFrac(1, 1)) then
                    if equalFrac(coeff, makeFrac(-1, 1)) then
                        io.write("-")
                    else
                        printFrac(coeff)
                    end
                end
            else
                if equalFrac(coeff, makeFrac(1, 1)) then
                    io.write(" + ")
                elseif equalFrac(coeff, makeFrac(-1, 1)) then
                    io.write(" - ")
                elseif lessFrac(makeFrac(0, 1), coeff) then
                    io.write(" + ")
                    printFrac(coeff)
                else
                    io.write(" - ")
                    printFrac(negateFrac(coeff))
                end
            end

            local pwr = p + 1 - j
            if pwr>1 then
                io.write("n^"..pwr)
            else
                io.write("n")
            end
        end
    end
    print()
    return nil
end

-- main
for i=0,9 do
    faulhaber(i)
end

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

Mathematica / Wolfram Language

ClearAll[Faulhaber]
Faulhaber[n_, 0] := n
Faulhaber[n_, p_] := n^(p + 1)/(p + 1) + 1/2 n^p + Sum[BernoulliB[k]/k! p!/(p - k + 1)! n^(p - k + 1), {k, 2, p}]
Table[{p, Faulhaber[n, p]}, {p, 0, 9}] // Grid

Output:

0	n
1	n/2+n^2/2
2	n/6+n^2/2+n^3/3
3	n^2/4+n^3/2+n^4/4
4	-(n/30)+n^3/3+n^4/2+n^5/5
5	-(n^2/12)+(5 n^4)/12+n^5/2+n^6/6
6	n/42-n^3/6+n^5/2+n^6/2+n^7/7
7	n^2/12-(7 n^4)/24+(7 n^6)/12+n^7/2+n^8/8
8	-(n/30)+(2 n^3)/9-(7 n^5)/15+(2 n^7)/3+n^8/2+n^9/9
9	-((3 n^2)/20)+n^4/2-(7 n^6)/10+(3 n^8)/4+n^9/2+n^10/10

Maxima

sum1(p):=sum(stirling2(p,k)*pochhammer(n-k+1,k+1)/(k+1),k,0,p)$
sum2(p):=sum((-1)^j*binomial(p+1,j)*bern(j)*n^(p-j+1),j,0,p)/(p+1)$

makelist(expand(sum1(p)-sum2(p)),p,1,10);
[0,0,0,0,0,0,0,0,0,0]

for p from 0 thru 9 do print(expand(sum2(p)));

Output:

n
n^2/2+n/2
n^3/3+n^2/2+n/6
n^4/4+n^3/2+n^2/4
n^5/5+n^4/2+n^3/3-n/30
n^6/6+n^5/2+(5*n^4)/12-n^2/12
n^7/7+n^6/2+n^5/2-n^3/6+n/42
n^8/8+n^7/2+(7*n^6)/12-(7*n^4)/24+n^2/12
n^9/9+n^8/2+(2*n^7)/3-(7*n^5)/15+(2*n^3)/9-n/30
n^10/10+n^9/2+(3*n^8)/4-(7*n^6)/10+n^4/2-(3*n^2)/20

Modula-2

Translation of: C#

MODULE Faulhaber;
FROM EXCEPTIONS IMPORT AllocateSource,ExceptionSource,GetMessage,RAISE;
FROM FormatString IMPORT FormatString;
FROM Terminal IMPORT WriteString,WriteLn,ReadChar;

VAR TextWinExSrc : ExceptionSource;

(* Helper Functions *)
PROCEDURE Abs(n : INTEGER) : INTEGER;
BEGIN
    IF n < 0 THEN
        RETURN -n
    END;
    RETURN n
END Abs;

PROCEDURE Binomial(n,k : INTEGER) : INTEGER;
VAR i,num,denom : INTEGER;
BEGIN
    IF (n < 0) OR (k < 0) OR (n < k) THEN
        RAISE(TextWinExSrc, 0, "Argument Exception.")
    END;
    IF (n = 0) OR (k = 0) THEN
        RETURN 1
    END;
    num := 1;
    FOR i:=k+1 TO n DO
        num := num * i
    END;
    denom := 1;
    FOR i:=2 TO n - k DO
        denom := denom * i
    END;
    RETURN num / denom
END Binomial;

PROCEDURE GCD(a,b : INTEGER) : INTEGER;
BEGIN
    IF b = 0 THEN
        RETURN a
    END;
    RETURN GCD(b, a MOD b)
END GCD;

PROCEDURE WriteInteger(n : INTEGER);
VAR buf : ARRAY[0..15] OF CHAR;
BEGIN
    FormatString("%i", buf, n);
    WriteString(buf)
END WriteInteger;

(* Fraction Handling *)
TYPE Frac = RECORD
    num,denom : INTEGER;
END;

PROCEDURE InitFrac(n,d : INTEGER) : Frac;
VAR nn,dd,g : INTEGER;
BEGIN
    IF d = 0 THEN
        RAISE(TextWinExSrc, 0, "The denominator must not be zero.")
    END;
    IF n = 0 THEN
        d := 1
    ELSIF d < 0 THEN
        n := -n;
        d := -d
    END;
    g := Abs(GCD(n, d));
    IF g > 1 THEN
        n := n / g;
        d := d / g
    END;
    RETURN Frac{n, d}
END InitFrac;

PROCEDURE EqualFrac(a,b : Frac) : BOOLEAN;
BEGIN
    RETURN (a.num = b.num) AND (a.denom = b.denom)
END EqualFrac;

PROCEDURE LessFrac(a,b : Frac) : BOOLEAN;
BEGIN
    RETURN a.num * b.denom < b.num * a.denom
END LessFrac;

PROCEDURE NegateFrac(f : Frac) : Frac;
BEGIN
    RETURN Frac{-f.num, f.denom}
END NegateFrac;

PROCEDURE SubFrac(lhs,rhs : Frac) : Frac;
BEGIN
    RETURN InitFrac(lhs.num * rhs.denom - lhs.denom * rhs.num, rhs.denom * lhs.denom)
END SubFrac;

PROCEDURE MultFrac(lhs,rhs : Frac) : Frac;
BEGIN
    RETURN InitFrac(lhs.num * rhs.num, lhs.denom * rhs.denom)
END MultFrac;

PROCEDURE Bernoulli(n : INTEGER) : Frac;
VAR
    a : ARRAY[0..15] OF Frac;
    i,j,m : INTEGER;
BEGIN
    IF n < 0 THEN
        RAISE(TextWinExSrc, 0, "n may not be negative or zero.")
    END;
    FOR m:=0 TO n DO
        a[m] := Frac{1, m + 1};
        FOR j:=m TO 1 BY -1 DO
            a[j-1] := MultFrac(SubFrac(a[j-1], a[j]), Frac{j, 1})
        END
    END;
    IF n # 1 THEN RETURN a[0] END;
    RETURN NegateFrac(a[0])
END Bernoulli;

PROCEDURE WriteFrac(f : Frac);
BEGIN
    WriteInteger(f.num);
    IF f.denom # 1 THEN
        WriteString("/");
        WriteInteger(f.denom)
    END
END WriteFrac;

(* Target *)
PROCEDURE Faulhaber(p : INTEGER);
VAR
    j,pwr,sign : INTEGER;
    q,coeff : Frac;
BEGIN
    WriteInteger(p);
    WriteString(" : ");
    q := InitFrac(1, p + 1);
    sign := -1;
    FOR j:=0 TO p DO
        sign := -1 * sign;
        coeff := MultFrac(MultFrac(MultFrac(q, Frac{sign, 1}), Frac{Binomial(p + 1, j), 1}), Bernoulli(j));
        IF EqualFrac(coeff, Frac{0, 1}) THEN CONTINUE END;
        IF j = 0 THEN
            IF NOT EqualFrac(coeff, Frac{1, 1}) THEN
                IF EqualFrac(coeff, Frac{-1, 1}) THEN
                    WriteString("-")
                ELSE
                    WriteFrac(coeff)
                END
            END
        ELSE
            IF EqualFrac(coeff, Frac{1, 1}) THEN
                WriteString(" + ")
            ELSIF EqualFrac(coeff, Frac{-1, 1}) THEN
                WriteString(" - ")
            ELSIF LessFrac(Frac{0, 1}, coeff) THEN
                WriteString(" + ");
                WriteFrac(coeff)
            ELSE
                WriteString(" - ");
                WriteFrac(NegateFrac(coeff))
            END
        END;
        pwr := p + 1 - j;
        IF pwr > 1 THEN
            WriteString("n^");
            WriteInteger(pwr)
        ELSE
            WriteString("n")
        END
    END;
    WriteLn
END Faulhaber;

(* Main *)
VAR i : INTEGER;
BEGIN
    FOR i:=0 TO 9 DO
        Faulhaber(i)
    END;
    ReadChar
END Faulhaber.

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

Nim

Translation of: Kotlin

import math, rationals

type
  Fraction = Rational[int]
  FaulhaberSequence = seq[Fraction]

const
  Zero = 0 // 1
  One = 1 // 1
  MinusOne = -1 // 1
  Powers = ["⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹"]

#---------------------------------------------------------------------------------------------------

func bernoulli(n: Natural): Fraction =
  ## Return nth Bernoulli coefficient.

  var a = newSeq[Fraction](n + 1)
  for m in 0..n:
    a[m] = 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.

  var a = 1 // (n + 1)
  var sign = -1
  for k in 0..n:
    sign = -sign
    result.add(a * sign * binom(n + 1, k) * bernoulli(k))

#---------------------------------------------------------------------------------------------------

func npower(k: Natural): string =
  ## Return the string representing "n" at power "k".

  if k == 0: return ""
  if k == 1: return "n"
  var k = k
  result = "n"
  while k != 0:
    result.insert(Powers[k mod 10], 1)
    k = k div 10

#---------------------------------------------------------------------------------------------------

func `$`(fs: FaulhaberSequence): string =
  ## Return the string representing a Faulhaber sequence.
  
  for i, coeff in fs:

    # Process coefficient.
    if coeff.num == 0: continue
    if i == 0:
      if coeff == MinusOne: result.add(" - ")
      elif coeff != One: result.add($coeff)
    else:
      if coeff == One: result.add(" + ")
      elif coeff == MinusOne: result.add(" - ")
      elif coeff > Zero: result.add(" + " & $coeff)
      else: result.add(" - " & $(-coeff))

    # Process power of "n".
    let pwr = fs.len - i
    result.add(npower(pwr))

#———————————————————————————————————————————————————————————————————————————————————————————————————

for n in 0..9:
  echo n, ": ", faulhaber(n)

Output:

0: n
1: 1/2n² + 1/2n
2: 1/3n³ + 1/2n² + 1/6n
3: 1/4n⁴ + 1/2n³ + 1/4n²
4: 1/5n⁵ + 1/2n⁴ + 1/3n³ - 1/30n
5: 1/6n⁶ + 1/2n⁵ + 5/12n⁴ - 1/12n²
6: 1/7n⁷ + 1/2n⁶ + 1/2n⁵ - 1/6n³ + 1/42n
7: 1/8n⁸ + 1/2n⁷ + 7/12n⁶ - 7/24n⁴ + 1/12n²
8: 1/9n⁹ + 1/2n⁸ + 2/3n⁷ - 7/15n⁵ + 2/9n³ - 1/30n
9: 1/10n¹⁰ + 1/2n⁹ + 3/4n⁸ - 7/10n⁶ + 1/2n⁴ - 3/20n²

PARI/GP

PARI/GP has 2 built in functions: bernfrac(n) for Bernoulli numbers and bernpol(n) for Bernoulli polynomials. Using Bernoulli polynomials in PARI/GP is more simple, clear and much faster. (See version #2).

Works with: PARI/GP version 2.9.1 and above

Version #1. Using Bernoulli numbers.

This version is using "Faulhaber's" formula based on Bernoulli numbers.
It's not worth using Bernoulli numbers in PARI/GP, because too much cleaning if you are avoiding "dirty" (but correct) result.
Note: Find ssubstr() function here on RC.

\\ Using "Faulhaber's" formula based on Bernoulli numbers. aev 2/7/17
\\ In str string replace all occurrences of the search string ssrch with the replacement string srepl. aev  3/8/16
sreplace(str,ssrch,srepl)={
  my(sn=#str,ssn=#ssrch,srn=#srepl,sres="",Vi,Vs,Vz,vin,vin1,vi,L=List(),tok,zi,js=1);
  if(sn==0,return("")); if(ssn==0||ssn>sn,return(str));
  \\ Vi - found ssrch indexes
  Vi=sfindalls(str,ssrch); vin=#Vi;
  if(vin==0,return(str));
  vin1=vin+1; Vi=Vec(Vi,vin1); Vi[vin1]=sn+1;
  for(i=1,vin1, vi=Vi[i];
  for(j=js,sn, \\print("ij:",i,"/",j,": ",sres);
    if(j!=vi, sres=concat(sres,ssubstr(str,j,1)),
              sres=concat(sres,srepl); js=j+ssn; break) 
  ); \\fend j
  ); \\fend i
  return(sres);
}
B(n)=(bernfrac(n));
Comb(n,k)={my(r=0); if(k<=n, r=n!/(n-k)!/k!); return(r)};
Faulhaber2(p)={
  my(s="",s1="",s2="",c1=0,c2=0);
  for(j=0,p, c1=(-1)^j*Comb(p+1,j)*B(j); c2=(p+1-j);
    s2="*n";
    if(c1==0, next);
    if(c2==1, s1="", s1=Str("^",c2));
    s=Str(s,c1,s2,s1,"+") );
  s=ssubstr(s,1,#s-1); s=sreplace(s,"1*n","n"); s=sreplace(s,"+-","-");
  if(p==0, s="n", s=Str("(",s,")/",p+1)); print(p,": ",s);
}
{\\ Testing:
for(i=0,9, Faulhaber2(i))}

Output:

0: n
1: (n^2+n)/2
2: (n^3+3/2*n^2+1/2*n)/3
3: (n^4+2*n^3+n^2)/4
4: (n^5+5/2*n^4+5/3*n^3-1/6*n)/5
5: (n^6+3*n^5+5/2*n^4-1/2*n^2)/6
6: (n^7+7/2*n^6+7/2*n^5-7/6*n^3+1/6*n)/7
7: (n^8+4*n^7+14/3*n^6-7/3*n^4+2/3*n^2)/8
8: (n^9+9/2*n^8+6*n^7-21/5*n^5+2*n^3-3/10*n)/9
9: (n^10+5*n^9+15/2*n^8-7*n^6+5*n^4-3/2*n^2)/10
time = 16 ms.

Version #2. Using Bernoulli polynomials.

This version is using the sums of pth powers formula from Bernoulli polynomials. It has small, simple and clear code, and produces instant result.

\\ Using a formula based on Bernoulli polynomials. aev 2/5/17
Faulhaber1(m)={
  my(B,B1,B2,Bn);
  Bn=bernpol(m+1);
  x=n+1; B1=eval(Bn); x=0; B2=eval(Bn);
  Bn=(B1-B2)/(m+1); if(m==0, Bn=Bn-1);
  print(m,": ",Bn);
}
{\\ Testing:
  for(i=0,9, Faulhaber1(i))}

Output:

0: n                          
1: 1/2*n^2 + 1/2*n
2: 1/3*n^3 + 1/2*n^2 + 1/6*n
3: 1/4*n^4 + 1/2*n^3 + 1/4*n^2
4: 1/5*n^5 + 1/2*n^4 + 1/3*n^3 - 1/30*n
5: 1/6*n^6 + 1/2*n^5 + 5/12*n^4 - 1/12*n^2
6: 1/7*n^7 + 1/2*n^6 + 1/2*n^5 - 1/6*n^3 + 1/42*n
7: 1/8*n^8 + 1/2*n^7 + 7/12*n^6 - 7/24*n^4 + 1/12*n^2
8: 1/9*n^9 + 1/2*n^8 + 2/3*n^7 - 7/15*n^5 + 2/9*n^3 - 1/30*n
9: 1/10*n^10 + 1/2*n^9 + 3/4*n^8 - 7/10*n^6 + 1/2*n^4 - 3/20*n^2
> ##
  ***   last result computed in 0 ms

Pascal

A console program that runs under Lazarus or Delphi. Does not make use of Bernoulli numbers.

program Faulhaber;

{$IFDEF FPC} // Lazarus
  {$MODE Delphi} // ensure Lazarus accepts Delphi-style code
  {$ASSERTIONS+} // by default, Lazarus does not compile 'Assert' statements
{$ELSE}     // Delphi
  {$APPTYPE CONSOLE}
{$ENDIF}

uses SysUtils;

type TRational = record
  Num, Den : integer; // where Den > 0 and Num, Den are coprime
end;

const
  ZERO : TRational = ( Num: 0; Den : 1);
  HALF : TRational = ( Num: 1; Den : 2);

// Construct rational a/b, assuming b > 0.
function Rational( const a, b : integer) : TRational;
var
  t, x, y : integer;
begin
  if b <= 0 then raise SysUtils.Exception.Create( 'Denominator must be > 0');
  // Find HCF of a and b (Euclid's algorithm) and cancel it out.
  x := Abs(a);
  y := b;
  while y <> 0 do begin
    t := x mod y;
    x := y;
    y := t;
  end;
  result.Num := a div x;
  result.Den := b div x
end;

function Prod( r, s : TRational) : TRational; // result := r*s
begin
  result := Rational( r.Num*s.Num, r.Den*s.Den);
end;

procedure DecRat( var r : TRational;
                const s : TRational); // r := r - s
begin
  r := Rational( r.Num*s.Den - s.Num*r.Den, r.Den * s.Den);
end;

// Write a term such as ' - (7/10)n^6' to the console.
procedure WriteTerm( coeff : TRational;
                     index : integer;
                     printPlus : boolean);
begin
  if Coeff.Num = 0 then exit;
  with coeff do begin
    if Num < 0 then Write(' - ')
    else if printPlus then Write(' + ');
    // Put brackets round a fractional coefficient
    if (Den > 1) then Write('(');
    // If coefficient is 1, don't write it
    if (Den > 1) or (Abs(Num) > 1) then Write( Abs(Num));
    // Write denominator if it's not 1
    if (Den > 1) then Write('/', Den, ')');
  end;
  Write('n');
  if index > 1 then Write('^', index);
end;

{-------------------------------------------------------------------------------
Main routine. Calculation of Faulhaber polynomials
  F_p(n) = 1^p + 2^p + ... + n^p,  p = 0, 1, ..., p_max
}
var
  p_max : integer;
  c : array of array of TRational;
  i, j, p : integer;
  coeff_of_n : TRational;
begin
  // User types program name, optionally followed by maximum power p (defaults to 9)
  if ParamCount = 0 then p_max := 9
                    else p_max := SysUtils.StrToInt( ParamStr(1));

  // c[p, i] is coefficient of n^i in the polynomial F_p(n).
  // Initialize all coefficients to 0.
  SetLength( c, p_max + 1, p_max + 2);
  for i := 0 to p_max do
    for j := 0 to p_max + 1 do
      c[i, j] := ZERO;

  c[0, 1] := Rational(1, 1); // F_0(n) = n, special case
  for p := 1 to p_max do begin
    // Initialize calculation of coefficient of n, needed if p is even.
    // If p is odd, still calculate it as a check on the working (should be 0).
    // Calculation uses the fact that F_p(1) = 1.
    coeff_of_n := Rational(1, 1);

    c[p, p+1] := Rational(1, p + 1);
    DecRat( coeff_of_n, c[p, p + 1]);
    c[p, p] := HALF;
    DecRat( coeff_of_n, c[p, p]);
    i := p - 1;
    while (i >= 2) do begin
      c[p, i] := Prod( Rational(p, i), c[p - 1, i - 1]);
      DecRat( coeff_of_n, c[p, i]);
      dec(i, 2);
    end;
    if i = 1 then // p is even
      c[p, 1] := coeff_of_n // = the Bernoulli number B_p
    else // p is odd
      Assert( coeff_of_n.Num = 0); // just checking
  end; // for p

  // Print the result
  for p := 0 to p_max do begin
    Write( 'F_', p, '(n) = ');
    for j := p + 1 downto 1 do WriteTerm( c[p, j], j, j <= p);
    WriteLn;
  end;
end.

Output:

F_0(n) = n
F_1(n) = (1/2)n^2 + (1/2)n
F_2(n) = (1/3)n^3 + (1/2)n^2 + (1/6)n
F_3(n) = (1/4)n^4 + (1/2)n^3 + (1/4)n^2
F_4(n) = (1/5)n^5 + (1/2)n^4 + (1/3)n^3 - (1/30)n
F_5(n) = (1/6)n^6 + (1/2)n^5 + (5/12)n^4 - (1/12)n^2
F_6(n) = (1/7)n^7 + (1/2)n^6 + (1/2)n^5 - (1/6)n^3 + (1/42)n
F_7(n) = (1/8)n^8 + (1/2)n^7 + (7/12)n^6 - (7/24)n^4 + (1/12)n^2
F_8(n) = (1/9)n^9 + (1/2)n^8 + (2/3)n^7 - (7/15)n^5 + (2/9)n^3 - (1/30)n
F_9(n) = (1/10)n^10 + (1/2)n^9 + (3/4)n^8 - (7/10)n^6 + (1/2)n^4 - (3/20)n^2

Perl

use 5.014;
use Math::Algebra::Symbols;

sub bernoulli_number {
    my ($n) = @_;

    return 0 if $n > 1 && $n % 2;

    my @A;
    for my $m (0 .. $n) {
        $A[$m] = symbols(1) / ($m + 1);

        for (my $j = $m ; $j > 0 ; $j--) {
            $A[$j - 1] = $j * ($A[$j - 1] - $A[$j]);
        }
    }

    return $A[0];
}

sub binomial {
    my ($n, $k) = @_;
    return 1 if $k == 0 || $n == $k;
    binomial($n - 1, $k - 1) + binomial($n - 1, $k);
}

sub faulhaber_s_formula {
    my ($p) = @_;

    my $formula = 0;
    for my $j (0 .. $p) {
        $formula += binomial($p + 1, $j)
                 *  bernoulli_number($j)
                 *  symbols('n')**($p + 1 - $j);
    }

    (symbols(1) / ($p + 1) * $formula)
        =~ s/\$n/n/gr =~ s/\*\*/^/gr =~ s/\*/ /gr;
}

foreach my $i (0 .. 9) {
    say "$i: ", faulhaber_s_formula($i);
}

Output:

0: n
1: 1/2 n+1/2 n^2
2: 1/6 n+1/2 n^2+1/3 n^3
3: 1/4 n^2+1/2 n^3+1/4 n^4
4: -1/30 n+1/3 n^3+1/2 n^4+1/5 n^5
5: -1/12 n^2+5/12 n^4+1/2 n^5+1/6 n^6
6: 1/42 n-1/6 n^3+1/2 n^5+1/2 n^6+1/7 n^7
7: 1/12 n^2-7/24 n^4+7/12 n^6+1/2 n^7+1/8 n^8
8: -1/30 n+2/9 n^3-7/15 n^5+2/3 n^7+1/2 n^8+1/9 n^9
9: -3/20 n^2+1/2 n^4-7/10 n^6+3/4 n^8+1/2 n^9+1/10 n^10

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
    integer 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
 
procedure faulhaber(integer p)
    string res = sprintf("%d : ", p)
    frac q = {1, p+1}
    for j=0 to p do
        frac bj = bernoulli(j)
        if frac_ne(bj,frac_zero) then
            frac coeff = frac_mul({binomial(p+1,j),p+1},bj)
            string s = frac_sprint(coeff)
            if j=0 then
                if s="1" then
                    s = ""
                end if
            else
                if s[1]='-' then
                    s[1..1] = " - "
                else
                    s[1..0] = " + "
                end if
            end if
            res &= s&"n"
            integer pwr = p+1-j
            if pwr>1 then
                res &= sprintf("^%d", pwr)
            end if
        end if
    end for
    printf(1,"%s\n",{res})
end procedure
 
for i=0 to 9 do
    faulhaber(i)
end for

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

Python

The following implementation does not use Bernoulli numbers, but Stirling numbers of the second kind, based on the relation: $m^{n}=\sum _{k=0}^{n}S_{n}^{k}(m)_{k}=\sum _{k=0}^{n}S_{n}^{k}k!{m \choose k}$ .

Then summing: $\sum _{j=0}^{m}j^{n}=\sum _{j=0}^{m}\sum _{k=0}^{n}S_{n}^{k}k!{j \choose k}=\sum _{k=0}^{n}S_{n}^{k}k!{m+1 \choose k+1}=\sum _{k=0}^{n}S_{n}^{k}{\frac {(m+1)_{k+1}}{k+1}}$ .

One has then to expand the product $(m+1)_{k+1}$ in order to get a polynomial in the variable $m$ . Also, for the sum of $j^{0}$ , the sum is too large by one (since we start at $j=0$ ), this has to be taken into account.

Note: a number of the formulae above are invisible to the majority of browsers, including Chrome, IE/Edge, Safari and Opera. They may (subject to the installation of necessary fonts) be visible to some Firefox installations.

from fractions import Fraction

def nextu(a):
    n = len(a)
    a.append(1)
    for i in range(n - 1, 0, -1):
        a[i] = i * a[i] + a[i - 1]
    return a

def nextv(a):
    n = len(a) - 1
    b = [(1 - n) * x for x in a]
    b.append(1)
    for i in range(n):
        b[i + 1] += a[i]
    return b

def sumpol(n):
    u = [0, 1]
    v = [[1], [1, 1]]
    yield [Fraction(0), Fraction(1)]
    for i in range(1, n):
        v.append(nextv(v[-1]))
        t = [0] * (i + 2)
        p = 1
        for j, r in enumerate(u):
            r = Fraction(r, j + 1)
            for k, s in enumerate(v[j + 1]):
                t[k] += r * s
        yield t
        u = nextu(u)

def polstr(a):
    s = ""
    q = False
    n = len(a) - 1
    for i, x in enumerate(reversed(a)):
        i = n - i
        if i < 2:
            m = "n" if i == 1 else ""
        else:
            m = "n^%d" % i
        c = str(abs(x))
        if i > 0:
            if c == "1":
                c = ""
            else:
                m = " " + m
        if x != 0:
            if q:
                t = " + " if x > 0 else " - "
                s += "%s%s%s" % (t, c, m)
            else:
                t = "" if x > 0 else "-"
                s = "%s%s%s" % (t, c, m)
                q = True
    if q:
        return s
    else:
        return "0"

for i, p in enumerate(sumpol(10)):
    print(i, ":", polstr(p))

Output:

0 : n
1 : 1/2 n^2 + 1/2 n
2 : 1/3 n^3 + 1/2 n^2 + 1/6 n
3 : 1/4 n^4 + 1/2 n^3 + 1/4 n^2
4 : 1/5 n^5 + 1/2 n^4 + 1/3 n^3 - 1/30 n
5 : 1/6 n^6 + 1/2 n^5 + 5/12 n^4 - 1/12 n^2
6 : 1/7 n^7 + 1/2 n^6 + 1/2 n^5 - 1/6 n^3 + 1/42 n
7 : 1/8 n^8 + 1/2 n^7 + 7/12 n^6 - 7/24 n^4 + 1/12 n^2
8 : 1/9 n^9 + 1/2 n^8 + 2/3 n^7 - 7/15 n^5 + 2/9 n^3 - 1/30 n
9 : 1/10 n^10 + 1/2 n^9 + 3/4 n^8 - 7/10 n^6 + 1/2 n^4 - 3/20 n^2

Racket

Racket will simplify rational numbers; if this code simplifies the expressions too much for your tastes (e.g. you like 1/1 * (n)) then tweak the simplify... clauses to taste.

#lang racket/base

(require racket/match
         racket/string
         math/number-theory)

(define simplify-arithmetic-expression
  (letrec ((s-a-e
            (match-lambda
              [(list (and op '+) l ... (list '+ m ...) r ...) (s-a-e `(,op ,@l ,@m ,@r))]
              [(list (and op '+) l ... (? number? n1) m ... (? number? n2) r ...) (s-a-e `(,op ,@l ,(+ n1 n2) ,@m ,@r))]
              [(list (and op '+) (app s-a-e l _) ... 0 (app s-a-e r _) ...) (s-a-e `(,op ,@l ,@r))]
              [(list (and op '+) (app s-a-e x _)) (values x #t)]
              [(list (and op '*) l ... (list '* m ...) r ...) (s-a-e `(,op ,@l ,@m ,@r))]
              [(list (and op '*) l ... (? number? n1) m ... (? number? n2) r ...) (s-a-e `(,op ,@l ,(* n1 n2) ,@m ,@r))]
              [(list (and op '*) (app s-a-e l _) ... 1 (app s-a-e r _) ...) (s-a-e `(,op ,@l ,@r))]
              [(list (and op '*) (app s-a-e l _) ... 0 (app s-a-e r _) ...) (values 0 #t)]
              [(list (and op '*) (app s-a-e x _)) (values x #t)]
              [(list 'expt (app s-a-e x x-simplified?) 1) (values x x-simplified?)]
              [(list op (app s-a-e a #f) ...) (values `(,op ,@a) #f)]
              [(list op (app s-a-e a _) ...) (s-a-e `(,op ,@a))]
              [e (values e #f)])))
    s-a-e))

(define (expression->infix-string e)
  (define (parenthesise-maybe s p?)
    (if p? (string-append "(" s ")") s))
  
  (letrec ((e->is
            (lambda (paren?)
              (match-lambda
                [(list (and op (or '+ '- '* '*)) (app (e->is #t) a p?) ...)
                 (define bits (map parenthesise-maybe a p?))
                 (define compound (string-join bits (format " ~a " op)))
                 (values (if paren? (string-append "(" compound ")") compound) #f)]
                [(list 'expt (app (e->is #t) x xp?) (app (e->is #t) n np?))
                 (values (format "~a^~a" (parenthesise-maybe x xp?) (parenthesise-maybe n np?)) #f)]                
                [(? number? (app number->string s)) (values s #f)]
                [(? symbol? (app symbol->string s)) (values s #f)]))))
    (define-values (str needs-parens?) ((e->is #f) e))
    str))

(define (faulhaber p)
  (define p+1 (add1 p))
  (define-values (simpler simplified?)
    (simplify-arithmetic-expression
     `(* ,(/ 1 p+1)
         (+ ,@(for/list ((j (in-range p+1)))
                `(* ,(* (expt -1 j)
                        (binomial p+1 j))
                    (* ,(bernoulli-number j)
                       (expt n ,(- p+1 j)))))))))
  simpler)

(for ((p (in-range 0 (add1 9))))
  (printf "f(~a) = ~a~%" p (expression->infix-string (faulhaber p))))

Output:

f(0) = n
f(1) = 1/2 * (n^2 + n)
f(2) = 1/3 * (n^3 + (3/2 * n^2) + (1/2 * n))
f(3) = 1/4 * (n^4 + (2 * n^3) + n^2)
f(4) = 1/5 * (n^5 + (5/2 * n^4) + (5/3 * n^3) + (-1/6 * n))
f(5) = 1/6 * (n^6 + (3 * n^5) + (5/2 * n^4) + (-1/2 * n^2))
f(6) = 1/7 * (n^7 + (7/2 * n^6) + (7/2 * n^5) + (-7/6 * n^3) + (1/6 * n))
f(7) = 1/8 * (n^8 + (4 * n^7) + (14/3 * n^6) + (-7/3 * n^4) + (2/3 * n^2))
f(8) = 1/9 * (n^9 + (9/2 * n^8) + (6 * n^7) + (-21/5 * n^5) + (2 * n^3) + (-3/10 * n))
f(9) = 1/10 * (n^10 + (5 * n^9) + (15/2 * n^8) + (-7 * n^6) + (5 * n^4) + (-3/2 * n^2))

Raku

(formerly Perl 6)

Works with: Rakudo version 2018.04.01

sub bernoulli_number($n) {

    return 1/2 if $n == 1;
    return 0/1 if $n % 2;

    my @A;
    for 0..$n -> $m {
        @A[$m] = 1 / ($m + 1);
        for $m, $m-1 ... 1 -> $j {
            @A[$j - 1] = $j * (@A[$j - 1] - @A[$j]);
        }
    }

    return @A[0];
}

sub binomial($n, $k) {
    $k == 0 || $n == $k ?? 1 !! binomial($n-1, $k-1) + binomial($n-1, $k);
}

sub faulhaber_s_formula($p) {

    my @formula = gather for 0..$p -> $j {
        take '('
            ~ join('/', (binomial($p+1, $j) * bernoulli_number($j)).Rat.nude)
            ~ ")*n^{$p+1 - $j}";
    }

    my $formula = join(' + ', @formula.grep({!m{'(0/1)*'}}));

    $formula .= subst(rx{ '(1/1)*' }, '', :g);
    $formula .= subst(rx{ '^1'» }, '', :g);

    "1/{$p+1} * ($formula)";
}

for 0..9 -> $p {
    say "f($p) = ", faulhaber_s_formula($p);
}

Output:

f(0) = 1/1 * (n)
f(1) = 1/2 * (n^2 + n)
f(2) = 1/3 * (n^3 + (3/2)*n^2 + (1/2)*n)
f(3) = 1/4 * (n^4 + (2/1)*n^3 + n^2)
f(4) = 1/5 * (n^5 + (5/2)*n^4 + (5/3)*n^3 + (-1/6)*n)
f(5) = 1/6 * (n^6 + (3/1)*n^5 + (5/2)*n^4 + (-1/2)*n^2)
f(6) = 1/7 * (n^7 + (7/2)*n^6 + (7/2)*n^5 + (-7/6)*n^3 + (1/6)*n)
f(7) = 1/8 * (n^8 + (4/1)*n^7 + (14/3)*n^6 + (-7/3)*n^4 + (2/3)*n^2)
f(8) = 1/9 * (n^9 + (9/2)*n^8 + (6/1)*n^7 + (-21/5)*n^5 + (2/1)*n^3 + (-3/10)*n)
f(9) = 1/10 * (n^10 + (5/1)*n^9 + (15/2)*n^8 + (-7/1)*n^6 + (5/1)*n^4 + (-3/2)*n^2)

RPL

Works with: HP version 49

≪ → p
  ≪ 0
     p 0 FOR m
        p 1 + m 1 + COMB
        p m - IBERNOULLI
        IF LASTARG 1 == THEN NEG END
        * EVAL + 'n' *
     -1 STEP
     p 1 + / EXPAN
≫ ≫ 'FAULH' STO

≪ P FAULH ≫ 'P' 0 9 1 SEQ

Output:

1: {'n' '1/2n^2+1/2n' '1/3n^3+1/2n^2+1/6n' '1/4n^4+1/2n^3+1/4n^2' '1/5n^5+1/2n^4+1/3n^3-1/30n' '1/6n^6+1/2n^5+5/12n^4-1/12n^2' '1/7n^7+1/2n^6+1/2n^5-1/6n^3+1/42n' '1/8n^8+1/2n^7+7/12n^6-7/24n^4+1/12n^2' '1/9n^9+1/2n^8+2/3n^7-7/15n^5+2/9n^3-1/30n' '1/10n^10+1/2n^9+3/4n^8-7/10n^6+1/2n^4-3/20n^2'}

Ruby

Translation of: C

def 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

    num = 1
    for i in k+1 .. n do
        num = num * i
    end

    denom = 1
    for i in 2 .. n-k do
        denom = denom * i
    end

    return num / denom
end

def bernoulli(n)
    if n < 0 then
        raise "n cannot be less than zero"
    end

    a = Array.new(16)
    for m in 0 .. n do
        a[m] = Rational(1, m + 1)
        for j in m.downto(1) do
            a[j-1] = (a[j-1] - a[j]) * Rational(j)
        end
    end

    if n != 1 then
        return a[0]
    end
    return -a[0]
end

def faulhaber(p)
    print("%d : " % [p])
    q = Rational(1, p + 1)
    sign = -1
    for j in 0 .. p do
        sign = -1 * sign
        coeff = q * Rational(sign) * Rational(binomial(p+1, j)) * bernoulli(j)
        if coeff == 0 then
            next
        end
        if j == 0 then
            if coeff != 1 then
                if coeff == -1 then
                    print "-"
                else
                    print coeff
                end
            end
        else
            if coeff == 1 then
                print " + "
            elsif coeff == -1 then
                print " - "
            elsif 0 < coeff then
                print " + "
                print coeff
            else
                print " - "
                print -coeff
            end
        end
        pwr = p + 1 - j
        if pwr > 1 then
            print "n^%d" % [pwr]
        else
            print "n"
        end
    end
    print "\n"
end

def main
    for i in 0 .. 9 do
        faulhaber(i)
    end
end

main()

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

Scala

Translation of: Java

import scala.math.Ordering.Implicits.infixOrderingOps

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 faulhaber(p: Int): Unit = {
    print(s"$p : ")
    val q = Frac(1, p + 1)
    var sign = -1
    for (j <- 0 to p) {
      sign *= -1
      val coeff = q * Frac(sign) * Frac(binomial(p + 1, j)) * bernoulli(j)
      if (Frac.ZERO != coeff) {
        if (j == 0) {
          if (Frac.ONE != coeff) {
            if (-Frac.ONE == coeff) {
              print('-')
            } else {
              print(coeff)
            }
          }
        } else {
          if (Frac.ONE == coeff) {
            print(" + ")
          } else if (-Frac.ONE == coeff) {
            print(" - ")
          } else if (coeff > Frac.ZERO) {
            print(s" + ${coeff}")
          } else {
            print(s" - ${-coeff}")
          }
        }
        val pwr = p + 1 - j
        if (pwr > 1) {
          print(s"n^${pwr}")
        } else {
          print('n')
        }
      }
    }
    println()
  }

  def main(args: Array[String]): Unit = {
    for (i <- 0 to 9) {
      faulhaber(i)
    }
  }
}

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

Sidef

func faulhaber_formula(p) {
    (p+1).of { |j|
        Poly(p - j + 1 => 1) * bernoulli(j) * binomial(p+1, j)
    }.sum / (p+1)
}

for p in (^10) {
    printf("%2d: %s\n", p, faulhaber_formula(p))
}

Output:

 0: x
 1: 1/2*x^2 + 1/2*x
 2: 1/3*x^3 + 1/2*x^2 + 1/6*x
 3: 1/4*x^4 + 1/2*x^3 + 1/4*x^2
 4: 1/5*x^5 + 1/2*x^4 + 1/3*x^3 - 1/30*x
 5: 1/6*x^6 + 1/2*x^5 + 5/12*x^4 - 1/12*x^2
 6: 1/7*x^7 + 1/2*x^6 + 1/2*x^5 - 1/6*x^3 + 1/42*x
 7: 1/8*x^8 + 1/2*x^7 + 7/12*x^6 - 7/24*x^4 + 1/12*x^2
 8: 1/9*x^9 + 1/2*x^8 + 2/3*x^7 - 7/15*x^5 + 2/9*x^3 - 1/30*x
 9: 1/10*x^10 + 1/2*x^9 + 3/4*x^8 - 7/10*x^6 + 1/2*x^4 - 3/20*x^2

Visual Basic .NET

Translation of: C#

Module Module1
    Function Gcd(a As Long, b As Long)
        If b = 0 Then
            Return a
        End If
        Return Gcd(b, a Mod b)
    End Function

    Class Frac
        ReadOnly num As Long
        ReadOnly denom As Long

        Public Shared ReadOnly ZERO As New Frac(0, 1)
        Public Shared ReadOnly ONE As 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")
            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

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

        Public Overloads Function GetHashCode() As Integer
            Dim hashCode = 1317992671
            hashCode = hashCode * -1521134295 + num.GetHashCode()
            hashCode = hashCode * -1521134295 + denom.GetHashCode()
            Return hashCode
        End Function

        Public Overloads Function ToString() As String
            If denom = 1 Then Return num.ToString()
            Return String.Format("{0}/{1}", num, 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")
        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 the first Bernoulli number
        If n <> 1 Then Return a(0)
        Return -a(0)
    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

    Sub Faulhaber(p As Integer)
        Console.Write("{0} : ", p)
        Dim q As New Frac(1, p + 1)
        Dim sign = -1
        For j = 0 To p
            sign *= -1
            Dim coeff = q * New Frac(sign, 1) * New Frac(Binomial(p + 1, j), 1) * Bernoulli(j)
            If Frac.ZERO = coeff Then Continue For
            If j = 0 Then
                If Frac.ONE <> coeff Then
                    If -Frac.ONE = coeff Then
                        Console.Write("-")
                    Else
                        Console.Write(coeff.ToString())
                    End If
                End If
            Else
                If Frac.ONE = coeff Then
                    Console.Write(" + ")
                ElseIf -Frac.ONE = coeff Then
                    Console.Write(" - ")
                ElseIf Frac.ZERO < coeff Then
                    Console.Write(" + {0}", coeff.ToString())
                Else
                    Console.Write(" - {0}", (-coeff).ToString())
                End If
            End If
            Dim pwr = p + 1 - j
            If pwr > 1 Then
                Console.Write("n^{0}", pwr)
            Else
                Console.Write("n")
            End If
        Next
        Console.WriteLine()
    End Sub

    Sub Main()
        For i = 0 To 9
            Faulhaber(i)
        Next
    End Sub
End Module

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

Wren

Translation of: Kotlin

Library: Wren-math

Library: Wren-rat

import "./math" for Int
import "./rat" for Rat

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] = Rat.new(1, m+1)
        var j = m
        while (j >= 1) {
            a[j-1] = (a[j-1] - a[j]) * Rat.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 faulhaber = Fn.new { |p|
    System.write("%(p) : ")
    var q = Rat.new(1, p+1)
    var sign = -1
    for (j in 0..p) {
        sign = sign * -1
        var b = Rat.new(binomial.call(p+1, j), 1)
        var coeff = q * Rat.new(sign, 1) * b * bernoulli.call(j)
        if (coeff != Rat.zero) {
            if (j == 0) {
                System.write((coeff == Rat.one) ? "" : (coeff == Rat.minusOne) ? "-" : "%(coeff)")
            } else {
                System.write((coeff == Rat.one) ? " + " : (coeff == Rat.minusOne) ? " - " :
                             (coeff > Rat.zero) ? " + %(coeff)" : " - %(-coeff)")
            }
            var pwr = p + 1 - j
            System.write((pwr > 1) ? "n^%(pwr)" : "n")            
        }
    }
    System.print()
}

for (i in 0..9) faulhaber.call(i)

Output:

0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

zkl

Library: GMP

GNU Multiple Precision Arithmetic Library

Uses code from the Bernoulli numbers task (copied here).

var [const] BN=Import("zklBigNum");	// libGMP (GNU MP Bignum Library)

fcn faulhaberFormula(p){  //-->(Rational,Rational...)
   [p..0,-1].pump(List(),'wrap(k){ B(k)*BN(p+1).binomial(k) })
   .apply('*(Rational(1,p+1)))
}

foreach p in (10){ 
   println("F(%d) --> %s".fmt(p,polyRatString(faulhaberFormula(p))))
}

class Rational{  // Weenie Rational class, can handle BigInts
   fcn init(_a,_b){ var a=_a, b=_b; normalize(); }
   fcn toString{ 
      if(b==1) a.toString()
      else     "%d/%d".fmt(a,b) 
   }
   var [proxy] isZero=fcn{ a==0   };
   fcn normalize{  // divide a and b by gcd
      g:= a.gcd(b);
      a/=g; b/=g;
      if(b<0){ a=-a; b=-b; } // denominator > 0
      self
   }
   fcn __opAdd(n){
      if(Rational.isChildOf(n)) self(a*n.b + b*n.a, b*n.b); // Rat + Rat
      else self(b*n + a, b);				    // Rat + Int
   }
   fcn __opSub(n){ self(a*n.b - b*n.a, b*n.b) }		    // Rat - Rat
   fcn __opMul(n){
      if(Rational.isChildOf(n)) self(a*n.a, b*n.b);	    // Rat * Rat
      else self(a*n, b);				    // Rat * Int
   }
   fcn __opDiv(n){ self(a*n.b,b*n.a) }			    // Rat / Rat
}

fcn B(N){	// calculate Bernoulli(n) --> Rational
   var A=List.createLong(100,0);  // aka static aka not thread safe
   foreach m in (N+1){
      A[m]=Rational(BN(1),BN(m+1));
      foreach j in ([m..1, -1]){ A[j-1]= (A[j-1] - A[j])*j; }
   }
   A[0]
}
fcn polyRatString(terms){ // (a1,a2...)-->"a1n + a2n^2 ..."
   str:=[1..].zipWith('wrap(n,a){ if(a.isZero) "" else "+ %sn^%s ".fmt(a,n) },
        terms)
   .pump(String)
   .replace(" 1n"," n").replace("n^1 ","n ").replace("+ -","- ");
   if(not str)     return(" ");  // all zeros
   if(str[0]=="+") str[1,*];     // leave leading space
   else            String("-",str[2,*]);
}

Output:

F(0) -->  n 
F(1) -->  1/2n + 1/2n^2 
F(2) -->  1/6n + 1/2n^2 + 1/3n^3 
F(3) -->  1/4n^2 + 1/2n^3 + 1/4n^4 
F(4) --> -1/30n + 1/3n^3 + 1/2n^4 + 1/5n^5 
F(5) --> -1/12n^2 + 5/12n^4 + 1/2n^5 + 1/6n^6 
F(6) -->  1/42n - 1/6n^3 + 1/2n^5 + 1/2n^6 + 1/7n^7 
F(7) -->  1/12n^2 - 7/24n^4 + 7/12n^6 + 1/2n^7 + 1/8n^8 
F(8) --> -1/30n + 2/9n^3 - 7/15n^5 + 2/3n^7 + 1/2n^8 + 1/9n^9 
F(9) --> -3/20n^2 + 1/2n^4 - 7/10n^6 + 3/4n^8 + 1/2n^9 + 1/10n^10