Attachment 194758 Details for Bug 229423 – whitespace in bessel functions

[patch] whitespace in bessel functions

z_bessel_whitespace.diff (text/plain), 8.57 KB, created by sgk on 2018-06-29 22:11:10 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: sgk

Created: 2018-06-29 22:11:10 UTC

Size: 8.57 KB

patch

obsolete

>Index: src/e_j0.c
>===================================================================
>--- src/e_j0.c	(revision 2101)
>+++ src/e_j0.c	(working copy)
>@@ -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 @@
>  * 	   (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 @@
> 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 @@
> 	 * 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,8 +292,8 @@
> 	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.
>  * We approximate pzero by
>Index: src/e_j1.c
>===================================================================
>--- src/e_j1.c	(revision 2101)
>+++ src/e_j1.c	(working copy)
>@@ -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 @@
>  * 	   (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 @@
> 	 * 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 @@
>                     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,8 +286,8 @@
>         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.
>  * We approximate pone by
>Index: src/e_j1f.c
>===================================================================
>--- src/e_j1f.c	(revision 2101)
>+++ src/e_j1f.c	(working copy)
>@@ -32,7 +32,7 @@
> 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 */
>Index: src/e_jn.c
>===================================================================
>--- src/e_jn.c	(revision 2101)
>+++ src/e_jn.c	(working copy)
>@@ -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 @@
>  * __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 @@
>  *	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 @@
> 	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 @@
> 	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 @@
> 		    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 @@
> 	    }
> 	} 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 @@
> 		}
> 	    } 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 @@
> 		 * 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 @@
> 	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
>      *		----------------------------------

Actions: View | Diff

Attachments on bug 229423: 194758