[seqfan] Re: A063747 sign of power series expansion of the Gamma function
Brendan McKay
Brendan.McKay at anu.edu.au
Sat May 20 09:36:48 CEST 2023
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/
More information about the SeqFan
mailing list