lib/libm/noieee_src/n_j1.c

   1 /*      $NetBSD: n_j1.c,v 1.7 2011/11/02 02:34:56 christos 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    = _HUGE,
 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
 153                 if (x != x)
 154                         return(x);
 155                 else
 156 #endif
 157                         return (copysign(x, zero));
 158         }
 159         y = fabs(x);
 160         if (y >= 2) {                   /* |x| >= 2.0 */
 161                 s = sin(y);
 162                 c = cos(y);
 163                 ss = -s-c;
 164                 cc = s-c;
 165                 if (y < .5*DBL_MAX) {   /* make sure y+y not overflow */
 166                     z = cos(y+y);
 167                     if ((s*c)<zero) cc = z/ss;
 168                     else            ss = z/cc;
 169                 }
 170         /*
 171          * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
 172          * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
 173          */
 174 #if !defined(__vax__) && !defined(tahoe)
 175                 if (y > two_129)         /* x > 2^129 */
 176                         z = (invsqrtpi*cc)/sqrt(y);
 177                 else
 178 #endif /* defined(__vax__) || defined(tahoe) */
 179                 {
 180                     u = pone(y); v = qone(y);
 181                     z = invsqrtpi*(u*cc-v*ss)/sqrt(y);
 182                 }
 183                 if (x < 0) return -z;
 184                 else     return  z;
 185         }
 186         if (y < 7.450580596923828125e-009) {    /* |x|<2**-27 */
 187             if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */
 188         }
 189         z = x*x;
 190         r =  z*(r00+z*(r01+z*(r02+z*r03)));
 191         s =  one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
 192         r *= x;
 193         return (x*0.5+r/s);
 194 }
 195
 196 static const double u0[5] = {
 197   -1.960570906462389484206891092512047539632e-0001,
 198    5.044387166398112572026169863174882070274e-0002,
 199   -1.912568958757635383926261729464141209569e-0003,
 200    2.352526005616105109577368905595045204577e-0005,
 201    -9.190991580398788465315411784276789663849e-0008,
 202 };
 203 static const double v0[5] = {
 204    1.991673182366499064031901734535479833387e-0002,
 205    2.025525810251351806268483867032781294682e-0004,
 206    1.356088010975162198085369545564475416398e-0006,
 207    6.227414523646214811803898435084697863445e-0009,
 208    1.665592462079920695971450872592458916421e-0011,
 209 };
 210
 211 double
 212 y1(double x)
 213 {
 214         double z, s, c, ss, cc, u, v;
 215     /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
 216         if (!finite(x)) {
 217 #if _IEEE
 218                 if (x < 0)
 219                         return(zero/zero);
 220                 else if (x > 0)
 221                         return (0);
 222                 else
 223                         return(x);
 224 #else
 225                 return (infnan(EDOM));
 226 #endif
 227         }
 228         if (x <= 0) {
 229 #if _IEEE
 230                 if (x == 0) return -one/zero;
 231 #endif
 232                 if(x == 0) return(infnan(-ERANGE));
 233 #if _IEEE
 234                 return (zero/zero);
 235 #else
 236                 return(infnan(EDOM));
 237 #endif
 238         }
 239         if (x >= 2) {                    /* |x| >= 2.0 */
 240                 s = sin(x);
 241                 c = cos(x);
 242                 ss = -s-c;
 243                 cc = s-c;
 244                 if (x < .5 * DBL_MAX) { /* make sure x+x not overflow */
 245                     z = cos(x+x);
 246                     if ((s*c)>zero) cc = z/ss;
 247                     else            ss = z/cc;
 248                 }
 249         /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
 250          * where x0 = x-3pi/4
 251          *      Better formula:
 252          *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
 253          *                      =  1/sqrt(2) * (sin(x) - cos(x))
 254          *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
 255          *                      = -1/sqrt(2) * (cos(x) + sin(x))
 256          * To avoid cancellation, use
 257          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 258          * to compute the worse one.
 259          */
 260 #if _IEEE
 261                 if (x>two_129) {
 262                         z = (invsqrtpi*ss)/sqrt(x);
 263                 } else
 264 #endif
 265                 {
 266                     u = pone(x); v = qone(x);
 267                     z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
 268                 }
 269                 return z;
 270         }
 271         if (x <= two_m54) {    /* x < 2**-54 */
 272             return (-tpi/x);
 273         }
 274         z = x*x;
 275         u = u0[0]+z*(u0[1]+z*(u0[2]+z*(u0[3]+z*u0[4])));
 276         v = one+z*(v0[0]+z*(v0[1]+z*(v0[2]+z*(v0[3]+z*v0[4]))));
 277         return (x*(u/v) + tpi*(j1(x)*log(x)-one/x));
 278 }
 279
 280 /* For x >= 8, the asymptotic expansions of pone is
 281  *      1 + 15/128 s^2 - 4725/2^15 s^4 - ...,   where s = 1/x.
 282  * We approximate pone by
 283  *      pone(x) = 1 + (R/S)
 284  * where  R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
 285  *        S = 1 + ps0*s^2 + ... + ps4*s^10
 286  * and
 287  *      | pone(x)-1-R/S | <= 2  ** ( -60.06)
 288  */
 289
 290 static const double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
 291    0.0,
 292    1.171874999999886486643746274751925399540e-0001,
 293    1.323948065930735690925827997575471527252e+0001,
 294    4.120518543073785433325860184116512799375e+0002,
 295    3.874745389139605254931106878336700275601e+0003,
 296    7.914479540318917214253998253147871806507e+0003,
 297 };
 298 static const double ps8[5] = {
 299    1.142073703756784104235066368252692471887e+0002,
 300    3.650930834208534511135396060708677099382e+0003,
 301    3.695620602690334708579444954937638371808e+0004,
 302    9.760279359349508334916300080109196824151e+0004,
 303    3.080427206278887984185421142572315054499e+0004,
 304 };
 305
 306 static const double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
 307    1.319905195562435287967533851581013807103e-0011,
 308    1.171874931906140985709584817065144884218e-0001,
 309    6.802751278684328781830052995333841452280e+0000,
 310    1.083081829901891089952869437126160568246e+0002,
 311    5.176361395331997166796512844100442096318e+0002,
 312    5.287152013633375676874794230748055786553e+0002,
 313 };
 314 static const double ps5[5] = {
 315    5.928059872211313557747989128353699746120e+0001,
 316    9.914014187336144114070148769222018425781e+0002,
 317    5.353266952914879348427003712029704477451e+0003,
 318    7.844690317495512717451367787640014588422e+0003,
 319    1.504046888103610723953792002716816255382e+0003,
 320 };
 321
 322 static const double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
 323    3.025039161373736032825049903408701962756e-0009,
 324    1.171868655672535980750284752227495879921e-0001,
 325    3.932977500333156527232725812363183251138e+0000,
 326    3.511940355916369600741054592597098912682e+0001,
 327    9.105501107507812029367749771053045219094e+0001,
 328    4.855906851973649494139275085628195457113e+0001,
 329 };
 330 static const double ps3[5] = {
 331    3.479130950012515114598605916318694946754e+0001,
 332    3.367624587478257581844639171605788622549e+0002,
 333    1.046871399757751279180649307467612538415e+0003,
 334    8.908113463982564638443204408234739237639e+0002,
 335    1.037879324396392739952487012284401031859e+0002,
 336 };
 337
 338 static const double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
 339    1.077108301068737449490056513753865482831e-0007,
 340    1.171762194626833490512746348050035171545e-0001,
 341    2.368514966676087902251125130227221462134e+0000,
 342    1.224261091482612280835153832574115951447e+0001,
 343    1.769397112716877301904532320376586509782e+0001,
 344    5.073523125888185399030700509321145995160e+0000,
 345 };
 346 static const double ps2[5] = {
 347    2.143648593638214170243114358933327983793e+0001,
 348    1.252902271684027493309211410842525120355e+0002,
 349    2.322764690571628159027850677565128301361e+0002,
 350    1.176793732871470939654351793502076106651e+0002,
 351    8.364638933716182492500902115164881195742e+0000,
 352 };
 353
 354 static double
 355 pone(double x)
 356 {
 357         const double *p,*q;
 358         double z,r,s;
 359         if (x >= 8.0)                      {p = pr8; q= ps8;}
 360         else if (x >= 4.54545211791992188) {p = pr5; q= ps5;}
 361         else if (x >= 2.85714149475097656) {p = pr3; q= ps3;}
 362         else /* if (x >= 2.0) */           {p = pr2; q= ps2;}
 363         z = one/(x*x);
 364         r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
 365         s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
 366         return (one + r/s);
 367 }
 368
 369
 370 /* For x >= 8, the asymptotic expansions of qone is
 371  *      3/8 s - 105/1024 s^3 - ..., where s = 1/x.
 372  * We approximate pone by
 373  *      qone(x) = s*(0.375 + (R/S))
 374  * where  R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
 375  *        S = 1 + qs1*s^2 + ... + qs6*s^12
 376  * and
 377  *      | qone(x)/s -0.375-R/S | <= 2  ** ( -61.13)
 378  */
 379
 380 static const double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
 381    0.0,
 382   -1.025390624999927207385863635575804210817e-0001,
 383   -1.627175345445899724355852152103771510209e+0001,
 384   -7.596017225139501519843072766973047217159e+0002,
 385   -1.184980667024295901645301570813228628541e+0004,
 386   -4.843851242857503225866761992518949647041e+0004,
 387 };
 388 static const double qs8[6] = {
 389    1.613953697007229231029079421446916397904e+0002,
 390    7.825385999233484705298782500926834217525e+0003,
 391    1.338753362872495800748094112937868089032e+0005,
 392    7.196577236832409151461363171617204036929e+0005,
 393    6.666012326177764020898162762642290294625e+0005,
 394   -2.944902643038346618211973470809456636830e+0005,
 395 };
 396
 397 static const double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
 398   -2.089799311417640889742251585097264715678e-0011,
 399   -1.025390502413754195402736294609692303708e-0001,
 400   -8.056448281239359746193011295417408828404e+0000,
 401   -1.836696074748883785606784430098756513222e+0002,
 402   -1.373193760655081612991329358017247355921e+0003,
 403   -2.612444404532156676659706427295870995743e+0003,
 404 };
 405 static const double qs5[6] = {
 406    8.127655013843357670881559763225310973118e+0001,
 407    1.991798734604859732508048816860471197220e+0003,
 408    1.746848519249089131627491835267411777366e+0004,
 409    4.985142709103522808438758919150738000353e+0004,
 410    2.794807516389181249227113445299675335543e+0004,
 411   -4.719183547951285076111596613593553911065e+0003,
 412 };
 413
 414 static const double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
 415   -5.078312264617665927595954813341838734288e-0009,
 416   -1.025378298208370901410560259001035577681e-0001,
 417   -4.610115811394734131557983832055607679242e+0000,
 418   -5.784722165627836421815348508816936196402e+0001,
 419   -2.282445407376317023842545937526967035712e+0002,
 420   -2.192101284789093123936441805496580237676e+0002,
 421 };
 422 static const double qs3[6] = {
 423    4.766515503237295155392317984171640809318e+0001,
 424    6.738651126766996691330687210949984203167e+0002,
 425    3.380152866795263466426219644231687474174e+0003,
 426    5.547729097207227642358288160210745890345e+0003,
 427    1.903119193388108072238947732674639066045e+0003,
 428   -1.352011914443073322978097159157678748982e+0002,
 429 };
 430
 431 static const double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
 432   -1.783817275109588656126772316921194887979e-0007,
 433   -1.025170426079855506812435356168903694433e-0001,
 434   -2.752205682781874520495702498875020485552e+0000,
 435   -1.966361626437037351076756351268110418862e+0001,
 436   -4.232531333728305108194363846333841480336e+0001,
 437   -2.137192117037040574661406572497288723430e+0001,
 438 };
 439 static const double qs2[6] = {
 440    2.953336290605238495019307530224241335502e+0001,
 441    2.529815499821905343698811319455305266409e+0002,
 442    7.575028348686454070022561120722815892346e+0002,
 443    7.393932053204672479746835719678434981599e+0002,
 444    1.559490033366661142496448853793707126179e+0002,
 445   -4.959498988226281813825263003231704397158e+0000,
 446 };
 447
 448 static double
 449 qone(double x)
 450 {
 451         const double *p,*q;
 452         double s,r,z;
 453         if (x >= 8.0)                      {p = qr8; q= qs8;}
 454         else if (x >= 4.54545211791992188) {p = qr5; q= qs5;}
 455         else if (x >= 2.85714149475097656) {p = qr3; q= qs3;}
 456         else /* if (x >= 2.0) */           {p = qr2; q= qs2;}
 457         z = one/(x*x);
 458         r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
 459         s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
 460         return (.375 + r/s)/x;
 461 }