Date: Wed, 11 Jul 2018 12:12:50 +0000 (UTC) From: Mark Johnston <markj@FreeBSD.org> To: src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-stable@freebsd.org, svn-src-stable-11@freebsd.org Subject: svn commit: r336196 - stable/11/lib/msun/src Message-ID: <201807111212.w6BCCox6074366@repo.freebsd.org>
next in thread | raw e-mail | index | archive | help
Author: markj Date: Wed Jul 11 12:12:49 2018 New Revision: 336196 URL: https://svnweb.freebsd.org/changeset/base/336196 Log: MFC r336089: Fix whitespace issues in bessel function routines. PR: 229423 Modified: stable/11/lib/msun/src/e_j0.c stable/11/lib/msun/src/e_j1.c stable/11/lib/msun/src/e_j1f.c stable/11/lib/msun/src/e_jn.c Directory Properties: stable/11/ (props changed) Modified: stable/11/lib/msun/src/e_j0.c ============================================================================== --- stable/11/lib/msun/src/e_j0.c Wed Jul 11 09:41:50 2018 (r336195) +++ stable/11/lib/msun/src/e_j0.c Wed Jul 11 12:12:49 2018 (r336196) @@ -1,4 +1,3 @@ - /* @(#)e_j0.c 1.3 95/01/18 */ /* * ==================================================== @@ -6,7 +5,7 @@ * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -33,20 +32,20 @@ __FBSDID("$FreeBSD$"); * (To avoid cancellation, use * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) * to compute the worse one.) - * + * * 3 Special cases * j0(nan)= nan * j0(0) = 1 * j0(inf) = 0 - * + * * Method -- y0(x): * 1. For x<2. - * Since + * Since * y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...) * therefore y0(x)-2/pi*j0(x)*ln(x) is an even function. * We use the following function to approximate y0, * y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2 - * where + * where * U(z) = u00 + u01*z + ... + u06*z^6 * V(z) = 1 + v01*z + ... + v04*z^4 * with absolute approximation error bounded by 2**-72. @@ -71,7 +70,7 @@ huge = 1e300, one = 1.0, invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */ tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ - /* R0/S0 on [0, 2.00] */ +/* R0/S0 on [0, 2.00] */ R02 = 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */ R03 = -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */ R04 = 1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */ @@ -157,7 +156,7 @@ __ieee754_y0(double x) * y0(Inf) = 0. * y0(-Inf) = NaN and raise invalid exception. */ - if(ix>=0x7ff00000) return vone/(x+x*x); + 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. */ @@ -293,7 +292,7 @@ pzero(double x) s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); return one+ r/s; } - + /* For x >= 8, the asymptotic expansions of qzero is * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. Modified: stable/11/lib/msun/src/e_j1.c ============================================================================== --- stable/11/lib/msun/src/e_j1.c Wed Jul 11 09:41:50 2018 (r336195) +++ stable/11/lib/msun/src/e_j1.c Wed Jul 11 12:12:49 2018 (r336196) @@ -1,4 +1,3 @@ - /* @(#)e_j1.c 1.3 95/01/18 */ /* * ==================================================== @@ -6,7 +5,7 @@ * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -34,16 +33,16 @@ __FBSDID("$FreeBSD$"); * (To avoid cancellation, use * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) * to compute the worse one.) - * + * * 3 Special cases * j1(nan)= nan * j1(0) = 0 * j1(inf) = 0 - * + * * Method -- y1(x): - * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN + * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN * 2. For x<2. - * Since + * Since * y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...) * therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function. * We use the following function to approximate y1, @@ -154,7 +153,7 @@ __ieee754_y1(double x) * y1(Inf) = 0. * y1(-Inf) = NaN and raise invalid exception. */ - if(ix>=0x7ff00000) return vone/(x+x*x); + 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. */ @@ -186,10 +185,10 @@ __ieee754_y1(double x) z = invsqrtpi*(u*ss+v*cc)/sqrt(x); } return z; - } + } if(ix<=0x3c900000) { /* x < 2**-54 */ return(-tpi/x); - } + } z = x*x; u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); @@ -287,7 +286,7 @@ pone(double x) s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); return one+ r/s; } - + /* For x >= 8, the asymptotic expansions of qone is * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. Modified: stable/11/lib/msun/src/e_j1f.c ============================================================================== --- stable/11/lib/msun/src/e_j1f.c Wed Jul 11 09:41:50 2018 (r336195) +++ stable/11/lib/msun/src/e_j1f.c Wed Jul 11 12:12:49 2018 (r336196) @@ -32,7 +32,7 @@ huge = 1e30, one = 1.0, invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */ tpi = 6.3661974669e-01, /* 0x3f22f983 */ - /* R0/S0 on [0,2] */ +/* R0/S0 on [0,2] */ r00 = -6.2500000000e-02, /* 0xbd800000 */ r01 = 1.4070566976e-03, /* 0x3ab86cfd */ r02 = -1.5995563444e-05, /* 0xb7862e36 */ Modified: stable/11/lib/msun/src/e_jn.c ============================================================================== --- stable/11/lib/msun/src/e_jn.c Wed Jul 11 09:41:50 2018 (r336195) +++ stable/11/lib/msun/src/e_jn.c Wed Jul 11 12:12:49 2018 (r336196) @@ -1,4 +1,3 @@ - /* @(#)e_jn.c 1.4 95/01/18 */ /* * ==================================================== @@ -6,7 +5,7 @@ * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -18,7 +17,7 @@ __FBSDID("$FreeBSD$"); * __ieee754_jn(n, x), __ieee754_yn(n, x) * floating point Bessel's function of the 1st and 2nd kind * of order n - * + * * Special cases: * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal; * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. @@ -37,7 +36,6 @@ __FBSDID("$FreeBSD$"); * yn(n,x) is similar in all respects, except * that forward recursion is used for all * values of n>1. - * */ #include "math.h" @@ -66,7 +64,7 @@ __ieee754_jn(int n, double x) ix = 0x7fffffff&hx; /* if J(n,NaN) is NaN */ if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x; - if(n<0){ + if(n<0){ n = -n; x = -x; hx ^= 0x80000000; @@ -77,14 +75,14 @@ __ieee754_jn(int n, double x) x = fabs(x); if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */ b = zero; - else if((double)n<=x) { + else if((double)n<=x) { /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ if(ix>=0x52D00000) { /* x > 2**302 */ - /* (x >> n**2) + /* (x >> n**2) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Let s=sin(x), c=cos(x), - * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then + * Let s=sin(x), c=cos(x), + * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2), then * * n sin(xn)*sqt2 cos(xn)*sqt2 * ---------------------------------- @@ -100,7 +98,7 @@ __ieee754_jn(int n, double x) case 3: temp = cos(x)-sin(x); break; } b = invsqrtpi*temp/sqrt(x); - } else { + } else { a = __ieee754_j0(x); b = __ieee754_j1(x); for(i=1;i<n;i++){ @@ -111,7 +109,7 @@ __ieee754_jn(int n, double x) } } else { if(ix<0x3e100000) { /* x < 2**-29 */ - /* x is tiny, return the first Taylor expansion of J(n,x) + /* x is tiny, return the first Taylor expansion of J(n,x) * J(n,x) = 1/n!*(x/2)^n - ... */ if(n>33) /* underflow */ @@ -126,14 +124,14 @@ __ieee754_jn(int n, double x) } } else { /* use backward recurrence */ - /* x x^2 x^2 + /* x x^2 x^2 * J(n,x)/J(n-1,x) = ---- ------ ------ ..... * 2n - 2(n+1) - 2(n+2) * - * 1 1 1 + * 1 1 1 * (for large x) = ---- ------ ------ ..... * 2n 2(n+1) 2(n+2) - * -- - ------ - ------ - + * -- - ------ - ------ - * x x x * * Let w = 2n/x and h=2/x, then the above quotient @@ -149,9 +147,9 @@ __ieee754_jn(int n, double x) * To determine how many terms needed, let * Q(0) = w, Q(1) = w(w+h) - 1, * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), - * When Q(k) > 1e4 good for single - * When Q(k) > 1e9 good for double - * When Q(k) > 1e17 good for quadruple + * When Q(k) > 1e4 good for single + * When Q(k) > 1e9 good for double + * When Q(k) > 1e17 good for quadruple */ /* determine k */ double t,v; @@ -237,11 +235,11 @@ __ieee754_yn(int n, double x) if(n==1) return(sign*__ieee754_y1(x)); if(ix==0x7ff00000) return zero; if(ix>=0x52D00000) { /* x > 2**302 */ - /* (x >> n**2) + /* (x >> n**2) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Let s=sin(x), c=cos(x), - * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then + * Let s=sin(x), c=cos(x), + * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2), then * * n sin(xn)*sqt2 cos(xn)*sqt2 * ----------------------------------
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201807111212.w6BCCox6074366>