/[jscoverage]/trunk/js/src/fdlibm/e_jn.c
ViewVC logotype

Annotation of /trunk/js/src/fdlibm/e_jn.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2 - (hide annotations)
Wed Aug 1 13:51:53 2007 UTC (14 years, 11 months ago) by siliconforks
File MIME type: text/plain
File size: 9092 byte(s)
Initial import.

1 siliconforks 2 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2     *
3     * ***** BEGIN LICENSE BLOCK *****
4     * Version: MPL 1.1/GPL 2.0/LGPL 2.1
5     *
6     * The contents of this file are subject to the Mozilla Public License Version
7     * 1.1 (the "License"); you may not use this file except in compliance with
8     * the License. You may obtain a copy of the License at
9     * http://www.mozilla.org/MPL/
10     *
11     * Software distributed under the License is distributed on an "AS IS" basis,
12     * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
13     * for the specific language governing rights and limitations under the
14     * License.
15     *
16     * The Original Code is Mozilla Communicator client code, released
17     * March 31, 1998.
18     *
19     * The Initial Developer of the Original Code is
20     * Sun Microsystems, Inc.
21     * Portions created by the Initial Developer are Copyright (C) 1998
22     * the Initial Developer. All Rights Reserved.
23     *
24     * Contributor(s):
25     *
26     * Alternatively, the contents of this file may be used under the terms of
27     * either of the GNU General Public License Version 2 or later (the "GPL"),
28     * or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
29     * in which case the provisions of the GPL or the LGPL are applicable instead
30     * of those above. If you wish to allow use of your version of this file only
31     * under the terms of either the GPL or the LGPL, and not to allow others to
32     * use your version of this file under the terms of the MPL, indicate your
33     * decision by deleting the provisions above and replace them with the notice
34     * and other provisions required by the GPL or the LGPL. If you do not delete
35     * the provisions above, a recipient may use your version of this file under
36     * the terms of any one of the MPL, the GPL or the LGPL.
37     *
38     * ***** END LICENSE BLOCK ***** */
39    
40     /* @(#)e_jn.c 1.4 95/01/18 */
41     /*
42     * ====================================================
43     * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
44     *
45     * Developed at SunSoft, a Sun Microsystems, Inc. business.
46     * Permission to use, copy, modify, and distribute this
47     * software is freely granted, provided that this notice
48     * is preserved.
49     * ====================================================
50     */
51    
52     /*
53     * __ieee754_jn(n, x), __ieee754_yn(n, x)
54     * floating point Bessel's function of the 1st and 2nd kind
55     * of order n
56     *
57     * Special cases:
58     * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
59     * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
60     * Note 2. About jn(n,x), yn(n,x)
61     * For n=0, j0(x) is called,
62     * for n=1, j1(x) is called,
63     * for n<x, forward recursion us used starting
64     * from values of j0(x) and j1(x).
65     * for n>x, a continued fraction approximation to
66     * j(n,x)/j(n-1,x) is evaluated and then backward
67     * recursion is used starting from a supposed value
68     * for j(n,x). The resulting value of j(0,x) is
69     * compared with the actual value to correct the
70     * supposed value of j(n,x).
71     *
72     * yn(n,x) is similar in all respects, except
73     * that forward recursion is used for all
74     * values of n>1.
75     *
76     */
77    
78     #include "fdlibm.h"
79    
80     #ifdef __STDC__
81     static const double
82     #else
83     static double
84     #endif
85     invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
86     two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
87     one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
88    
89     static double zero = 0.00000000000000000000e+00;
90    
91     #ifdef __STDC__
92     double __ieee754_jn(int n, double x)
93     #else
94     double __ieee754_jn(n,x)
95     int n; double x;
96     #endif
97     {
98     fd_twoints u;
99     int i,hx,ix,lx, sgn;
100     double a, b, temp, di;
101     double z, w;
102    
103     /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
104     * Thus, J(-n,x) = J(n,-x)
105     */
106     u.d = x;
107     hx = __HI(u);
108     ix = 0x7fffffff&hx;
109     lx = __LO(u);
110     /* if J(n,NaN) is NaN */
111     if((ix|((unsigned)(lx|-lx))>>31)>0x7ff00000) return x+x;
112     if(n<0){
113     n = -n;
114     x = -x;
115     hx ^= 0x80000000;
116     }
117     if(n==0) return(__ieee754_j0(x));
118     if(n==1) return(__ieee754_j1(x));
119     sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
120     x = fd_fabs(x);
121     if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */
122     b = zero;
123     else if((double)n<=x) {
124     /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
125     if(ix>=0x52D00000) { /* x > 2**302 */
126     /* (x >> n**2)
127     * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
128     * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
129     * Let s=sin(x), c=cos(x),
130     * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
131     *
132     * n sin(xn)*sqt2 cos(xn)*sqt2
133     * ----------------------------------
134     * 0 s-c c+s
135     * 1 -s-c -c+s
136     * 2 -s+c -c-s
137     * 3 s+c c-s
138     */
139     switch(n&3) {
140     case 0: temp = fd_cos(x)+fd_sin(x); break;
141     case 1: temp = -fd_cos(x)+fd_sin(x); break;
142     case 2: temp = -fd_cos(x)-fd_sin(x); break;
143     case 3: temp = fd_cos(x)-fd_sin(x); break;
144     }
145     b = invsqrtpi*temp/fd_sqrt(x);
146     } else {
147     a = __ieee754_j0(x);
148     b = __ieee754_j1(x);
149     for(i=1;i<n;i++){
150     temp = b;
151     b = b*((double)(i+i)/x) - a; /* avoid underflow */
152     a = temp;
153     }
154     }
155     } else {
156     if(ix<0x3e100000) { /* x < 2**-29 */
157     /* x is tiny, return the first Taylor expansion of J(n,x)
158     * J(n,x) = 1/n!*(x/2)^n - ...
159     */
160     if(n>33) /* underflow */
161     b = zero;
162     else {
163     temp = x*0.5; b = temp;
164     for (a=one,i=2;i<=n;i++) {
165     a *= (double)i; /* a = n! */
166     b *= temp; /* b = (x/2)^n */
167     }
168     b = b/a;
169     }
170     } else {
171     /* use backward recurrence */
172     /* x x^2 x^2
173     * J(n,x)/J(n-1,x) = ---- ------ ------ .....
174     * 2n - 2(n+1) - 2(n+2)
175     *
176     * 1 1 1
177     * (for large x) = ---- ------ ------ .....
178     * 2n 2(n+1) 2(n+2)
179     * -- - ------ - ------ -
180     * x x x
181     *
182     * Let w = 2n/x and h=2/x, then the above quotient
183     * is equal to the continued fraction:
184     * 1
185     * = -----------------------
186     * 1
187     * w - -----------------
188     * 1
189     * w+h - ---------
190     * w+2h - ...
191     *
192     * To determine how many terms needed, let
193     * Q(0) = w, Q(1) = w(w+h) - 1,
194     * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
195     * When Q(k) > 1e4 good for single
196     * When Q(k) > 1e9 good for double
197     * When Q(k) > 1e17 good for quadruple
198     */
199     /* determine k */
200     double t,v;
201     double q0,q1,h,tmp; int k,m;
202     w = (n+n)/(double)x; h = 2.0/(double)x;
203     q0 = w; z = w+h; q1 = w*z - 1.0; k=1;
204     while(q1<1.0e9) {
205     k += 1; z += h;
206     tmp = z*q1 - q0;
207     q0 = q1;
208     q1 = tmp;
209     }
210     m = n+n;
211     for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
212     a = t;
213     b = one;
214     /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
215     * Hence, if n*(log(2n/x)) > ...
216     * single 8.8722839355e+01
217     * double 7.09782712893383973096e+02
218     * long double 1.1356523406294143949491931077970765006170e+04
219     * then recurrent value may overflow and the result is
220     * likely underflow to zero
221     */
222     tmp = n;
223     v = two/x;
224     tmp = tmp*__ieee754_log(fd_fabs(v*tmp));
225     if(tmp<7.09782712893383973096e+02) {
226     for(i=n-1,di=(double)(i+i);i>0;i--){
227     temp = b;
228     b *= di;
229     b = b/x - a;
230     a = temp;
231     di -= two;
232     }
233     } else {
234     for(i=n-1,di=(double)(i+i);i>0;i--){
235     temp = b;
236     b *= di;
237     b = b/x - a;
238     a = temp;
239     di -= two;
240     /* scale b to avoid spurious overflow */
241     if(b>1e100) {
242     a /= b;
243     t /= b;
244     b = one;
245     }
246     }
247     }
248     b = (t*__ieee754_j0(x)/b);
249     }
250     }
251     if(sgn==1) return -b; else return b;
252     }
253    
254     #ifdef __STDC__
255     double __ieee754_yn(int n, double x)
256     #else
257     double __ieee754_yn(n,x)
258     int n; double x;
259     #endif
260     {
261     fd_twoints u;
262     int i,hx,ix,lx;
263     int sign;
264     double a, b, temp;
265    
266     u.d = x;
267     hx = __HI(u);
268     ix = 0x7fffffff&hx;
269     lx = __LO(u);
270     /* if Y(n,NaN) is NaN */
271     if((ix|((unsigned)(lx|-lx))>>31)>0x7ff00000) return x+x;
272     if((ix|lx)==0) return -one/zero;
273     if(hx<0) return zero/zero;
274     sign = 1;
275     if(n<0){
276     n = -n;
277     sign = 1 - ((n&1)<<1);
278     }
279     if(n==0) return(__ieee754_y0(x));
280     if(n==1) return(sign*__ieee754_y1(x));
281     if(ix==0x7ff00000) return zero;
282     if(ix>=0x52D00000) { /* x > 2**302 */
283     /* (x >> n**2)
284     * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
285     * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
286     * Let s=sin(x), c=cos(x),
287     * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
288     *
289     * n sin(xn)*sqt2 cos(xn)*sqt2
290     * ----------------------------------
291     * 0 s-c c+s
292     * 1 -s-c -c+s
293     * 2 -s+c -c-s
294     * 3 s+c c-s
295     */
296     switch(n&3) {
297     case 0: temp = fd_sin(x)-fd_cos(x); break;
298     case 1: temp = -fd_sin(x)-fd_cos(x); break;
299     case 2: temp = -fd_sin(x)+fd_cos(x); break;
300     case 3: temp = fd_sin(x)+fd_cos(x); break;
301     }
302     b = invsqrtpi*temp/fd_sqrt(x);
303     } else {
304     a = __ieee754_y0(x);
305     b = __ieee754_y1(x);
306     /* quit if b is -inf */
307     u.d = b;
308     for(i=1;i<n&&(__HI(u) != 0xfff00000);i++){
309     temp = b;
310     b = ((double)(i+i)/x)*b - a;
311     a = temp;
312     }
313     }
314     if(sign>0) return b; else return -b;
315     }

  ViewVC Help
Powered by ViewVC 1.1.24