1
/* -*- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
3
* ***** BEGIN LICENSE BLOCK *****
4
* Version: MPL 1.1/GPL 2.0/LGPL 2.1
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/
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
16
* The Original Code is Mozilla Communicator client code, released
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.
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.
38
* ***** END LICENSE BLOCK ***** */
40
/* @(#)e_jn.c 1.4 95/01/18 */
42
* ====================================================
43
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
49
* ====================================================
53
* __ieee754_jn(n, x), __ieee754_yn(n, x)
54
* floating point Bessel's function of the 1st and 2nd kind
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).
72
* yn(n,x) is similar in all respects, except
73
* that forward recursion is used for all
85
invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
86
two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
87
one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
89
static double zero = 0.00000000000000000000e+00;
92
double __ieee754_jn(int n, double x)
94
double __ieee754_jn(n,x)
100
double a, b, temp, di;
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)
110
/* if J(n,NaN) is NaN */
111
if((ix|((unsigned)(lx|-lx))>>31)>0x7ff00000) return x+x;
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) */
121
if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */
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 */
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
132
* n sin(xn)*sqt2 cos(xn)*sqt2
133
* ----------------------------------
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;
145
b = invsqrtpi*temp/fd_sqrt(x);
151
b = b*((double)(i+i)/x) - a; /* avoid underflow */
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 - ...
160
if(n>33) /* underflow */
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 */
171
/* use backward recurrence */
173
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
174
* 2n - 2(n+1) - 2(n+2)
177
* (for large x) = ---- ------ ------ .....
179
* -- - ------ - ------ -
182
* Let w = 2n/x and h=2/x, then the above quotient
183
* is equal to the continued fraction:
185
* = -----------------------
187
* w - -----------------
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
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;
211
for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
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
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--){
234
for(i=n-1,di=(double)(i+i);i>0;i--){
240
/* scale b to avoid spurious overflow */
248
b = (t*__ieee754_j0(x)/b);
251
if(sgn==1) return -b; else return b;
255
double __ieee754_yn(int n, double x)
257
double __ieee754_yn(n,x)
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;
277
sign = 1 - ((n&1)<<1);
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 */
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
289
* n sin(xn)*sqt2 cos(xn)*sqt2
290
* ----------------------------------
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;
302
b = invsqrtpi*temp/fd_sqrt(x);
306
/* quit if b is -inf */
308
for(i=1;i<n&&(__HI(u) != 0xfff00000);i++){
310
b = ((double)(i+i)/x)*b - a;
314
if(sign>0) return b; else return -b;
added
removed
src/js/fdlibm/e_jn.c
