According to POSIX.1-2008, the Bessel functions of second kind

should raise a divide-by-zero floating point exception for x = +-0
and an invalid floating point exception for x < 0 including x = -Inf.
Update the code to raise the exception and update the documentation
with hopefully better description of the behavior.

Reviewed by:	bde (code only)

7 files changed, 77 insertions, 47 deletions

diff --git a/lib/msun/man/j0.3 b/lib/msun/man/j0.3
index 91849da..7e1b790 100644
--- a/lib/msun/man/j0.3
+++ b/lib/msun/man/j0.3
@@ -28,7 +28,7 @@
 .\"     from: @(#)j0.3	6.7 (Berkeley) 4/19/91
 .\" $FreeBSD$
 .\"
-.Dd February 18, 2008
+.Dd March 10, 2015
 .Dt J0 3
 .Os
 .Sh NAME
@@ -77,24 +77,17 @@
 The functions
 .Fn j0 ,
 .Fn j0f ,
-.Fn j1
+.Fn j1 ,
 and
 .Fn j1f
-compute the
-.Em Bessel function of the first kind of the order
-0 and the
-.Em order
-1, respectively,
-for the
-real value
+compute the Bessel function of the first kind of orders
+0 and 1 for the real value
 .Fa x ;
 the functions
 .Fn jn
 and
 .Fn jnf
-compute the
-.Em Bessel function of the first kind of the integer
-.Em order
+compute the Bessel function of the first kind of the integer order
 .Fa n
 for the real value
 .Fa x .
@@ -105,13 +98,8 @@ The functions
 .Fn y1 ,
 and
 .Fn y1f
-compute the linearly independent
-.Em Bessel function of the second kind of the order
-0 and the
-.Em order
-1, respectively,
-for the
-positive
+compute the linearly independent Bessel function of the second kind 
+of orders 0 and 1 for the positive
 .Em real
 value
 .Fa x ;
@@ -119,9 +107,7 @@ the functions
 .Fn yn
 and
 .Fn ynf
-compute the
-.Em Bessel function of the second kind for the integer
-.Em order
+compute the Bessel function of the second kind for the integer order
 .Fa n
 for the positive
 .Em real
@@ -141,11 +127,20 @@ and
 .Fn ynf .
 If
 .Fa x
