[seqfan] Re: A063747 sign of power series expansion of the Gamma function
Brendan McKay
Brendan.McKay at anu.edu.au
Sat May 20 09:30:55 CEST 2023
Concerning rounding error, the following is not rigorous but could be
made rigorous.
Since all the numbers involved in the calculation of the coefficients of
Gamma(1-x) are less than 2 in absolute value, we can upper-bound the
floating error by the fixed-point error. For some epsilon round about
10^(number of digits) the fixed-point (absolute, not relative) error in
a calculation is about at most N*epsilon where N is the number of
arithmetic operations. This says that the absolute error in the
coefficients of Gamma(1-x) is at most about 3*n^2*epsilon.
In calculating the coefficients of 1/Gamma(1-x) there is another factor
of n^2 due to the recursion. However, there is also the possibility
that the intermediate values in the calculation of b[*] from a[*] blow
up even though the final answer is small. (Actually it doesn't happen,
but..) To compensate we can multiply by an additional maximum over i of
SUM[ a[j] |b[i-j]|, j=1..i ] / SUM[ a[j] b[i-j], j=1..i ] which is about
3e19.
So the absolute error in the coefficients of 1/Gamma(1-x) can't be more
than about 10^20*n^4*epsilon.
For n=100 and 200 digits, this gives a maximum absolute error of
10^(-172). Since none of the calculated values are less than 10^(-108),
this shows that the signs are all correct with a large margin of
safety. As I said, this is a crude calculation, but a more careful
calculation won't be much different.
Here are the first 201 numbers:
[1, -1, -1, 1, 1, 1, -1, -1, -1, 1,
1, 1, -1, -1, -1, -1, 1, 1, 1, -1,
-1, -1, -1, 1, 1, 1, 1, -1, -1, -1,
1, 1, 1, 1, -1, -1, -1, -1, 1, 1,
1, 1, 1, -1, -1, -1, -1, 1, 1, 1,
1, -1, -1, -1, -1, 1, 1, 1, 1, 1,
-1, -1, -1, -1, 1, 1, 1, 1, -1, -1,
-1, -1, -1, 1, 1, 1, 1, 1, -1, -1,
-1, -1, 1, 1, 1, 1, 1, -1, -1, -1,
-1, 1, 1, 1, 1, 1, -1, -1, -1, -1,
-1, 1, 1, 1, 1, 1, -1, -1, -1, -1,
-1, 1, 1, 1, 1, -1, -1, -1, -1, -1,
1, 1, 1, 1, 1, -1, -1, -1, -1, -1,
1, 1, 1, 1, 1, -1, -1, -1, -1, -1,
1, 1, 1, 1, 1, -1, -1, -1, -1, -1,
1, 1, 1, 1, 1, -1, -1, -1, -1, -1,
1, 1, 1, 1, 1, 1, -1, -1, -1, -1,
-1, 1, 1, 1, 1, 1, -1, -1, -1, -1,
-1, 1, 1, 1, 1, 1, -1, -1, -1, -1,
-1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1]
Brendan.
On 20/5/2023 4:43 pm, Brendan McKay via SeqFan wrote:
> Here is a quick way to calculate these numbers.
>
> First we find the Taylor expansion Gamma(1-x) = SUM[ a[i] x^i,
> i=0..infinity ].
>
> Define g(1,k) = gamma^k/k! and g(m,k) = Zeta(m)^k/(k! m^k) for m>1,
> where gamma is the Euler constant and Zeta is the Riemann zeta function.
>
> The coefficient a[i] equals f(i,i) where f is defined recursively:
> f(0,0)=1, f(n,0)=0 if n>0, f(n,m)=SUM[ g(m,k)*f(n-km,m-1),
> k=0..floor(n/m) ] .
>
> Both g and f should be memoised in floating point. There is no
> numerical issue from cancellation since all the terms are positive.
> The coefficients a[i] rapidly approach 1.
>
> Finally, take the reciprocal 1/Gamma(1-x) = SUM[ b[i] x^i,
> i=0..infinity ]. Clearly b[0]=1 and b[i]=-SUM[ a[j] b[i-j], j=1..i ]
> for i>0. In this case there is cancellation and some analysis is
> needed to be sure that the correct sign is achieved.
>
> Maple takes 1.7 seconds to get 101 terms using 200-digit arithmetic
> and 18 seconds to repeat the calculation using 3000-digit arithmetic.
> The answer is the same:
>
> [1, -1, -1, 1, 1, 1, -1, -1, -1, 1,
> 1, 1, -1, -1, -1, -1, 1, 1, 1, -1,
> -1, -1, -1, 1, 1, 1, 1, -1, -1, -1,
> 1, 1, 1, 1, -1, -1, -1, -1, 1, 1,
> 1, 1, 1, -1, -1, -1, -1, 1, 1, 1,
> 1, -1, -1, -1, -1, 1, 1, 1, 1, 1,
> -1, -1, -1, -1, 1, 1, 1, 1, -1, -1,
> -1, -1, -1, 1, 1, 1, 1, 1, -1, -1,
> -1, -1, 1, 1, 1, 1, 1, -1, -1, -1,
> -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1]
>
> The coefficients b[i] range in absolute value from 1 down to about
> 9e-108. I'm pretty sure that this combined with the fact that a[i] is
> always between 1/2 and 1 will allow a rigorous bound on the floating
> point error. See my next email, probably.
>
> Brendan.
>
> On 20/5/2023 6:59 am, Sean A. Irvine wrote:
>> Agreed.
>>
>> I'm not sure what our usual conventions are with regard to Pari programs
>> where precision is an issue.
>> It would seem to me that adding a "\p 1000" line to the program would at
>> least make it clear that choosing an appropriate precision is important
>> here.
>>
>> Sean.
>>
>>
>> On Sat, 20 May 2023 at 08:52, Hugo Pfoertner<yae9911 at gmail.com> wrote:
>>
>>> Yes, this is a delicate matter. We sometimes have to rely on the
>>> users of
>>> the programs to choose a level of accuracy appropriate to the
>>> problem. It
>>> is also sometimes not obvious that floating-point arithmetic is
>>> being used
>>> in the background. And then one can hardly assume that everyone has
>>> half a
>>> century of experience in numerical computing.
>>>
>>> On Fri, May 19, 2023 at 10:33 PM Sean A. Irvine<sairvin at gmail.com>
>>> wrote:
>>>
>>>> I agree with those values.
>>>> But (at least for me), if you run Pari without the "\p" you get the
>>>> wrong
>>>> result.
>>>>
>>>> Sean.
>>>>
>>>>
>>>> On Sat, 20 May 2023 at 08:28, Hugo Pfoertner<yae9911 at gmail.com>
>>>> wrote:
>>>>
>>>>> With my(x='x+O('x^150)); apply(sign, Vec(1/(gamma(1-x)))) and \p1000
>>> and
>>>>> \p10000 I get identical results for the first 100 terms.
>>>>> [1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1,
>>>>> -1,
>>>> -1,
>>>>> -1, -1, 1, 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1,
>>>>> 1, 1,
>>>> 1,
>>>>> -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1,
>>>>> -1,
>>>> -1,
>>>>> 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1,
>>>>> 1, 1,
>>>> 1,
>>>>> 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1] wrong?
>>>>>
>>>>> On Fri, May 19, 2023 at 10:19 PM Sean A. Irvine<sairvin at gmail.com>
>>>> wrote:
>>>>>> Because it is still wrong unless you set the precision higher.
>>>>>>
>>>>>> On Sat, 20 May 2023 at 08:10, Hugo Pfoertner<yae9911 at gmail.com>
>>>> wrote:
>>>>>>> Why not use the result of Michel Marcus's new PARI program? Is
>>> there
>>>>> any
>>>>>>> doubt that PARI can calculate e.g. 100 terms correctly? Then we
>>> have
>>>> a
>>>>>> very
>>>>>>> fundamental problem.
>>>>>>>
>>>>>>>
>>>>>>> On Fri, May 19, 2023 at 9:34 PM Sean A. Irvine<sairvin at gmail.com>
>>>>>> wrote:
>>>>>>>> Yes, I will make an update.
>>>>>>>>
>>>>>>>> Thank you Jean-François and Brendan for taking the time to
>>> verify.
>>>>>>>> Sean.
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> On Sat, 20 May 2023 at 07:04, Neil Sloane<njasloane at gmail.com>
>>>>> wrote:
>>>>>>>>> Well, if a(46) is wrong (which seems to be the case), could
>>>> someone
>>>>>>>> please
>>>>>>>>> make the correction? And probably all the terms after a(46)
>>>> should
>>>>>> be
>>>>>>>>> deleted, and we should give it keyword "more".
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> Best regards
>>>>>>>>> Neil
>>>>>>>>>
>>>>>>>>> Neil J. A. Sloane, Chairman, OEIS Foundation.
>>>>>>>>> Also Visiting Scientist, Math. Dept., Rutgers University,
>>>>>>>>> Email:njasloane at gmail.com
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> On Fri, May 19, 2023 at 3:14 AM Jean-François Alcover <
>>>>>>>>> jf.alcover at gmail.com>
>>>>>>>>> wrote:
>>>>>>>>>
>>>>>>>>>> Mma's SeriesCoefficient[1/Gamma[1 - x], {x, 0, n}] // Sign
>>>> gives
>>>>> -1
>>>>>>> for
>>>>>>>>>> n=46 (after 15 mn !)
>>>>>>>>>>
>>>>>>>>>> Le ven. 19 mai 2023 à 06:52, Brendan McKay via SeqFan <
>>>>>>>>>> seqfan at list.seqfan.eu>
>>>>>>>>>> a écrit :
>>>>>>>>>>
>>>>>>>>>>> My typo, it is e-38. Oops. B/
>>>>>>>>>>>
>>>>>>>>>>> On 19/5/2023 1:11 pm, Allan Wechsler wrote:
>>>>>>>>>>>> One of you said "e-38" and the other said "e-58". Do you
>>>>> really
>>>>>>>> have
>>>>>>>>>>>> 20 orders of magnitude of disagreement?
>>>>>>>>>>>>
>>>>>>>>>>>> On Thu, May 18, 2023 at 11:07 PM Brendan McKay via SeqFan
>>>>>>>>>>>> <seqfan at list.seqfan.eu> wrote:
>>>>>>>>>>>>
>>>>>>>>>>>> I also get -0.444..e-58 * x^46 in Maple using 60
>>>> digits.
>>>>>>> The
>>>>>>>>>>>> calculation
>>>>>>>>>>>> is numerically unstable though.
>>>>>>>>>>>>
>>>>>>>>>>>> B/
>>>>>>>>>>>>
>>>>>>>>>>>> On 19/5/2023 12:31 pm, Sean A. Irvine wrote:
>>>>>>>>>>>> > Hi,
>>>>>>>>>>>> >
>>>>>>>>>>>> > I believe the terms in A063747 are incorrect for
>>>> n>=46.
>>>>>>>>>>>> >
>>>>>>>>>>>> >https://oeis.org/A063747
>>>>>>>>>>>> >
>>>>>>>>>>>> > My attempts to reproduce the calculation in Maple
>>> and
>>>>>> using
>>>>>>>>>>>> constructible
>>>>>>>>>>>> > reals in Java both yield a coefficient, like
>>>>>>> -0.44458...e-38
>>>>>>>> *
>>>>>>>>>>>> x^46, from
>>>>>>>>>>>> > which a(46) should be -1 rather than 1.
>>>>>>>>>>>> >
>>>>>>>>>>>> > Is someone able to make an independent check?
>>>>>>>>>>>> >
>>>>>>>>>>>> > I was unable to run the existing Pari program
>>> beyond
>>>>>> n=15,
>>>>>>>> but
>>>>>>>>> I
>>>>>>>>>>>> assume
>>>>>>>>>>>> > there is some way to tell Pari to try harder?
>>>>>>>>>>>> >
>>>>>>>>>>>> > (I previously emailed the author of the sequence,
>>> but
>>>>>>>> received
>>>>>>>>>>>> no response.)
>>>>>>>>>>>> >
>>>>>>>>>>>> > Sean.
>>>>>>>>>>>> >
>>>>>>>>>>>> > --
>>>>>>>>>>>> > Seqfan Mailing list -http://list.seqfan.eu/
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> --
>>>>>>>>>>>> Seqfan Mailing list -http://list.seqfan.eu/
>>>>>>>>>>>>
>>>>>>>>>>> --
>>>>>>>>>>> Seqfan Mailing list -http://list.seqfan.eu/
>>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>> Seqfan Mailing list -http://list.seqfan.eu/
>>>>>>>>>>
>>>>>>>>> --
>>>>>>>>> Seqfan Mailing list -http://list.seqfan.eu/
>>>>>>>>>
>>>>>>>> --
>>>>>>>> Seqfan Mailing list -http://list.seqfan.eu/
>>>>>>>>
>>>>>>> --
>>>>>>> Seqfan Mailing list -http://list.seqfan.eu/
>>>>>>>
>>>>>> --
>>>>>> Seqfan Mailing list -http://list.seqfan.eu/
>>>>>>
>>>>> --
>>>>> Seqfan Mailing list -http://list.seqfan.eu/
>>>>>
>>>> --
>>>> Seqfan Mailing list -http://list.seqfan.eu/
>>>>
>>> --
>>> Seqfan Mailing list -http://list.seqfan.eu/
>>>
>> --
>> Seqfan Mailing list -http://list.seqfan.eu/
>
> --
> Seqfan Mailing list - http://list.seqfan.eu/
More information about the SeqFan
mailing list