[seqfan] Re: A063747 sign of power series expansion of the Gamma function
Brendan McKay
Brendan.McKay at anu.edu.au
Sat May 20 11:37:36 CEST 2023
I did it will interval arithmetic and got the same 201 terms.
It should be checked using a program other than Maple due to the
possibility of errors in the numerical evaluation of gamma and Zeta(n).
B/
On 20/5/2023 5:36 pm, Brendan McKay via SeqFan wrote:
> Of course epsilon is about 10^(-number of digits) not 10^(number of
> digits). B/
>
> On 20/5/2023 5:30 pm, Brendan McKay via SeqFan wrote:
>> 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/
>>
>> --
>> Seqfan Mailing list - http://list.seqfan.eu/
>
>
> --
> Seqfan Mailing list - http://list.seqfan.eu/
More information about the SeqFan
mailing list