Mercurial > hg > CbC > CbC_gcc

comparison libquadmath/math/jnq.c @ 68:561a7518be6b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update gcc-4.6
author Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
date Sun, 21 Aug 2011 07:07:55 +0900
parents
children 04ced10e8804
comparison
equal deleted inserted replaced
67:f6334be47118 68:561a7518be6b
1 /*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 */
11
12 /* Modifications for 128-bit long double are
13 Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
14 and are incorporated herein by permission of the author. The author
15 reserves the right to distribute this material elsewhere under different
16 copying permissions. These modifications are distributed here under
17 the following terms:
18
19 This library is free software; you can redistribute it and/or
20 modify it under the terms of the GNU Lesser General Public
21 License as published by the Free Software Foundation; either
22 version 2.1 of the License, or (at your option) any later version.
23
24 This library is distributed in the hope that it will be useful,
25 but WITHOUT ANY WARRANTY; without even the implied warranty of
26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 Lesser General Public License for more details.
28
29 You should have received a copy of the GNU Lesser General Public
30 License along with this library; if not, write to the Free Software
31 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
32
33 /*
34 * __ieee754_jn(n, x), __ieee754_yn(n, x)
35 * floating point Bessel's function of the 1st and 2nd kind
36 * of order n
37 *
38 * Special cases:
39 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
40 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
41 * Note 2. About jn(n,x), yn(n,x)
42 * For n=0, j0(x) is called,
43 * for n=1, j1(x) is called,
44 * for n<x, forward recursion us used starting
45 * from values of j0(x) and j1(x).
46 * for n>x, a continued fraction approximation to
47 * j(n,x)/j(n-1,x) is evaluated and then backward
48 * recursion is used starting from a supposed value
49 * for j(n,x). The resulting value of j(0,x) is
50 * compared with the actual value to correct the
51 * supposed value of j(n,x).
52 *
53 * yn(n,x) is similar in all respects, except
54 * that forward recursion is used for all
55 * values of n>1.
56 *
57 */
58
59 #include "quadmath-imp.h"
60
61 static const __float128
62 invsqrtpi = 5.6418958354775628694807945156077258584405E-1Q,
63 two = 2.0e0Q,
64 one = 1.0e0Q,
65 zero = 0.0Q;
66
67
68 __float128
69 jnq (int n, __float128 x)
70 {
71 uint32_t se;
72 int32_t i, ix, sgn;
73 __float128 a, b, temp, di;
74 __float128 z, w;
75 ieee854_float128 u;
76
77
78 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
79 * Thus, J(-n,x) = J(n,-x)
80 */
81
82 u.value = x;
83 se = u.words32.w0;
84 ix = se & 0x7fffffff;
85
86 /* if J(n,NaN) is NaN */
87 if (ix >= 0x7fff0000)
88 {
89 if ((u.words32.w0 & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3)
90 return x + x;
91 }
92
93 if (n < 0)
94 {
95 n = -n;
96 x = -x;
97 se ^= 0x80000000;
98 }
99 if (n == 0)
100 return (j0q (x));
101 if (n == 1)
102 return (j1q (x));
103 sgn = (n & 1) & (se >> 31); /* even n -- 0, odd n -- sign(x) */
104 x = fabsq (x);
105
106 if (x == 0.0Q || ix >= 0x7fff0000) /* if x is 0 or inf */
107 b = zero;
108 else if ((__float128) n <= x)
109 {
110 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
111 if (ix >= 0x412D0000)
112 { /* x > 2**302 */
113
114 /* ??? Could use an expansion for large x here. */
115
116 /* (x >> n**2)
117 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
118 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
119 * Let s=sin(x), c=cos(x),
120 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
121 *
122 * n sin(xn)*sqt2 cos(xn)*sqt2
123 * ----------------------------------
124 * 0 s-c c+s
125 * 1 -s-c -c+s
126 * 2 -s+c -c-s
127 * 3 s+c c-s
128 */
129 __float128 s;
130 __float128 c;
131 sincosq (x, &s, &c);
132 switch (n & 3)
133 {
134 case 0:
135 temp = c + s;
136 break;
137 case 1:
138 temp = -c + s;
139 break;
140 case 2:
141 temp = -c - s;
142 break;
143 case 3:
144 temp = c - s;
145 break;
146 }
147 b = invsqrtpi * temp / sqrtq (x);
148 }
149 else
150 {
151 a = j0q (x);
152 b = j1q (x);
153 for (i = 1; i < n; i++)
154 {
155 temp = b;
156 b = b * ((__float128) (i + i) / x) - a; /* avoid underflow */
157 a = temp;
158 }
159 }
160 }
161 else
162 {
163 if (ix < 0x3fc60000)
164 { /* x < 2**-57 */
165 /* x is tiny, return the first Taylor expansion of J(n,x)
166 * J(n,x) = 1/n!*(x/2)^n - ...
167 */
168 if (n >= 400) /* underflow, result < 10^-4952 */
169 b = zero;
170 else
171 {
172 temp = x * 0.5;
173 b = temp;
174 for (a = one, i = 2; i <= n; i++)
175 {
176 a *= (__float128) i; /* a = n! */
177 b *= temp; /* b = (x/2)^n */
178 }
179 b = b / a;
180 }
181 }
182 else
183 {
184 /* use backward recurrence */
185 /* x x^2 x^2
186 * J(n,x)/J(n-1,x) = ---- ------ ------ .....
187 * 2n - 2(n+1) - 2(n+2)
188 *
189 * 1 1 1
190 * (for large x) = ---- ------ ------ .....
191 * 2n 2(n+1) 2(n+2)
192 * -- - ------ - ------ -
193 * x x x
194 *
195 * Let w = 2n/x and h=2/x, then the above quotient
196 * is equal to the continued fraction:
197 * 1
198 * = -----------------------
199 * 1
200 * w - -----------------
201 * 1
202 * w+h - ---------
203 * w+2h - ...
204 *
205 * To determine how many terms needed, let
206 * Q(0) = w, Q(1) = w(w+h) - 1,
207 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
208 * When Q(k) > 1e4 good for single
209 * When Q(k) > 1e9 good for double
210 * When Q(k) > 1e17 good for quadruple
211 */
212 /* determine k */
213 __float128 t, v;
214 __float128 q0, q1, h, tmp;
215 int32_t k, m;
216 w = (n + n) / (__float128) x;
217 h = 2.0Q / (__float128) x;
218 q0 = w;
219 z = w + h;
220 q1 = w * z - 1.0Q;
221 k = 1;
222 while (q1 < 1.0e17Q)
223 {
224 k += 1;
225 z += h;
226 tmp = z * q1 - q0;
227 q0 = q1;
228 q1 = tmp;
229 }
230 m = n + n;
231 for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
232 t = one / (i / x - t);
233 a = t;
234 b = one;
235 /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
236 * Hence, if n*(log(2n/x)) > ...
237 * single 8.8722839355e+01
238 * double 7.09782712893383973096e+02
239 * __float128 1.1356523406294143949491931077970765006170e+04
240 * then recurrent value may overflow and the result is
241 * likely underflow to zero
242 */
243 tmp = n;
244 v = two / x;
245 tmp = tmp * logq (fabsq (v * tmp));
246
247 if (tmp < 1.1356523406294143949491931077970765006170e+04Q)
248 {
249 for (i = n - 1, di = (__float128) (i + i); i > 0; i--)
250 {
251 temp = b;
252 b *= di;
253 b = b / x - a;
254 a = temp;
255 di -= two;
256 }
257 }
258 else
259 {
260 for (i = n - 1, di = (__float128) (i + i); i > 0; i--)
261 {
262 temp = b;
263 b *= di;
264 b = b / x - a;
265 a = temp;
266 di -= two;
267 /* scale b to avoid spurious overflow */
268 if (b > 1e100Q)
269 {
270 a /= b;
271 t /= b;
272 b = one;
273 }
274 }
275 }
276 b = (t * j0q (x) / b);
277 }
278 }
279 if (sgn == 1)
280 return -b;
281 else
282 return b;
283 }
284
285 __float128
286 ynq (int n, __float128 x)
287 {
288 uint32_t se;
289 int32_t i, ix;
290 int32_t sign;
291 __float128 a, b, temp;
292 ieee854_float128 u;
293
294 u.value = x;
295 se = u.words32.w0;
296 ix = se & 0x7fffffff;
297
298 /* if Y(n,NaN) is NaN */
299 if (ix >= 0x7fff0000)
300 {
301 if ((u.words32.w0 & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3)
302 return x + x;
303 }
304 if (x <= 0.0Q)
305 {
306 if (x == 0.0Q)
307 return -HUGE_VALQ + x;
308 if (se & 0x80000000)
309 return zero / (zero * x);
310 }
311 sign = 1;
312 if (n < 0)
313 {
314 n = -n;
315 sign = 1 - ((n & 1) << 1);
316 }
317 if (n == 0)
318 return (y0q (x));
319 if (n == 1)
320 return (sign * y1q (x));
321 if (ix >= 0x7fff0000)
322 return zero;
323 if (ix >= 0x412D0000)
324 { /* x > 2**302 */
325
326 /* ??? See comment above on the possible futility of this. */
327
328 /* (x >> n**2)
329 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
330 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
331 * Let s=sin(x), c=cos(x),
332 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
333 *
334 * n sin(xn)*sqt2 cos(xn)*sqt2
335 * ----------------------------------
336 * 0 s-c c+s
337 * 1 -s-c -c+s
338 * 2 -s+c -c-s
339 * 3 s+c c-s
340 */
341 __float128 s;
342 __float128 c;
343 sincosq (x, &s, &c);
344 switch (n & 3)
345 {
346 case 0:
347 temp = s - c;
348 break;
349 case 1:
350 temp = -s - c;
351 break;
352 case 2:
353 temp = -s + c;
354 break;
355 case 3:
356 temp = s + c;
357 break;
358 }
359 b = invsqrtpi * temp / sqrtq (x);
360 }
361 else
362 {
363 a = y0q (x);
364 b = y1q (x);
365 /* quit if b is -inf */
366 u.value = b;
367 se = u.words32.w0 & 0xffff0000;
368 for (i = 1; i < n && se != 0xffff0000; i++)
369 {
370 temp = b;
371 b = ((__float128) (i + i) / x) * b - a;
372 u.value = b;
373 se = u.words32.w0 & 0xffff0000;
374 a = temp;
375 }
376 }
377 if (sign > 0)
378 return b;
379 else
380 return -b;
381 }