[seqfan] Re: A063747 sign of power series expansion of the Gamma function
Hugo Pfoertner
yae9911 at gmail.com
Sat May 20 13:33:40 CEST 2023
With \p339 and Michel Marcus's PARI
my(x='x+O('x^201)); apply(sign, Vec(1/(gamma(1-x))))
confirms the reported 201 terms. 339 is obtained from round(n*log(n)/Pi)
suggested as a rough estimate of required (decimal) localprec in the pink
box discussion.
All programs provided within the sequence should contain an indication of
the required precision as a function of n and not just an accuracy value
that can only be understood with reference to this SeqFan discussion.
On Sat, May 20, 2023 at 11:38 AM Brendan McKay via SeqFan <
seqfan at list.seqfan.eu> wrote:
> 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/
>
>
> --
> Seqfan Mailing list - http://list.seqfan.eu/
>
More information about the SeqFan
mailing list