Bug report (and possible fix): regularized incomplete gamma function in specfunc/gamma_inc.c -- continued fraction error

[gsl at boonstra dot org]

Hi all


   The continued fraction series for the regularized incomplete gamma
function appears to malfunction for certain values of a > x.  See the
value table appended below for details.


   By using the P series instead in this range, I believe the bug is
fixed.  The diff that works for me is included below.


    - Brian



*** /tmp/gamma_incCVS.c Thu Feb 10 14:00:33 2005
--- /tmp/gamma_incBRIAN.c       Thu Feb 10 14:01:34 2005
***************
*** 515,537 ****
      }
    }
    else {
-     if(0.8*a < x) {
-       /* Continued fraction again. The convergence
-        * is a little slower here, but that is fine.
-        * We have to trade that off against the slow
-        * convergence of the series, which is the
-        * only other option.
-        */
-       return gamma_inc_Q_CF(a, x, result);
-     }
-     else {
        gsl_sf_result P;
        int stat_P = gamma_inc_P_series(a, x, &P);
        result->val  = 1.0 - P.val;
        result->err  = P.err;
        result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
        return stat_P;
-     }
    }
  }

--- 515,526 ----




This is a table of a and the result.val of gsl_sf_gamma_inc_Q_e( a, 800 )
for values of a between 870 and 1020.  As you can see, once the (0.8*a<x)
condition is violated and the P series is used, the values snap back to
their correct levels.


