15 Jan
2012
15 Jan
'12
1:37 a.m.
On Thu, Jan 12, 2012 at 1:49 PM, Bill Gosper <billgosper@gmail.com> wrote:
> 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
>
> We did. Here's a 2x2 giving psi[0](x+1)+gamma at 6 bits/whack, at a cost
slightly higher than
twice the two bits/whack. Product k>=0
{{((k + 1)*(-x + k + 1)*(x + k + 1)^2)/(64*(k + 3/2)^2*(x/2 + k + 2)*((x +
3)/2 + k)),
(21*x*((x^2/7 + (2*k + 10/7)*x + (102*k + 39)/21/3) + k^2))/(4*(x +
1)*(x + 2))},
{0, 1}}
(Upper right element).
Equivalent sum:
(21*Sum[((-1)^k*
k!*(((3*x^2 + 30*x + 102*k + 39)/21 + 2*k*x)/3 + k^2)*(x + k)!^2)/
(2^(4*(k - 1))*(k + 1/2)!^2*(x/2 + k + 1)*(x + 2*(k - 1) + 3)!*(x - k
- 1)!),
{k, 0, Infinity}]*Pi)/512
This can be expressed as a sum of three 5F4s or one 7f6 with square root
parameters.
Mma gives NINE 5f4s and 6f5s!
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.
>
> http://www.tweedledum.com/rwg/munafaccf.png
>
>
>> 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?)
>
Robert explained it very clearly. And w/o a spoiler warning!
--rwg
> (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
>
>
>
>