[Scipy-svn] r2041 - in trunk/Lib/special: cephes tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Jul 5 14:27:59 EDT 2006
Author: cookedm
Date: 2006-07-05 13:27:34 -0500 (Wed, 05 Jul 2006)
New Revision: 2041
Modified:
trunk/Lib/special/cephes/airy.c
trunk/Lib/special/cephes/hyp2f1.c
trunk/Lib/special/cephes/hyperg.c
trunk/Lib/special/cephes/i0.c
trunk/Lib/special/cephes/i1.c
trunk/Lib/special/cephes/iv.c
trunk/Lib/special/cephes/j0.c
trunk/Lib/special/cephes/j1.c
trunk/Lib/special/cephes/jn.c
trunk/Lib/special/cephes/jv.c
trunk/Lib/special/cephes/k0.c
trunk/Lib/special/cephes/k1.c
trunk/Lib/special/cephes/kn.c
trunk/Lib/special/cephes/mconf.h
trunk/Lib/special/cephes/psi.c
trunk/Lib/special/cephes/struve.c
trunk/Lib/special/cephes/yn.c
trunk/Lib/special/tests/test_basic.py
Log:
upgrade bessel functions (and related) in cephes to cephes 2.8
Modified: trunk/Lib/special/cephes/airy.c
===================================================================
--- trunk/Lib/special/cephes/airy.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/airy.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -51,9 +51,8 @@
/* airy.c */
/*
-Cephes Math Library Release 2.1: January, 1989
-Copyright 1984, 1987, 1989 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
@@ -823,7 +822,15 @@
};
#endif
-#ifndef ANSIPROT
+#ifdef ANSIPROT
+extern double fabs ( double );
+extern double exp ( double );
+extern double sqrt ( double );
+extern double polevl ( double, void *, int );
+extern double p1evl ( double, void *, int );
+extern double sin ( double );
+extern double cos ( double );
+#else
double fabs(), exp(), sqrt();
double polevl(), p1evl(), sin(), cos();
#endif
Modified: trunk/Lib/special/cephes/hyp2f1.c
===================================================================
--- trunk/Lib/special/cephes/hyp2f1.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/hyp2f1.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -31,11 +31,10 @@
* Transformation for x < -0.5
* Psi function expansion if x > 0.5 and c - a - b integer
* Conditionally, a recurrence on c to make c-a-b > 0
- *
- * x < -1 AMS 15.3.7 transformation applied (Travis Oliphant).
- * valid for b,a,c,(b-a) != integer and (c-a),(c-b) != negative integer
- *
*
+ * x < -1 AMS 15.3.7 transformation applied (Travis Oliphant)
+ * valid for b,a,c,(b-a) != integer and (c-a),(c-b) != negative integer
+ *
* x >= 1 is rejected (unless special cases are present)
*
* The parameters a, b, c are considered to be integer
@@ -64,9 +63,8 @@
/*
-Cephes Math Library Release 2.2: November, 1992
-Copyright 1984, 1987, 1992 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier
*/
@@ -94,13 +92,22 @@
#define ETHRESH 1.0e-12
-#ifndef ANSIPROT
-double fabs(), pow(), round(), Gamma(), log(), exp(), psi();
-static double hyt2f1();
-static double hys2f1();
-#else
+#ifdef ANSIPROT
+extern double fabs ( double );
+extern double pow ( double, double );
+extern double round ( double );
+extern double gamma ( double );
+extern double log ( double );
+extern double exp ( double );
+extern double psi ( double );
static double hyt2f1(double, double, double, double, double *);
static double hys2f1(double, double, double, double, double *);
+double hyp2f1(double, double, double, double);
+#else
+double fabs(), pow(), round(), gamma(), log(), exp(), psi();
+static double hyt2f1();
+static double hys2f1();
+double hyp2f1();
#endif
extern double MAXNUM, MACHEP, NAN;
@@ -110,7 +117,6 @@
double d, d1, d2, e;
double p, q, r, s, y, ax;
double ia, ib, ic, id, err;
-double t1;
int flag, i, aid;
err = 0.0;
@@ -120,7 +126,6 @@
ia = round(a); /* nearest integer to a */
ib = round(b);
-
if( a <= 0 )
{
if( fabs(a-ia) < EPS ) /* a is a negative integer */
@@ -133,23 +138,19 @@
flag |= 2;
}
-
-if( fabs(ax-1.0) > EPS ) /* |x| != 1.0 */
- {
- if( fabs(b-c) < EPS ) /* b = c */
- {
+if (fabs(ax-1.0) > EPS) { /* |x| != 1.0 */
+ if( fabs(b-c) < EPS ) { /* b = c */
y = pow( s, -a ); /* s to the -a power */
goto hypdon;
- }
-
- if( fabs(a-c) < EPS ) /* a = c */
- {
+ }
+ if( fabs(a-c) < EPS ) { /* a = c */
y = pow( s, -b ); /* s to the -b power */
goto hypdon;
- }
+ }
}
+
if( c <= 0.0 )
{
ic = round(c); /* nearest integer to c */
@@ -167,7 +168,22 @@
if( flag ) /* function is a polynomial */
goto hypok;
+if (x < -1.0) {
+ double t1;
+ r = -x;
+ p = hyp2f1(a, 1-c+a, 1-b+a, 1.0/x);
+ q = hyp2f1(b, 1-c+b, 1-b+a, 1.0/x);
+ p *= pow(r, -a);
+ q *= pow(r, -b);
+ t1 = gamma(c);
+ s = t1*gamma(b-a)/(gamma(b)*gamma(c-a));
+ y = t1*gamma(a-b)/(gamma(a)*gamma(c-b));
+ return s*p + y*q;
+}
+if( ax > 1.0 ) /* series diverges */
+ goto hypdiv;
+
p = c - a;
ia = round(p); /* nearest integer to c-a */
if( (ia <= 0.0) && (fabs(p-ia) < EPS) ) /* negative int c - a */
@@ -182,7 +198,6 @@
id = round(d); /* nearest integer to d */
q = fabs(d-id);
-
/* Thanks to Christian Burger <BURGER at DMRHRZ11.HRZ.Uni-Marburg.DE>
* for reporting a bug here. */
if( fabs(ax-1.0) < EPS ) /* |x| == 1.0 */
@@ -198,7 +213,7 @@
}
if( d <= 0.0 )
goto hypdiv;
- y = Gamma(c)*Gamma(d)/(Gamma(p)*Gamma(r));
+ y = gamma(c)*gamma(d)/(gamma(p)*gamma(r));
goto hypdon;
}
@@ -238,11 +253,6 @@
if( flag & 12 )
goto hypf; /* negative integer c-a or c-b */
-
-if( ax > 1.0 ) /* series diverges */
- goto hypdiv;
-
-
hypok:
y = hyt2f1( a, b, c, x, &err );
@@ -264,27 +274,15 @@
/* The alarm exit */
hypdiv:
-
-/* Added by Travis Oliphant */
-
-if ((x < -1) & (fabs(b-a-(ib-ia)) < EPS)) { /* Handle negative values of x */
- r = -x;
- p = hyp2f1(a, 1-c+a, 1-b+a, 1.0/x);
- q = hyp2f1(b, 1-c+b, 1-a+b, 1.0/x);
- p *= pow (r, -a);
- q *= pow (r, -b);
- t1 = Gamma(c);
- s = t1*Gamma(b-a)/(Gamma(b)*Gamma(c-a));
- y = t1*Gamma(a-b)/(Gamma(a)*Gamma(c-b));
- return s*p + y*q;
-}
-
mtherr( "hyp2f1", OVERFLOW );
return( MAXNUM );
}
+
+
+
/* Apply transformations for |x| near 1
* then call the power series
*/
@@ -322,9 +320,9 @@
goto done;
/* If power series fails, then apply AMS55 #15.3.6 */
q = hys2f1( a, b, 1.0-d, s, &err );
- q *= Gamma(d) /(Gamma(c-a) * Gamma(c-b));
+ q *= gamma(d) /(gamma(c-a) * gamma(c-b));
r = pow(s,d) * hys2f1( c-a, c-b, d+1.0, s, &err1 );
- r *= Gamma(-d)/(Gamma(a) * Gamma(b));
+ r *= gamma(-d)/(gamma(a) * gamma(b));
y = q + r;
q = fabs(q); /* estimate cancellation error */
@@ -333,7 +331,7 @@
r = q;
err += err1 + (MACHEP*r)/y;
- y *= Gamma(c);
+ y *= gamma(c);
goto done;
}
else
@@ -358,9 +356,9 @@
/* sum for t = 0 */
y = psi(1.0) + psi(1.0+e) - psi(a+d1) - psi(b+d1) - ax;
- y /= Gamma(e+1.0);
+ y /= gamma(e+1.0);
- p = (a+d1) * (b+d1) * s / Gamma(e+2.0); /* Poch for t=1 */
+ p = (a+d1) * (b+d1) * s / gamma(e+2.0); /* Poch for t=1 */
t = 1.0;
do
{
@@ -377,7 +375,7 @@
if( id == 0.0 )
{
- y *= Gamma(c)/(Gamma(a)*Gamma(b));
+ y *= gamma(c)/(gamma(a)*gamma(b));
goto psidon;
}
@@ -397,10 +395,10 @@
y1 += p;
}
nosum:
- p = Gamma(c);
- y1 *= Gamma(e) * p / (Gamma(a+d1) * Gamma(b+d1));
+ p = gamma(c);
+ y1 *= gamma(e) * p / (gamma(a+d1) * gamma(b+d1));
- y *= p / (Gamma(a+d2) * Gamma(b+d2));
+ y *= p / (gamma(a+d2) * gamma(b+d2));
if( (aid & 1) != 0 )
y = -y;
Modified: trunk/Lib/special/cephes/hyperg.c
===================================================================
--- trunk/Lib/special/cephes/hyperg.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/hyperg.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -60,21 +60,27 @@
/*
-Cephes Math Library Release 2.1: November, 1988
-Copyright 1984, 1987, 1988 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
-#ifndef ANSIPROT
-double exp(), fabs();
-static double hy1f1p();
-static double hy1f1a();
-#else
+#ifdef ANSIPROT
+extern double exp ( double );
+extern double log ( double );
+extern double gamma ( double );
+extern double lgam ( double );
+extern double fabs ( double );
+double hyp2f0 ( double, double, double, int, double * );
static double hy1f1p(double, double, double, double *);
static double hy1f1a(double, double, double, double *);
-double hyp2f0( double,double,double,int,double*);
+double hyperg (double, double, double);
+#else
+double exp(), log(), gamma(), lgam(), fabs(), hyp2f0();
+static double hy1f1p();
+static double hy1f1a();
+double hyperg();
#endif
extern double MAXNUM, MACHEP;
@@ -125,10 +131,7 @@
double *err;
{
double n, a0, sum, t, u, temp;
-#ifndef ANSIPROT
-double fabs();
-#endif
-double an, bn, maxt, pcanc;
+double an, bn, maxt;
/* set up for power series summation */
@@ -139,8 +142,7 @@
n = 1.0;
t = 1.0;
maxt = 0.0;
-pcanc = 1.0; /* estimate 100% error if problem */
-*err = pcanc;
+*err = 1.0;
while( t > MACHEP )
@@ -160,7 +162,8 @@
temp = fabs(u);
if( (temp > 1.0 ) && (maxt > (MAXNUM/temp)) )
{
- goto blowup;
+ *err = 1.0; /* blowup: estimate 100% error */
+ return sum;
}
a0 *= u;
@@ -186,12 +189,8 @@
if( sum != 0.0 )
maxt /= fabs(sum);
maxt *= MACHEP; /* this way avoids multiply overflow */
-pcanc = fabs( MACHEP * n + maxt );
+*err = fabs( MACHEP * n + maxt );
-blowup:
-
-*err = pcanc;
-
return( sum );
}
@@ -218,9 +217,6 @@
double *err;
{
double h1, h2, t, u, temp, acanc, asum, err1, err2;
-#ifndef ANSIPROT
-double exp(), log(), Gamma(), lgam(), fabs(), hyp2f0();
-#endif
if( x == 0 )
{
@@ -241,14 +237,14 @@
h1 = hyp2f0( a, a-b+1, -1.0/x, 1, &err1 );
-temp = exp(u) / Gamma(b-a);
+temp = exp(u) / gamma(b-a);
h1 *= temp;
err1 *= temp;
h2 = hyp2f0( b-a, 1.0-a, 1.0/x, 2, &err2 );
if( a < 0 )
- temp = exp(t) / Gamma(a);
+ temp = exp(t) / gamma(a);
else
temp = exp( t - lgam(a) );
@@ -265,7 +261,7 @@
if( b < 0 )
{
- temp = Gamma(b);
+ temp = gamma(b);
asum *= temp;
acanc *= fabs(temp);
}
@@ -291,9 +287,6 @@
int type; /* determines what converging factor to use */
double *err;
{
-#ifndef ANSIPROT
-double fabs();
-#endif
double a0, alast, t, tlast, maxt;
double n, an, bn, u, sum, temp;
Modified: trunk/Lib/special/cephes/i0.c
===================================================================
--- trunk/Lib/special/cephes/i0.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/i0.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -70,13 +70,11 @@
/*
-Cephes Math Library Release 2.0: April, 1987
-Copyright 1984, 1987 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
-#include <stdio.h>
/* Chebyshev coefficients for exp(-x) I0(x)
* in the interval [0,8].
@@ -353,7 +351,11 @@
};
#endif
-#ifndef ANSIPROT
+#ifdef ANSIPROT
+extern double chbevl ( double, void *, int );
+extern double exp ( double );
+extern double sqrt ( double );
+#else
double chbevl(), exp(), sqrt();
#endif
Modified: trunk/Lib/special/cephes/i1.c
===================================================================
--- trunk/Lib/special/cephes/i1.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/i1.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -71,9 +71,8 @@
/*
-Cephes Math Library Release 2.0: March, 1987
-Copyright 1985, 1987 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1985, 1987, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
@@ -351,7 +350,12 @@
#endif
/* i1.c */
-#ifndef ANSIPROT
+#ifdef ANSIPROT
+extern double chbevl ( double, void *, int );
+extern double exp ( double );
+extern double sqrt ( double );
+extern double fabs ( double );
+#else
double chbevl(), exp(), sqrt(), fabs();
#endif
Modified: trunk/Lib/special/cephes/iv.c
===================================================================
--- trunk/Lib/special/cephes/iv.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/iv.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -22,7 +22,7 @@
* function, according to the formula
*
* v -x
- * Iv(x) = (x/2) e hyperg( v+0.5, 2v+1, 2x ) / Gamma(v+1)
+ * Iv(x) = (x/2) e hyperg( v+0.5, 2v+1, 2x ) / gamma(v+1)
*
* If v is a negative integer, then v is replaced by -v.
*
@@ -47,17 +47,23 @@
/*
-Cephes Math Library Release 2.1: November, 1988
-Copyright 1984, 1987, 1988 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
-#ifndef ANSIPROT
-double hyperg(), exp(), Gamma(), log(), fabs(), floor();
+#ifdef ANSIPROT
+extern double hyperg ( double, double, double );
+extern double exp ( double );
+extern double gamma ( double );
+extern double log ( double );
+extern double fabs ( double );
+extern double floor ( double );
+#else
+double hyperg(), exp(), gamma(), log(), fabs(), floor();
#endif
-extern double MACHEP, MAXNUM, NAN;
+extern double MACHEP, MAXNUM;
double iv( v, x )
double v, x;
@@ -82,7 +88,7 @@
if( t != v )
{
mtherr( "iv", DOMAIN );
- return( NAN );
+ return( 0.0 );
}
if( v != 2.0 * floor(v/2.0) )
sign = -1;
@@ -104,7 +110,7 @@
ax = fabs(x);
t = v * log( 0.5 * ax ) - x;
-t = sign * exp(t) / Gamma( v + 1.0 );
+t = sign * exp(t) / gamma( v + 1.0 );
ax = v + 0.5;
return( t * hyperg( ax, 2.0 * ax, 2.0 * x ) );
}
Modified: trunk/Lib/special/cephes/j0.c
===================================================================
--- trunk/Lib/special/cephes/j0.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/j0.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -84,9 +84,8 @@
*/
/*
-Cephes Math Library Release 2.1: January, 1989
-Copyright 1984, 1987, 1989 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
*/
/* Note: all coefficients satisfy the relative error criterion
@@ -459,8 +458,17 @@
};
#endif
-#ifndef ANSIPROT
-double j0(), polevl(), p1evl(), log(), sin(), cos(), sqrt();
+#ifdef ANSIPROT
+extern double polevl ( double, void *, int );
+extern double p1evl ( double, void *, int );
+extern double log ( double );
+extern double sin ( double );
+extern double cos ( double );
+extern double sqrt ( double );
+double j0 ( double );
+#else
+double polevl(), p1evl(), log(), sin(), cos(), sqrt();
+double j0();
#endif
extern double TWOOPI, SQ2OPI, PIO4;
Modified: trunk/Lib/special/cephes/j1.c
===================================================================
--- trunk/Lib/special/cephes/j1.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/j1.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -73,9 +73,8 @@
/*
-Cephes Math Library Release 2.1: January, 1989
-Copyright 1984, 1987, 1989 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
*/
/*
@@ -445,8 +444,17 @@
#define Z2 (*(double *)DZ2)
#endif
-#ifndef ANSIPROT
-double j1(), polevl(), p1evl(), log(), sin(), cos(), sqrt();
+#ifdef ANSIPROT
+extern double polevl ( double, void *, int );
+extern double p1evl ( double, void *, int );
+extern double log ( double );
+extern double sin ( double );
+extern double cos ( double );
+extern double sqrt ( double );
+double j1 ( double );
+#else
+double polevl(), p1evl(), log(), sin(), cos(), sqrt();
+double j1();
#endif
extern double TWOOPI, THPIO4, SQ2OPI;
@@ -457,7 +465,7 @@
w = x;
if( x < 0 )
- return -j1(-x);
+ return -j1(-x);
if( w <= 5.0 )
{
Modified: trunk/Lib/special/cephes/jn.c
===================================================================
--- trunk/Lib/special/cephes/jn.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/jn.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -42,12 +42,15 @@
*/
/* jn.c
-Cephes Math Library Release 2.0: April, 1987
-Copyright 1984, 1987 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
-#ifndef ANSIPROT
+#ifdef ANSIPROT
+extern double fabs ( double );
+extern double j0 ( double );
+extern double j1 ( double );
+#else
double fabs(), j0(), j1();
#endif
extern double MACHEP;
@@ -86,7 +89,7 @@
double y = x*x;
return sign * 0.125 * y * (1 - y / 12.);
} else {
- return( sign * (2.0 * j1(x) / x - j0(x)) );
+ return( sign * (2.0 * j1(x) / x - j0(x)) );
}
}
Modified: trunk/Lib/special/cephes/jv.c
===================================================================
--- trunk/Lib/special/cephes/jv.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/jv.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -42,9 +42,8 @@
/*
-Cephes Math Library Release 2.2: July, 1992
-Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
*/
@@ -57,21 +56,37 @@
#define MAXGAM 171.624376956302725
#endif
-#ifndef ANSIPROT
-double fabs(), floor(), frexp(), polevl(), j0(), j1(), sqrt(), cbrt();
-double exp(), log(), sin(), cos(), acos(), pow(), Gamma(), lgam();
-static double recur(), jvs(), hankel(), jnx(), jnt();
-int airy();
-#else
+#ifdef ANSIPROT
+extern int airy ( double, double *, double *, double *, double * );
+extern double fabs ( double );
+extern double floor ( double );
+extern double frexp ( double, int * );
+extern double polevl ( double, void *, int );
+extern double j0 ( double );
+extern double j1 ( double );
+extern double sqrt ( double );
+extern double cbrt ( double );
+extern double exp ( double );
+extern double log ( double );
+extern double sin ( double );
+extern double cos ( double );
+extern double acos ( double );
+extern double pow ( double, double );
+extern double gamma ( double );
+extern double lgam ( double );
static double recur(double *, double, double *, int);
static double jvs(double, double);
static double hankel(double, double);
static double jnx(double, double);
static double jnt(double, double);
-extern int airy ( double x, double *ai, double *aip, double *bi, double *bip );
+#else
+int airy();
+double fabs(), floor(), frexp(), polevl(), j0(), j1(), sqrt(), cbrt();
+double exp(), log(), sin(), cos(), acos(), pow(), gamma(), lgam();
+static double recur(), jvs(), hankel(), jnx(), jnt();
#endif
-extern double MAXNUM, MACHEP, MINLOG, MAXLOG, NAN;
+extern double MAXNUM, MACHEP, MINLOG, MAXLOG;
#define BIG 1.44115188075855872E+17
double jv( n, x )
@@ -109,7 +124,8 @@
if( (x < 0.0) && (y != an) )
{
mtherr( "Jv", DOMAIN );
- return (NAN);
+ y = 0.0;
+ goto done;
}
y = fabs(x);
@@ -259,7 +275,8 @@
double *newn;
int cancel;
{
-double pkm2, pkm1, pk, pkp1, qkm2, qkm1;
+double pkm2, pkm1, pk, qkm2, qkm1;
+/* double pkp1; */
double k, ans, qk, xk, yk, r, t, kf;
static double big = BIG;
int nflag, ctr;
@@ -357,7 +374,7 @@
do
{
pkm2 = (pkm1 * r - pk * x) / x;
- pkp1 = pk;
+ /* pkp1 = pk; */
pk = pkm1;
pkm1 = pkm2;
r -= 2.0;
@@ -436,9 +453,9 @@
&& (n > 0.0)
&& (n < (MAXGAM-1.0)) )
{
- t = pow( 0.5*x, n ) / Gamma( n + 1.0 );
+ t = pow( 0.5*x, n ) / gamma( n + 1.0 );
#if DEBUG
-printf( "pow(.5*x, %.4e)/Gamma(n+1)=%.5e\n", n, t );
+printf( "pow(.5*x, %.4e)/gamma(n+1)=%.5e\n", n, t );
#endif
y *= t;
}
Modified: trunk/Lib/special/cephes/k0.c
===================================================================
--- trunk/Lib/special/cephes/k0.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/k0.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -70,9 +70,8 @@
*/
/*
-Cephes Math Library Release 2.0: April, 1987
-Copyright 1984, 1987 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
@@ -274,7 +273,13 @@
#endif
/* k0.c */
-#ifndef ANSIPROT
+#ifdef ANSIPROT
+extern double chbevl ( double, void *, int );
+extern double exp ( double );
+extern double i0 ( double );
+extern double log ( double );
+extern double sqrt ( double );
+#else
double chbevl(), exp(), i0(), log(), sqrt();
#endif
extern double PI;
Modified: trunk/Lib/special/cephes/k1.c
===================================================================
--- trunk/Lib/special/cephes/k1.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/k1.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -70,9 +70,8 @@
*/
/*
-Cephes Math Library Release 2.0: April, 1987
-Copyright 1984, 1987 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
@@ -277,7 +276,13 @@
};
#endif
-#ifndef ANSIPROT
+#ifdef ANSIPROT
+extern double chbevl ( double, void *, int );
+extern double exp ( double );
+extern double i1 ( double );
+extern double log ( double );
+extern double sqrt ( double );
+#else
double chbevl(), exp(), i1(), log(), sqrt();
#endif
extern double PI;
Modified: trunk/Lib/special/cephes/kn.c
===================================================================
--- trunk/Lib/special/cephes/kn.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/kn.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -37,10 +37,8 @@
/*
-Cephes Math Library Release 2.1: November, 1988
-Copyright 1984, 1987, 1988 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
*/
@@ -83,7 +81,12 @@
#define EUL 5.772156649015328606065e-1
#define MAXFAC 31
-#ifndef ANSIPROT
+#ifdef ANSIPROT
+extern double fabs ( double );
+extern double exp ( double );
+extern double log ( double );
+extern double sqrt ( double );
+#else
double fabs(), exp(), log(), sqrt();
#endif
extern double MACHEP, MAXNUM, MAXLOG, PI;
Modified: trunk/Lib/special/cephes/mconf.h
===================================================================
--- trunk/Lib/special/cephes/mconf.h 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/mconf.h 2006-07-05 18:27:34 UTC (rev 2041)
@@ -199,3 +199,4 @@
/* Variable for error reporting. See mtherr.c. */
extern int merror;
+#define gamma Gamma
Modified: trunk/Lib/special/cephes/psi.c
===================================================================
--- trunk/Lib/special/cephes/psi.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/psi.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -1,6 +1,6 @@
/* psi.c
*
- * Psi (diGamma) function
+ * Psi (digamma) function
*
*
* SYNOPSIS:
@@ -16,7 +16,7 @@
* psi(x) = -- ln | (x)
* dx
*
- * is the logarithmic derivative of the Gamma function.
+ * is the logarithmic derivative of the gamma function.
* For integer x,
* n-1
* -
@@ -52,9 +52,8 @@
*/
/*
-Cephes Math Library Release 2.2: July, 1992
-Copyright 1984, 1987, 1992 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
@@ -109,7 +108,12 @@
#define EUL 0.57721566490153286061
-#ifndef ANSIPROT
+#ifdef ANSIPROT
+extern double floor ( double );
+extern double log ( double );
+extern double tan ( double );
+extern double polevl ( double, void *, int );
+#else
double floor(), log(), tan(), polevl();
#endif
extern double PI, MAXNUM;
Modified: trunk/Lib/special/cephes/struve.c
===================================================================
--- trunk/Lib/special/cephes/struve.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/struve.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -31,32 +31,28 @@
/*
-Cephes Math Library Release 2.1: January, 1989
-Copyright 1984, 1987, 1989 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.81: June, 2000
+Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
*/
-
-#import "mconf.h"
-
-#define ANSIPROT
+#include "mconf.h"
#define DEBUG 0
-#ifndef ANSIPROT
-double Gamma(), pow(), sqrt(), yn(), yv(), jv(), fabs(), floor();
-double sin(), cos();
+#ifdef ANSIPROT
+extern double gamma ( double );
+extern double pow ( double, double );
+extern double sqrt ( double );
+extern double yn ( int, double );
+extern double jv ( double, double );
+extern double fabs ( double );
+extern double floor ( double );
+extern double sin ( double );
+extern double cos ( double );
+double yv ( double, double );
+double onef2 (double, double, double, double, double * );
+double threef0 (double, double, double, double, double * );
#else
-double onef2( double,double,double,double,double*);
-double threef0( double,double,double,double,double*);
-extern double fabs(double);
-extern double floor(double);
-extern double sqrt(double);
-extern double sin ( double x );
-extern double cos ( double x );
-extern double pow ( double x, double y );
-extern double Gamma ( double x );
-extern double jv ( double n, double x );
-extern double yn ( int n, double x );
-double struve(double, double);
-double yv(double, double);
+double gamma(), pow(), sqrt(), yn(), yv(), jv(), fabs(), floor();
+double sin(), cos();
+double onef2(), threef0();
#endif
static double stop = 1.37e-17;
extern double MACHEP, INFINITY;
@@ -225,12 +221,14 @@
double onef2err, threef0err;
if (x == 0.0) {
- if ((v>-1) || ((floor(v)-v)==0.5)) return 0.0;
- if (v<-1) {
- if ((int)(floor(0.5-v)-1) % 2) return -INFINITY;
- else return INFINITY;
- }
- return 2.0/PI;
+ if (v > -1) {
+ return 0.0;
+ } else if (v < -1) {
+ if ((int)(floor(0.5-v)-1) % 2) return -INFINITY;
+ else return INFINITY;
+ } else {
+ return 2.0/PI;
+ }
}
f = floor(v);
@@ -271,13 +269,13 @@
if( onef2err <= threef0err )
{
- g = Gamma( v + 1.5 );
+ g = gamma( v + 1.5 );
y = y * h * t / ( 0.5 * f * g );
return(y);
}
else
{
- g = Gamma( v + 0.5 );
+ g = gamma( v + 0.5 );
ya = ya * h / ( f * g );
ya = ya + yv( v, x );
return(ya);
Modified: trunk/Lib/special/cephes/yn.c
===================================================================
--- trunk/Lib/special/cephes/yn.c 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/cephes/yn.c 2006-07-05 18:27:34 UTC (rev 2041)
@@ -48,13 +48,16 @@
*/
/*
-Cephes Math Library Release 2.1: December, 1988
-Copyright 1984, 1987 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
-#ifndef ANSIPROT
+#ifdef ANSIPROT
+extern double y0 ( double );
+extern double y1 ( double );
+extern double log ( double );
+#else
double y0(), y1(), log();
#endif
extern double MAXNUM, MAXLOG;
Modified: trunk/Lib/special/tests/test_basic.py
===================================================================
--- trunk/Lib/special/tests/test_basic.py 2006-07-05 15:08:01 UTC (rev 2040)
+++ trunk/Lib/special/tests/test_basic.py 2006-07-05 18:27:34 UTC (rev 2041)
@@ -1155,11 +1155,23 @@
pass
class test_hyp2f1(ScipyTestCase):
-
def check_hyp2f1(self):
- hyp = hyp2f1(1,1,2,.5)
- hrl = -(1/.5)*log(1-.5)
- assert_almost_equal(hyp,hrl,8)
+ # a collection of special cases taken from AMS 55
+ values = [[0.5, 1, 1.5, 0.2**2, 0.5/0.2*log((1+0.2)/(1-0.2))],
+ [0.5, 1, 1.5, -0.2**2, 1./0.2*arctan(0.2)],
+ [1, 1, 2, 0.2, -1/0.2*log(1-0.2)],
+ [3, 3.5, 1.5, 0.2**2,
+ 0.5/0.2/(-5)*((1+0.2)**(-5)-(1-0.2)**(-5))],
+ [-3, 3, 0.5, sin(0.2)**2, cos(2*3*0.2)],
+ [3, 4, 8, 1, gamma(8)*gamma(8-4-3)/gamma(8-3)/gamma(8-4)],
+ [3, 2, 3-2+1, -1, 1./2**3*sqrt(pi)*
+ gamma(1+3-2)/gamma(1+0.5*3-2)/gamma(0.5+0.5*3)],
+ [4, 0.5+4, 5./6+4, 1./9, (0.75)**4*sqrt(pi)*
+ gamma(5./6+2./3*4)/gamma(0.5+4./3)*gamma(5./6+4./3)],
+ ]
+ for i, (a, b, c, x, v) in enumerate(values):
+ cv = hyp2f1(a, b, c, x)
+ assert_almost_equal(cv, v, 8, err_msg='test #%d' % i)
class test_hyp3f0(ScipyTestCase):
@@ -1182,11 +1194,19 @@
assert_array_almost_equal(hypu,hprl,12)
class test_i0(ScipyTestCase):
-
def check_i0(self):
- oiz = i0(.1)
- oizr = iv(0,.1)
- assert_almost_equal(oiz,oizr,8)
+ values = [[0.0, 1.0],
+ [1e-10, 1.0],
+ [0.1, 0.9071009258],
+ [0.5, 0.6450352706],
+ [1.0, 0.4657596077],
+ [2.5, 0.2700464416],
+ [5.0, 0.1835408126],
+ [20.0, 0.0897803119],
+ ]
+ for i, (x, v) in enumerate(values):
+ cv = i0(x) * exp(-x)
+ assert_almost_equal(cv, v, 8, err_msg='test #%d' % i)
class test_i0e(ScipyTestCase):
@@ -1198,11 +1218,18 @@
class test_i1(ScipyTestCase):
def check_i1(self):
+ values = [[0.0, 0.0],
+ [1e-10, 0.4999999999500000e-10],
+ [0.1, 0.0452984468],
+ [0.5, 0.1564208032],
+ [1.0, 0.2079104154],
+ [5.0, 0.1639722669],
+ [20.0, 0.0875062222],
+ ]
+ for i, (x, v) in enumerate(values):
+ cv = i1(x) * exp(-x)
+ assert_almost_equal(cv, v, 8, err_msg='test #%d' % i)
- oi1 = i1(.1)
- oi1r = iv(1,.1)
- assert_almost_equal(oi1,oi1r,8)
-
class test_i1e(ScipyTestCase):
def check_i1e(self):
More information about the Scipy-svn
mailing list