[math-fun] Peter Luschny's gamma function approximation-formulas collection
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. http://gosper.org/munafacmu34.png I'll bet we can do better with a different shape to the rational function. --rwg
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
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
>
>
>
>
Alternatively, here's a straight 18 parameter quad z!<http://www.tweedledum.com/rwg/farquad18.png>on |z|<1/4, which suffices. --rwg 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. http://gosper.org/munafacmu34.png I'll bet we can do better with a different shape to the rational function.
--rwg
On Sun, Jan 15, 2012 at 3:59 PM, Bill Gosper <billgosper@gmail.com> wrote:
Alternatively, here's a straight 18 parameter quad z!<http://www.tweedledum.com/rwg/farquad18.png>on |z|<1/4, which suffices. --rwg
ACK, no it doesn't! Nevermind. --rwg
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. http://gosper.org/munafacmu34.png I'll bet we can do better with a different shape to the rational function.
--rwg
On Sun, Jan 15, 2012 at 4:06 PM, Bill Gosper <billgosper@gmail.com> wrote:
On Sun, Jan 15, 2012 at 3:59 PM, Bill Gosper <billgosper@gmail.com> wrote:
Alternatively, here's a straight 18 parameter quad z!<http://www.tweedledum.com/rwg/farquad18.png>on |z|<1/4, which suffices. --rwg
ACK, no it doesn't! Nevermind. --rwg
But it does when combined with this 15 parameter quad z!<http://www.tweedledum.com/rwg/farquad2.png>on 1/4 < z < 3/4. --rwg (All these coefficients could stand another Remez iteration or so, but it's unlikely to affect even the lsb of the outputs.)
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. http://gosper.org/munafacmu34.png I'll bet we can do better with a different shape to the rational function.
--rwg
participants (1)
-
Bill Gosper