(28 48 49 50) Bernoulli Numbers

[Albert Chan]

Hi, John Keith

I was also surprised that O(m) to O(√m) does not translate to much speed gain.
Yes, O(√m) version eventually is faster, but not by much.

The reason is D = 2*(2^m-1) factors between 2 and m+1 matched d factors, with few false positives.
This make O(m) code very efficient. (it skipped D nonfactors)

m = 2 to 52 step 2, there are only 3 m's that have false positives. (bold)

m = 24:
d = 2*3*5*7*13
D = 2*3²*5*7*13*17*241

m = 40:
d = 2*3*5*11*41
D = 2*3*5²*11*17*31*41*61681

m = 50:
d = 2*3*11
D = 2*3*11*31*251*601*1801*4051

---

Instead of half q's, we may use a few q's, to reduce search space.
With 1 q, it already beat O(√m) version, all the way to m=52 (IEEE double limit)

Code:

function min_d(m) -- B(even m>=2) reduced denominator
    local n, d = 2^m-1, 1
    for p = 3, m/2+1, 2 do            
        if n%p ~= 0 then continue end
        repeat n = n/p until n%p ~= 0
        if m%(p-1)==0 then d = d*p end
    end
    if n%(m+1)==0 then d = d*(m+1) end -- q for p=2
    return 2*d
end