Deleted Added
sdiff udiff text old ( 279240 ) new ( 279493 )
full compact
e_j1.c (279240) e_j1.c (279493)
1
2/* @(#)e_j1.c 1.3 95/01/18 */
3/*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunSoft, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14#include <sys/cdefs.h>
1
2/* @(#)e_j1.c 1.3 95/01/18 */
3/*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunSoft, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14#include <sys/cdefs.h>
15__FBSDID("$FreeBSD: head/lib/msun/src/e_j1.c 279240 2015-02-24 16:45:16Z pfg $");
15__FBSDID("$FreeBSD: head/lib/msun/src/e_j1.c 279493 2015-03-01 20:32:47Z kargl $");
16
17/* __ieee754_j1(x), __ieee754_y1(x)
18 * Bessel function of the first and second kinds of order zero.
19 * Method -- j1(x):
20 * 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
21 * 2. Reduce x to |x| since j1(x)=-j1(-x), and
22 * for x in (0,2)
23 * j1(x) = x/2 + x*z*R0/S0, where z = x*x;

--- 233 unchanged lines hidden (view full) ---

257static const double ps2[5] = {
258 2.14364859363821409488e+01, /* 0x40356FBD, 0x8AD5ECDC */
259 1.25290227168402751090e+02, /* 0x405F5293, 0x14F92CD5 */
260 2.32276469057162813669e+02, /* 0x406D08D8, 0xD5A2DBD9 */
261 1.17679373287147100768e+02, /* 0x405D6B7A, 0xDA1884A9 */
262 8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */
263};
264
16
17/* __ieee754_j1(x), __ieee754_y1(x)
18 * Bessel function of the first and second kinds of order zero.
19 * Method -- j1(x):
20 * 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
21 * 2. Reduce x to |x| since j1(x)=-j1(-x), and
22 * for x in (0,2)
23 * j1(x) = x/2 + x*z*R0/S0, where z = x*x;

--- 233 unchanged lines hidden (view full) ---

257static const double ps2[5] = {
258 2.14364859363821409488e+01, /* 0x40356FBD, 0x8AD5ECDC */
259 1.25290227168402751090e+02, /* 0x405F5293, 0x14F92CD5 */
260 2.32276469057162813669e+02, /* 0x406D08D8, 0xD5A2DBD9 */
261 1.17679373287147100768e+02, /* 0x405D6B7A, 0xDA1884A9 */
262 8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */
263};
264
265 static double pone(double x)
265static __inline double
266pone(double x)
266{
267 const double *p,*q;
268 double z,r,s;
269 int32_t ix;
270 GET_HIGH_WORD(ix,x);
271 ix &= 0x7fffffff;
272 if(ix>=0x40200000) {p = pr8; q= ps8;}
273 else if(ix>=0x40122E8B){p = pr5; q= ps5;}

--- 79 unchanged lines hidden (view full) ---

353 2.95333629060523854548e+01, /* 0x403D888A, 0x78AE64FF */
354 2.52981549982190529136e+02, /* 0x406F9F68, 0xDB821CBA */
355 7.57502834868645436472e+02, /* 0x4087AC05, 0xCE49A0F7 */
356 7.39393205320467245656e+02, /* 0x40871B25, 0x48D4C029 */
357 1.55949003336666123687e+02, /* 0x40637E5E, 0x3C3ED8D4 */
358 -4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */
359};
360
267{
268 const double *p,*q;
269 double z,r,s;
270 int32_t ix;
271 GET_HIGH_WORD(ix,x);
272 ix &= 0x7fffffff;
273 if(ix>=0x40200000) {p = pr8; q= ps8;}
274 else if(ix>=0x40122E8B){p = pr5; q= ps5;}

--- 79 unchanged lines hidden (view full) ---

354 2.95333629060523854548e+01, /* 0x403D888A, 0x78AE64FF */
355 2.52981549982190529136e+02, /* 0x406F9F68, 0xDB821CBA */
356 7.57502834868645436472e+02, /* 0x4087AC05, 0xCE49A0F7 */
357 7.39393205320467245656e+02, /* 0x40871B25, 0x48D4C029 */
358 1.55949003336666123687e+02, /* 0x40637E5E, 0x3C3ED8D4 */
359 -4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */
360};
361
361 static double qone(double x)
362static __inline double
363qone(double x)
362{
363 const double *p,*q;
364 double s,r,z;
365 int32_t ix;
366 GET_HIGH_WORD(ix,x);
367 ix &= 0x7fffffff;
368 if(ix>=0x40200000) {p = qr8; q= qs8;}
369 else if(ix>=0x40122E8B){p = qr5; q= qs5;}
370 else if(ix>=0x4006DB6D){p = qr3; q= qs3;}
371 else {p = qr2; q= qs2;} /* ix>=0x40000000 */
372 z = one/(x*x);
373 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
374 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
375 return (.375 + r/s)/x;
376}
364{
365 const double *p,*q;
366 double s,r,z;
367 int32_t ix;
368 GET_HIGH_WORD(ix,x);
369 ix &= 0x7fffffff;
370 if(ix>=0x40200000) {p = qr8; q= qs8;}
371 else if(ix>=0x40122E8B){p = qr5; q= qs5;}
372 else if(ix>=0x4006DB6D){p = qr3; q= qs3;}
373 else {p = qr2; q= qs2;} /* ix>=0x40000000 */
374 z = one/(x*x);
375 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
376 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
377 return (.375 + r/s)/x;
378}