Fixes incorporated during OctConf 2018.  Closing report.

Dear all,
I corrected a bug in the .m file and during the sprint code at the OctConf I worked with Carne to make a good patch. I attach the new changeset

(file #40115)

Dear all,
I just saw that there is a problem with my previous files (they are missing), now I hope that is all ok.

(file #39689)

I don't see a patch or changeset, just new files.

Changing the status to in progress.

We need a changeset with a proper commit message and tests.

What is the next step for getting this committed?

Nir and Marco, is this still "in progress", or do you want to push for it to be included in 4.2.0?

Bear in mind that review resources are particularly tight at this time.

Excellent, Nir. I see that you do not remove converged values in __gammainc_lentz as you did for gammainc_s (jj subindex). Isn't it dangerous? For instance, with your previous file #37329, if I do

gammainc(rand(4),rand(4))

quite often I have no convergence. The problem seems to me that at a certain iteration, some values have converged, but at the next iteration their error is a bit above the tolerance and other converged, and so on. However, I cannot reproduce this behavior with your new file #37342.

Yes, a simple scalar C++ implementation of that loop improves the speed quite a bit.
(after mkoctfile __gammainv_lentz.cc )

(file #37341, file #37342)

Yes, unfortunately the new implementation can be very slow:

m = 5;
n = 100;
a = rand (m, n);
x = rand (m, n);
tic
for i = 1:n
  y = gammainc (x(:, i), a(:, i));
endfor
t = toc #old: 0.4 s, new: 19 s

Profiling shows that the main bottleneck is the ngammainc_l subfunction, where the while loop can run thousands of times before the tolerance is reached. I streamlined this loop a bit (attached), but it's seemingly impossible to relax the tolerance without losing accuracy (tests start to fail).

(file #37329)

So, I'm uploading the latest version with

• name changed to gammainc (I feel very close to the final implemetation)
• vectorized D (thanks Nir for the idea)
• D(a,x) changed into D(x,a) (as gammaninc(x,a))
• updated test with smaller tolerances

Marco

(file #37328)

I think bug #40362 is triggered by

gammainc(600,600)

Current ngammainc handles without problems those arguments.
If you download file #37225:, call it gammainc.m, rename the inside function to gammainc, and put the file in your path, I think you should be able to overcome bug #40362. Can you please check that?

Must be compiler differences, since lgamma is a C++  built-in. I suggest incorporating D into ngammainc and then seeing we can agree on the ngammainc tolerances.

What do you get for the intermediate values in D(10, 10)? I have

> format hex
> log(10)

  40026bb1bbb55516

> gammaln(11)

  402e357590954d15

> D(10, 10)

 3fc0039b0be96dba

I get 3 tolerances not met with the new D on Octave 4.0.2:

*** test
 assert(D(10,10),0.12511003572113329898476497894772544708420990097708600,-7*eps)
!!!!! test failed
ASSERT errors for:  assert (D (10, 10),0.12511003572113329898476497894772544708420990097708600,-7 * eps)

  Location  |  Observed  |  Expected  |  Reason
     ()        0.12511      0.12511      Rel err 3.3277e-15 exceeds tol 1.5543e-15

*** test
 assert(D(10,-10),6.069903492836947893782490257188118487200593203967e7,-8*eps)
!!!!! test failed
ASSERT errors for:  assert (D (10, -10),6.069903492836947893782490257188118487200593203967e7,-8 * eps)

  Location  |  Observed  |  Expected  |  Reason
     ()      60699034.928 60699034.928   Rel err 5.1553e-15 exceeds tol 1.7764e-15

*** test
 assert(D(9,-10),-6.069903492836947893782490257188118487200593203967e7,-8*eps)
!!!!! test failed
ASSERT errors for:  assert (D (9, -10),-6.069903492836947893782490257188118487200593203967e7,-8 * eps)

  Location  |  Observed  |  Expected  |  Reason
     ()      -60699034.928 -60699034.928   Rel err 5.1553e-15 exceeds tol 1.7764e-15

The the handling of arrays is really nice -- not a single for loop!

Here is an updated version that copies most of the documentation from the existing C++ implementation.

I hadn't realised previously that the current implementation is in C++.  That presumably will make it much faster than the current one.  Would you propose to replace the existing C++ function by this one, or would you want to call the C++ function for tail=="lower"?  Or should this be translated this into C++, now that you have the algorithms selected?

(file #37225)

Hmmm...  Before this can be applied, it will definitely need to handle vectors, even if only by loops.

There are many special cases in this code.  I'd imagine that at least some of them can be vectorised efficiently.  The iterative components are probably the most expensive, and they seem to be vectorisable, although iteration would need to continue until all had converged, and so the actual FLOP count could be much higher than currently.

I'll add it to my to-do list, but will need reminding...

Yes, I release my contribution under the GPL

Thanks, Marco.
I added the exponential sum algorithm for small integer a and small x, which gives near-exact results.
Also, I found a more recent (2012) discussion involving Temme at http://oai.cwi.nl/oai/asset/20433/20433B.pdf Their Fortran code is posted at http://personales.unican.es/gila/incgam.zip , unfortunately with no license.

(file #37114)

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)

In case it helps with development, I've just added gammainc to Symbolic.  Here's a useful recipe for experimenting:

>> xs = sym (20); x = double (xs); a = 1;
>> A = gammainc (x, a, 'upper')
A =    2.0612e-09
>> B = double (gammainc (xs, a, 'upper'))
B =    2.0612e-09
>> relerr = (A-B)/abs(B)
relerr =   -2.0034e-08

I also put some 'xtest' BIST into '@sym/gammainc' that should pass once this bug is fixed.

I'm worried about the state of special functions in Octave :(

Thanks Nir,

but that code is quite inaccurate, try, for instance, gamain(2,3).
I started to implement, but how many corner cases! I used Numerical Recipes in Fortran77 (chapter 6.2). It is possible to use either a series or a continued fraction. I include my result, with some tests (the last still fails). Matlab is generally better for large input values, unless completely wrong. scaledupper and scaledlower options are not finished, yet. You can test it in Matlab (no endif...). The reference values in the tests were taken by wolframalpha. Don't know how to continue: modification of slatec functions (reliable, but forget x<0) or fixing the m-file.

(file #37075)

I had almost finished to implement the options scaledupper and scaledlower when I realized that Matlab supports negative x, while slatec (and thus Octave) does not. On the other hand, gammainc(-10,10) gives a complex result, while wolframalpha a real result (Matlab result's real part) and the usual series definition agrees with wolframalpha. gammainc(-10,10.1) gives the correct result, in Matlab and as the usual series. Now I think that a plain m file should be used (like expint, for instance).

Marco: I haven't looked at your patch yet, but your description sounds like the right direction to get this fixed in Octave (adding more Fortran subroutines along the same lines as xgammainc and xsgammainc).

I attach xcgmainc.f, a modification of slatec dgamic, which computes the complementary incomplete Gamma functions. The modifications are

1) transformation into a subroutine (like xgmainc.f by jwe)
2) division by Gamma(A) (not present in dgamic)

It requires d9gmic (not in liboctave/cruft/slatec-fn). When invoked with call xcgmainc(6.0d0,60.0d0,result) it gives result=6.1802e-20. Am I on the right way for the incorporation into Octave? Because Matlab allows also 'scaledupper' and 'scaledlower', which I would like to implement. Of course, the single version is still missing.

(file #37056)

The problem is really in
http://hg.savannah.gnu.org/hgweb/octave/file/@/libinterp/corefcn/gammainc.cc
There, gammainc(a, x, 'upper') is computed as 1 - gammainc(a, x), which expectedly evaluates to 0 if gammainc(a, x) is within eps of 1.

If we want more accurate results, we need a separate function to compute the complement.