1 2/* @(#)e_j0.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_j0.c 275518 2014-12-05 19:00:55Z kargl $"); |
16 17/* __ieee754_j0(x), __ieee754_y0(x) 18 * Bessel function of the first and second kinds of order zero. 19 * Method -- j0(x): 20 * 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ... 21 * 2. Reduce x to |x| since j0(x)=j0(-x), and 22 * for x in (0,2) 23 * j0(x) = 1-z/4+ z^2*R0/S0, where z = x*x; --- 86 unchanged lines hidden (view full) --- 110 u = pzero(x); v = qzero(x); 111 z = invsqrtpi*(u*cc-v*ss)/sqrt(x); 112 } 113 return z; 114 } 115 if(ix<0x3f200000) { /* |x| < 2**-13 */ 116 if(huge+x>one) { /* raise inexact if x != 0 */ 117 if(ix<0x3e400000) return one; /* |x|<2**-27 */ |
118 else return one - x*x/4; |
119 } 120 } 121 z = x*x; 122 r = z*(R02+z*(R03+z*(R04+z*R05))); 123 s = one+z*(S01+z*(S02+z*(S03+z*S04))); 124 if(ix < 0x3FF00000) { /* |x| < 1.00 */ 125 return one + z*(-0.25+(r/s)); 126 } else { --- 255 unchanged lines hidden --- |