(28 48 49 50) Bernoulli Numbers

[Albert Chan]

Here is lua function for zigzag numbers, build 1 zigzag at a time

Code:

function zigzag(t)      -- Euler zigzag numbers
    t = t or {[0] = 1}  -- t[0], t[-0] = T(0), T(1)
    local i = #t + 1
    t[i] = 0; t[-i] = 0 -- make room
    for j=2-i, i    do t[j]=t[j]+t[j-1] end -- zig
    for j=i-1,-i,-1 do t[j]=t[j]+t[j+1] end -- zag
    return t            -- t[i], t[-i] = T(next even), T(next odd)
end

lua> m = 16
lua> L = m/2-1          -- 1(excluded) to m-1 step 2, zigzags to reach = ((m-1)-1)/2 = m/2-1
lua> t = {[0] = 1}     -- t[0], t[-0] = T(0), T(1)
lua> for k=1,L do t=zigzag(t) end
lua> t[-#t]                -- = T(m-1)
1903757312

lua> p = 2^m
lua> n = (-1)^L * m * t[-#t] * 2/p
lua> d = 2*(p-1)       -- B(m) denominator = 2*odd
lua> n, d, n/d
-929569      131070      -7.092156862745098

---

Array t can be extended, since biggest element, T(m-1), does not overflow yet.

lua> m = 100
lua> t[-#t] = require'mapm'.new(t[-#t]) -- pollute array with mapm number
lua> L2 = m/2-1
lua> for k=1,L2-L do t=zigzag(t) end
lua> t[-#t] -- = T(m-1)
4.560851661680111182104382953145169718558E+136

-- Note: mapm +/-/* is exact, above holds *exact* T(m-1)
-- mapm.gcd use efficient binary gcd algorithm. I don't bother with denominator = 2*odd

lua> ;mapm.places(136) -- T(m-1) exponent more than enough precision
lua> n = (-1)^L2 * m * t[-#t]
lua> d = mapm.pow(2, m)
lua> d = d*(d-1)
lua> g = mapm.gcd(n, d)
lua> n, d = n/g, d/g
lua> n:tofixed(0) .. ' / ' .. d:tofixed(0)      -- B(m) reduced fraction

-94598037819122125295227433069493721872702841533066936133385696204311395415197247​711 / 33330