Mon 15 May 2017 08:18:53 PM UTC, comment #48:
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)

Fri 24 Mar 2017 02:12:48 PM UTC, comment #47:
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)

Tue 07 Mar 2017 06:05:52 PM UTC, comment #46:
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)

Tue 07 Feb 2017 01:17:11 PM UTC, comment #45:
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)

Sun 05 Feb 2017 11:33:08 AM UTC, comment #44:
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.

Tue 31 Jan 2017 01:48:44 PM UTC, comment #43:
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.

Wed 25 Jan 2017 07:25:47 PM UTC, comment #42:
Changing status to "patch reviewed"
Maintainers: please see comment #41

Thu 17 Nov 2016 10:58:42 AM UTC, comment #41:
What is the next step for getting this committed?

Sun 03 Jul 2016 09:23:59 AM UTC, comment #40:
For me it is ready for pushing.
Marco

Sat 02 Jul 2016 06:14:17 AM UTC, comment #39:
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.

Thu 02 Jun 2016 12:38:10 PM UTC, comment #38:
yes, actually __gammainc_lentz currently deals with scalar inputs, so there is only one convergence at a time to worry about

Thu 02 Jun 2016 07:51:54 AM UTC, comment #37:
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.

Wed 01 Jun 2016 07:10:20 PM UTC, comment #36:
*should be mkoctfile __gammainc_lentz.cc

Wed 01 Jun 2016 07:08:02 PM UTC, comment #35:
Yes, a simple scalar C++ implementation of that loop improves the speed quite a bit.
(after mkoctfile __gammainv_lentz.cc )
(file #37341, file #37342)

Wed 01 Jun 2016 03:42:38 PM UTC, comment #34:
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)

Wed 01 Jun 2016 02:01:03 PM UTC, comment #33:
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)

Wed 01 Jun 2016 01:01:30 PM UTC, comment #32:
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?

Wed 01 Jun 2016 12:53:09 PM UTC, comment #31:
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)

Wed 01 Jun 2016 06:30:00 AM UTC, comment #30:
Yes, in octave 3.6.2 as well as 3.8.0 this fixes the reported problems with bug #40362. Thank you.

Tue 31 May 2016 10:04:01 AM UTC, comment #29:
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?

Tue 31 May 2016 07:55:33 AM UTC, comment #28:
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?

Thu 26 May 2016 04:43:15 PM UTC, comment #27:
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.

Tue 24 May 2016 06:35:55 PM UTC, comment #26:
ans = 40026bb1bbb55516
ans = 402e357590954d16
ans = 3fc0039b0be96db2
So, gammaln(11) (and therefore D(10,10)) differs in the last digit.

Tue 24 May 2016 05:06:56 PM UTC, comment #25:
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

Tue 24 May 2016 04:58:36 PM UTC, comment #24:
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.

Tue 24 May 2016 03:55:26 PM UTC, comment #23:
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

Mon 23 May 2016 07:00:18 AM UTC, comment #22:
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)

Sun 22 May 2016 11:13:40 AM UTC, comment #21:
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)

Fri 20 May 2016 04:46:05 PM UTC, comment #20:
Thanks Andrew and Marco. I added a license notice and now provide for handling of array x,a inputs.
(file #37213)

Fri 20 May 2016 06:44:33 AM UTC, comment #19:
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...

Fri 20 May 2016 06:32:01 AM UTC, comment #18:
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

Thu 19 May 2016 01:50:11 PM UTC, comment #17:
Yes, I release my contribution under the GPL

Thu 19 May 2016 11:59:27 AM UTC, comment #16:
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)

Mon 09 May 2016 04:18:20 PM UTC, comment #15:
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)

Mon 09 May 2016 01:52:26 PM UTC, comment #14:
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)

Fri 06 May 2016 04:57:58 PM UTC, comment #13:
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  alog(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^200200^200*GammaRegularized[200, 200] = 17.39844385537915051351229001995819645096209733412625912131
y2 = newgammainc(200,200,'scaledlower'); #fractional error of 2E14 from Gamma(201)e^200200^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  alog(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)

Fri 06 May 2016 06:54:57 AM UTC, comment #12:
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^aexp(x)/Gamma(a+1) is discussed. I already notice that that standard way exp(alog(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!

Fri 06 May 2016 05:10:10 AM UTC, comment #11:
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 :(

Thu 05 May 2016 08:14:04 PM UTC, comment #10:
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 ~1000eps. 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)

Thu 05 May 2016 12:47:36 PM UTC, comment #9:
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)

Wed 04 May 2016 04:25:52 PM UTC, comment #8:
There is possibly useful LGPL code in different languages at http://people.sc.fsu.edu/~jburkardt/cpp_src/asa032/asa032.html

Wed 04 May 2016 02:23:13 PM UTC, comment #7:
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).

Wed 04 May 2016 06:54:07 AM UTC, comment #6:
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.

Tue 03 May 2016 10:46:55 PM UTC, comment #5:
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).

Tue 03 May 2016 01:02:44 PM UTC, comment #4:
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]

Tue 03 May 2016 10:56:45 AM UTC, comment #3:
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)

Mon 02 May 2016 06:54:36 AM UTC, comment #2:
Yes Nir,
the function dgamic.f should probabily be used. I can give a look.

Thu 28 Apr 2016 06:11:03 PM UTC, comment #1:
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.

Thu 28 Apr 2016 06:02:53 PM UTC, original submission:
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.
