[Scipy-svn] r3144 - trunk/Lib/special/cephes

[scipy-svn@scip...]

Author: cookedm
Date: 2007-07-03 12:41:27 -0500 (Tue, 03 Jul 2007)
New Revision: 3144

Modified:
   trunk/Lib/special/cephes/j0.c
   trunk/Lib/special/cephes/j1.c
   trunk/Lib/special/cephes/k0.c
   trunk/Lib/special/cephes/k1.c
   trunk/Lib/special/cephes/kn.c
   trunk/Lib/special/cephes/yn.c
Log:
special: More complete handling of domains/singularities for real argument to Bessel functions


Modified: trunk/Lib/special/cephes/j0.c
===================================================================
--- trunk/Lib/special/cephes/j0.c	2007-07-03 11:29:01 UTC (rev 3143)
+++ trunk/Lib/special/cephes/j0.c	2007-07-03 17:41:27 UTC (rev 3144)
@@ -513,7 +513,7 @@
 #define PIO4 .78539816339744830962
 #define SQ2OPI .79788456080286535588
 */
-extern double MAXNUM;
+extern double INFINITY, NAN;
 
 double y0(x)
 double x;
@@ -522,11 +522,13 @@
 
 if( x <= 5.0 )
 	{
-	if( x <= 0.0 )
-		{
-		mtherr( "y0", DOMAIN );
-		return( -MAXNUM );
-		}
+	if (x == 0.0) {
+		mtherr("y0", SING);
+		return -INFINITY;
+	} else if (x < 0.0) {
+		mtherr("y0", DOMAIN);
+		return NAN;
+	}
 	z = x * x;
 	w = polevl( z, YP, 7) / p1evl( z, YQ, 7 );
 	w += TWOOPI * log(x) * j0(x);

Modified: trunk/Lib/special/cephes/j1.c
===================================================================
--- trunk/Lib/special/cephes/j1.c	2007-07-03 11:29:01 UTC (rev 3143)
+++ trunk/Lib/special/cephes/j1.c	2007-07-03 17:41:27 UTC (rev 3144)
@@ -485,7 +485,7 @@
 }
 
 
-extern double MAXNUM;
+extern double INFINITY, NAN;
 
 double y1(x)
 double x;
@@ -494,11 +494,13 @@
 
 if( x <= 5.0 )
 	{
-	if( x <= 0.0 )
-		{
-		mtherr( "y1", DOMAIN );
-		return( -MAXNUM );
-		}
+	if (x == 0.0) {
+		mtherr("y1", SING);
+		return -INFINITY;
+	} else if (x <= 0.0) {
+		mtherr("y1", DOMAIN);
+		return NAN;
+	}
 	z = x * x;
 	w = x * (polevl( z, YP, 5 ) / p1evl( z, YQ, 8 ));
 	w += TWOOPI * ( j1(x) * log(x)  -  1.0/x );

Modified: trunk/Lib/special/cephes/k0.c
===================================================================
--- trunk/Lib/special/cephes/k0.c	2007-07-03 11:29:01 UTC (rev 3143)
+++ trunk/Lib/special/cephes/k0.c	2007-07-03 17:41:27 UTC (rev 3144)
@@ -283,18 +283,20 @@
 double chbevl(), exp(), i0(), log(), sqrt();
 #endif
 extern double PI;
-extern double MAXNUM;
+extern double INFINITY, NAN;
 
 double k0(x)
 double x;
 {
 double y, z;
 
-if( x <= 0.0 )
-	{
-	mtherr( "k0", DOMAIN );
-	return( MAXNUM );
-	}
+if (x == 0.0) {
+	mtherr("k0", SING);
+	return INFINITY;
+} else if (x < 0.0) {
+	mtherr("k0", DOMAIN);
+	return NAN;
+}
 
 if( x <= 2.0 )
 	{
@@ -315,11 +317,13 @@
 {
 double y;
 
-if( x <= 0.0 )
-	{
+if (x == 0.0) {
+	mtherr("k0e", SING);
+	return INFINITY;
+} else if (x < 0.0) {
 	mtherr( "k0e", DOMAIN );
-	return( MAXNUM );
-	}
+	return NAN;
+}
 
 if( x <= 2.0 )
 	{

Modified: trunk/Lib/special/cephes/k1.c
===================================================================
--- trunk/Lib/special/cephes/k1.c	2007-07-03 11:29:01 UTC (rev 3143)
+++ trunk/Lib/special/cephes/k1.c	2007-07-03 17:41:27 UTC (rev 3144)
@@ -286,19 +286,21 @@
 double chbevl(), exp(), i1(), log(), sqrt();
 #endif
 extern double PI;
-extern double MINLOG, MAXNUM;
+extern double MINLOG, INFINITY, NAN;
 
 double k1(x)
 double x;
 {
 double y, z;
 
+if (x == 0.0) {
+	mtherr("k1", SING);
+	return INFINITY;
+} else if (x < 0.0) {
+	mtherr("k1", DOMAIN);
+	return NAN;
+}
 z = 0.5 * x;
-if( z <= 0.0 )
-	{
-	mtherr( "k1", DOMAIN );
-	return( MAXNUM );
-	}
 
 if( x <= 2.0 )
 	{
@@ -318,11 +320,13 @@
 {
 double y;
 
-if( x <= 0.0 )
-	{
-	mtherr( "k1e", DOMAIN );
-	return( MAXNUM );
-	}
+if (x == 0.0) {
+	mtherr("k1e", SING);
+	return INFINITY;
+} else if (x < 0.0) {
+	mtherr("k1e", DOMAIN);
+	return NAN;
+}
 
 if( x <= 2.0 )
 	{

Modified: trunk/Lib/special/cephes/kn.c
===================================================================
--- trunk/Lib/special/cephes/kn.c	2007-07-03 11:29:01 UTC (rev 3143)
+++ trunk/Lib/special/cephes/kn.c	2007-07-03 17:41:27 UTC (rev 3144)
@@ -89,7 +89,7 @@
 #else
 double fabs(), exp(), log(), sqrt();
 #endif
-extern double MACHEP, MAXNUM, MAXLOG, PI;
+extern double MACHEP, MAXNUM, MAXLOG, PI, INFINITY, NAN;
 
 double kn( nn, x )
 int nn;
@@ -111,14 +111,15 @@
 	return( MAXNUM );
 	}
 
-if( x <= 0.0 )
-	{
-	if( x < 0.0 )
-		mtherr( "kn", DOMAIN );
-	else
-		mtherr( "kn", SING );
-	return( MAXNUM );
+if(x <= 0.0) {
+	if( x < 0.0 ) {
+		mtherr("kn", DOMAIN);
+                return NAN;
+	} else {
+		mtherr("kn", SING);
+		return INFINITY;
 	}
+}
 
 
 if( x > 9.55 )

Modified: trunk/Lib/special/cephes/yn.c
===================================================================
--- trunk/Lib/special/cephes/yn.c	2007-07-03 11:29:01 UTC (rev 3143)
+++ trunk/Lib/special/cephes/yn.c	2007-07-03 17:41:27 UTC (rev 3144)
@@ -91,17 +91,12 @@
 
 /* test for overflow */
 if (x == 0.0) {
-    return -INFINITY;
+	mtherr("yn", SING);
+	return -INFINITY;
+} else if (x < 0.0) {
+	mtherr("yn", DOMAIN);
+        return NAN;
 }
-else if( x < 0.0 )
-	{
-	mtherr( "yn", SING );
-#ifdef NANS
-        return (NAN);
-#else
-	return( -MAXNUM );
-#endif
-	}
 
 /* forward recurrence on n */