Fixes incorporated during OctConf 2018. Closing report.

Dear all, I corrected a bug in the .m file and added a description about the "scaled" option in the help. Then I remake the patch to work with the new organization of the folder (liboctave/cruft has been renamed in liboctave/external). You can find the patch attached or you can pull it from the gammainc bookmark on my repository
https://bitbucket.org/M_Ginesi/octave
(file #40696)

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 repost the patch because the old one stopped working (mercurial couldn't remove gammainc.cc since a line about the copyright was modified). Now it work with the last version (I'm on revision 23250).
(file #39928)

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)

Dear all,
I have made some minor modifications to gammainc.m and __gammainc_lentz.cc (files #39669 and #39670 )in such a way to adapt to the class of the input (as the actual gammainc do: if one between x and a is single, then the output is single as well).
Then I tried to make a patch (file #39668 ). I was able to build octave with this, but it's the first time I try to make a patch in which there is more to do than to modify a single file, so I'm not sure it is "wellmade". So if everyone has some suggestion, I would like to help in making this patch definitive.

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.

Changing status to "patch reviewed"
Maintainers: please see comment #41

What is the next step for getting this committed?

For me it is ready for pushing.
Marco

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.

yes, actually __gammainc_lentz currently deals with scalar inputs, so there is only one convergence at a time to worry about

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
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.

*should be mkoctfile __gammainc_lentz.cc

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

First of all, the trick with
is missing in ngammainc_l. Now included.
Second, ngammainc_l is used for x>0, therefore always real results. What about to write an oct file just for ngammainc_l? it should not be difficult.
(file #37335)

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)

Well done on fixing bug #40362 as well.
I'm still concerned about the loss of performance relative to the current C++/Fortran implementation. Do you think it won't be an issue, or is it worth working out which cases the builtin version can handle and calling it where possible?

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)

Yes, in octave 3.6.2 as well as 3.8.0 this fixes the reported problems with bug #40362. Thank you.

I think bug #40362 is triggered by
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?

I did not take the time to follow the whole discussion through, but do you propose to replace the lowlevel computation of the inverse gamma function? If so, does this fix bug #40362?

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

ans = 40026bb1bbb55516
ans = 402e357590954d16
ans = 3fc0039b0be96db2
So, gammaln(11) (and therefore D(10,10)) differs in the last digit.

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

Weird. It passes all tests on my 4.0.0 and 4.1.0+. Failure for D(10,10) means there is something in exp, log or gammaln. The new code is not used in this case.

I get 3 tolerances not met with the new D on Octave 4.0.2:
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.3277e15 exceeds tol 1.5543e15
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.1553e15 exceeds tol 1.7764e15
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.1553e15 exceeds tol 1.7764e15

I include my new version of the subfunction D. It allows to reduce a lot the tolerances for the tests in ngammainc. It not vectorized, yet.
(file #37232)

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)

Thanks Andrew and Marco. I added a license notice and now provide for handling of array x,a inputs.
(file #37213)

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 todo list, but will need reminding...

I release my contribution under the GPL, too.
The only functionality missing with respect to the current Octave implementation (and Matlab as well) is the possibility to work with a and x vectors. I do not see an easy way to vectorize: in expint.m, for instance, there is a simple

Yes, I release my contribution under the GPL

Thanks to you both for your great work.
I've tidied the formatting to comply with Octave's coding standards, but haven't yet touched the docstring. The code has no licence. Are you both happy for it to be GPL'd?
FYI, the changes I've made to comply with the standards are:
 endif/endwhile/endfunction instead of end
 double quotes instead of single quotes
 ## instead of % to introduce fullline comments
 space before function call parenthesis
 spacing around operators
 ; after the "assert" statements
 lines wrapped after 72~76 columns.
In addition
 strcmp > strcmpi
 Use of +=, *= etc when possible, since they are actually faster than "n = n + 1" etc.
 Removed what seemed to be a duplicated "aj = a;" before the comment
since there was also one after that comment.
(file #37206)

Thanks, Marco.
I added the exponential sum algorithm for small integer a and small x, which gives nearexact 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)

Hi,
I include my final proposal. It has corner cases first, then a series and finally a continued fraction expansion. All tests gammainc(200,200,tail) are passed with at most relative error 78*eps. I used all the ideas discussed here (thanks Nir). Two points are still missing
1) accurate implementation of D(a,x) (seems hard to implement)
2) Matlab is exact for a integer. I found a close form (http://mathworld.wolfram.com/IncompleteGammaFunction.html, formula (2)), but can't implement without over/underflow.
(file #37113)

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 9E15 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 2E15 from Gamma(201)*e^200*200^200*GammaRegularized[200, 200] = 17.39844385537915051351229001995819645096209733412625912131
y2 = newgammainc(200,200,'scaledlower'); #fractional error of 2E14 from Gamma(201)*e^200*200^200*(1  GammaRegularized[200, 200]) = 18.06540667677922164306585933961844442075718103589675543606
y3 = newgammainc(200,200,'upper'); #fractional error of 5E14 from GammaRegularized[200, 200] = 0.490596581992763674972174540425472956845027831400859498457
y4 = newgammainc(200,200,'lower'); #fractional error of 5E14 from 1  GammaRegularized[200, 200] = 0.509403418007236325027825459574527043154972168599140501542
scale2 = exp(x  a*log(x) + gammaln(a)); #fractional error of 5E14 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 1E14
y4 = newgammainc(200,200,'lower'); #fractional error of 1E14
(file #37091)

Hi Nir. If I understand correctly, you (Gautschi) can avoid overflow/underflow by computing the series for the scaled version and scaling at the end. The Legendre continued fraction is instead formula (6.2.6) of Numerical Recipes. I used (6.2.7) (claimed to be faster), but only for abs(x) > a+1. Anyway, I found this interesting paper
https://www.researchgate.net/publication/231929375_A_Set_of_Algorithms_for_the_Incomplete_Gamma_Functions
In section 5.1 it is clear that either the series or the continued fraction method is used. But it is interesting the first part of the paper, in which the accurate computation of D(a,x)=x^a*exp(x)/Gamma(a+1) is discussed. I already notice that that standard way exp(a*log(x)xgammaln(a+1)) is not always accurate (and may explain the loss of 34 digits).
To summarize, I would guess (or hope) that the series/continued fraction method + accurate computation of D(a,x) would manage all the cases (instead of currently two series and two continued fractions).
By the way, "improvement of special functions": what a nice project for next GSOC!

In case it helps with development, I've just added gammainc to Symbolic. Here's a useful recipe for experimenting:
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 :(

Hi Marco I added two more methods to your mfile to deal with positive a and x without overflow. Now all the tests pass, albeit at fractional tolerances of ~1000*eps. I also streamlined the processing of scaled* options.
Based on our results so far, I think we should first implement gammainc as an mfile, since the available Fortran programs don't work for all cases and are difficult to debug. Once we have a working implementation, we can look at increasing efficiency by calling C++ or Fortran for some branches.
(file #37083)

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 mfile.
(file #37075)

There is possibly useful LGPL code in different languages at http://people.sc.fsu.edu/~jburkardt/cpp_src/asa032/asa032.html

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).

Is there a reason why the file xgmainc.f (by jwe) has inside the subroutine xgammainc? How should I call the files and the subroutines (I think I will need two more subroutines, for the options "scaledlower" and "scaledupper")? Then, I will need some advice in order to incorporate it into Octave.

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).

Yes, it's the correct answer! [note the different orderings of a and x  xcgmainc(6.0d0,60.0d0,result) = gammainc(60, 6, "upper") = GammaRegularized[6, 60] in Wolfram Alpha]

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/slatecfn). When invoked with call xcgmainc(6.0d0,60.0d0,result) it gives result=6.1802e20. 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)

Yes Nir,
the function dgamic.f should probabily be used. I can give a look.

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.

Consider this:
In Matlab:
>> gammainc(60,6,'upper')
ans =
6.1802e20
In Octave 4.0.0 (same for 3.8.2):
>> gammainc(60,6,'upper')
ans = 0
This causes large problems, as norminv(6.1802e20) = 9.0659 (the expected result), whereas norminv(0) = Inf, and subsequent code crashes.
