Threaded Mode | Linear Mode

John Keith · 09-05-2023, 05:57 PM

(09-02-2023 05:44 PM)Albert Chan Wrote: We may be able to correct for tiny errors. B(m) largest denominator is 2*(2^m-1).

d = 2*(2^12-1) = 8190 (largest d)
d * -0.25312 = -2073.0528            → B(12) = -2073/8190 = -691/2730

Denominator likely smaller, but required work. d = product(p, for (p-1) | even m)
This is why d is divisible by 6. 2-1=1, 3-1=2, both divides even m.

d = 2*3*5*7*13 = 2730
d * -0.25312 = -691.0176             → B(12) = -691/2730

Amazingly, we get correct B(m = 12), without correction (algorithm based on above link)
If m divisible by 4, B(m) is negative. (positive otherwise)

lua> m = 12
lua> d = 2*(2^m-1)
lua> b = 4*d/(d+2) * fac(m)/pi^m
lua> 1+b, d
2073.4899880820685      8190
lua> _ + b*2^-m
2073.995967083065

--> B(12) = -2073/8190 = -691/2730

lua> m = 16
lua> d = 2*(2^m-1)
lua> b = 4*d/(d+2) * fac(m)/pi^m
lua> 1+b, d
929555.7943024995      131070
lua> _ + b*2^-m
929569.9781830278
lua> _ + b*3^-m
929569.9997771184

--> B(16) = -929569/131070 = -3617/510

Rounding errors might cause off-by-1 error with B(m) numerator.
It may be safer to drop numbers by 1/2, equivalent to rounding without 1+

Code:

gcd = require'factor'.gcd fac = fn'x: tgamma(x+1)' roughB = fn'm: I.real(-2*fac(m)/(2*pi*I)^m) * (1+2^-m+3^-m)' ratioB = fn'm,n,d: d=2*(2^m-1); n=round(roughB(m)*d); m=gcd(n,d); n/m, d/m'

This is quite interesting, I haven't seen this method before. Can it be modified to avoid complex numbers so that it can be used on older and simpler calculators?

Additionally, I have tested the method that I referred to above using zigzag (actually tangent) numbers, and the accuracy is far better than the method using Stirling numbers and factorials. On the HP-48 (12 digits) the results are correct to all 12 digits up to B(18) and to 11 or 12 digits up to B(40), which is as far as I tested. Programs are here. Also see this paper for further information and algorithms.