← Back to branch summary

~mc.../inkscape/inkscape

« back to all changes in this revision

Viewing changes to src/extension/script/js/fdlibm/e_jn.c

  • Committer: mental
  • Date: 2006-01-16 02:36:01 UTC
  • Revision ID: mental@users.sourceforge.net-20060116023601-wkr0h7edl5veyudq
moving trunk for module inkscape

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/extension/script/js/fdlibm/e_jn.c
 
1
/* -*- 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
}