This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: AW: hypergeometric function
- From: Brian Gough <bjg at network-theory dot co dot uk>
- To: <s dot koch at seh-zuelpich dot de>
- Cc: gsl-discuss at sources dot redhat dot com
- Date: Wed, 5 Dec 2001 19:38:42 +0000 (GMT)
- Subject: Re: AW: hypergeometric function
- References: <15373.15152.311318.362697@debian><NEBBIBGCPKOPEIODNFDIEEIHCAAA.s.koch@seh-zuelpich.de>
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;