Re: AW: hypergeometric function

[Brian Gough <bjg at network-theory dot co dot uk>]

Stefan Koch writes:
 > With the following function I produced the output shown in the attached file
 > Hyper.jpg. You can see the discontinuity clearly. Please note that I
 > multiply the output of the hypergeometric function by a factor of 100.
 > I used gsl 1.0 from GnuWin32, the precompiled version (dll). I got the
 > dicontinuity with the free download version from Network Theory Ltd. too.

Ok, I see it now at x=13..14.  I didn't see it on the original scale.

Thanks for the bug report.  I've found the problem, so a revised
version will be available in the next release or you can apply the
patch below and rebuild from source to fix it.

regards
Brian Gough

Index: specfunc/hyperg_1F1.c
===================================================================
RCS file: /cvs/gsl/gsl/specfunc/hyperg_1F1.c,v
retrieving revision 1.67
retrieving revision 1.69
diff -u -r1.67 -r1.69
--- hyperg_1F1.c	2001/11/19 22:32:10	1.67
+++ hyperg_1F1.c	2001/12/05 19:34:41	1.69
@@ -1084,13 +1084,28 @@
 	Ma0p1b = (b*(a0+x)*Ma0b + x*(a0-b)*Ma0bp1)/(a0*b);
       }
 
-      Mnm1 = Ma0b;
-      Mn   = Ma0p1b;
-      for(n=a0+1; n<a; n++) {
-    	Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
-    	Mnm1 = Mn;
-    	Mn   = Mnp1;
-      }
+      /* Initialise the recurrence correctly BJG */
+
+      if (a0 >= a)
+        { 
+          Mn = Ma0b;
+        }
+      else if (a0 + 1>= a)
+        {
+          Mn = Ma0p1b;
+        }
+      else
+        {
+          Mnm1 = Ma0b;
+          Mn   = Ma0p1b;
+
+          for(n=a0+1; n<a; n++) {
+            Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
+            Mnm1 = Mn;
+            Mn   = Mnp1;
+          }
+        }
+
       result->val  = Mn;
       result->err  = (fabs(x) + 1.0) * GSL_DBL_EPSILON *  fabs(Mn);
       result->err *= fabs(b-a)+1.0;