1234567891011121314151617181920212223242526272829303132333435363738394041 |
- #!/usr/bin/ruby
- # The formula for calculating the sum of consecutive
- # numbers raised to a given power, such as:
- # 1^p + 2^p + 3^p + ... + n^p
- # where p is a positive integer.
- # See also:
- # https://en.wikipedia.org/wiki/Faulhaber%27s_formula
- # The Faulhaber's formula
- # See: https://en.wikipedia.org/wiki/Faulhaber%27s_formula
- func faulhaber_formula (n, p) {
- sum(0..p, {|k|
- binomial(p + 1, k) * bernoulli(k) * n**(p + 1 - k)
- }) / (p+1)
- }
- # Alternate expression using Bernoulli polynomials
- # See: https://en.wikipedia.org/wiki/Faulhaber%27s_formula#Alternate_expressions
- func bernoulli_polynomials (n, x) {
- sum(0..n, {|k|
- binomial(n, k) * bernoulli(n - k) * x**k
- })
- }
- func faulhaber_formula_2 (n, p) {
- (bernoulli_polynomials(p + 1, n) - bernoulli(p + 1)) / (p + 1)
- }
- for p in (0..20) {
- var n = 100.irand
- var r = faulhaber_formula(n, p)
- assert_eq(r, faulhaber_formula_2(n, p))
- assert_eq(r, n.faulhaber_sum(p)) # built-in
- printf("faulhaber_#{p}(%2d) = %s\n", n, r)
- }
|