Faulhaber's formula

From Rosetta Code
Faulhaber's formula is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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.

See also


EchoLisp[edit]

 
(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


GAP[edit]

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

J[edit]

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

Maxima[edit]

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

PARI/GP[edit]

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.[edit]

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.[edit]

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

Perl[edit]

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

Perl 6[edit]

Works with: rakudo version 2015.12
use experimental :cached;
 
sub bernoulli_number($n) is cached {
 
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) is cached {
$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)

Python[edit]

The following implementation does not use Bernoulli numbers, but Stirling numbers of the second kind, based on the relation: .

Then summing: .

One has then to expand the product in order to get a polynomial in the variable . Also, for the sum of , the sum is too large by one (since we start at ), 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 fronts) be visible to Firefox.


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[edit]

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

Sidef[edit]

func faulhaber_s_formula(p) {
 
var formula = gather {
{ |j|
take "(#{binomial(p+1, j) * j.bernfrac -> as_rat})*n^#{p+1 - j}"
} << 0..p
}
 
formula.grep! { !.contains('(0)*') }.join!(' + ')
 
formula -= /\(1\)\*/g
formula -= /\^1\b/g
formula.gsub!(/\(([^+]*?)\)/, { _ })
 
"1/#{p + 1} * (#{formula})"
}
 
{ |p|
printf("%2d: %s\n", p, faulhaber_s_formula(p))
} << ^10
Output:
 0: 1/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)

By not simplifying the formulas, we can have a much cleaner code:

func faulhaber_s_formula(p) {
"1/#{p + 1} * (" + gather {
{ |j|
take "#{binomial(p+1, j) * j.bernfrac -> as_rat}*n^#{p+1 - j}"
} << 0..p
}.join(' + ') + ")"
}
 
{ |p|
printf("%2d: %s\n", p, faulhaber_s_formula(p))
} << ^10
Output:
 0: 1/1 * (1*n^1)
 1: 1/2 * (1*n^2 + 1*n^1)
 2: 1/3 * (1*n^3 + 3/2*n^2 + 1/2*n^1)
 3: 1/4 * (1*n^4 + 2*n^3 + 1*n^2 + 0*n^1)
 4: 1/5 * (1*n^5 + 5/2*n^4 + 5/3*n^3 + 0*n^2 + -1/6*n^1)
 5: 1/6 * (1*n^6 + 3*n^5 + 5/2*n^4 + 0*n^3 + -1/2*n^2 + 0*n^1)
 6: 1/7 * (1*n^7 + 7/2*n^6 + 7/2*n^5 + 0*n^4 + -7/6*n^3 + 0*n^2 + 1/6*n^1)
 7: 1/8 * (1*n^8 + 4*n^7 + 14/3*n^6 + 0*n^5 + -7/3*n^4 + 0*n^3 + 2/3*n^2 + 0*n^1)
 8: 1/9 * (1*n^9 + 9/2*n^8 + 6*n^7 + 0*n^6 + -21/5*n^5 + 0*n^4 + 2*n^3 + 0*n^2 + -3/10*n^1)
 9: 1/10 * (1*n^10 + 5*n^9 + 15/2*n^8 + 0*n^7 + -7*n^6 + 0*n^5 + 5*n^4 + 0*n^3 + -3/2*n^2 + 0*n^1)

zkl[edit]

Uses GMP (Gnu Multi Precision library) and 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