n_j0.c revision 1.2
1 /* $NetBSD: n_j0.c,v 1.2 1997/10/20 14:12:44 ragge Exp $ */ 2 /*- 3 * Copyright (c) 1992, 1993 4 * The Regents of the University of California. All rights reserved. 5 * 6 * Redistribution and use in source and binary forms, with or without 7 * modification, are permitted provided that the following conditions 8 * are met: 9 * 1. Redistributions of source code must retain the above copyright 10 * notice, this list of conditions and the following disclaimer. 11 * 2. Redistributions in binary form must reproduce the above copyright 12 * notice, this list of conditions and the following disclaimer in the 13 * documentation and/or other materials provided with the distribution. 14 * 3. All advertising materials mentioning features or use of this software 15 * must display the following acknowledgement: 16 * This product includes software developed by the University of 17 * California, Berkeley and its contributors. 18 * 4. Neither the name of the University nor the names of its contributors 19 * may be used to endorse or promote products derived from this software 20 * without specific prior written permission. 21 * 22 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 23 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 24 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 25 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 26 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 27 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 28 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 29 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 30 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 31 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 32 * SUCH DAMAGE. 33 */ 34 35 #ifndef lint 36 #if 0 37 static char sccsid[] = "@(#)j0.c 8.2 (Berkeley) 11/30/93"; 38 #endif 39 #endif /* not lint */ 40 41 /* 42 * 16 December 1992 43 * Minor modifications by Peter McIlroy to adapt non-IEEE architecture. 44 */ 45 46 /* 47 * ==================================================== 48 * Copyright (C) 1992 by Sun Microsystems, Inc. 49 * 50 * Developed at SunPro, a Sun Microsystems, Inc. business. 51 * Permission to use, copy, modify, and distribute this 52 * software is freely granted, provided that this notice 53 * is preserved. 54 * ==================================================== 55 * 56 * ******************* WARNING ******************** 57 * This is an alpha version of SunPro's FDLIBM (Freely 58 * Distributable Math Library) for IEEE double precision 59 * arithmetic. FDLIBM is a basic math library written 60 * in C that runs on machines that conform to IEEE 61 * Standard 754/854. This alpha version is distributed 62 * for testing purpose. Those who use this software 63 * should report any bugs to 64 * 65 * fdlibm-comments (at) sunpro.eng.sun.com 66 * 67 * -- K.C. Ng, Oct 12, 1992 68 * ************************************************ 69 */ 70 71 /* double j0(double x), y0(double x) 72 * Bessel function of the first and second kinds of order zero. 73 * Method -- j0(x): 74 * 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ... 75 * 2. Reduce x to |x| since j0(x)=j0(-x), and 76 * for x in (0,2) 77 * j0(x) = 1-z/4+ z^2*R0/S0, where z = x*x; 78 * (precision: |j0-1+z/4-z^2R0/S0 |<2**-63.67 ) 79 * for x in (2,inf) 80 * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) 81 * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) 82 * as follow: 83 * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) 84 * = 1/sqrt(2) * (cos(x) + sin(x)) 85 * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4) 86 * = 1/sqrt(2) * (sin(x) - cos(x)) 87 * (To avoid cancellation, use 88 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) 89 * to compute the worse one.) 90 * 91 * 3 Special cases 92 * j0(nan)= nan 93 * j0(0) = 1 94 * j0(inf) = 0 95 * 96 * Method -- y0(x): 97 * 1. For x<2. 98 * Since 99 * y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...) 100 * therefore y0(x)-2/pi*j0(x)*ln(x) is an even function. 101 * We use the following function to approximate y0, 102 * y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2 103 * where 104 * U(z) = u0 + u1*z + ... + u6*z^6 105 * V(z) = 1 + v1*z + ... + v4*z^4 106 * with absolute approximation error bounded by 2**-72. 107 * Note: For tiny x, U/V = u0 and j0(x)~1, hence 108 * y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27) 109 * 2. For x>=2. 110 * y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0)) 111 * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) 112 * by the method mentioned above. 113 * 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0. 114 */ 115 116 #include "mathimpl.h" 117 #include <float.h> 118 #include <errno.h> 119 120 #if defined(vax) || defined(tahoe) 121 #define _IEEE 0 122 #else 123 #define _IEEE 1 124 #define infnan(x) (0.0) 125 #endif 126 127 static double pzero __P((double)), qzero __P((double)); 128 129 static double 130 huge = 1e300, 131 zero = 0.0, 132 one = 1.0, 133 invsqrtpi= 5.641895835477562869480794515607725858441e-0001, 134 tpi = 0.636619772367581343075535053490057448, 135 /* R0/S0 on [0, 2.00] */ 136 r02 = 1.562499999999999408594634421055018003102e-0002, 137 r03 = -1.899792942388547334476601771991800712355e-0004, 138 r04 = 1.829540495327006565964161150603950916854e-0006, 139 r05 = -4.618326885321032060803075217804816988758e-0009, 140 s01 = 1.561910294648900170180789369288114642057e-0002, 141 s02 = 1.169267846633374484918570613449245536323e-0004, 142 s03 = 5.135465502073181376284426245689510134134e-0007, 143 s04 = 1.166140033337900097836930825478674320464e-0009; 144 145 double 146 j0(x) 147 double x; 148 { 149 double z, s,c,ss,cc,r,u,v; 150 151 if (!finite(x)) 152 if (_IEEE) return one/(x*x); 153 else return (0); 154 x = fabs(x); 155 if (x >= 2.0) { /* |x| >= 2.0 */ 156 s = sin(x); 157 c = cos(x); 158 ss = s-c; 159 cc = s+c; 160 if (x < .5 * DBL_MAX) { /* make sure x+x not overflow */ 161 z = -cos(x+x); 162 if ((s*c)<zero) cc = z/ss; 163 else ss = z/cc; 164 } 165 /* 166 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) 167 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) 168 */ 169 if (_IEEE && x> 6.80564733841876927e+38) /* 2^129 */ 170 z = (invsqrtpi*cc)/sqrt(x); 171 else { 172 u = pzero(x); v = qzero(x); 173 z = invsqrtpi*(u*cc-v*ss)/sqrt(x); 174 } 175 return z; 176 } 177 if (x < 1.220703125e-004) { /* |x| < 2**-13 */ 178 if (huge+x > one) { /* raise inexact if x != 0 */ 179 if (x < 7.450580596923828125e-009) /* |x|<2**-27 */ 180 return one; 181 else return (one - 0.25*x*x); 182 } 183 } 184 z = x*x; 185 r = z*(r02+z*(r03+z*(r04+z*r05))); 186 s = one+z*(s01+z*(s02+z*(s03+z*s04))); 187 if (x < one) { /* |x| < 1.00 */ 188 return (one + z*(-0.25+(r/s))); 189 } else { 190 u = 0.5*x; 191 return ((one+u)*(one-u)+z*(r/s)); 192 } 193 } 194 195 static double 196 u00 = -7.380429510868722527422411862872999615628e-0002, 197 u01 = 1.766664525091811069896442906220827182707e-0001, 198 u02 = -1.381856719455968955440002438182885835344e-0002, 199 u03 = 3.474534320936836562092566861515617053954e-0004, 200 u04 = -3.814070537243641752631729276103284491172e-0006, 201 u05 = 1.955901370350229170025509706510038090009e-0008, 202 u06 = -3.982051941321034108350630097330144576337e-0011, 203 v01 = 1.273048348341237002944554656529224780561e-0002, 204 v02 = 7.600686273503532807462101309675806839635e-0005, 205 v03 = 2.591508518404578033173189144579208685163e-0007, 206 v04 = 4.411103113326754838596529339004302243157e-0010; 207 208 double 209 y0(x) 210 double x; 211 { 212 double z, s, c, ss, cc, u, v; 213 /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ 214 if (!finite(x)) 215 if (_IEEE) 216 return (one/(x+x*x)); 217 else 218 return (0); 219 if (x == 0) 220 if (_IEEE) return (-one/zero); 221 else return(infnan(-ERANGE)); 222 if (x<0) 223 if (_IEEE) return (zero/zero); 224 else return (infnan(EDOM)); 225 if (x >= 2.00) { /* |x| >= 2.0 */ 226 /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) 227 * where x0 = x-pi/4 228 * Better formula: 229 * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) 230 * = 1/sqrt(2) * (sin(x) + cos(x)) 231 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) 232 * = 1/sqrt(2) * (sin(x) - cos(x)) 233 * To avoid cancellation, use 234 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) 235 * to compute the worse one. 236 */ 237 s = sin(x); 238 c = cos(x); 239 ss = s-c; 240 cc = s+c; 241 /* 242 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) 243 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) 244 */ 245 if (x < .5 * DBL_MAX) { /* make sure x+x not overflow */ 246 z = -cos(x+x); 247 if ((s*c)<zero) cc = z/ss; 248 else ss = z/cc; 249 } 250 if (_IEEE && x > 6.80564733841876927e+38) /* > 2^129 */ 251 z = (invsqrtpi*ss)/sqrt(x); 252 else { 253 u = pzero(x); v = qzero(x); 254 z = invsqrtpi*(u*ss+v*cc)/sqrt(x); 255 } 256 return z; 257 } 258 if (x <= 7.450580596923828125e-009) { /* x < 2**-27 */ 259 return (u00 + tpi*log(x)); 260 } 261 z = x*x; 262 u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); 263 v = one+z*(v01+z*(v02+z*(v03+z*v04))); 264 return (u/v + tpi*(j0(x)*log(x))); 265 } 266 267 /* The asymptotic expansions of pzero is 268 * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. 269 * For x >= 2, We approximate pzero by 270 * pzero(x) = 1 + (R/S) 271 * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 272 * S = 1 + ps0*s^2 + ... + ps4*s^10 273 * and 274 * | pzero(x)-1-R/S | <= 2 ** ( -60.26) 275 */ 276 static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 277 0.0, 278 -7.031249999999003994151563066182798210142e-0002, 279 -8.081670412753498508883963849859423939871e+0000, 280 -2.570631056797048755890526455854482662510e+0002, 281 -2.485216410094288379417154382189125598962e+0003, 282 -5.253043804907295692946647153614119665649e+0003, 283 }; 284 static double ps8[5] = { 285 1.165343646196681758075176077627332052048e+0002, 286 3.833744753641218451213253490882686307027e+0003, 287 4.059785726484725470626341023967186966531e+0004, 288 1.167529725643759169416844015694440325519e+0005, 289 4.762772841467309430100106254805711722972e+0004, 290 }; 291 292 static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 293 -1.141254646918944974922813501362824060117e-0011, 294 -7.031249408735992804117367183001996028304e-0002, 295 -4.159610644705877925119684455252125760478e+0000, 296 -6.767476522651671942610538094335912346253e+0001, 297 -3.312312996491729755731871867397057689078e+0002, 298 -3.464333883656048910814187305901796723256e+0002, 299 }; 300 static double ps5[5] = { 301 6.075393826923003305967637195319271932944e+0001, 302 1.051252305957045869801410979087427910437e+0003, 303 5.978970943338558182743915287887408780344e+0003, 304 9.625445143577745335793221135208591603029e+0003, 305 2.406058159229391070820491174867406875471e+0003, 306 }; 307 308 static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 309 -2.547046017719519317420607587742992297519e-0009, 310 -7.031196163814817199050629727406231152464e-0002, 311 -2.409032215495295917537157371488126555072e+0000, 312 -2.196597747348830936268718293366935843223e+0001, 313 -5.807917047017375458527187341817239891940e+0001, 314 -3.144794705948885090518775074177485744176e+0001, 315 }; 316 static double ps3[5] = { 317 3.585603380552097167919946472266854507059e+0001, 318 3.615139830503038919981567245265266294189e+0002, 319 1.193607837921115243628631691509851364715e+0003, 320 1.127996798569074250675414186814529958010e+0003, 321 1.735809308133357510239737333055228118910e+0002, 322 }; 323 324 static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 325 -8.875343330325263874525704514800809730145e-0008, 326 -7.030309954836247756556445443331044338352e-0002, 327 -1.450738467809529910662233622603401167409e+0000, 328 -7.635696138235277739186371273434739292491e+0000, 329 -1.119316688603567398846655082201614524650e+0001, 330 -3.233645793513353260006821113608134669030e+0000, 331 }; 332 static double ps2[5] = { 333 2.222029975320888079364901247548798910952e+0001, 334 1.362067942182152109590340823043813120940e+0002, 335 2.704702786580835044524562897256790293238e+0002, 336 1.538753942083203315263554770476850028583e+0002, 337 1.465761769482561965099880599279699314477e+0001, 338 }; 339 340 static double pzero(x) 341 double x; 342 { 343 double *p,*q,z,r,s; 344 if (x >= 8.00) {p = pr8; q= ps8;} 345 else if (x >= 4.54545211791992188) {p = pr5; q= ps5;} 346 else if (x >= 2.85714149475097656) {p = pr3; q= ps3;} 347 else /* if (x >= 2.00) */ {p = pr2; q= ps2;} 348 z = one/(x*x); 349 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 350 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); 351 return one+ r/s; 352 } 353 354 355 /* For x >= 8, the asymptotic expansions of qzero is 356 * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. 357 * We approximate pzero by 358 * qzero(x) = s*(-1.25 + (R/S)) 359 * where R = qr0 + qr1*s^2 + qr2*s^4 + ... + qr5*s^10 360 * S = 1 + qs0*s^2 + ... + qs5*s^12 361 * and 362 * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) 363 */ 364 static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 365 0.0, 366 7.324218749999350414479738504551775297096e-0002, 367 1.176820646822526933903301695932765232456e+0001, 368 5.576733802564018422407734683549251364365e+0002, 369 8.859197207564685717547076568608235802317e+0003, 370 3.701462677768878501173055581933725704809e+0004, 371 }; 372 static double qs8[6] = { 373 1.637760268956898345680262381842235272369e+0002, 374 8.098344946564498460163123708054674227492e+0003, 375 1.425382914191204905277585267143216379136e+0005, 376 8.033092571195144136565231198526081387047e+0005, 377 8.405015798190605130722042369969184811488e+0005, 378 -3.438992935378666373204500729736454421006e+0005, 379 }; 380 381 static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 382 1.840859635945155400568380711372759921179e-0011, 383 7.324217666126847411304688081129741939255e-0002, 384 5.835635089620569401157245917610984757296e+0000, 385 1.351115772864498375785526599119895942361e+0002, 386 1.027243765961641042977177679021711341529e+0003, 387 1.989977858646053872589042328678602481924e+0003, 388 }; 389 static double qs5[6] = { 390 8.277661022365377058749454444343415524509e+0001, 391 2.077814164213929827140178285401017305309e+0003, 392 1.884728877857180787101956800212453218179e+0004, 393 5.675111228949473657576693406600265778689e+0004, 394 3.597675384251145011342454247417399490174e+0004, 395 -5.354342756019447546671440667961399442388e+0003, 396 }; 397 398 static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 399 4.377410140897386263955149197672576223054e-0009, 400 7.324111800429115152536250525131924283018e-0002, 401 3.344231375161707158666412987337679317358e+0000, 402 4.262184407454126175974453269277100206290e+0001, 403 1.708080913405656078640701512007621675724e+0002, 404 1.667339486966511691019925923456050558293e+0002, 405 }; 406 static double qs3[6] = { 407 4.875887297245871932865584382810260676713e+0001, 408 7.096892210566060535416958362640184894280e+0002, 409 3.704148226201113687434290319905207398682e+0003, 410 6.460425167525689088321109036469797462086e+0003, 411 2.516333689203689683999196167394889715078e+0003, 412 -1.492474518361563818275130131510339371048e+0002, 413 }; 414 415 static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 416 1.504444448869832780257436041633206366087e-0007, 417 7.322342659630792930894554535717104926902e-0002, 418 1.998191740938159956838594407540292600331e+0000, 419 1.449560293478857407645853071687125850962e+0001, 420 3.166623175047815297062638132537957315395e+0001, 421 1.625270757109292688799540258329430963726e+0001, 422 }; 423 static double qs2[6] = { 424 3.036558483552191922522729838478169383969e+0001, 425 2.693481186080498724211751445725708524507e+0002, 426 8.447837575953201460013136756723746023736e+0002, 427 8.829358451124885811233995083187666981299e+0002, 428 2.126663885117988324180482985363624996652e+0002, 429 -5.310954938826669402431816125780738924463e+0000, 430 }; 431 432 static double qzero(x) 433 double x; 434 { 435 double *p,*q, s,r,z; 436 if (x >= 8.00) {p = qr8; q= qs8;} 437 else if (x >= 4.54545211791992188) {p = qr5; q= qs5;} 438 else if (x >= 2.85714149475097656) {p = qr3; q= qs3;} 439 else /* if (x >= 2.00) */ {p = qr2; q= qs2;} 440 z = one/(x*x); 441 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 442 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); 443 return (-.125 + r/s)/x; 444 } 445
Indexes created Sun Oct 12 09:09:55 GMT 2025