Faulhaber's formula: Difference between revisions
Line 222: | Line 222: | ||
</pre> |
</pre> |
||
=={{header|Python}}== |
== {{header|Python}} == |
||
<lang python>from fractions import Fraction |
<lang python>from fractions import Fraction |
||
Line 240: | Line 240: | ||
return b |
return b |
||
def |
def sumpol(n): |
||
u = [0, 1] |
u = [0, 1] |
||
v = [[1], [1, 1]] |
v = [[1], [1, 1]] |
||
yield [Fraction(0), Fraction(1)] |
|||
for i in range(1, n): |
for i in range(1, n): |
||
v.append(nextv(v[-1])) |
v.append(nextv(v[-1])) |
||
Line 253: | Line 253: | ||
for k, s in enumerate(v[j + 1]): |
for k, s in enumerate(v[j + 1]): |
||
t[k] += r * s |
t[k] += r * s |
||
yield t |
|||
u = nextu(u) |
u = nextu(u) |
||
return w |
|||
def polstr(a): |
def polstr(a): |
||
Line 287: | Line 286: | ||
return "0" |
return "0" |
||
for i, p in enumerate( |
for i, p in enumerate(sumpol(10)): |
||
print(i, ":", polstr(p)) |
print(i, ":", polstr(p)) |
||
Line 300: | Line 299: | ||
8 : 1/9 n^9 + 1/2 n^8 + 2/3 n^7 - 7/15 n^5 + 2/9 n^3 - 1/30 n |
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</lang> |
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</lang> |
||
=={{header|Racket}}== |
=={{header|Racket}}== |
||
Revision as of 16:19, 5 May 2016
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
<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
<lang python>from fractions import Fraction
def nextu(a):
n = len(a) a.append(a[-1] * n) 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): p *= j + 1 r = Fraction(r, p) 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 q: if x > 0: s += " + %s%s" % (c, m) elif x < 0: s += " - %s%s" % (c, m) elif x > 0: s = "%s%s" % (c, m) q = True elif x < 0: s = "-%s%s" % (c, m) q = True if q: return s else: return "0"
for i, p in enumerate(sumpol(10)):
print(i, ":", polstr(p))
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</lang>
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)