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 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 |
265static __inline double 266pone(double x) |
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 |
362static __inline double 363qone(double x) |
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} |