[Octave-bug-tracker] [bug #47800] gammainc(x, a, "upper") rounds down to zero if output is below eps

[Nir Krakauer]

Follow-up Comment #13, bug #47800 (project octave):

Yes, the scale factor seems to have more error than the series sum. We could
adopt Temme's idea, but then we will need to accurately compute his "tempered
gamma". Failing that, I was able to improve the error slightly:

#with previously sent program:
x = 200; a = 200;
scale = exp(x - a*log(x) + gammaln(a+1)); #fractional error of 9E-15 compared
to Wolfram Alpha's result of Gamma(201)*e^200*200^-200 =
35.46385053215837215657814935957664087171927837002301455738
y = newgammainc(200,200,'scaledupper'); #essentially the series sum --
fractional error of 2E-15 from Gamma(201)*e^200*200^-200*GammaRegularized[200,
200] = 17.39844385537915051351229001995819645096209733412625912131
y2 = newgammainc(200,200,'scaledlower'); #fractional error of 2E-14 from
Gamma(201)*e^200*200^-200*(1 - GammaRegularized[200, 200]) =
18.06540667677922164306585933961844442075718103589675543606
y3 = newgammainc(200,200,'upper'); #fractional error of 5E-14 from
GammaRegularized[200, 200] =
0.490596581992763674972174540425472956845027831400859498457
y4 = newgammainc(200,200,'lower'); #fractional error of 5E-14 from 1 -
GammaRegularized[200, 200] =
0.509403418007236325027825459574527043154972168599140501542
scale2 = exp(x - a*log(x) + gammaln(a)); #fractional error of 5E-14 compared
to Wolfram Alpha's result of Gamma(200)*e^200*200^-200 =
0.177319252660791860782890746797883204358596391850115072786

#with new attached version, where scale instead of scale2 is used to compute
the unscaled variants:
y3 = newgammainc(200,200,'upper'); #fractional error of 1E-14
y4 = newgammainc(200,200,'lower'); #fractional error of 1E-14

(file #37091)
    _______________________________________________________

Additional Item Attachment:

File name: newgammainc.m                  Size:5 KB


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?47800>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/