faulhaber_s_polynomials.sf

#!/usr/bin/ruby

# See also:
#   https://en.wikipedia.org/wiki/Faulhaber%27s_formula
#   https://rosettacode.org/wiki/Faulhaber%27s_formula

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 (0..20) {
    printf("%2d: %s\n", p, faulhaber_formula(p))
}

__END__
 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
10: 1/11*x^11 + 1/2*x^10 + 5/6*x^9 - x^7 + x^5 - 1/2*x^3 + 5/66*x
11: 1/12*x^12 + 1/2*x^11 + 11/12*x^10 - 11/8*x^8 + 11/6*x^6 - 11/8*x^4 + 5/12*x^2
12: 1/13*x^13 + 1/2*x^12 + x^11 - 11/6*x^9 + 22/7*x^7 - 33/10*x^5 + 5/3*x^3 - 691/2730*x
13: 1/14*x^14 + 1/2*x^13 + 13/12*x^12 - 143/60*x^10 + 143/28*x^8 - 143/20*x^6 + 65/12*x^4 - 691/420*x^2
14: 1/15*x^15 + 1/2*x^14 + 7/6*x^13 - 91/30*x^11 + 143/18*x^9 - 143/10*x^7 + 91/6*x^5 - 691/90*x^3 + 7/6*x
15: 1/16*x^16 + 1/2*x^15 + 5/4*x^14 - 91/24*x^12 + 143/12*x^10 - 429/16*x^8 + 455/12*x^6 - 691/24*x^4 + 35/4*x^2
16: 1/17*x^17 + 1/2*x^16 + 4/3*x^15 - 14/3*x^13 + 52/3*x^11 - 143/3*x^9 + 260/3*x^7 - 1382/15*x^5 + 140/3*x^3 - 3617/510*x
17: 1/18*x^18 + 1/2*x^17 + 17/12*x^16 - 17/3*x^14 + 221/9*x^12 - 2431/30*x^10 + 1105/6*x^8 - 11747/45*x^6 + 595/3*x^4 - 3617/60*x^2
18: 1/19*x^19 + 1/2*x^18 + 3/2*x^17 - 34/5*x^15 + 34*x^13 - 663/5*x^11 + 1105/3*x^9 - 23494/35*x^7 + 714*x^5 - 3617/10*x^3 + 43867/798*x
19: 1/20*x^20 + 1/2*x^19 + 19/12*x^18 - 323/40*x^16 + 323/7*x^14 - 4199/20*x^12 + 4199/6*x^10 - 223193/140*x^8 + 2261*x^6 - 68723/40*x^4 + 43867/84*x^2
20: 1/21*x^21 + 1/2*x^20 + 5/3*x^19 - 19/2*x^17 + 1292/21*x^15 - 323*x^13 + 41990/33*x^11 - 223193/63*x^9 + 6460*x^7 - 68723/10*x^5 + 219335/63*x^3 - 174611/330*x