Attachment #194758 for bug #229423

View | Details | Raw Unified | Return to bug 229423
Collapse All | Expand All

Lines 1-4 Link Here

(-)src/e_j0.c (-9 / +8 lines)
1
2	/* @(#)e_j0.c 1.3 95/01/18 */	1	/* @(#)e_j0.c 1.3 95/01/18 */
3	/*	2	/*
4	* ====================================================	3	* ====================================================
Lines 6-12 Link Here
6	*	5	*
7	* Developed at SunSoft, a Sun Microsystems, Inc. business.	6	* Developed at SunSoft, a Sun Microsystems, Inc. business.
8	* Permission to use, copy, modify, and distribute this	7	* Permission to use, copy, modify, and distribute this
9	* software is freely granted, provided that this notice	8	* software is freely granted, provided that this notice
10	* is preserved.	9	* is preserved.
11	* ====================================================	10	* ====================================================
12	*/	11	*/
Lines 33-52 Link Here
33	* (To avoid cancellation, use	32	* (To avoid cancellation, use
34	* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))	33	* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
35	* to compute the worse one.)	34	* to compute the worse one.)
36	*	35	*
37	* 3 Special cases	36	* 3 Special cases
38	* j0(nan)= nan	37	* j0(nan)= nan
39	* j0(0) = 1	38	* j0(0) = 1
40	* j0(inf) = 0	39	* j0(inf) = 0
41	*	40	*
42	* Method -- y0(x):	41	* Method -- y0(x):
43	* 1. For x<2.	42	* 1. For x<2.
44	* Since	43	* Since
45	* y0(x) = 2/pi(j0(x)(ln(x/2)+Euler) + x^2/4 - ...)	44	* y0(x) = 2/pi(j0(x)(ln(x/2)+Euler) + x^2/4 - ...)
46	* therefore y0(x)-2/pij0(x)ln(x) is an even function.	45	* therefore y0(x)-2/pij0(x)ln(x) is an even function.
47	* We use the following function to approximate y0,	46	* We use the following function to approximate y0,
48	* y0(x) = U(z)/V(z) + (2/pi)(j0(x)ln(x)), z= x^2	47	* y0(x) = U(z)/V(z) + (2/pi)(j0(x)ln(x)), z= x^2
49	* where	48	* where
50	* U(z) = u00 + u01z + ... + u06z^6	49	* U(z) = u00 + u01z + ... + u06z^6
51	* V(z) = 1 + v01z + ... + v04z^4	50	* V(z) = 1 + v01z + ... + v04z^4
52	* with absolute approximation error bounded by 2**-72.	51	* with absolute approximation error bounded by 2**-72.
Lines 71-77 Link Here
71	one = 1.0,	70	one = 1.0,
72	invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */	71	invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
73	tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */	72	tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
74	/* R0/S0 on [0, 2.00] */	73	/* R0/S0 on [0, 2.00] */
75	R02 = 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */	74	R02 = 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */
76	R03 = -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */	75	R03 = -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */
77	R04 = 1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */	76	R04 = 1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */
Lines 157-163 Link Here
157	* y0(Inf) = 0.	156	* y0(Inf) = 0.
158	* y0(-Inf) = NaN and raise invalid exception.	157	* y0(-Inf) = NaN and raise invalid exception.
159	*/	158	*/
160	if(ix>=0x7ff00000) return vone/(x+x*x);	159	if(ix>=0x7ff00000) return vone/(x+x*x);
161	/* y0(+-0) = -inf and raise divide-by-zero exception. */	160	/* y0(+-0) = -inf and raise divide-by-zero exception. */
162	if((ix\|lx)==0) return -one/vzero;	161	if((ix\|lx)==0) return -one/vzero;
163	/* y0(x<0) = NaN and raise invalid exception. */	162	/* y0(x<0) = NaN and raise invalid exception. */
Lines 293-300 Link Here
293	s = one+z(q[0]+z(q[1]+z(q[2]+z(q[3]+z*q[4]))));	292	s = one+z(q[0]+z(q[1]+z(q[2]+z(q[3]+z*q[4]))));
294	return one+ r/s;	293	return one+ r/s;
295	}	294	}
296
297		295
		296
