1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21 /*
22 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
23 */
24 /*
25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
26 * Use is subject to license terms.
27 */
28
29 #pragma weak __j0f = j0f
30 #pragma weak __j1f = j1f
31 #pragma weak __jnf = jnf
32 #pragma weak __y0f = y0f
33 #pragma weak __y1f = y1f
34 #pragma weak __ynf = ynf
35
36 #include "libm.h"
37 #include <float.h>
38
39 #if defined(__i386) && !defined(__amd64)
40 extern int __swapRP(int);
41 #endif
42
43 static const float
44 zerof = 0.0f,
45 onef = 1.0f;
46
47 static const double C[] = {
48 0.0,
49 -0.125,
50 0.25,
51 0.375,
52 0.5,
53 1.0,
54 2.0,
55 8.0,
56 0.5641895835477562869480794515607725858441, /* 1/sqrt(pi) */
57 0.636619772367581343075535053490057448, /* 2/pi */
58 1.0e9,
59 };
60
61 #define zero C[0]
62 #define neighth C[1]
63 #define quarter C[2]
64 #define three8 C[3]
65 #define half C[4]
85 0.5068271181053546392490184353e4,
86 0.1829817905472769960535671664e3,
87 1.0,
88 -0.1731210995701068539185611951e3, /* qr */
89 -0.5522559165936166961235240613e3,
90 -0.5604935606637346590614529613e3,
91 -0.2200430300226009379477365011e3,
92 -0.323869355375648849771296746e2,
93 -0.14294979207907956223499258e1,
94 -0.834690374102384988158918e-2,
95 0.1107975037248683865326709645e5, /* qs */
96 0.3544581680627082674651471873e5,
97 0.3619118937918394132179019059e5,
98 0.1439895563565398007471485822e5,
99 0.2190277023344363955930226234e4,
100 0.106695157020407986137501682e3,
101 1.0,
102 };
103
104 #define pr Cj0y0
105 #define ps (Cj0y0+7)
106 #define qr (Cj0y0+14)
107 #define qs (Cj0y0+21)
108
109 static const double Cj0[] = {
110 -2.500000000000003622131880894830476755537e-0001, /* r0 */
111 1.095597547334830263234433855932375353303e-0002,
112 -1.819734750463320921799187258987098087697e-0004,
113 9.977001946806131657544212501069893930846e-0007,
114 1.0, /* s0 */
115 1.867609810662950169966782360588199673741e-0002,
116 1.590389206181565490878430827706972074208e-0004,
117 6.520867386742583632375520147714499522721e-0007,
118 9.999999999999999942156495584397047660949e-0001, /* r1 */
119 -2.389887722731319130476839836908143731281e-0001,
120 1.293359476138939027791270393439493640570e-0002,
121 -2.770985642343140122168852400228563364082e-0004,
122 2.905241575772067678086738389169625218912e-0006,
123 -1.636846356264052597969042009265043251279e-0008,
124 5.072306160724884775085431059052611737827e-0011,
125 -8.187060730684066824228914775146536139112e-0014,
126 5.422219326959949863954297860723723423842e-0017,
127 1.0, /* s1 */
128 1.101122772686807702762104741932076228349e-0002,
129 6.140169310641649223411427764669143978228e-0005,
130 2.292035877515152097976946119293215705250e-0007,
131 6.356910426504644334558832036362219583789e-0010,
132 1.366626326900219555045096999553948891401e-0012,
133 2.280399586866739522891837985560481180088e-0015,
134 2.801559820648939665270492520004836611187e-0018,
135 2.073101088320349159764410261466350732968e-0021,
136 };
137
138 #define r0 Cj0
139 #define s0 (Cj0+4)
140 #define r1 (Cj0+8)
141 #define s1 (Cj0+17)
142
143 static const double Cy0[] = {
144 -7.380429510868722526754723020704317641941e-0002, /* u0 */
145 1.772607102684869924301459663049874294814e-0001,
146 -1.524370666542713828604078090970799356306e-0002,
147 4.650819100693891757143771557629924591915e-0004,
148 -7.125768872339528975036316108718239946022e-0006,
149 6.411017001656104598327565004771515257146e-0008,
150 -3.694275157433032553021246812379258781665e-0010,
151 1.434364544206266624252820889648445263842e-0012,
152 -3.852064731859936455895036286874139896861e-0015,
153 7.182052899726138381739945881914874579696e-0018,
154 -9.060556574619677567323741194079797987200e-0021,
155 7.124435467408860515265552217131230511455e-0024,
156 -2.709726774636397615328813121715432044771e-0027,
157 1.0, /* v0 */
158 4.678678931512549002587702477349214886475e-0003,
159 9.486828955529948534822800829497565178985e-0006,
160 1.001495929158861646659010844136682454906e-0008,
161 4.725338116256021660204443235685358593611e-0012,
162 };
163
164 #define u0 Cy0
165 #define v0 (Cy0+13)
166
167 static const double Cj1y1[] = {
168 -0.4435757816794127857114720794e7, /* pr0 */
169 -0.9942246505077641195658377899e7,
170 -0.6603373248364939109255245434e7,
171 -0.1523529351181137383255105722e7,
172 -0.1098240554345934672737413139e6,
173 -0.1611616644324610116477412898e4,
174 -0.4435757816794127856828016962e7, /* ps0 */
175 -0.9934124389934585658967556309e7,
176 -0.6585339479723087072826915069e7,
177 -0.1511809506634160881644546358e7,
178 -0.1072638599110382011903063867e6,
179 -0.1455009440190496182453565068e4,
180 0.3322091340985722351859704442e5, /* qr0 */
181 0.8514516067533570196555001171e5,
182 0.6617883658127083517939992166e5,
183 0.1849426287322386679652009819e5,
184 0.1706375429020768002061283546e4,
185 0.3526513384663603218592175580e2,
186 0.7087128194102874357377502472e6, /* qs0 */
187 0.1819458042243997298924553839e7,
188 0.1419460669603720892855755253e7,
189 0.4002944358226697511708610813e6,
190 0.3789022974577220264142952256e5,
191 0.8638367769604990967475517183e3,
192 };
193
194 #define pr0 Cj1y1
195 #define ps0 (Cj1y1+6)
196 #define qr0 (Cj1y1+12)
197 #define qs0 (Cj1y1+18)
198
199 static const double Cj1[] = {
200 -6.250000000000002203053200981413218949548e-0002, /* a0 */
201 1.600998455640072901321605101981501263762e-0003,
202 -1.963888815948313758552511884390162864930e-0005,
203 8.263917341093549759781339713418201620998e-0008,
204 1.0e0, /* b0 */
205 1.605069137643004242395356851797873766927e-0002,
206 1.149454623251299996428500249509098499383e-0004,
207 3.849701673735260970379681807910852327825e-0007,
208 4.999999999999999995517408894340485471724e-0001,
209 -6.003825028120475684835384519945468075423e-0002,
210 2.301719899263321828388344461995355419832e-0003,
211 -4.208494869238892934859525221654040304068e-0005,
212 4.377745135188837783031540029700282443388e-0007,
213 -2.854106755678624335145364226735677754179e-0009,
214 1.234002865443952024332943901323798413689e-0011,
215 -3.645498437039791058951273508838177134310e-0014,
216 7.404320596071797459925377103787837414422e-0017,
217 -1.009457448277522275262808398517024439084e-0019,
218 8.520158355824819796968771418801019930585e-0023,
219 -3.458159926081163274483854614601091361424e-0026,
220 1.0e0, /* b1 */
221 4.923499437590484879081138588998986303306e-0003,
222 1.054389489212184156499666953501976688452e-0005,
223 1.180768373106166527048240364872043816050e-0008,
224 5.942665743476099355323245707680648588540e-0012,
225 };
226
227 #define a0 Cj1
228 #define b0 (Cj1+4)
229 #define a1 (Cj1+8)
230 #define b1 (Cj1+20)
231
232 static const double Cy1[] = {
233 -1.960570906462389461018983259589655961560e-0001, /* c0 */
234 4.931824118350661953459180060007970291139e-0002,
235 -1.626975871565393656845930125424683008677e-0003,
236 1.359657517926394132692884168082224258360e-0005,
237 1.0e0, /* d0 */
238 2.565807214838390835108224713630901653793e-0002,
239 3.374175208978404268650522752520906231508e-0004,
240 2.840368571306070719539936935220728843177e-0006,
241 1.396387402048998277638900944415752207592e-0008,
242 -1.960570906462389473336339614647555351626e-0001, /* c1 */
243 5.336268030335074494231369159933012844735e-0002,
244 -2.684137504382748094149184541866332033280e-0003,
245 5.737671618979185736981543498580051903060e-0005,
246 -6.642696350686335339171171785557663224892e-0007,
247 4.692417922568160354012347591960362101664e-0009,
248 -2.161728635907789319335231338621412258355e-0011,
249 6.727353419738316107197644431844194668702e-0014,
250 -1.427502986803861372125234355906790573422e-0016,
251 2.020392498726806769468143219616642940371e-0019,
252 -1.761371948595104156753045457888272716340e-0022,
253 7.352828391941157905175042420249225115816e-0026,
254 1.0e0, /* d1 */
255 5.029187436727947764916247076102283399442e-0003,
256 1.102693095808242775074856548927801750627e-0005,
257 1.268035774543174837829534603830227216291e-0008,
258 6.579416271766610825192542295821308730206e-0012,
259 };
260
261 #define c0 Cy1
262 #define d0 (Cy1+4)
263 #define c1 (Cy1+9)
264 #define d1 (Cy1+21)
265
266
267 /* core of j0f computation; assumes fx is finite */
268 static double
269 __k_j0f(float fx)
270 {
271 double x, z, s, c, ss, cc, r, t, p0, q0;
272 int ix, i;
273
274 ix = *(int *)&fx & ~0x80000000;
275 x = fabs((double)fx);
276 if (ix > 0x41000000) {
277 /* x > 8; see comments in j0.c */
278 s = sin(x);
279 c = cos(x);
280 if (signbit(s) != signbit(c)) {
281 ss = s - c;
282 cc = -cos(x + x) / ss;
283 } else {
284 cc = s + c;
285 ss = -cos(x + x) / cc;
286 }
287 if (ix > 0x501502f9) {
288 /* x > 1.0e10 */
289 p0 = one;
290 q0 = neighth / x;
291 } else {
292 t = eight / x;
293 z = t * t;
294 p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] +
295 z * (pr[4] + z * (pr[5] + z * pr[6])))))) /
296 (ps[0] + z * (ps[1] + z * (ps[2] + z * (ps[3] +
297 z * (ps[4] + z * (ps[5] + z))))));
298 q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
299 z * (qr[4] + z * (qr[5] + z * qr[6])))))) /
300 (qs[0] + z * (qs[1] + z * (qs[2] + z * (qs[3] +
301 z * (qs[4] + z * (qs[5] + z))))))) * t;
302 }
303 return (isqrtpi * (p0 * cc - q0 * ss) / sqrt(x));
304 }
305 if (ix <= 0x3727c5ac) {
306 /* x <= 1.0e-5 */
307 if (ix <= 0x219392ef) /* x <= 1.0e-18 */
308 return (one - x);
309 return (one - x * x * quarter);
310 }
311 z = x * x;
312 if (ix <= 0x3fa3d70a) {
313 /* x <= 1.28 */
314 r = r0[0] + z * (r0[1] + z * (r0[2] + z * r0[3]));
315 s = s0[0] + z * (s0[1] + z * (s0[2] + z * s0[3]));
316 return (one + z * (r / s));
317 }
318 r = r1[8];
319 s = s1[8];
320 for (i = 7; i >= 0; i--) {
321 r = r * z + r1[i];
322 s = s * z + s1[i];
323 }
324 return (r / s);
325 }
326
327 float
328 j0f(float fx)
329 {
330 float f;
331 int ix;
332 #if defined(__i386) && !defined(__amd64)
333 int rp;
334 #endif
335
336 ix = *(int *)&fx & ~0x80000000;
337 if (ix >= 0x7f800000) { /* nan or inf */
338 if (ix > 0x7f800000)
339 return (fx * fx);
340 return (zerof);
341 }
342
343 #if defined(__i386) && !defined(__amd64)
344 rp = __swapRP(fp_extended);
345 #endif
346 f = (float)__k_j0f(fx);
347 #if defined(__i386) && !defined(__amd64)
348 if (rp != fp_extended)
349 (void) __swapRP(rp);
350 #endif
351 return (f);
352 }
353
354 /* core of y0f computation; assumes fx is finite and positive */
355 static double
356 __k_y0f(float fx)
357 {
358 double x, z, s, c, ss, cc, t, p0, q0, u, v;
359 int ix, i;
360
361 ix = *(int *)&fx;
362 x = (double)fx;
363 if (ix > 0x41000000) {
364 /* x > 8; see comments in j0.c */
365 s = sin(x);
366 c = cos(x);
367 if (signbit(s) != signbit(c)) {
368 ss = s - c;
369 cc = -cos(x + x) / ss;
370 } else {
371 cc = s + c;
372 ss = -cos(x + x) / cc;
373 }
374 if (ix > 0x501502f9) {
375 /* x > 1.0e10 */
376 p0 = one;
377 q0 = neighth / x;
378 } else {
379 t = eight / x;
380 z = t * t;
381 p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] +
382 z * (pr[4] + z * (pr[5] + z * pr[6])))))) /
383 (ps[0] + z * (ps[1] + z * (ps[2] + z * (ps[3] +
384 z * (ps[4] + z * (ps[5] + z))))));
385 q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
386 z * (qr[4] + z * (qr[5] + z * qr[6])))))) /
387 (qs[0] + z * (qs[1] + z * (qs[2] + z * (qs[3] +
388 z * (qs[4] + z * (qs[5] + z))))))) * t;
389 }
390 return (isqrtpi * (p0 * ss + q0 * cc) / sqrt(x));
391 }
392 if (ix <= 0x219392ef) /* x <= 1.0e-18 */
393 return (u0[0] + tpi * log(x));
394 z = x * x;
395 u = u0[12];
396 for (i = 11; i >= 0; i--)
397 u = u * z + u0[i];
398 v = v0[0] + z * (v0[1] + z * (v0[2] + z * (v0[3] + z * v0[4])));
399 return (u / v + tpi * (__k_j0f(fx) * log(x)));
400 }
401
402 float
403 y0f(float fx)
404 {
405 float f;
406 int ix;
407 #if defined(__i386) && !defined(__amd64)
408 int rp;
409 #endif
410
411 ix = *(int *)&fx;
412 if ((ix & ~0x80000000) > 0x7f800000) /* nan */
413 return (fx * fx);
414 if (ix <= 0) { /* zero or negative */
415 if ((ix << 1) == 0)
416 return (-onef / zerof);
417 return (zerof / zerof);
418 }
419 if (ix == 0x7f800000) /* +inf */
420 return (zerof);
421
422 #if defined(__i386) && !defined(__amd64)
423 rp = __swapRP(fp_extended);
424 #endif
425 f = (float)__k_y0f(fx);
426 #if defined(__i386) && !defined(__amd64)
427 if (rp != fp_extended)
428 (void) __swapRP(rp);
429 #endif
430 return (f);
431 }
432
433 /* core of j1f computation; assumes fx is finite */
434 static double
435 __k_j1f(float fx)
436 {
437 double x, z, s, c, ss, cc, r, t, p1, q1;
438 int i, ix, sgn;
439
440 ix = *(int *)&fx;
441 sgn = (unsigned)ix >> 31;
442 ix &= ~0x80000000;
443 x = fabs((double)fx);
444 if (ix > 0x41000000) {
445 /* x > 8; see comments in j1.c */
446 s = sin(x);
447 c = cos(x);
448 if (signbit(s) != signbit(c)) {
449 cc = s - c;
450 ss = cos(x + x) / cc;
451 } else {
452 ss = -s - c;
453 cc = cos(x + x) / ss;
454 }
455 if (ix > 0x501502f9) {
456 /* x > 1.0e10 */
457 p1 = one;
458 q1 = three8 / x;
459 } else {
460 t = eight / x;
461 z = t * t;
462 p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z *
463 (pr0[3] + z * (pr0[4] + z * pr0[5]))))) /
464 (ps0[0] + z * (ps0[1] + z * (ps0[2] + z *
465 (ps0[3] + z * (ps0[4] + z * (ps0[5] + z))))));
466 q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
467 (qr0[3] + z * (qr0[4] + z * qr0[5]))))) /
468 (qs0[0] + z * (qs0[1] + z * (qs0[2] + z *
469 (qs0[3] + z * (qs0[4] + z * (qs0[5] + z))))))) * t;
470 }
471 t = isqrtpi * (p1 * cc - q1 * ss) / sqrt(x);
472 return ((sgn)? -t : t);
473 }
474 if (ix <= 0x3727c5ac) {
475 /* x <= 1.0e-5 */
476 if (ix <= 0x219392ef) /* x <= 1.0e-18 */
477 t = half * x;
478 else
479 t = x * (half + neighth * x * x);
480 return ((sgn)? -t : t);
481 }
482 z = x * x;
483 if (ix < 0x3fa3d70a) {
484 /* x < 1.28 */
485 r = a0[0] + z * (a0[1] + z * (a0[2] + z * a0[3]));
486 s = b0[0] + z * (b0[1] + z * (b0[2] + z * b0[3]));
487 t = x * half + x * (z * (r / s));
488 } else {
489 r = a1[11];
490 for (i = 10; i >= 0; i--)
491 r = r * z + a1[i];
492 s = b1[0] + z * (b1[1] + z * (b1[2] + z * (b1[3] + z * b1[4])));
493 t = x * (r / s);
494 }
495 return ((sgn)? -t : t);
496 }
497
498 float
499 j1f(float fx)
500 {
501 float f;
502 int ix;
503 #if defined(__i386) && !defined(__amd64)
504 int rp;
505 #endif
506
507 ix = *(int *)&fx & ~0x80000000;
508 if (ix >= 0x7f800000) /* nan or inf */
509 return (onef / fx);
510
511 #if defined(__i386) && !defined(__amd64)
512 rp = __swapRP(fp_extended);
513 #endif
514 f = (float)__k_j1f(fx);
515 #if defined(__i386) && !defined(__amd64)
516 if (rp != fp_extended)
517 (void) __swapRP(rp);
518 #endif
519 return (f);
520 }
521
522 /* core of y1f computation; assumes fx is finite and positive */
523 static double
524 __k_y1f(float fx)
525 {
526 double x, z, s, c, ss, cc, u, v, p1, q1, t;
527 int i, ix;
528
529 ix = *(int *)&fx;
530 x = (double)fx;
531 if (ix > 0x41000000) {
532 /* x > 8; see comments in j1.c */
533 s = sin(x);
534 c = cos(x);
535 if (signbit(s) != signbit(c)) {
536 cc = s - c;
537 ss = cos(x + x) / cc;
538 } else {
539 ss = -s - c;
540 cc = cos(x + x) / ss;
541 }
542 if (ix > 0x501502f9) {
543 /* x > 1.0e10 */
544 p1 = one;
545 q1 = three8 / x;
546 } else {
547 t = eight / x;
548 z = t * t;
549 p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z *
550 (pr0[3] + z * (pr0[4] + z * pr0[5]))))) /
551 (ps0[0] + z * (ps0[1] + z * (ps0[2] + z *
552 (ps0[3] + z * (ps0[4] + z * (ps0[5] + z))))));
553 q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
554 (qr0[3] + z * (qr0[4] + z * qr0[5]))))) /
555 (qs0[0] + z * (qs0[1] + z * (qs0[2] + z *
556 (qs0[3] + z * (qs0[4] + z * (qs0[5] + z))))))) * t;
557 }
558 return (isqrtpi * (p1 * ss + q1 * cc) / sqrt(x));
559 }
560 if (ix <= 0x219392ef) /* x <= 1.0e-18 */
561 return (-tpi / x);
562 z = x * x;
563 if (ix < 0x3fa3d70a) {
564 /* x < 1.28 */
565 u = c0[0] + z * (c0[1] + z * (c0[2] + z * c0[3]));
566 v = d0[0] + z * (d0[1] + z * (d0[2] + z * (d0[3] + z * d0[4])));
567 } else {
568 u = c1[11];
569 for (i = 10; i >= 0; i--)
570 u = u * z + c1[i];
571 v = d1[0] + z * (d1[1] + z * (d1[2] + z * (d1[3] + z * d1[4])));
572 }
573 return (x * (u / v) + tpi * (__k_j1f(fx) * log(x) - one / x));
574 }
575
576 float
577 y1f(float fx)
578 {
579 float f;
580 int ix;
581 #if defined(__i386) && !defined(__amd64)
582 int rp;
583 #endif
584
585 ix = *(int *)&fx;
586 if ((ix & ~0x80000000) > 0x7f800000) /* nan */
587 return (fx * fx);
588 if (ix <= 0) { /* zero or negative */
589 if ((ix << 1) == 0)
590 return (-onef / zerof);
591 return (zerof / zerof);
592 }
593 if (ix == 0x7f800000) /* +inf */
594 return (zerof);
595
596 #if defined(__i386) && !defined(__amd64)
597 rp = __swapRP(fp_extended);
598 #endif
599 f = (float)__k_y1f(fx);
600 #if defined(__i386) && !defined(__amd64)
601 if (rp != fp_extended)
602 (void) __swapRP(rp);
603 #endif
604 return (f);
605 }
606
607 float
608 jnf(int n, float fx)
609 {
610 double a, b, temp, x, z, w, t, q0, q1, h;
611 float f;
612 int i, ix, sgn, m, k;
613 #if defined(__i386) && !defined(__amd64)
614 int rp;
615 #endif
616
617 if (n < 0) {
618 n = -n;
619 fx = -fx;
620 }
621 if (n == 0)
622 return (j0f(fx));
623 if (n == 1)
624 return (j1f(fx));
625
626 ix = *(int *)&fx;
627 sgn = (n & 1)? ((unsigned)ix >> 31) : 0;
628 ix &= ~0x80000000;
629 if (ix >= 0x7f800000) { /* nan or inf */
630 if (ix > 0x7f800000)
631 return (fx * fx);
632 return ((sgn)? -zerof : zerof);
633 }
634 if ((ix << 1) == 0)
635 return ((sgn)? -zerof : zerof);
636
637 #if defined(__i386) && !defined(__amd64)
638 rp = __swapRP(fp_extended);
639 #endif
640 fx = fabsf(fx);
641 x = (double)fx;
642 if ((double)n <= x) {
643 /* safe to use J(n+1,x) = 2n/x * J(n,x) - J(n-1,x) */
644 a = __k_j0f(fx);
645 b = __k_j1f(fx);
646 for (i = 1; i < n; i++) {
647 temp = b;
648 b = b * ((double)(i + i) / x) - a;
649 a = temp;
650 }
651 f = (float)b;
652 #if defined(__i386) && !defined(__amd64)
653 if (rp != fp_extended)
654 (void) __swapRP(rp);
655 #endif
656 return ((sgn)? -f : f);
657 }
658 if (ix < 0x3089705f) {
659 /* x < 1.0e-9; use J(n,x) = 1/n! * (x / 2)^n */
660 if (n > 6)
661 n = 6; /* result underflows to zero for n >= 6 */
662 b = t = half * x;
663 a = one;
664 for (i = 2; i <= n; i++) {
665 b *= t;
666 a *= (double)i;
667 }
668 b /= a;
669 } else {
670 /*
671 * Use the backward recurrence:
672 *
673 * x x^2 x^2
674 * J(n,x)/J(n-1,x) = ---- - ------ - ------ .....
675 * 2n 2(n+1) 2(n+2)
676 *
677 * Let w = 2n/x and h = 2/x. Then the above quotient
678 * is equal to the continued fraction:
679 * 1
680 * = -----------------------
681 * 1
682 * w - -----------------
683 * 1
684 * w+h - ---------
685 * w+2h - ...
686 *
687 * To determine how many terms are needed, run the
688 * recurrence
689 *
690 * Q(0) = w,
691 * Q(1) = w(w+h) - 1,
692 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2).
693 *
694 * Then when Q(k) > 1e4, k is large enough for single
695 * precision.
696 */
697 /* XXX NOT DONE - rework this */
698 w = (n + n) / x;
699 h = two / x;
700 q0 = w;
701 z = w + h;
702 q1 = w * z - one;
703 k = 1;
704 while (q1 < big) {
705 k++;
706 z += h;
707 temp = z * q1 - q0;
708 q0 = q1;
709 q1 = temp;
710 }
711 m = n + n;
712 t = zero;
713 for (i = (n + k) << 1; i >= m; i -= 2)
714 t = one / ((double)i / x - t);
715 a = t;
716 b = one;
717 /*
718 * estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
719 * hence, if n*(log(2n/x)) > ...
720 * single 8.8722839355e+01
721 * double 7.09782712893383973096e+02
722 * then recurrent value may overflow and the result is
723 * likely underflow to zero
724 */
725 temp = (double)n;
726 temp *= log((two / x) * temp);
727 if (temp < 7.09782712893383973096e+02) {
728 for (i = n - 1; i > 0; i--) {
729 temp = b;
730 b = b * ((double)(i + i) / x) - a;
731 a = temp;
732 }
733 } else {
734 for (i = n - 1; i > 0; i--) {
735 temp = b;
736 b = b * ((double)(i + i) / x) - a;
737 a = temp;
738 if (b > 1.0e100) {
739 a /= b;
740 t /= b;
741 b = one;
742 }
743 }
744 }
745 b = (t * __k_j0f(fx) / b);
746 }
747 f = (float)b;
748 #if defined(__i386) && !defined(__amd64)
749 if (rp != fp_extended)
750 (void) __swapRP(rp);
751 #endif
752 return ((sgn)? -f : f);
753 }
754
755 float
756 ynf(int n, float fx)
757 {
758 double a, b, temp, x;
759 float f;
760 int i, sign, ix;
761 #if defined(__i386) && !defined(__amd64)
762 int rp;
763 #endif
764
765 sign = 0;
766 if (n < 0) {
767 n = -n;
768 if (n & 1)
769 sign = 1;
770 }
771 if (n == 0)
772 return (y0f(fx));
773 if (n == 1)
774 return ((sign)? -y1f(fx) : y1f(fx));
775
776 ix = *(int *)&fx;
777 if ((ix & ~0x80000000) > 0x7f800000) /* nan */
778 return (fx * fx);
779 if (ix <= 0) { /* zero or negative */
780 if ((ix << 1) == 0)
781 return (-onef / zerof);
782 return (zerof / zerof);
783 }
784 if (ix == 0x7f800000) /* +inf */
785 return (zerof);
786
787 #if defined(__i386) && !defined(__amd64)
788 rp = __swapRP(fp_extended);
789 #endif
790 a = __k_y0f(fx);
791 b = __k_y1f(fx);
792 x = (double)fx;
793 for (i = 1; i < n; i++) {
794 temp = b;
795 b *= (double)(i + i) / x;
796 if (b <= -DBL_MAX)
797 break;
798 b -= a;
799 a = temp;
800 }
801 f = (float)b;
802 #if defined(__i386) && !defined(__amd64)
803 if (rp != fp_extended)
804 (void) __swapRP(rp);
805 #endif
806 return ((sign)? -f : f);
807 }
|
1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25
26 /*
27 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __j0f = j0f
32 #pragma weak __j1f = j1f
33 #pragma weak __jnf = jnf
34 #pragma weak __y0f = y0f
35 #pragma weak __y1f = y1f
36 #pragma weak __ynf = ynf
37
38 #include "libm.h"
39 #include <float.h>
40
41 #if defined(__i386) && !defined(__amd64)
42 extern int __swapRP(int);
43 #endif
44
45 static const float zerof = 0.0f, onef = 1.0f;
46
47 static const double C[] = {
48 0.0,
49 -0.125,
50 0.25,
51 0.375,
52 0.5,
53 1.0,
54 2.0,
55 8.0,
56 0.5641895835477562869480794515607725858441, /* 1/sqrt(pi) */
57 0.636619772367581343075535053490057448, /* 2/pi */
58 1.0e9,
59 };
60
61 #define zero C[0]
62 #define neighth C[1]
63 #define quarter C[2]
64 #define three8 C[3]
65 #define half C[4]
85 0.5068271181053546392490184353e4,
86 0.1829817905472769960535671664e3,
87 1.0,
88 -0.1731210995701068539185611951e3, /* qr */
89 -0.5522559165936166961235240613e3,
90 -0.5604935606637346590614529613e3,
91 -0.2200430300226009379477365011e3,
92 -0.323869355375648849771296746e2,
93 -0.14294979207907956223499258e1,
94 -0.834690374102384988158918e-2,
95 0.1107975037248683865326709645e5, /* qs */
96 0.3544581680627082674651471873e5,
97 0.3619118937918394132179019059e5,
98 0.1439895563565398007471485822e5,
99 0.2190277023344363955930226234e4,
100 0.106695157020407986137501682e3,
101 1.0,
102 };
103
104 #define pr Cj0y0
105 #define ps (Cj0y0 + 7)
106 #define qr (Cj0y0 + 14)
107 #define qs (Cj0y0 + 21)
108
109 static const double Cj0[] = {
110 -2.500000000000003622131880894830476755537e-0001, /* r0 */
111 1.095597547334830263234433855932375353303e-0002,
112 -1.819734750463320921799187258987098087697e-0004,
113 9.977001946806131657544212501069893930846e-0007,
114 1.0, /* s0 */
115 1.867609810662950169966782360588199673741e-0002,
116 1.590389206181565490878430827706972074208e-0004,
117 6.520867386742583632375520147714499522721e-0007,
118 9.999999999999999942156495584397047660949e-0001, /* r1 */
119 -2.389887722731319130476839836908143731281e-0001,
120 1.293359476138939027791270393439493640570e-0002,
121 -2.770985642343140122168852400228563364082e-0004,
122 2.905241575772067678086738389169625218912e-0006,
123 -1.636846356264052597969042009265043251279e-0008,
124 5.072306160724884775085431059052611737827e-0011,
125 -8.187060730684066824228914775146536139112e-0014,
126 5.422219326959949863954297860723723423842e-0017,
127 1.0, /* s1 */
128 1.101122772686807702762104741932076228349e-0002,
129 6.140169310641649223411427764669143978228e-0005,
130 2.292035877515152097976946119293215705250e-0007,
131 6.356910426504644334558832036362219583789e-0010,
132 1.366626326900219555045096999553948891401e-0012,
133 2.280399586866739522891837985560481180088e-0015,
134 2.801559820648939665270492520004836611187e-0018,
135 2.073101088320349159764410261466350732968e-0021,
136 };
137
138 #define r0 Cj0
139 #define s0 (Cj0 + 4)
140 #define r1 (Cj0 + 8)
141 #define s1 (Cj0 + 17)
142
143 static const double Cy0[] = {
144 -7.380429510868722526754723020704317641941e-0002, /* u0 */
145 1.772607102684869924301459663049874294814e-0001,
146 -1.524370666542713828604078090970799356306e-0002,
147 4.650819100693891757143771557629924591915e-0004,
148 -7.125768872339528975036316108718239946022e-0006,
149 6.411017001656104598327565004771515257146e-0008,
150 -3.694275157433032553021246812379258781665e-0010,
151 1.434364544206266624252820889648445263842e-0012,
152 -3.852064731859936455895036286874139896861e-0015,
153 7.182052899726138381739945881914874579696e-0018,
154 -9.060556574619677567323741194079797987200e-0021,
155 7.124435467408860515265552217131230511455e-0024,
156 -2.709726774636397615328813121715432044771e-0027,
157 1.0, /* v0 */
158 4.678678931512549002587702477349214886475e-0003,
159 9.486828955529948534822800829497565178985e-0006,
160 1.001495929158861646659010844136682454906e-0008,
161 4.725338116256021660204443235685358593611e-0012,
162 };
163
164 #define u0 Cy0
165 #define v0 (Cy0 + 13)
166
167 static const double Cj1y1[] = {
168 -0.4435757816794127857114720794e7, /* pr0 */
169 -0.9942246505077641195658377899e7,
170 -0.6603373248364939109255245434e7,
171 -0.1523529351181137383255105722e7,
172 -0.1098240554345934672737413139e6,
173 -0.1611616644324610116477412898e4,
174 -0.4435757816794127856828016962e7, /* ps0 */
175 -0.9934124389934585658967556309e7,
176 -0.6585339479723087072826915069e7,
177 -0.1511809506634160881644546358e7,
178 -0.1072638599110382011903063867e6,
179 -0.1455009440190496182453565068e4,
180 0.3322091340985722351859704442e5, /* qr0 */
181 0.8514516067533570196555001171e5,
182 0.6617883658127083517939992166e5,
183 0.1849426287322386679652009819e5,
184 0.1706375429020768002061283546e4,
185 0.3526513384663603218592175580e2,
186 0.7087128194102874357377502472e6, /* qs0 */
187 0.1819458042243997298924553839e7,
188 0.1419460669603720892855755253e7,
189 0.4002944358226697511708610813e6,
190 0.3789022974577220264142952256e5,
191 0.8638367769604990967475517183e3,
192 };
193
194 #define pr0 Cj1y1
195 #define ps0 (Cj1y1 + 6)
196 #define qr0 (Cj1y1 + 12)
197 #define qs0 (Cj1y1 + 18)
198
199 static const double Cj1[] = {
200 -6.250000000000002203053200981413218949548e-0002, /* a0 */
201 1.600998455640072901321605101981501263762e-0003,
202 -1.963888815948313758552511884390162864930e-0005,
203 8.263917341093549759781339713418201620998e-0008,
204 1.0e0, /* b0 */
205 1.605069137643004242395356851797873766927e-0002,
206 1.149454623251299996428500249509098499383e-0004,
207 3.849701673735260970379681807910852327825e-0007,
208 4.999999999999999995517408894340485471724e-0001,
209 -6.003825028120475684835384519945468075423e-0002,
210 2.301719899263321828388344461995355419832e-0003,
211 -4.208494869238892934859525221654040304068e-0005,
212 4.377745135188837783031540029700282443388e-0007,
213 -2.854106755678624335145364226735677754179e-0009,
214 1.234002865443952024332943901323798413689e-0011,
215 -3.645498437039791058951273508838177134310e-0014,
216 7.404320596071797459925377103787837414422e-0017,
217 -1.009457448277522275262808398517024439084e-0019,
218 8.520158355824819796968771418801019930585e-0023,
219 -3.458159926081163274483854614601091361424e-0026,
220 1.0e0, /* b1 */
221 4.923499437590484879081138588998986303306e-0003,
222 1.054389489212184156499666953501976688452e-0005,
223 1.180768373106166527048240364872043816050e-0008,
224 5.942665743476099355323245707680648588540e-0012,
225 };
226
227 #define a0 Cj1
228 #define b0 (Cj1 + 4)
229 #define a1 (Cj1 + 8)
230 #define b1 (Cj1 + 20)
231
232 static const double Cy1[] = {
233 -1.960570906462389461018983259589655961560e-0001, /* c0 */
234 4.931824118350661953459180060007970291139e-0002,
235 -1.626975871565393656845930125424683008677e-0003,
236 1.359657517926394132692884168082224258360e-0005,
237 1.0e0, /* d0 */
238 2.565807214838390835108224713630901653793e-0002,
239 3.374175208978404268650522752520906231508e-0004,
240 2.840368571306070719539936935220728843177e-0006,
241 1.396387402048998277638900944415752207592e-0008,
242 -1.960570906462389473336339614647555351626e-0001, /* c1 */
243 5.336268030335074494231369159933012844735e-0002,
244 -2.684137504382748094149184541866332033280e-0003,
245 5.737671618979185736981543498580051903060e-0005,
246 -6.642696350686335339171171785557663224892e-0007,
247 4.692417922568160354012347591960362101664e-0009,
248 -2.161728635907789319335231338621412258355e-0011,
249 6.727353419738316107197644431844194668702e-0014,
250 -1.427502986803861372125234355906790573422e-0016,
251 2.020392498726806769468143219616642940371e-0019,
252 -1.761371948595104156753045457888272716340e-0022,
253 7.352828391941157905175042420249225115816e-0026,
254 1.0e0, /* d1 */
255 5.029187436727947764916247076102283399442e-0003,
256 1.102693095808242775074856548927801750627e-0005,
257 1.268035774543174837829534603830227216291e-0008,
258 6.579416271766610825192542295821308730206e-0012,
259 };
260
261 #define c0 Cy1
262 #define d0 (Cy1 + 4)
263 #define c1 (Cy1 + 9)
264 #define d1 (Cy1 + 21)
265
266 /* core of j0f computation; assumes fx is finite */
267 static double
268 __k_j0f(float fx)
269 {
270 double x, z, s, c, ss, cc, r, t, p0, q0;
271 int ix, i;
272
273 ix = *(int *)&fx & ~0x80000000;
274 x = fabs((double)fx);
275
276 if (ix > 0x41000000) {
277 /* x > 8; see comments in j0.c */
278 s = sin(x);
279 c = cos(x);
280
281 if (signbit(s) != signbit(c)) {
282 ss = s - c;
283 cc = -cos(x + x) / ss;
284 } else {
285 cc = s + c;
286 ss = -cos(x + x) / cc;
287 }
288
289 if (ix > 0x501502f9) {
290 /* x > 1.0e10 */
291 p0 = one;
292 q0 = neighth / x;
293 } else {
294 t = eight / x;
295 z = t * t;
296 p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] + z *
297 (pr[4] + z * (pr[5] + z * pr[6])))))) / (ps[0] + z *
298 (ps[1] + z * (ps[2] + z * (ps[3] + z * (ps[4] + z *
299 (ps[5] + z))))));
300 q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
301 z * (qr[4] + z * (qr[5] + z * qr[6])))))) / (qs[0] +
302 z * (qs[1] + z * (qs[2] + z * (qs[3] + z * (qs[4] +
303 z * (qs[5] + z))))))) * t;
304 }
305
306 return (isqrtpi * (p0 * cc - q0 * ss) / sqrt(x));
307 }
308
309 if (ix <= 0x3727c5ac) {
310 /* x <= 1.0e-5 */
311 if (ix <= 0x219392ef) /* x <= 1.0e-18 */
312 return (one - x);
313
314 return (one - x * x * quarter);
315 }
316
317 z = x * x;
318
319 if (ix <= 0x3fa3d70a) {
320 /* x <= 1.28 */
321 r = r0[0] + z * (r0[1] + z * (r0[2] + z * r0[3]));
322 s = s0[0] + z * (s0[1] + z * (s0[2] + z * s0[3]));
323 return (one + z * (r / s));
324 }
325
326 r = r1[8];
327 s = s1[8];
328
329 for (i = 7; i >= 0; i--) {
330 r = r * z + r1[i];
331 s = s * z + s1[i];
332 }
333
334 return (r / s);
335 }
336
337 float
338 j0f(float fx)
339 {
340 float f;
341 int ix;
342
343 #if defined(__i386) && !defined(__amd64)
344 int rp;
345 #endif
346
347 ix = *(int *)&fx & ~0x80000000;
348
349 if (ix >= 0x7f800000) { /* nan or inf */
350 if (ix > 0x7f800000)
351 return (fx * fx);
352
353 return (zerof);
354 }
355
356 #if defined(__i386) && !defined(__amd64)
357 rp = __swapRP(fp_extended);
358 #endif
359 f = (float)__k_j0f(fx);
360 #if defined(__i386) && !defined(__amd64)
361 if (rp != fp_extended)
362 (void) __swapRP(rp);
363 #endif
364 return (f);
365 }
366
367 /* core of y0f computation; assumes fx is finite and positive */
368 static double
369 __k_y0f(float fx)
370 {
371 double x, z, s, c, ss, cc, t, p0, q0, u, v;
372 int ix, i;
373
374 ix = *(int *)&fx;
375 x = (double)fx;
376
377 if (ix > 0x41000000) {
378 /* x > 8; see comments in j0.c */
379 s = sin(x);
380 c = cos(x);
381
382 if (signbit(s) != signbit(c)) {
383 ss = s - c;
384 cc = -cos(x + x) / ss;
385 } else {
386 cc = s + c;
387 ss = -cos(x + x) / cc;
388 }
389
390 if (ix > 0x501502f9) {
391 /* x > 1.0e10 */
392 p0 = one;
393 q0 = neighth / x;
394 } else {
395 t = eight / x;
396 z = t * t;
397 p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] + z *
398 (pr[4] + z * (pr[5] + z * pr[6])))))) / (ps[0] + z *
399 (ps[1] + z * (ps[2] + z * (ps[3] + z * (ps[4] + z *
400 (ps[5] + z))))));
401 q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
402 z * (qr[4] + z * (qr[5] + z * qr[6])))))) / (qs[0] +
403 z * (qs[1] + z * (qs[2] + z * (qs[3] + z * (qs[4] +
404 z * (qs[5] + z))))))) * t;
405 }
406
407 return (isqrtpi * (p0 * ss + q0 * cc) / sqrt(x));
408 }
409
410 if (ix <= 0x219392ef) /* x <= 1.0e-18 */
411 return (u0[0] + tpi * log(x));
412
413 z = x * x;
414 u = u0[12];
415
416 for (i = 11; i >= 0; i--)
417 u = u * z + u0[i];
418
419 v = v0[0] + z * (v0[1] + z * (v0[2] + z * (v0[3] + z * v0[4])));
420 return (u / v + tpi * (__k_j0f(fx) * log(x)));
421 }
422
423 float
424 y0f(float fx)
425 {
426 float f;
427 int ix;
428
429 #if defined(__i386) && !defined(__amd64)
430 int rp;
431 #endif
432
433 ix = *(int *)&fx;
434
435 if ((ix & ~0x80000000) > 0x7f800000) /* nan */
436 return (fx * fx);
437
438 if (ix <= 0) { /* zero or negative */
439 if ((ix << 1) == 0)
440 return (-onef / zerof);
441
442 return (zerof / zerof);
443 }
444
445 if (ix == 0x7f800000) /* +inf */
446 return (zerof);
447
448 #if defined(__i386) && !defined(__amd64)
449 rp = __swapRP(fp_extended);
450 #endif
451 f = (float)__k_y0f(fx);
452 #if defined(__i386) && !defined(__amd64)
453 if (rp != fp_extended)
454 (void) __swapRP(rp);
455 #endif
456 return (f);
457 }
458
459 /* core of j1f computation; assumes fx is finite */
460 static double
461 __k_j1f(float fx)
462 {
463 double x, z, s, c, ss, cc, r, t, p1, q1;
464 int i, ix, sgn;
465
466 ix = *(int *)&fx;
467 sgn = (unsigned)ix >> 31;
468 ix &= ~0x80000000;
469 x = fabs((double)fx);
470
471 if (ix > 0x41000000) {
472 /* x > 8; see comments in j1.c */
473 s = sin(x);
474 c = cos(x);
475
476 if (signbit(s) != signbit(c)) {
477 cc = s - c;
478 ss = cos(x + x) / cc;
479 } else {
480 ss = -s - c;
481 cc = cos(x + x) / ss;
482 }
483
484 if (ix > 0x501502f9) {
485 /* x > 1.0e10 */
486 p1 = one;
487 q1 = three8 / x;
488 } else {
489 t = eight / x;
490 z = t * t;
491 p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z * (pr0[3] +
492 z * (pr0[4] + z * pr0[5]))))) / (ps0[0] + z *
493 (ps0[1] + z * (ps0[2] + z * (ps0[3] + z * (ps0[4] +
494 z * (ps0[5] + z))))));
495 q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
496 (qr0[3] + z * (qr0[4] + z * qr0[5]))))) / (qs0[0] +
497 z * (qs0[1] + z * (qs0[2] + z * (qs0[3] + z *
498 (qs0[4] + z * (qs0[5] + z))))))) * t;
499 }
500
501 t = isqrtpi * (p1 * cc - q1 * ss) / sqrt(x);
502 return ((sgn) ? -t : t);
503 }
504
505 if (ix <= 0x3727c5ac) {
506 /* x <= 1.0e-5 */
507 if (ix <= 0x219392ef) /* x <= 1.0e-18 */
508 t = half * x;
509 else
510 t = x * (half + neighth * x * x);
511
512 return ((sgn) ? -t : t);
513 }
514
515 z = x * x;
516
517 if (ix < 0x3fa3d70a) {
518 /* x < 1.28 */
519 r = a0[0] + z * (a0[1] + z * (a0[2] + z * a0[3]));
520 s = b0[0] + z * (b0[1] + z * (b0[2] + z * b0[3]));
521 t = x * half + x * (z * (r / s));
522 } else {
523 r = a1[11];
524
525 for (i = 10; i >= 0; i--)
526 r = r * z + a1[i];
527
528 s = b1[0] + z * (b1[1] + z * (b1[2] + z * (b1[3] + z * b1[4])));
529 t = x * (r / s);
530 }
531
532 return ((sgn) ? -t : t);
533 }
534
535 float
536 j1f(float fx)
537 {
538 float f;
539 int ix;
540
541 #if defined(__i386) && !defined(__amd64)
542 int rp;
543 #endif
544
545 ix = *(int *)&fx & ~0x80000000;
546
547 if (ix >= 0x7f800000) /* nan or inf */
548 return (onef / fx);
549
550 #if defined(__i386) && !defined(__amd64)
551 rp = __swapRP(fp_extended);
552 #endif
553 f = (float)__k_j1f(fx);
554 #if defined(__i386) && !defined(__amd64)
555 if (rp != fp_extended)
556 (void) __swapRP(rp);
557 #endif
558 return (f);
559 }
560
561 /* core of y1f computation; assumes fx is finite and positive */
562 static double
563 __k_y1f(float fx)
564 {
565 double x, z, s, c, ss, cc, u, v, p1, q1, t;
566 int i, ix;
567
568 ix = *(int *)&fx;
569 x = (double)fx;
570
571 if (ix > 0x41000000) {
572 /* x > 8; see comments in j1.c */
573 s = sin(x);
574 c = cos(x);
575
576 if (signbit(s) != signbit(c)) {
577 cc = s - c;
578 ss = cos(x + x) / cc;
579 } else {
580 ss = -s - c;
581 cc = cos(x + x) / ss;
582 }
583
584 if (ix > 0x501502f9) {
585 /* x > 1.0e10 */
586 p1 = one;
587 q1 = three8 / x;
588 } else {
589 t = eight / x;
590 z = t * t;
591 p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z * (pr0[3] +
592 z * (pr0[4] + z * pr0[5]))))) / (ps0[0] + z *
593 (ps0[1] + z * (ps0[2] + z * (ps0[3] + z * (ps0[4] +
594 z * (ps0[5] + z))))));
595 q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
596 (qr0[3] + z * (qr0[4] + z * qr0[5]))))) / (qs0[0] +
597 z * (qs0[1] + z * (qs0[2] + z * (qs0[3] + z *
598 (qs0[4] + z * (qs0[5] + z))))))) * t;
599 }
600
601 return (isqrtpi * (p1 * ss + q1 * cc) / sqrt(x));
602 }
603
604 if (ix <= 0x219392ef) /* x <= 1.0e-18 */
605 return (-tpi / x);
606
607 z = x * x;
608
609 if (ix < 0x3fa3d70a) {
610 /* x < 1.28 */
611 u = c0[0] + z * (c0[1] + z * (c0[2] + z * c0[3]));
612 v = d0[0] + z * (d0[1] + z * (d0[2] + z * (d0[3] + z * d0[4])));
613 } else {
614 u = c1[11];
615
616 for (i = 10; i >= 0; i--)
617 u = u * z + c1[i];
618
619 v = d1[0] + z * (d1[1] + z * (d1[2] + z * (d1[3] + z * d1[4])));
620 }
621
622 return (x * (u / v) + tpi * (__k_j1f(fx) * log(x) - one / x));
623 }
624
625 float
626 y1f(float fx)
627 {
628 float f;
629 int ix;
630
631 #if defined(__i386) && !defined(__amd64)
632 int rp;
633 #endif
634
635 ix = *(int *)&fx;
636
637 if ((ix & ~0x80000000) > 0x7f800000) /* nan */
638 return (fx * fx);
639
640 if (ix <= 0) { /* zero or negative */
641 if ((ix << 1) == 0)
642 return (-onef / zerof);
643
644 return (zerof / zerof);
645 }
646
647 if (ix == 0x7f800000) /* +inf */
648 return (zerof);
649
650 #if defined(__i386) && !defined(__amd64)
651 rp = __swapRP(fp_extended);
652 #endif
653 f = (float)__k_y1f(fx);
654 #if defined(__i386) && !defined(__amd64)
655 if (rp != fp_extended)
656 (void) __swapRP(rp);
657 #endif
658 return (f);
659 }
660
661 float
662 jnf(int n, float fx)
663 {
664 double a, b, temp, x, z, w, t, q0, q1, h;
665 float f;
666 int i, ix, sgn, m, k;
667
668 #if defined(__i386) && !defined(__amd64)
669 int rp;
670 #endif
671
672 if (n < 0) {
673 n = -n;
674 fx = -fx;
675 }
676
677 if (n == 0)
678 return (j0f(fx));
679
680 if (n == 1)
681 return (j1f(fx));
682
683 ix = *(int *)&fx;
684 sgn = (n & 1) ? ((unsigned)ix >> 31) : 0;
685 ix &= ~0x80000000;
686
687 if (ix >= 0x7f800000) { /* nan or inf */
688 if (ix > 0x7f800000)
689 return (fx * fx);
690
691 return ((sgn) ? -zerof : zerof);
692 }
693
694 if ((ix << 1) == 0)
695 return ((sgn) ? -zerof : zerof);
696
697 #if defined(__i386) && !defined(__amd64)
698 rp = __swapRP(fp_extended);
699 #endif
700 fx = fabsf(fx);
701 x = (double)fx;
702
703 if ((double)n <= x) {
704 /* safe to use J(n+1,x) = 2n/x * J(n,x) - J(n-1,x) */
705 a = __k_j0f(fx);
706 b = __k_j1f(fx);
707
708 for (i = 1; i < n; i++) {
709 temp = b;
710 b = b * ((double)(i + i) / x) - a;
711 a = temp;
712 }
713
714 f = (float)b;
715 #if defined(__i386) && !defined(__amd64)
716 if (rp != fp_extended)
717 (void) __swapRP(rp);
718 #endif
719 return ((sgn) ? -f : f);
720 }
721
722 if (ix < 0x3089705f) {
723 /* x < 1.0e-9; use J(n,x) = 1/n! * (x / 2)^n */
724 if (n > 6)
725 n = 6; /* result underflows to zero for n >= 6 */
726
727 b = t = half * x;
728 a = one;
729
730 for (i = 2; i <= n; i++) {
731 b *= t;
732 a *= (double)i;
733 }
734
735 b /= a;
736 } else {
737 /* BEGIN CSTYLED */
738 /*
739 * Use the backward recurrence:
740 *
741 * x x^2 x^2
742 * J(n,x)/J(n-1,x) = ---- - ------ - ------ .....
743 * 2n 2(n+1) 2(n+2)
744 *
745 * Let w = 2n/x and h = 2/x. Then the above quotient
746 * is equal to the continued fraction:
747 * 1
748 * = -----------------------
749 * 1
750 * w - -----------------
751 * 1
752 * w+h - ---------
753 * w+2h - ...
754 *
755 * To determine how many terms are needed, run the
756 * recurrence
757 *
758 * Q(0) = w,
759 * Q(1) = w(w+h) - 1,
760 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2).
761 *
762 * Then when Q(k) > 1e4, k is large enough for single
763 * precision.
764 */
765 /* END CSTYLED */
766 /* XXX NOT DONE - rework this */
767 w = (n + n) / x;
768 h = two / x;
769 q0 = w;
770 z = w + h;
771 q1 = w * z - one;
772 k = 1;
773
774 while (q1 < big) {
775 k++;
776 z += h;
777 temp = z * q1 - q0;
778 q0 = q1;
779 q1 = temp;
780 }
781
782 m = n + n;
783 t = zero;
784
785 for (i = (n + k) << 1; i >= m; i -= 2)
786 t = one / ((double)i / x - t);
787
788 a = t;
789 b = one;
790
791 /*
792 * estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
793 * hence, if n*(log(2n/x)) > ...
794 * single 8.8722839355e+01
795 * double 7.09782712893383973096e+02
796 * then recurrent value may overflow and the result is
797 * likely underflow to zero
798 */
799 temp = (double)n;
800 temp *= log((two / x) * temp);
801
802 if (temp < 7.09782712893383973096e+02) {
803 for (i = n - 1; i > 0; i--) {
804 temp = b;
805 b = b * ((double)(i + i) / x) - a;
806 a = temp;
807 }
808 } else {
809 for (i = n - 1; i > 0; i--) {
810 temp = b;
811 b = b * ((double)(i + i) / x) - a;
812 a = temp;
813
814 if (b > 1.0e100) {
815 a /= b;
816 t /= b;
817 b = one;
818 }
819 }
820 }
821
822 b = (t * __k_j0f(fx) / b);
823 }
824
825 f = (float)b;
826 #if defined(__i386) && !defined(__amd64)
827 if (rp != fp_extended)
828 (void) __swapRP(rp);
829 #endif
830 return ((sgn) ? -f : f);
831 }
832
833 float
834 ynf(int n, float fx)
835 {
836 double a, b, temp, x;
837 float f;
838 int i, sign, ix;
839
840 #if defined(__i386) && !defined(__amd64)
841 int rp;
842 #endif
843
844 sign = 0;
845
846 if (n < 0) {
847 n = -n;
848
849 if (n & 1)
850 sign = 1;
851 }
852
853 if (n == 0)
854 return (y0f(fx));
855
856 if (n == 1)
857 return ((sign) ? -y1f(fx) : y1f(fx));
858
859 ix = *(int *)&fx;
860
861 if ((ix & ~0x80000000) > 0x7f800000) /* nan */
862 return (fx * fx);
863
864 if (ix <= 0) { /* zero or negative */
865 if ((ix << 1) == 0)
866 return (-onef / zerof);
867
868 return (zerof / zerof);
869 }
870
871 if (ix == 0x7f800000) /* +inf */
872 return (zerof);
873
874 #if defined(__i386) && !defined(__amd64)
875 rp = __swapRP(fp_extended);
876 #endif
877 a = __k_y0f(fx);
878 b = __k_y1f(fx);
879 x = (double)fx;
880
881 for (i = 1; i < n; i++) {
882 temp = b;
883 b *= (double)(i + i) / x;
884
885 if (b <= -DBL_MAX)
886 break;
887
888 b -= a;
889 a = temp;
890 }
891
892 f = (float)b;
893 #if defined(__i386) && !defined(__amd64)
894 if (rp != fp_extended)
895 (void) __swapRP(rp);
896 #endif
897 return ((sign) ? -f : f);
898 }
|