Faulhaber's formula

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

• Bernoulli numbers and binomial coefficients on Wikipedia.

EchoLisp

<lang scheme> (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) ))) </lang>

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

J

Implementation:

<lang J>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.

)</lang>

Task example:

<lang J> 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.</lang>

Double checking our work:

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

Perl

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

}</lang>

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

Works with: rakudo version 2015.12

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

}</lang>

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

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)_{k+1}$ in order to get a polynomial in the variable $m+1$. 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.

<lang python>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 c == "1" and i > 0:
           c = ""
           m = m.strip()
       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))</lang>

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

</lang>

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

<lang ruby>func bernoulli_number((1)) { 1/2 } func bernoulli_number({.is_odd}) { 0/1 }

func bernoulli_number(n) is cached {

   var a = [];
   0.to(n).each { |m|
       a[m] = 1/(m + 1)
       m.downto(1).each { |j|
           a[j-1] = j*(a[j-1] - a[j])
       }
   }

   return a[0]

}

func faulhaber_s_formula(p) {

   var formula = gather {
       0.to(p).each { |j|
           take "(#{binomial(p+1, j) * bernoulli_number(j) -> as_rat})*n^#{p+1 - j}"
       }
   }

   formula.grep! { !.contains('(0)*') }.join!(' + ')

   formula -= /\(1\)\*/g
   formula -= /\^1\b/g

   "1/#{p + 1} * (#{formula})"

}

0.to(9).each { |p|

   printf("%2d: %s\n", p, faulhaber_s_formula(p))

}</lang>

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