298	/* For x >= 8, the asymptotic expansions of qzero is	297	/* For x >= 8, the asymptotic expansions of qzero is
299	* -1/8 s + 75/1024 s^3 - ..., where s = 1/x.	298	* -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
300	* We approximate pzero by	299	* We approximate pzero by

Lines 1-4 Link Here

(-)src/e_j1.c (-10 / +9 lines)
1
2	/* @(#)e_j1.c 1.3 95/01/18 */	1	/* @(#)e_j1.c 1.3 95/01/18 */
3	/*	2	/*
4	* ====================================================	3	* ====================================================
Lines 6-12 Link Here
6	*	5	*
7	* Developed at SunSoft, a Sun Microsystems, Inc. business.	6	* Developed at SunSoft, a Sun Microsystems, Inc. business.
8	* Permission to use, copy, modify, and distribute this	7	* Permission to use, copy, modify, and distribute this
9	* software is freely granted, provided that this notice	8	* software is freely granted, provided that this notice
10	* is preserved.	9	* is preserved.
11	* ====================================================	10	* ====================================================
12	*/	11	*/
Lines 34-49 Link Here
34	* (To avoid cancellation, use	33	* (To avoid cancellation, use
35	* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))	34	* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
36	* to compute the worse one.)	35	* to compute the worse one.)
37	*	36	*
38	* 3 Special cases	37	* 3 Special cases
39	* j1(nan)= nan	38	* j1(nan)= nan
40	* j1(0) = 0	39	* j1(0) = 0
41	* j1(inf) = 0	40	* j1(inf) = 0
42	*	41	*
43	* Method -- y1(x):	42	* Method -- y1(x):
44	* 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN	43	* 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
45	* 2. For x<2.	44	* 2. For x<2.
46	* Since	45	* Since
47	* y1(x) = 2/pi(j1(x)(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)	46	* y1(x) = 2/pi(j1(x)(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
48	* therefore y1(x)-2/pij1(x)ln(x)-1/x is an odd function.	47	* therefore y1(x)-2/pij1(x)ln(x)-1/x is an odd function.
49	* We use the following function to approximate y1,	48	* We use the following function to approximate y1,
Lines 154-160 Link Here
154	* y1(Inf) = 0.	153	* y1(Inf) = 0.
155	* y1(-Inf) = NaN and raise invalid exception.	154	* y1(-Inf) = NaN and raise invalid exception.
156	*/	155	*/
157	if(ix>=0x7ff00000) return vone/(x+x*x);	156	if(ix>=0x7ff00000) return vone/(x+x*x);
158	/* y1(+-0) = -inf and raise divide-by-zero exception. */	157	/* y1(+-0) = -inf and raise divide-by-zero exception. */
159	if((ix\|lx)==0) return -one/vzero;	158	if((ix\|lx)==0) return -one/vzero;
160	/* y1(x<0) = NaN and raise invalid exception. */	159	/* y1(x<0) = NaN and raise invalid exception. */
Lines 186-195 Link Here
186	z = invsqrtpi(uss+v*cc)/sqrt(x);	185	z = invsqrtpi(uss+v*cc)/sqrt(x);
187	}	186	}
188	return z;	187	return z;
189	}	188	}
190	if(ix<=0x3c900000) { /* x < 2*-54 /	189	if(ix<=0x3c900000) { /* x < 2*-54 /
191	return(-tpi/x);	190	return(-tpi/x);
192	}	191	}
193	z = x*x;	192	z = x*x;
194	u = U0[0]+z(U0[1]+z(U0[2]+z(U0[3]+zU0[4])));	193	u = U0[0]+z(U0[1]+z(U0[2]+z(U0[3]+zU0[4])));
195	v = one+z(V0[0]+z(V0[1]+z(V0[2]+z(V0[3]+z*V0[4]))));	194	v = one+z(V0[0]+z(V0[1]+z(V0[2]+z(V0[3]+z*V0[4]))));
Lines 287-294 Link Here
287	s = one+z(q[0]+z(q[1]+z(q[2]+z(q[3]+z*q[4]))));	286	s = one+z(q[0]+z(q[1]+z(q[2]+z(q[3]+z*q[4]))));
288	return one+ r/s;	287	return one+ r/s;
289	}	288	}
290
291		289
		290
292	/* For x >= 8, the asymptotic expansions of qone is	291	/* For x >= 8, the asymptotic expansions of qone is
293	* 3/8 s - 105/1024 s^3 - ..., where s = 1/x.	292	* 3/8 s - 105/1024 s^3 - ..., where s = 1/x.
294	* We approximate pone by	293	* We approximate pone by

Lines 32-38 Link Here

(-)src/e_j1f.c (-1 / +1 lines)
32	one = 1.0,	32	one = 1.0,
33	invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */	33	invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */
34	tpi = 6.3661974669e-01, /* 0x3f22f983 */	34	tpi = 6.3661974669e-01, /* 0x3f22f983 */
35	/* R0/S0 on [0,2] */	35	/* R0/S0 on [0,2] */
36	r00 = -6.2500000000e-02, /* 0xbd800000 */	36	r00 = -6.2500000000e-02, /* 0xbd800000 */
37	r01 = 1.4070566976e-03, /* 0x3ab86cfd */	37	r01 = 1.4070566976e-03, /* 0x3ab86cfd */
38	r02 = -1.5995563444e-05, /* 0xb7862e36 */	38	r02 = -1.5995563444e-05, /* 0xb7862e36 */

Lines 1-4 Link Here

(-)src/e_jn.c (-20 / +18 lines)
1
2	/* @(#)e_jn.c 1.4 95/01/18 */	1	/* @(#)e_jn.c 1.4 95/01/18 */
3	/*	2	/*
4	* ====================================================	3	* ====================================================
Lines 6-12 Link Here
6	*	5	*
7	* Developed at SunSoft, a Sun Microsystems, Inc. business.	6	* Developed at SunSoft, a Sun Microsystems, Inc. business.
8	* Permission to use, copy, modify, and distribute this	7	* Permission to use, copy, modify, and distribute this
9	* software is freely granted, provided that this notice	8	* software is freely granted, provided that this notice
10	* is preserved.	9	* is preserved.
11	* ====================================================	10	* ====================================================
12	*/	11	*/
Lines 18-24 Link Here
18	* __ieee754_jn(n, x), __ieee754_yn(n, x)	17	* __ieee754_jn(n, x), __ieee754_yn(n, x)
19	* floating point Bessel's function of the 1st and 2nd kind	18	* floating point Bessel's function of the 1st and 2nd kind
20	* of order n	19	* of order n
21	*	20	*
22	* Special cases:	21	* Special cases:
23	* y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;	22	* y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
24	* y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.	23	* y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
Lines 37-43 Link Here
37	* yn(n,x) is similar in all respects, except	36	* yn(n,x) is similar in all respects, except
38	* that forward recursion is used for all	37	* that forward recursion is used for all
39	* values of n>1.	38	* values of n>1.
40	*
41	*/	39	*/
42		40
43	#include "math.h"	41	#include "math.h"
Lines 66-72 Link Here
66	ix = 0x7fffffff&hx;	64	ix = 0x7fffffff&hx;
67	/* if J(n,NaN) is NaN */	65	/* if J(n,NaN) is NaN */
68	if((ix\|((u_int32_t)(lx\|-lx))>>31)>0x7ff00000) return x+x;	66	if((ix\|((u_int32_t)(lx\|-lx))>>31)>0x7ff00000) return x+x;
69	if(n<0){	67	if(n<0){
70	n = -n;	68	n = -n;
71	x = -x;	69	x = -x;
72	hx ^= 0x80000000;	70	hx ^= 0x80000000;
Lines 77-90 Link Here
77	x = fabs(x);	75	x = fabs(x);
78	if((ix\|lx)==0\|\|ix>=0x7ff00000) /* if x is 0 or inf */	76	if((ix\|lx)==0\|\|ix>=0x7ff00000) /* if x is 0 or inf */
79	b = zero;	77	b = zero;
80	else if((double)n<=x) {	78	else if((double)n<=x) {
81	/* Safe to use J(n+1,x)=2n/x J(n,x)-J(n-1,x) /	79	/* Safe to use J(n+1,x)=2n/x J(n,x)-J(n-1,x) /
82	if(ix>=0x52D00000) { /* x > 2*302 /	80	if(ix>=0x52D00000) { /* x > 2*302 /
83	/* (x >> n**2)	81	/* (x >> n**2)
84	* Jn(x) = cos(x-(2n+1)pi/4)sqrt(2/x*pi)	82	* Jn(x) = cos(x-(2n+1)pi/4)sqrt(2/x*pi)
85	* Yn(x) = sin(x-(2n+1)pi/4)sqrt(2/x*pi)	83	* Yn(x) = sin(x-(2n+1)pi/4)sqrt(2/x*pi)
86	* Let s=sin(x), c=cos(x),	84	* Let s=sin(x), c=cos(x),
87	* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then	85	* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2), then
88	*	86	*
89	* n sin(xn)sqt2 cos(xn)sqt2	87	* n sin(xn)sqt2 cos(xn)sqt2
90	* ----------------------------------	88	* ----------------------------------
Lines 100-106 Link Here
100	case 3: temp = cos(x)-sin(x); break;	98	case 3: temp = cos(x)-sin(x); break;
101	}	99	}
102	b = invsqrtpi*temp/sqrt(x);	100	b = invsqrtpi*temp/sqrt(x);
103	} else {	101	} else {
104	a = __ieee754_j0(x);	102	a = __ieee754_j0(x);
105	b = __ieee754_j1(x);	103	b = __ieee754_j1(x);
106	for(i=1;i<n;i++){	104	for(i=1;i<n;i++){
Lines 111-117 Link Here
111	}	109	}
112	} else {	110	} else {
113	if(ix<0x3e100000) { /* x < 2*-29 /	111	if(ix<0x3e100000) { /* x < 2*-29 /
114	/* x is tiny, return the first Taylor expansion of J(n,x)	112	/* x is tiny, return the first Taylor expansion of J(n,x)
115	* J(n,x) = 1/n!*(x/2)^n - ...	113	* J(n,x) = 1/n!*(x/2)^n - ...
116	*/	114	*/
117	if(n>33) /* underflow */	115	if(n>33) /* underflow */
Lines 126-139 Link Here
126	}	124	}
127	} else {	125	} else {
128	/* use backward recurrence */	126	/* use backward recurrence */
129	/* x x^2 x^2	127	/* x x^2 x^2
130	* J(n,x)/J(n-1,x) = ---- ------ ------ .....	128	* J(n,x)/J(n-1,x) = ---- ------ ------ .....
131	* 2n - 2(n+1) - 2(n+2)	129	* 2n - 2(n+1) - 2(n+2)
132	*	130	*
133	* 1 1 1	131	* 1 1 1
134	* (for large x) = ---- ------ ------ .....	132	* (for large x) = ---- ------ ------ .....
135	* 2n 2(n+1) 2(n+2)	133	* 2n 2(n+1) 2(n+2)
136	* -- - ------ - ------ -	134	* -- - ------ - ------ -
137	* x x x	135	* x x x
138	*	136	*
139	* Let w = 2n/x and h=2/x, then the above quotient	137	* Let w = 2n/x and h=2/x, then the above quotient
Lines 149-157 Link Here
149	* To determine how many terms needed, let	147	* To determine how many terms needed, let
150	* Q(0) = w, Q(1) = w(w+h) - 1,	148	* Q(0) = w, Q(1) = w(w+h) - 1,
151	* Q(k) = (w+kh)Q(k-1) - Q(k-2),	149	* Q(k) = (w+kh)Q(k-1) - Q(k-2),
152	* When Q(k) > 1e4 good for single	150	* When Q(k) > 1e4 good for single
153	* When Q(k) > 1e9 good for double	151	* When Q(k) > 1e9 good for double
154	* When Q(k) > 1e17 good for quadruple	152	* When Q(k) > 1e17 good for quadruple
155	*/	153	*/
156	/* determine k */	154	/* determine k */
157	double t,v;	155	double t,v;
Lines 237-247 Link Here
237	if(n==1) return(sign*__ieee754_y1(x));	235	if(n==1) return(sign*__ieee754_y1(x));
238	if(ix==0x7ff00000) return zero;	236	if(ix==0x7ff00000) return zero;
239	if(ix>=0x52D00000) { /* x > 2*302 /	237	if(ix>=0x52D00000) { /* x > 2*302 /
240	/* (x >> n**2)	238	/* (x >> n**2)
241	* Jn(x) = cos(x-(2n+1)pi/4)sqrt(2/x*pi)	239	* Jn(x) = cos(x-(2n+1)pi/4)sqrt(2/x*pi)
242	* Yn(x) = sin(x-(2n+1)pi/4)sqrt(2/x*pi)	240	* Yn(x) = sin(x-(2n+1)pi/4)sqrt(2/x*pi)
243	* Let s=sin(x), c=cos(x),	241	* Let s=sin(x), c=cos(x),
244	* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then	242	* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2), then
245	*	243	*
246	* n sin(xn)sqt2 cos(xn)sqt2	244	* n sin(xn)sqt2 cos(xn)sqt2
247	* ----------------------------------	245	* ----------------------------------

Return to bug 229423