d/freebsd-nq
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

Browse Source 
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)
This commit is contained in:
Steve Kargl 2015-03-10 17:10:54 +00:00
parent 25023a8a88
commit 186f620727
7 changed files with 77 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
lib/msun/man/j0.3
View File

@ -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 0 or a sufficiently small positive number, these routines
is \*(Pm0, these routines
will generate a divide-by-zero exception and return -\*(If.
If
.Fa x
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

16
lib/msun/src/e_j0.c
View File

@ -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

13
lib/msun/src/e_j0f.c
View File

@ -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

16
lib/msun/src/e_j1.c
View File

@ -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);

13
lib/msun/src/e_j1f.c
View File

@ -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);

10
lib/msun/src/e_jn.c
View File

@ -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;

11
lib/msun/src/e_jnf.c
View File

@ -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;