1 /* $NetBSD: n_j1.c,v 1.6 2003/08/07 16:44:51 agc 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. Neither the name of the University nor the names of its contributors 15 * may be used to endorse or promote products derived from this software 16 * without specific prior written permission. 17 * 18 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 19 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 20 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 21 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 22 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 23 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 24 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 25 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 26 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 27 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 28 * SUCH DAMAGE. 29 */ 30 31 #ifndef lint 32 #if 0 33 static char sccsid[] = "@(#)j1.c 8.2 (Berkeley) 11/30/93"; 34 #endif 35 #endif /* not lint */ 36 37 /* 38 * 16 December 1992 39 * Minor modifications by Peter McIlroy to adapt non-IEEE architecture. 40 */ 41 42 /* 43 * ==================================================== 44 * Copyright (C) 1992 by Sun Microsystems, Inc. 45 * 46 * Developed at SunPro, a Sun Microsystems, Inc. business. 47 * Permission to use, copy, modify, and distribute this 48 * software is freely granted, provided that this notice 49 * is preserved. 50 * ==================================================== 51 * 52 * ******************* WARNING ******************** 53 * This is an alpha version of SunPro's FDLIBM (Freely 54 * Distributable Math Library) for IEEE double precision 55 * arithmetic. FDLIBM is a basic math library written 56 * in C that runs on machines that conform to IEEE 57 * Standard 754/854. This alpha version is distributed 58 * for testing purpose. Those who use this software 59 * should report any bugs to 60 * 61 * fdlibm-comments@sunpro.eng.sun.com 62 * 63 * -- K.C. Ng, Oct 12, 1992 64 * ************************************************ 65 */ 66 67 /* double j1(double x), y1(double x) 68 * Bessel function of the first and second kinds of order zero. 69 * Method -- j1(x): 70 * 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ... 71 * 2. Reduce x to |x| since j1(x)=-j1(-x), and 72 * for x in (0,2) 73 * j1(x) = x/2 + x*z*R0/S0, where z = x*x; 74 * (precision: |j1/x - 1/2 - R0/S0 |<2**-61.51 ) 75 * for x in (2,inf) 76 * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1)) 77 * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) 78 * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) 79 * as follows: 80 * cos(x1) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) 81 * = 1/sqrt(2) * (sin(x) - cos(x)) 82 * sin(x1) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) 83 * = -1/sqrt(2) * (sin(x) + cos(x)) 84 * (To avoid cancellation, use 85 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) 86 * to compute the worse one.) 87 * 88 * 3 Special cases 89 * j1(nan)= nan 90 * j1(0) = 0 91 * j1(inf) = 0 92 * 93 * Method -- y1(x): 94 * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN 95 * 2. For x<2. 96 * Since 97 * y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...) 98 * therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function. 99 * We use the following function to approximate y1, 100 * y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2 101 * where for x in [0,2] (abs err less than 2**-65.89) 102 * U(z) = u0 + u1*z + ... + u4*z^4 103 * V(z) = 1 + v1*z + ... + v5*z^5 104 * Note: For tiny x, 1/x dominate y1 and hence 105 * y1(tiny) = -2/pi/tiny, (choose tiny<2**-54) 106 * 3. For x>=2. 107 * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) 108 * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) 109 * by method mentioned above. 110 */ 111 112 #include "mathimpl.h" 113 #include <float.h> 114 #include <errno.h> 115 116 #if defined(__vax__) || defined(tahoe) 117 #define _IEEE 0 118 #else 119 #define _IEEE 1 120 #define infnan(x) (0.0) 121 #endif 122 123 static double pone (double), qone (double); 124 125 static const double 126 huge = 1e300, 127 zero = 0.0, 128 one = 1.0, 129 invsqrtpi= 5.641895835477562869480794515607725858441e-0001, 130 tpi = 0.636619772367581343075535053490057448, 131 132 /* R0/S0 on [0,2] */ 133 r00 = -6.250000000000000020842322918309200910191e-0002, 134 r01 = 1.407056669551897148204830386691427791200e-0003, 135 r02 = -1.599556310840356073980727783817809847071e-0005, 136 r03 = 4.967279996095844750387702652791615403527e-0008, 137 s01 = 1.915375995383634614394860200531091839635e-0002, 138 s02 = 1.859467855886309024045655476348872850396e-0004, 139 s03 = 1.177184640426236767593432585906758230822e-0006, 140 s04 = 5.046362570762170559046714468225101016915e-0009, 141 s05 = 1.235422744261379203512624973117299248281e-0011; 142 143 #define two_129 6.80564733841876926e+038 /* 2^129 */ 144 #define two_m54 5.55111512312578270e-017 /* 2^-54 */ 145 146 double 147 j1(double x) 148 { 149 double z, s,c,ss,cc,r,u,v,y; 150 y = fabs(x); 151 if (!finite(x)) { /* Inf or NaN */ 152 if (_IEEE && x != x) 153 return(x); 154 else 155 return (copysign(x, zero)); 156 } 157 y = fabs(x); 158 if (y >= 2) { /* |x| >= 2.0 */ 159 s = sin(y); 160 c = cos(y); 161 ss = -s-c; 162 cc = s-c; 163 if (y < .5*DBL_MAX) { /* make sure y+y not overflow */ 164 z = cos(y+y); 165 if ((s*c)<zero) cc = z/ss; 166 else ss = z/cc; 167 } 168 /* 169 * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) 170 * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) 171 */ 172 #if !defined(__vax__) && !defined(tahoe) 173 if (y > two_129) /* x > 2^129 */ 174 z = (invsqrtpi*cc)/sqrt(y); 175 else 176 #endif /* defined(__vax__) || defined(tahoe) */ 177 { 178 u = pone(y); v = qone(y); 179 z = invsqrtpi*(u*cc-v*ss)/sqrt(y); 180 } 181 if (x < 0) return -z; 182 else return z; 183 } 184 if (y < 7.450580596923828125e-009) { /* |x|<2**-27 */ 185 if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */ 186 } 187 z = x*x; 188 r = z*(r00+z*(r01+z*(r02+z*r03))); 189 s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); 190 r *= x; 191 return (x*0.5+r/s); 192 } 193 194 static const double u0[5] = { 195 -1.960570906462389484206891092512047539632e-0001, 196 5.044387166398112572026169863174882070274e-0002, 197 -1.912568958757635383926261729464141209569e-0003, 198 2.352526005616105109577368905595045204577e-0005, 199 -9.190991580398788465315411784276789663849e-0008, 200 }; 201 static const double v0[5] = { 202 1.991673182366499064031901734535479833387e-0002, 203 2.025525810251351806268483867032781294682e-0004, 204 1.356088010975162198085369545564475416398e-0006, 205 6.227414523646214811803898435084697863445e-0009, 206 1.665592462079920695971450872592458916421e-0011, 207 }; 208 209 double 210 y1(double x) 211 { 212 double z, s, c, ss, cc, u, v; 213 /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ 214 if (!finite(x)) { 215 if (!_IEEE) return (infnan(EDOM)); 216 else if (x < 0) 217 return(zero/zero); 218 else if (x > 0) 219 return (0); 220 else 221 return(x); 222 } 223 if (x <= 0) { 224 if (_IEEE && x == 0) return -one/zero; 225 else if(x == 0) return(infnan(-ERANGE)); 226 else if(_IEEE) return (zero/zero); 227 else return(infnan(EDOM)); 228 } 229 if (x >= 2) { /* |x| >= 2.0 */ 230 s = sin(x); 231 c = cos(x); 232 ss = -s-c; 233 cc = s-c; 234 if (x < .5 * DBL_MAX) { /* make sure x+x not overflow */ 235 z = cos(x+x); 236 if ((s*c)>zero) cc = z/ss; 237 else ss = z/cc; 238 } 239 /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0)) 240 * where x0 = x-3pi/4 241 * Better formula: 242 * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) 243 * = 1/sqrt(2) * (sin(x) - cos(x)) 244 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) 245 * = -1/sqrt(2) * (cos(x) + sin(x)) 246 * To avoid cancellation, use 247 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) 248 * to compute the worse one. 249 */ 250 if (_IEEE && x>two_129) { 251 z = (invsqrtpi*ss)/sqrt(x); 252 } else { 253 u = pone(x); v = qone(x); 254 z = invsqrtpi*(u*ss+v*cc)/sqrt(x); 255 } 256 return z; 257 } 258 if (x <= two_m54) { /* x < 2**-54 */ 259 return (-tpi/x); 260 } 261 z = x*x; 262 u = u0[0]+z*(u0[1]+z*(u0[2]+z*(u0[3]+z*u0[4]))); 263 v = one+z*(v0[0]+z*(v0[1]+z*(v0[2]+z*(v0[3]+z*v0[4])))); 264 return (x*(u/v) + tpi*(j1(x)*log(x)-one/x)); 265 } 266 267 /* For x >= 8, the asymptotic expansions of pone is 268 * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x. 269 * We approximate pone by 270 * pone(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 * | pone(x)-1-R/S | <= 2 ** ( -60.06) 275 */ 276 277 static const double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 278 0.0, 279 1.171874999999886486643746274751925399540e-0001, 280 1.323948065930735690925827997575471527252e+0001, 281 4.120518543073785433325860184116512799375e+0002, 282 3.874745389139605254931106878336700275601e+0003, 283 7.914479540318917214253998253147871806507e+0003, 284 }; 285 static const double ps8[5] = { 286 1.142073703756784104235066368252692471887e+0002, 287 3.650930834208534511135396060708677099382e+0003, 288 3.695620602690334708579444954937638371808e+0004, 289 9.760279359349508334916300080109196824151e+0004, 290 3.080427206278887984185421142572315054499e+0004, 291 }; 292 293 static const double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 294 1.319905195562435287967533851581013807103e-0011, 295 1.171874931906140985709584817065144884218e-0001, 296 6.802751278684328781830052995333841452280e+0000, 297 1.083081829901891089952869437126160568246e+0002, 298 5.176361395331997166796512844100442096318e+0002, 299 5.287152013633375676874794230748055786553e+0002, 300 }; 301 static const double ps5[5] = { 302 5.928059872211313557747989128353699746120e+0001, 303 9.914014187336144114070148769222018425781e+0002, 304 5.353266952914879348427003712029704477451e+0003, 305 7.844690317495512717451367787640014588422e+0003, 306 1.504046888103610723953792002716816255382e+0003, 307 }; 308 309 static const double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 310 3.025039161373736032825049903408701962756e-0009, 311 1.171868655672535980750284752227495879921e-0001, 312 3.932977500333156527232725812363183251138e+0000, 313 3.511940355916369600741054592597098912682e+0001, 314 9.105501107507812029367749771053045219094e+0001, 315 4.855906851973649494139275085628195457113e+0001, 316 }; 317 static const double ps3[5] = { 318 3.479130950012515114598605916318694946754e+0001, 319 3.367624587478257581844639171605788622549e+0002, 320 1.046871399757751279180649307467612538415e+0003, 321 8.908113463982564638443204408234739237639e+0002, 322 1.037879324396392739952487012284401031859e+0002, 323 }; 324 325 static const double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 326 1.077108301068737449490056513753865482831e-0007, 327 1.171762194626833490512746348050035171545e-0001, 328 2.368514966676087902251125130227221462134e+0000, 329 1.224261091482612280835153832574115951447e+0001, 330 1.769397112716877301904532320376586509782e+0001, 331 5.073523125888185399030700509321145995160e+0000, 332 }; 333 static const double ps2[5] = { 334 2.143648593638214170243114358933327983793e+0001, 335 1.252902271684027493309211410842525120355e+0002, 336 2.322764690571628159027850677565128301361e+0002, 337 1.176793732871470939654351793502076106651e+0002, 338 8.364638933716182492500902115164881195742e+0000, 339 }; 340 341 static double 342 pone(double x) 343 { 344 const double *p,*q; 345 double z,r,s; 346 if (x >= 8.0) {p = pr8; q= ps8;} 347 else if (x >= 4.54545211791992188) {p = pr5; q= ps5;} 348 else if (x >= 2.85714149475097656) {p = pr3; q= ps3;} 349 else /* if (x >= 2.0) */ {p = pr2; q= ps2;} 350 z = one/(x*x); 351 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 352 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); 353 return (one + r/s); 354 } 355 356 357 /* For x >= 8, the asymptotic expansions of qone is 358 * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. 359 * We approximate pone by 360 * qone(x) = s*(0.375 + (R/S)) 361 * where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10 362 * S = 1 + qs1*s^2 + ... + qs6*s^12 363 * and 364 * | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13) 365 */ 366 367 static const double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 368 0.0, 369 -1.025390624999927207385863635575804210817e-0001, 370 -1.627175345445899724355852152103771510209e+0001, 371 -7.596017225139501519843072766973047217159e+0002, 372 -1.184980667024295901645301570813228628541e+0004, 373 -4.843851242857503225866761992518949647041e+0004, 374 }; 375 static const double qs8[6] = { 376 1.613953697007229231029079421446916397904e+0002, 377 7.825385999233484705298782500926834217525e+0003, 378 1.338753362872495800748094112937868089032e+0005, 379 7.196577236832409151461363171617204036929e+0005, 380 6.666012326177764020898162762642290294625e+0005, 381 -2.944902643038346618211973470809456636830e+0005, 382 }; 383 384 static const double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 385 -2.089799311417640889742251585097264715678e-0011, 386 -1.025390502413754195402736294609692303708e-0001, 387 -8.056448281239359746193011295417408828404e+0000, 388 -1.836696074748883785606784430098756513222e+0002, 389 -1.373193760655081612991329358017247355921e+0003, 390 -2.612444404532156676659706427295870995743e+0003, 391 }; 392 static const double qs5[6] = { 393 8.127655013843357670881559763225310973118e+0001, 394 1.991798734604859732508048816860471197220e+0003, 395 1.746848519249089131627491835267411777366e+0004, 396 4.985142709103522808438758919150738000353e+0004, 397 2.794807516389181249227113445299675335543e+0004, 398 -4.719183547951285076111596613593553911065e+0003, 399 }; 400 401 static const double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 402 -5.078312264617665927595954813341838734288e-0009, 403 -1.025378298208370901410560259001035577681e-0001, 404 -4.610115811394734131557983832055607679242e+0000, 405 -5.784722165627836421815348508816936196402e+0001, 406 -2.282445407376317023842545937526967035712e+0002, 407 -2.192101284789093123936441805496580237676e+0002, 408 }; 409 static const double qs3[6] = { 410 4.766515503237295155392317984171640809318e+0001, 411 6.738651126766996691330687210949984203167e+0002, 412 3.380152866795263466426219644231687474174e+0003, 413 5.547729097207227642358288160210745890345e+0003, 414 1.903119193388108072238947732674639066045e+0003, 415 -1.352011914443073322978097159157678748982e+0002, 416 }; 417 418 static const double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 419 -1.783817275109588656126772316921194887979e-0007, 420 -1.025170426079855506812435356168903694433e-0001, 421 -2.752205682781874520495702498875020485552e+0000, 422 -1.966361626437037351076756351268110418862e+0001, 423 -4.232531333728305108194363846333841480336e+0001, 424 -2.137192117037040574661406572497288723430e+0001, 425 }; 426 static const double qs2[6] = { 427 2.953336290605238495019307530224241335502e+0001, 428 2.529815499821905343698811319455305266409e+0002, 429 7.575028348686454070022561120722815892346e+0002, 430 7.393932053204672479746835719678434981599e+0002, 431 1.559490033366661142496448853793707126179e+0002, 432 -4.959498988226281813825263003231704397158e+0000, 433 }; 434 435 static double 436 qone(double x) 437 { 438 const double *p,*q; 439 double s,r,z; 440 if (x >= 8.0) {p = qr8; q= qs8;} 441 else if (x >= 4.54545211791992188) {p = qr5; q= qs5;} 442 else if (x >= 2.85714149475097656) {p = qr3; q= qs3;} 443 else /* if (x >= 2.0) */ {p = qr2; q= qs2;} 444 z = one/(x*x); 445 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 446 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); 447 return (.375 + r/s)/x; 448 } 449