Date: Tue, 10 Mar 2015 17:10:55 +0000 (UTC) From: Steve Kargl <kargl@FreeBSD.org> To: src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-head@freebsd.org Subject: svn commit: r279856 - in head/lib/msun: man src Message-ID: <201503101710.t2AHAt0V028926@svn.freebsd.org>
next in thread | raw e-mail | index | archive | help
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;
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201503101710.t2AHAt0V028926>