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)

6 files changed, 57 insertions, 22 deletions

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;