{
{ 870 ,0.992422523945 } ,
{ 871 ,0.993111511866 } ,
{ 872 ,0.993744891076 } ,
{ 873 ,0.99432516717 } ,
{ 874 ,0.994856120805 } ,
{ 875 ,0.995343420979 } ,
{ 876 ,0.995788288413 } ,
{ 877 ,0.996195268676 } ,
{ 878 ,0.996574107188 } ,
{ 879 ,0.99691039566 } ,
{ 880 ,0.997231656484 } ,
{ 881 ,0.997500651595 } ,
{ 882 ,0.997803095401 } ,
{ 883 ,0.99797563731 } ,
{ 884 ,0.998116146641 } ,
{ 885 ,0.998605913831 } ,
{ 886 ,0.998785694801 } ,
{ 887 ,0.998326871442 } ,
{ 888 ,0.998340407608 } ,
{ 889 ,0.997814135102 } ,
{ 890 ,0.999657677983 } ,
{ 891 ,1.00007721129 } ,
{ 892 ,0.996092642105 } ,
{ 893 ,1.00064883702 } ,
{ 894 ,0.999618771111 } ,
{ 895 ,0.976656195765 } ,
{ 896 ,0.974468555553 } ,
{ 897 ,1.00573119105 } ,
{ 898 ,0.892945232274 } ,
{ 899 ,1.04096589162 } ,
{ 900 ,0.764266553047 } ,
{ 901 ,0.846149495572 } ,
{ 902 ,0.697324914548 } ,
{ 903 ,0.940849132517 } ,
{ 904 ,0.711076195463 } ,
{ 905 ,0.520921931192 } ,
{ 906 ,-10.9540823686 } ,
{ 907 ,-0.146534882293 } ,
{ 908 ,0.0882262678587 } ,
{ 909 ,0.139608307551 } ,
{ 910 ,0.030659237637 } ,
{ 911 ,-0.101356846252 } ,
{ 912 ,-0.0168428760234 } ,
{ 913 ,-0.0329260052005 } ,
{ 914 ,0.0044909367932 } ,
{ 915 ,-0.00155232057479 } ,
{ 916 ,0.00372447545372 } ,
{ 917 ,-0.00268357345391 } ,
{ 918 ,0.000606832396998 } ,
{ 919 ,-0.000329803783915 } ,
{ 920 ,-0.00114976625823 } ,
{ 921 ,0.000311719031751 } ,
{ 922 ,8.98670094573 10^-005 } ,
{ 923 ,7.44390786713 10^-005 } ,
{ 924 ,1.70636846254 10^-005 } ,
{ 925 ,5.0441924117 10^-005 } ,
{ 926 ,2.88034533718 10^-006 } ,
{ 927 ,-6.28445403526 10^-006 } ,
{ 928 ,-5.38447011763 10^-006 } ,
{ 929 ,-4.60862191278 10^-006 } ,
{ 930 ,-3.94050720462 10^-006 } ,
{ 931 ,-3.36578487502 10^-006 } ,
{ 932 ,-2.87193110093 10^-006 } ,
{ 933 ,-2.44802228196 10^-006 } ,
{ 934 ,-2.08454205134 10^-006 } ,
{ 935 ,-1.77320994802 10^-006 } ,
{ 936 ,-1.50682953877 10^-006 } ,
{ 937 ,-1.2791539753 10^-006 } ,
{ 938 ,-1.08476715484 10^-006 } ,
{ 939 ,-9.18978821825 10^-007 } ,
{ 940 ,-7.7773210573 10^-007 } ,
{ 941 ,-6.57522134585 10^-007 } ,
{ 942 ,-5.55324497055 10^-007 } ,
{ 943 ,-4.68532448197 10^-007 } ,
{ 944 ,-3.94901865921 10^-007 } ,
{ 945 ,-3.32503067382 10^-007 } ,
{ 946 ,-2.79678687667 10^-007 } ,
{ 947 ,-2.3500690778 10^-007 } ,
{ 948 ,-1.97269395733 10^-007 } ,
{ 949 ,-1.65423394005 10^-007 } ,
{ 950 ,-1.3857744945 10^-007 } ,
{ 951 ,-1.15970338246 10^-007 } ,
{ 952 ,-9.69527893893 10^-008 } ,
{ 953 ,-8.0971655896 10^-008 } ,
{ 954 ,-6.75562237924 10^-008 } ,
{ 955 ,-5.63063855504 10^-008 } ,
{ 956 ,-4.68824372849 10^-008 } ,
{ 957 ,-3.89962881087 10^-008 } ,
{ 958 ,-3.2403895903 10^-008 } ,
{ 959 ,-2.68987667123 10^-008 } ,
{ 960 ,-2.23063753124 10^-008 } ,
{ 961 ,-1.84793824792 10^-008 } ,
{ 962 ,-1.52935403703 10^-008 } ,
{ 963 ,-1.26441914233 10^-008 } ,
{ 964 ,-1.0443278487 10^-008 } ,
{ 965 ,-8.61679472239 10^-009 } ,
{ 966 ,-7.10261129063 10^-009 } ,
{ 967 ,-5.84862915253 10^-009 } ,
{ 968 ,-4.81120856257 10^-009 } ,
{ 969 ,-3.95383617418 10^-009 } ,
{ 970 ,-3.24599519306 10^-009 } ,
{ 971 ,-2.66220881687 10^-009 } ,
{ 972 ,-2.18123136996 10^-009 } ,
{ 973 ,-1.78536515851 10^-009 } ,
{ 974 ,-1.45988420357 10^-009 } ,
{ 975 ,-1.19254871696 10^-009 } ,
{ 976 ,-9.73196522826 10^-010 } ,
{ 977 ,-7.93399642385 10^-010 } ,
{ 978 ,-6.46175994122 10^-010 } ,
{ 979 ,-5.25747651985 10^-010 } ,
{ 980 ,-4.2733838364 10^-010 } ,
{ 981 ,-3.47004287032 10^-010 } ,
{ 982 ,-2.81492281643 10^-010 } ,
{ 983 ,-2.28122012487 10^-010 } ,
{ 984 ,-1.84687408901 10^-010 } ,
{ 985 ,-1.49374723059 10^-010 } ,
{ 986 ,-1.20694369076 10^-010 } ,
{ 987 ,-9.74243049918 10^-011 } ,
{ 988 ,-7.85630575073 10^-011 } ,
{ 989 ,-6.32907923378 10^-011 } ,
{ 990 ,-5.09370894706 10^-011 } ,
{ 991 ,-4.09542992749 10^-011 } ,
{ 992 ,-3.28955382693 10^-011 } ,
{ 993 ,-2.63965374583 10^-011 } ,
{ 994 ,-2.11606858482 10^-011 } ,
{ 995 ,-1.69467207697 10^-011 } ,
{ 996 ,-1.35586081438 10^-011 } ,
{ 997 ,-1.08372325406 10^-011 } ,
{ 998 ,-8.65358110375 10^-012 } ,
{ 999 ,-6.90315911028 10^-012 } ,
{ 1000 ,0.999999999994 } ,
{ 1001 ,0.999999999996 } ,
{ 1002 ,0.999999999997 } ,
{ 1003 ,0.999999999997 } ,
{ 1004 ,0.999999999998 } ,
{ 1005 ,0.999999999998 } ,
{ 1006 ,0.999999999999 } ,
{ 1007 ,0.999999999999 } ,
{ 1008 ,0.999999999999 } ,
{ 1009 ,0.999999999999 } ,
{ 1010 ,0.999999999999 } ,
{ 1011 ,1 } ,
{ 1012 ,1 } ,
{ 1013 ,1 } ,
{ 1014 ,1 } ,
{ 1015 ,1 } ,
{ 1016 ,1 } ,
{ 1017 ,1 } ,
{ 1018 ,1 } ,
{ 1019 ,1 }
}