svn commit: r279856 - in head/lib/msun: man src
Steve Kargl
kargl at FreeBSD.org
Tue Mar 10 17:10:57 UTC 2015
Author: kargl
Date: Tue Mar 10 17:10:54 2015
New Revision: 279856
URL: https://svnweb.freebsd.org/changeset/base/279856
Log:
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)
Modified:
head/lib/msun/man/j0.3
head/lib/msun/src/e_j0.c
head/lib/msun/src/e_j0f.c
head/lib/msun/src/e_j1.c
head/lib/msun/src/e_j1f.c
head/lib/msun/src/e_jn.c
head/lib/msun/src/e_jnf.c
Modified: head/lib/msun/man/j0.3
==============================================================================
--- head/lib/msun/man/j0.3 Tue Mar 10 17:04:11 2015 (r279855)
+++ head/lib/msun/man/j0.3 Tue Mar 10 17:10:54 2015 (r279856)
@@ -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
Modified: head/lib/msun/src/e_j0.c
==============================================================================
--- head/lib/msun/src/e_j0.c Tue Mar 10 17:04:11 2015 (r279855)
+++ head/lib/msun/src/e_j0.c Tue Mar 10 17:10:54 2015 (r279856)
@@ -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
Modified: head/lib/msun/src/e_j0f.c
==============================================================================
--- head/lib/msun/src/e_j0f.c Tue Mar 10 17:04:11 2015 (r279855)
+++ head/lib/msun/src/e_j0f.c Tue Mar 10 17:10:54 2015 (r279856)
@@ -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
Modified: head/lib/msun/src/e_j1.c
==============================================================================
--- head/lib/msun/src/e_j1.c Tue Mar 10 17:04:11 2015 (r279855)
+++ head/lib/msun/src/e_j1.c Tue Mar 10 17:10:54 2015 (r279856)
@@ -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);
Modified: head/lib/msun/src/e_j1f.c
==============================================================================
--- head/lib/msun/src/e_j1f.c Tue Mar 10 17:04:11 2015 (r279855)
+++ head/lib/msun/src/e_j1f.c Tue Mar 10 17:10:54 2015 (r279856)
@@ -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);
Modified: head/lib/msun/src/e_jn.c
==============================================================================
--- head/lib/msun/src/e_jn.c Tue Mar 10 17:04:11 2015 (r279855)
+++ head/lib/msun/src/e_jn.c Tue Mar 10 17:10:54 2015 (r279856)
@@ -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;
Modified: head/lib/msun/src/e_jnf.c
==============================================================================
--- head/lib/msun/src/e_jnf.c Tue Mar 10 17:04:11 2015 (r279855)
+++ head/lib/msun/src/e_jnf.c Tue Mar 10 17:10:54 2015 (r279856)
@@ -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;
More information about the svn-src-all
mailing list