On Tue, Jan 10, 2012 at 5:27 PM, Bill Gosper <billgosper@gmail.com> wrote:
wds>
http://www.luschny.de/math/factorial/approx/SimpleCases.html
Luschny describes efforts like Gosper's latest as "pimp my formula" :) ------- Ok, here's 34.6 digits on z>12, nine BFs.
OOPS, ten! But morally less--a1 and b1 need only be 29D,...,a5,b5 only 21D.
http://gosper.org/munafacmu34.png I'll bet we can do better with a different shape to the rational function.
--rwg
I.e., there are >= 19 other ten parameter rational functions to try, some of which are probably superior.
Here's something I just sent to Robert: On Thu, Jan 12, 2012 at 2:11 AM, Bill Gosper <billgosper@gmail.com> wrote:
Hi Robert. Here's a 2 bit/whack matrix product for digamma. You obviously don't need matrices nor even vectors to use it, especially if you run k backwards.
(c81) 'PRODUCT(MATRIX([(-X+K+1)/(2*(2*K+3)),X*(X+3*K+3)/(2*(K+1)*(X+K+1))],[0,1]),K,0,INF) = MATRIX([0,RHS(D73)],[0,1]);
inf /===\ [ - x + k + 1 x (x + 3 k + 3) ] | | [ ----------- --------------------- ] [ 0 psi (x + 1) + %gamma ] (d81) | | [ 2 (2 k + 3) 2 (k + 1) (x + k + 1) ] = [ 0 ] | | [ ] [ ] k = 0 [ 0 1 ] [ 0 1 ]
(c82) SUBST([X = 1/2,NOUNIFY(PRODUCT) = LAMBDA([XP,V,LO,HI],PRUD(XP,V,LO,29))],%);
[ 1 5228238559942243328529758531923102907210593 ] [ -------------------- ------------------------------------------- ] (d82) [ 70328211781017665536 8519130717916716392166434494153620203765760 ] = [ ] [ 0 1 ]
[ 0 2 - 2 log(2) ]
[ ] [ 0 1 ]
(c83) DFLOAT(%);
[ 1.42190448850558d-20 0.61370563888011d0 ] (d83) [ ] = [ 0.0d0 1.0d0 ]
[ 0.0d0 0.61370563888011d0 ]
[ ] [ 0.0d0 1.0d0 ]
I can do a lot better than 2 bits/whack, but the whacks get more expensive. It may be like grouping terms pairwise and claiming twice the convergence. Julian expressed interest in reviewing matrix path invariance, so maybe we can investigate this Saturday. --Bill
On Wed, Jan 11, 2012 at 10:37 PM, Bill Gosper <billgosper@gmail.com>wrote:
I was lengthening the CF a term at a time, heading toward another ten parameter, quad precision Remez, when I got stuck--FindRoot can't find a6! I guess it can't stand the thought of making that ridiculous plot label any taller.
Robert, I think I can simply provide a decently convergent series for
digamma. Will that do? --Bill
This is interesting. At two bits/whack, you need about 52 whacks for quad: (c96) BFLOAT(SUBST([X = 1/2,NOUNIFY(PRODUCT) = LAMBDA([XP,V,LO,HI],PRUD(XP,V, LO,52))],'PRODUCT(MATRIX([(-X+K+1)/(2*(2*K+3)),X*(X+3*K+3)/(2*(K+1)*(X+K+1))],[0,1]),K,0,INF) = MATRIX([0,PSI[0](X+1)+%GAMMA],[0,1])),35); (d96) matrix([1.15195809757741209902413633949005b-34, 6.1370563888010938116553575708364686b-1], [0.0b0, 1.0b0]) = [ 0.0b0 6.1370563888010938116553575708364687b-1 ] [ ] [ 0.0b0 1.0b0 ] But offsetting the arg to be roughly 52/2 (cost: ~26 adds) reduces the whacks to 37: (c97) BFLOAT(SUBST([X = 26+1/2,NOUNIFY(PRODUCT) = LAMBDA([XP,V,LO,HI],PRUD(XP,V,LO,37))],'PRODUCT(MATRIX([(-X+K+1)/(2*(2*K+3)),X*(X+3*K+3)/(2*(K+1)*(X+K+1))],[0,1]),K,0,INF) = MATRIX([0,PSI[0](X+1)+%GAMMA],[0,1])),35); (d97) matrix([1.2351588921475039695714809452331839b-35, 3.8731096731165662046581399543667894b0], [0.0b0, 1.0b0]) = [ 0.0b0 3.8731096731165662046581399543667894b0 ] [ ] [ 0.0b0 1.0b0 ] (OK kids, why?) (I believe this zmod1 = 1/2 is worst case.) A completely other approach would be the log-derivative of the ten parameter factorial scheme: psi0[z+1] ~ (g[z] + 810*z^5*(1 + z - Coth[1/z] + z*Log[z^3*Sinh[1/z]]))/(1620*z^6) with g[z] ~z f'[z]-5f[z], where f[z] was the rational Remez fudgefunction. The ten terms plus 12 divides may sound a lot cheaper than 37 whacks plus 26 adds, but those are cheap whacks. --Bill