[seqfan] Re: A063747 sign of power series expansion of the Gamma function

Brendan McKay Brendan.McKay at anu.edu.au
Sat May 20 08:43:38 CEST 2023 

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/

More information about the SeqFan mailing list