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

Haskell[edit]

Bernouilli polynomials[edit]

import Data.Ratio ((%), numerator, denominator)
import Data.List (intercalate, transpose)
import Control.Arrow ((&&&), (***))
import Data.Char (isSpace)
import Data.Monoid ((<>))
 
-- FAULHABER -------------------------------------------------------------------
faulhaber :: [[Rational]]
faulhaber =
tail $
scanl
(\rs n ->
let xs = zipWith ((*) . (n %)) [2 ..] rs
in 1 - sum xs : xs)
[]
[0 ..]
 
-- EXPRESSION STRINGS ----------------------------------------------------------
polynomials :: [[(String, String)]]
polynomials = fmap ((ratioPower =<<) . reverse . flip zip [1 ..]) faulhaber
 
-- 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 (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 if null s
then []
else [(sn, sr)]
 
-- 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
, if blank a || head a == '-'
then " "
else " + "
, 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 case l of
0 -> 0
_ ->
case head xs of
'-' -> l - 1
_ -> l
 
-- TEST ------------------------------------------------------------------------
main :: IO ()
main = (putStrLn . unlines . expressionTable . take 10) polynomials
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[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

Kotlin[edit]

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

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