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

   Bernoulli numbers.
   evaluate binomial coefficients.


See also

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


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