Faulhaber's formula
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
GAP
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.
<lang gap>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</lang>
Haskell
Bernouilli polynomials
<lang Haskell>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</lang>
- 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
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>
Kotlin
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. <lang scala>// 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)
}</lang>
- 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
<lang maxima>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)));</lang>
- 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
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).
Version #1. Using Bernoulli numbers.
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.
<lang parigp>
\\ 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))} </lang>
- 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.
This version is using the sums of pth powers formula from Bernoulli polynomials. It has small, simple and clear code, and produces instant result. <lang parigp> \\ 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))}
</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 > ## *** last result computed in 0 ms
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
<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: .
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.
<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 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))</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 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</lang>
- 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: <lang ruby>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</lang>
- 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
Uses GMP (Gnu Multi Precision library) and code from the Bernoulli numbers task (copied here). <lang zkl>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)))
}</lang> <lang zkl>foreach p in (10){
println("F(%d) --> %s".fmt(p,polyRatString(faulhaberFormula(p))))
}</lang> <lang zkl>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
}</lang> <lang zkl>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,*]);
}</lang>
- 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