#2173 gamma_incomplete(1/2, -24+10*%i) not accurate

Raymond Toy

Maxima says gamma_incomplete(1/2,-20+10*%i),float is

4.93408396056858e+9 %i + 1.91200542673871e+9

Wolfram/Alpha says the value is

4.9341163e+9 + 1.911933e+9

Maxima's answer is only accurate to about 5 digits.


  • Raymond Toy

    Raymond Toy - 2011-03-24

    Setting *gamma-imag* to .1d0 (so the region where we use the continued fraction is enlarged), improves the result to

    4.93411631550489e+9 %i + 1.911933769523828e+9

    Wolfram/Alpha gives

    4.9341163155048952691479351508925718965852838271239234... × 10^9 i +
    1.9119337695238284027493471541730604272204040299891632... × 10^9

    Essentially full accuracy.

  • Raymond Toy

    Raymond Toy - 2011-03-24
    • summary: gamma_incomplete(1/2, -2410%i) not accurate --> gamma_incomplete(1/2, -24+10*%i) not accurate
  • Dan Gildea

    Dan Gildea - 2012-01-28

    Maxima 5.26.0_26_gc4216e7 http://maxima.sourceforge.net
    using Lisp CMU Common Lisp 20a (20A Unicode)
    (%i1) gamma_incomplete(1/2,-20+10*%i),float;
    (%o1) 9.902410782885888e+7 %i + 3.4166055039087765e+7

  • Raymond Toy

    Raymond Toy - 2012-01-31

    Sorry. The description text is wrong. The subject line has the correct parameters for the gamma function:

    gamma_incomplete(1/2, -24+10*%i), float;
    4.934116315504894e+9 %i + 1.911933769523827e+9

    Wolfram alpha says
    1.911933769523828402749347154173060427220404029989163... × 10^9 +
    4.934116315504895269147935150892571896585283827123923... × 10^9 i

    So it looks like we are accurate to full precision now.

    Marking as pending/worksforme.

  • Raymond Toy

    Raymond Toy - 2012-01-31
    • status: open --> pending
  • Dan Gildea

    Dan Gildea - 2012-01-31

    (%i6) i : integrate(sin(x)^2/x,x);
    (%o6) (2*log(x)+gamma_incomplete(0,2*%i*x)+gamma_incomplete(0,-2*%i*x))/4

    I think that the indefinite integral above is correct.
    However, if we evaluate it numerically:

    (%i7) i,x:10,numer,expand;
    (%o7) 1.129082636318599
    (%i8) i,x:100,numer,expand;
    (%o8) 1.167036722342424e+67

    we get something that does not match the result of numerical integration:

    (%i12) quad_qags(sin(x)^2/x,x,10,100);
    (%o12) [1.175691679966213,1.71048903216243e-11,651,0]

    Not sure if this is a branch cut issue or what, but the results for
    gamma_incomplete here do not match wolfram alpha.

  • Raymond Toy

    Raymond Toy - 2012-01-31
    • status: pending --> open
  • Raymond Toy

    Raymond Toy - 2012-01-31

    Yes, this is a known bug with gamma_incomplete: it's totally wrong for large arguments. There's already a bug filed for that, I think.

  • Raymond Toy

    Raymond Toy - 2012-08-16
    • status: open --> closed
  • Raymond Toy

    Raymond Toy - 2012-08-16

    This was fixed sometime ago.


Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

No, thanks