(12C+) Bernoulli Number

[Albert Chan]

(07-28-2019 12:02 AM)Albert Chan Wrote:  1 64 729 4096 15625 46656 117649 // value of 1^6 to 7^6
63 665 3367 11529 31031 70993     // forward differences
602 2702 8162 19502 39962
2100 5460 11340 20460
3360 5880 9120
2520 3240
720

Sum of k^6 formula = \(1\binom{n}{1}+63\binom{n}{2}+602\binom{n}{3}+2100\binom{n}{4}+3360\binom{n}{5}​+2520\binom{n}{6}+720\binom{n}{7}\)

B(6) = Linear term coefficient = 1/1 - 63/2 + 602/3 - 2100/4 + 3360/5 - 2520/6 + 720/7 = 1/42

Note: above Σk^6 formula is from 0 to n = B(7,n)/7 + n^6

We may use symmetry, and do 0^6 to 3^6, then flip it, to get values (-3)^6 .. (+3)^6
This keep table numbers smaller, thus smaller forward difference numbers:

\(x^6 =
729\binom{x+3}{0}
-665\binom{x+3}{1}
+602\binom{x+3}{2}
-540\binom{x+3}{3}
+480\binom{x+3}{4}
-360\binom{x+3}{5}
+720\binom{x+3}{6}\)

Integrate (actually sum, but look like integration)

\(\large \frac{B^7}{7} \normalsize =
729\binom{x+3}{1}
-665\binom{x+3}{2}
+602\binom{x+3}{3}
-540\binom{x+3}{4}
+480\binom{x+3}{5}
-360\binom{x+3}{6}
+720\binom{x+3}{7} \) + C

Differentiate, evaluated at x=0, we get B(6) = RHS constant term.

First half of terms is not as easy to evaluate, because there was no factor x
(with a factor x, derivative at x=0 make the term goes away!)

But, it may still worth it, to keep forward difference table entries small.
Note: for B(m), n = floor(m/2), k = 1 .. n

\(\displaystyle \binom{x+n}{k} =
\left(\frac{x+n}{k}\right) \binom{x+n-1}{k-1}
\)

\(\displaystyle \binom{x+n}{k}' \bigg\rvert_{x=0} =
\left(\frac{1}{k}\right) \binom{n-1}{k-1} +
\left(\frac{n}{k}\right) \binom{x+n-1}{k-1}'\bigg\rvert_{x=0}
\)

I had managed to undo recursion, into iterative sums (written prove welcome!)

\(\displaystyle \binom{x+n}{k}' \bigg\rvert_{x=0}
\;=\; \sum_{j=1}^{k} \frac{(-1)^{j+1}}{j}\; \binom{n}{k-j}
\;=\; \binom{n}{k}\;\sum_{j=n-(k-1)}^{n} \frac{1}{j}
\)

Note that coefficient formula good for half terms (the other half is easy, see code)

Code:

from gmpy2 import mpq, mpz, xmpz
from itertools import izip, count

def bernoulli(m):
    if m == 0: return 1
    lst = [0, 1] + [ mpz(i)**m for i in range(2, (m+3)//2) ]
    lst = [(-x if m&1 else x) for x in lst[-1-(m&1):0:-1]] + lst
    lst = map(xmpz, lst)        # speedup fdiff_inplace
    fdiff_inplace(lst)          # forward diff table
    return sum(c*d for (c,d) in izip(bernoulli_coef(m//2), lst))

def fdiff_inplace(f):
    n = len(f) - 1
    for i in xrange(n):         # n forward diff(s)
        for j in xrange(n,i,-1):
            f[j] -= f[j-1]

def bernoulli_coef(n, mpq=mpq):
    'diff(comb(x+n,k)), k=1,2,3,..., evaluate at x=0'
    c, s, N = mpz(1), 0, n+1
    for j in xrange(n, 0, -1):  # tricky half
        c = c * j // (N-j)
        s += mpq(1, j)
        yield c * s
    s = mpq(1, N)
    for c in count(1):          # easy half
        yield s
        s *= mpq(-c, N+c)

I did not optimize away B(odd m), just to see if it matches ... it does! B(1) = -1/2, B(odd>1) = 0

>>> for m in range(2,31,2): print m, bernoulli(m)
...
2      1/6
4      -1/30
6      1/42
8      -1/30
10     5/66
12     -691/2730
14     7/6
16     -3617/510
18     43867/798
20     -174611/330
22     854513/138
24     -236364091/2730
26     8553103/6
28     -23749461029L/870
30     8615841276005L/14322