-is negative, these routines will generate an invalid exception and
-return \*(Na.
+is negative, including -\*(If, these routines will generate an invalid
+exception and return \*(Na.
+If
+.Fa x
+is \*(Pm0, these routines
+will generate a divide-by-zero exception and return -\*(If.
 If
 .Fa x
-is 0 or a sufficiently small positive number, these routines
+is a sufficiently small positive number, then 
+.Fn y1 ,
+.Fn y1f ,
+.Fn yn ,
+and
+.Fn ynf 
 will generate an overflow exception and return -\*(If.
 .Sh SEE ALSO
 .Xr math 3
diff --git a/lib/msun/src/e_j0.c b/lib/msun/src/e_j0.c
index 9e269e7..365ffe5 100644
--- a/lib/msun/src/e_j0.c
+++ b/lib/msun/src/e_j0.c
@@ -64,6 +64,8 @@ __FBSDID("$FreeBSD$");
 
 static double pzero(double), qzero(double);
 
+static const volatile double vone = 1, vzero = 0;
+
 static const double
 huge 	= 1e300,
 one	= 1.0,
@@ -150,10 +152,16 @@ __ieee754_y0(double x)
 
 	EXTRACT_WORDS(hx,lx,x);
         ix = 0x7fffffff&hx;
-    /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
-	if(ix>=0x7ff00000) return  one/(x+x*x); 
-        if((ix|lx)==0) return -one/zero;
-        if(hx<0) return zero/zero;
+	/*
+	 * y0(NaN) = NaN.
+	 * y0(Inf) = 0.
+	 * y0(-Inf) = NaN and raise invalid exception.
+	 */
+	if(ix>=0x7ff00000) return vone/(x+x*x); 
+	/* y0(+-0) = -inf and raise divide-by-zero exception. */
+	if((ix|lx)==0) return -one/vzero;
+	/* y0(x<0) = NaN and raise invalid exception. */
+	if(hx<0) return vzero/vzero;
         if(ix >= 0x40000000) {  /* |x| >= 2.0 */
         /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
          * where x0 = x-pi/4
diff --git a/lib/msun/src/e_j0f.c b/lib/msun/src/e_j0f.c
index 0479957..e86ed64 100644
--- a/lib/msun/src/e_j0f.c
+++ b/lib/msun/src/e_j0f.c
@@ -16,11 +16,17 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
+/*
+ * See e_j0.c for complete comments.
+ */
+
 #include "math.h"
 #include "math_private.h"
 
 static float pzerof(float), qzerof(float);
 
+static const volatile float vone = 1,  vzero = 0;
+
 static const float
 huge 	= 1e30,
 one	= 1.0,
@@ -107,10 +113,9 @@ __ieee754_y0f(float x)
 
 	GET_FLOAT_WORD(hx,x);
         ix = 0x7fffffff&hx;
-    /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
-	if(ix>=0x7f800000) return  one/(x+x*x);
-        if(ix==0) return -one/zero;
-        if(hx<0) return zero/zero;
+	if(ix>=0x7f800000) return  vone/(x+x*x);
+	if(ix==0) return -one/vzero;
+	if(hx<0) return vzero/vzero;
         if(ix >= 0x40000000) {  /* |x| >= 2.0 */
         /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
          * where x0 = x-pi/4
diff --git a/lib/msun/src/e_j1.c b/lib/msun/src/e_j1.c
index 2dc0ba1..4922058 100644
--- a/lib/msun/src/e_j1.c
+++ b/lib/msun/src/e_j1.c
@@ -64,6 +64,8 @@ __FBSDID("$FreeBSD$");
 
 static double pone(double), qone(double);
 
+static const volatile double vone = 1, vzero = 0;
+
 static const double
 huge    = 1e300,
 one	= 1.0,
@@ -147,10 +149,16 @@ __ieee754_y1(double x)
 
 	EXTRACT_WORDS(hx,lx,x);
         ix = 0x7fffffff&hx;
-    /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
-	if(ix>=0x7ff00000) return  one/(x+x*x); 
-        if((ix|lx)==0) return -one/zero;
-        if(hx<0) return zero/zero;
+	/*
+	 * y1(NaN) = NaN.
+	 * y1(Inf) = 0.
+	 * y1(-Inf) = NaN and raise invalid exception.
+	 */
+	if(ix>=0x7ff00000) return  vone/(x+x*x); 
+	/* y1(+-0) = -inf and raise divide-by-zero exception. */
+        if((ix|lx)==0) return -one/vzero;
+	/* y1(x<0) = NaN and raise invalid exception. */
+        if(hx<0) return vzero/vzero;
         if(ix >= 0x40000000) {  /* |x| >= 2.0 */
                 s = sin(x);
                 c = cos(x);
diff --git a/lib/msun/src/e_j1f.c b/lib/msun/src/e_j1f.c
index 6077a69..c39c548 100644
--- a/lib/msun/src/e_j1f.c
+++ b/lib/msun/src/e_j1f.c
@@ -16,11 +16,17 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
+/*
+ * See e_j1.c for complete comments.
+ */
+
 #include "math.h"
 #include "math_private.h"
 
 static float ponef(float), qonef(float);
 
+static const volatile float vone = 1, vzero = 0;
+
 static const float
 huge    = 1e30,
 one	= 1.0,
@@ -104,10 +110,9 @@ __ieee754_y1f(float x)
 
 	GET_FLOAT_WORD(hx,x);
         ix = 0x7fffffff&hx;
-    /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
-	if(ix>=0x7f800000) return  one/(x+x*x);
-        if(ix==0) return -one/zero;
-        if(hx<0) return zero/zero;
+	if(ix>=0x7f800000) return  vone/(x+x*x);
+	if(ix==0) return -one/vzero;
+	if(hx<0) return vzero/vzero;
         if(ix >= 0x40000000) {  /* |x| >= 2.0 */
                 s = sinf(x);
                 c = cosf(x);
diff --git a/lib/msun/src/e_jn.c b/lib/msun/src/e_jn.c
index 8b0bc62..eefd4ff 100644
--- a/lib/msun/src/e_jn.c
+++ b/lib/msun/src/e_jn.c
@@ -43,6 +43,8 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
+static const volatile double vone = 1, vzero = 0;
+
 static const double
 invsqrtpi=  5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
 two   =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
@@ -220,10 +222,12 @@ __ieee754_yn(int n, double x)
 
 	EXTRACT_WORDS(hx,lx,x);
 	ix = 0x7fffffff&hx;
-    /* if Y(n,NaN) is NaN */
+	/* yn(n,NaN) = NaN */
 	if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
-	if((ix|lx)==0) return -one/zero;
-	if(hx<0) return zero/zero;
+	/* yn(n,+-0) = -inf and raise divide-by-zero exception. */
+	if((ix|lx)==0) return -one/vzero;
+	/* yn(n,x<0) = NaN and raise invalid exception. */
+	if(hx<0) return vzero/vzero;
 	sign = 1;
 	if(n<0){
 		n = -n;
diff --git a/lib/msun/src/e_jnf.c b/lib/msun/src/e_jnf.c
index f564aec..965feeb 100644
--- a/lib/msun/src/e_jnf.c
+++ b/lib/msun/src/e_jnf.c
@@ -16,9 +16,15 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
+/*
+ * See e_jn.c for complete comments.
+ */
+
 #include "math.h"
 #include "math_private.h"
 
+static const volatile float vone = 1, vzero = 0;
+
 static const float
 two   =  2.0000000000e+00, /* 0x40000000 */
 one   =  1.0000000000e+00; /* 0x3F800000 */
@@ -172,10 +178,9 @@ __ieee754_ynf(int n, float x)
 
 	GET_FLOAT_WORD(hx,x);
 	ix = 0x7fffffff&hx;
-    /* if Y(n,NaN) is NaN */
 	if(ix>0x7f800000) return x+x;
-	if(ix==0) return -one/zero;
-	if(hx<0) return zero/zero;
+	if(ix==0) return -one/vzero;
+	if(hx<0) return vzero/vzero;
 	sign = 1;
 	if(n<0){
 		n = -n;