lib/libm/noieee_src/n_j0.c

   1 /*      $NetBSD: n_j0.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[] = "@(#)j0.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 j0(double x), y0(double x)
  68  * Bessel function of the first and second kinds of order zero.
  69  * Method -- j0(x):
  70  *      1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
  71  *      2. Reduce x to |x| since j0(x)=j0(-x),  and
  72  *         for x in (0,2)
  73  *              j0(x) = 1-z/4+ z^2*R0/S0,  where z = x*x;
  74  *         (precision:  |j0-1+z/4-z^2R0/S0 |<2**-63.67 )
  75  *         for x in (2,inf)
  76  *              j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
  77  *         where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
  78  *         as follow:
  79  *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
  80  *                      = 1/sqrt(2) * (cos(x) + sin(x))
  81  *              sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
  82  *                      = 1/sqrt(2) * (sin(x) - cos(x))
  83  *         (To avoid cancellation, use
  84  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  85  *          to compute the worse one.)
  86  *
  87  *      3 Special cases
  88  *              j0(nan)= nan
  89  *              j0(0) = 1
  90  *              j0(inf) = 0
  91  *
  92  * Method -- y0(x):
  93  *      1. For x<2.
  94  *         Since
  95  *              y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
  96  *         therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
  97  *         We use the following function to approximate y0,
  98  *              y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
  99  *         where
 100  *              U(z) = u0 + u1*z + ... + u6*z^6
 101  *              V(z) = 1  + v1*z + ... + v4*z^4
 102  *         with absolute approximation error bounded by 2**-72.
 103  *         Note: For tiny x, U/V = u0 and j0(x)~1, hence
 104  *              y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
 105  *      2. For x>=2.
 106  *              y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
 107  *         where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
 108  *         by the method mentioned above.
 109  *      3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
 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 pzero (double), qzero (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                 /* R0/S0 on [0, 2.00] */
 132 r02 =   1.562499999999999408594634421055018003102e-0002,
 133 r03 =  -1.899792942388547334476601771991800712355e-0004,
 134 r04 =   1.829540495327006565964161150603950916854e-0006,
 135 r05 =  -4.618326885321032060803075217804816988758e-0009,
 136 s01 =   1.561910294648900170180789369288114642057e-0002,
 137 s02 =   1.169267846633374484918570613449245536323e-0004,
 138 s03 =   5.135465502073181376284426245689510134134e-0007,
 139 s04 =   1.166140033337900097836930825478674320464e-0009;
 140
 141 double
 142 j0(double x)
 143 {
 144         double z, s,c,ss,cc,r,u,v;
 145
 146         if (!finite(x)) {
 147 #if _IEEE
 148                 return one/(x*x);
 149 #else
 150                 return (0);
 151 #endif
 152         }
 153         x = fabs(x);
 154         if (x >= 2.0) { /* |x| >= 2.0 */
 155                 s = sin(x);
 156                 c = cos(x);
 157                 ss = s-c;
 158                 cc = s+c;
 159                 if (x < .5 * DBL_MAX) {  /* make sure x+x not overflow */
 160                     z = -cos(x+x);
 161                     if ((s*c)<zero) cc = z/ss;
 162                     else            ss = z/cc;
 163                 }
 164         /*
 165          * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 166          * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 167          */
 168 #if _IEEE
 169                 if (x > 6.80564733841876927e+38) /* 2^129 */
 170                         z = (invsqrtpi*cc)/sqrt(x);
 171                 else
 172 #endif
 173                 {
 174                     u = pzero(x); v = qzero(x);
 175                     z = invsqrtpi*(u*cc-v*ss)/sqrt(x);
 176                 }
 177                 return z;
 178         }
 179         if (x < 1.220703125e-004) {                /* |x| < 2**-13 */
 180             if (huge+x > one) {                    /* raise inexact if x != 0 */
 181                 if (x < 7.450580596923828125e-009) /* |x|<2**-27 */
 182                         return one;
 183                 else return (one - 0.25*x*x);
 184             }
 185         }
 186         z = x*x;
 187         r =  z*(r02+z*(r03+z*(r04+z*r05)));
 188         s =  one+z*(s01+z*(s02+z*(s03+z*s04)));
 189         if (x < one) {                  /* |x| < 1.00 */
 190             return (one + z*(-0.25+(r/s)));
 191         } else {
 192             u = 0.5*x;
 193             return ((one+u)*(one-u)+z*(r/s));
 194         }
 195 }
 196
 197 static const double
 198 u00 =  -7.380429510868722527422411862872999615628e-0002,
 199 u01 =   1.766664525091811069896442906220827182707e-0001,
 200 u02 =  -1.381856719455968955440002438182885835344e-0002,
 201 u03 =   3.474534320936836562092566861515617053954e-0004,
 202 u04 =  -3.814070537243641752631729276103284491172e-0006,
 203 u05 =   1.955901370350229170025509706510038090009e-0008,
 204 u06 =  -3.982051941321034108350630097330144576337e-0011,
 205 v01 =   1.273048348341237002944554656529224780561e-0002,
 206 v02 =   7.600686273503532807462101309675806839635e-0005,
 207 v03 =   2.591508518404578033173189144579208685163e-0007,
 208 v04 =   4.411103113326754838596529339004302243157e-0010;
 209
 210 double
 211 y0(double x)
 212 {
 213         double z, s, c, ss, cc, u, v;
 214     /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
 215         if (!finite(x)) {
 216 #if _IEEE
 217                         return (one/(x+x*x));
 218 #else
 219                         return (0);
 220 #endif
 221         }
 222         if (x == 0) {
 223 #if _IEEE
 224                 return (-one/zero);
 225 #else
 226                 return(infnan(-ERANGE));
 227 #endif
 228         }
 229         if (x<0) {
 230 #if _IEEE
 231                 return (zero/zero);
 232 #else
 233                 return (infnan(EDOM));
 234 #endif
 235         }
 236         if (x >= 2.00) {        /* |x| >= 2.0 */
 237         /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
 238          * where x0 = x-pi/4
 239          *      Better formula:
 240          *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
 241          *                      =  1/sqrt(2) * (sin(x) + cos(x))
 242          *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
 243          *                      =  1/sqrt(2) * (sin(x) - cos(x))
 244          * To avoid cancellation, use
 245          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 246          * to compute the worse one.
 247          */
 248                 s = sin(x);
 249                 c = cos(x);
 250                 ss = s-c;
 251                 cc = s+c;
 252         /*
 253          * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 254          * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 255          */
 256                 if (x < .5 * DBL_MAX) {  /* make sure x+x not overflow */
 257                     z = -cos(x+x);
 258                     if ((s*c)<zero) cc = z/ss;
 259                     else            ss = z/cc;
 260                 }
 261 #if _IEEE
 262                 if (x > 6.80564733841876927e+38) /* > 2^129 */
 263                         z = (invsqrtpi*ss)/sqrt(x);
 264                 else
 265 #endif
 266                 {
 267                     u = pzero(x); v = qzero(x);
 268                     z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
 269                 }
 270                 return z;
 271         }
 272         if (x <= 7.450580596923828125e-009) {           /* x < 2**-27 */
 273             return (u00 + tpi*log(x));
 274         }
 275         z = x*x;
 276         u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
 277         v = one+z*(v01+z*(v02+z*(v03+z*v04)));
 278         return (u/v + tpi*(j0(x)*log(x)));
 279 }
 280
 281 /* The asymptotic expansions of pzero is
 282  *      1 - 9/128 s^2 + 11025/98304 s^4 - ...,  where s = 1/x.
 283  * For x >= 2, We approximate pzero by
 284  *      pzero(x) = 1 + (R/S)
 285  * where  R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
 286  *        S = 1 + ps0*s^2 + ... + ps4*s^10
 287  * and
 288  *      | pzero(x)-1-R/S | <= 2  ** ( -60.26)
 289  */
 290 static const double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
 291    0.0,
 292   -7.031249999999003994151563066182798210142e-0002,
 293   -8.081670412753498508883963849859423939871e+0000,
 294   -2.570631056797048755890526455854482662510e+0002,
 295   -2.485216410094288379417154382189125598962e+0003,
 296   -5.253043804907295692946647153614119665649e+0003,
 297 };
 298 static const double ps8[5] = {
 299    1.165343646196681758075176077627332052048e+0002,
 300    3.833744753641218451213253490882686307027e+0003,
 301    4.059785726484725470626341023967186966531e+0004,
 302    1.167529725643759169416844015694440325519e+0005,
 303    4.762772841467309430100106254805711722972e+0004,
 304 };
 305
 306 static const double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
 307   -1.141254646918944974922813501362824060117e-0011,
 308   -7.031249408735992804117367183001996028304e-0002,
 309   -4.159610644705877925119684455252125760478e+0000,
 310   -6.767476522651671942610538094335912346253e+0001,
 311   -3.312312996491729755731871867397057689078e+0002,
 312   -3.464333883656048910814187305901796723256e+0002,
 313 };
 314 static const double ps5[5] = {
 315    6.075393826923003305967637195319271932944e+0001,
 316    1.051252305957045869801410979087427910437e+0003,
 317    5.978970943338558182743915287887408780344e+0003,
 318    9.625445143577745335793221135208591603029e+0003,
 319    2.406058159229391070820491174867406875471e+0003,
 320 };
 321
 322 static const double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
 323   -2.547046017719519317420607587742992297519e-0009,
 324   -7.031196163814817199050629727406231152464e-0002,
 325   -2.409032215495295917537157371488126555072e+0000,
 326   -2.196597747348830936268718293366935843223e+0001,
 327   -5.807917047017375458527187341817239891940e+0001,
 328   -3.144794705948885090518775074177485744176e+0001,
 329 };
 330 static const double ps3[5] = {
 331    3.585603380552097167919946472266854507059e+0001,
 332    3.615139830503038919981567245265266294189e+0002,
 333    1.193607837921115243628631691509851364715e+0003,
 334    1.127996798569074250675414186814529958010e+0003,
 335    1.735809308133357510239737333055228118910e+0002,
 336 };
 337
 338 static const double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
 339   -8.875343330325263874525704514800809730145e-0008,
 340   -7.030309954836247756556445443331044338352e-0002,
 341   -1.450738467809529910662233622603401167409e+0000,
 342   -7.635696138235277739186371273434739292491e+0000,
 343   -1.119316688603567398846655082201614524650e+0001,
 344   -3.233645793513353260006821113608134669030e+0000,
 345 };
 346 static const double ps2[5] = {
 347    2.222029975320888079364901247548798910952e+0001,
 348    1.362067942182152109590340823043813120940e+0002,
 349    2.704702786580835044524562897256790293238e+0002,
 350    1.538753942083203315263554770476850028583e+0002,
 351    1.465761769482561965099880599279699314477e+0001,
 352 };
 353
 354 static double
 355 pzero(double x)
 356 {
 357         const double *p,*q;
 358         double z,r,s;
 359         if (x >= 8.00)                     {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.00) */          {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 qzero is
 371  *      -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
 372  * We approximate pzero by
 373  *      qzero(x) = s*(-1.25 + (R/S))
 374  * where  R = qr0 + qr1*s^2 + qr2*s^4 + ... + qr5*s^10
 375  *        S = 1 + qs0*s^2 + ... + qs5*s^12
 376  * and
 377  *      | qzero(x)/s +1.25-R/S | <= 2  ** ( -61.22)
 378  */
 379 static const double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
 380    0.0,
 381    7.324218749999350414479738504551775297096e-0002,
 382    1.176820646822526933903301695932765232456e+0001,
 383    5.576733802564018422407734683549251364365e+0002,
 384    8.859197207564685717547076568608235802317e+0003,
 385    3.701462677768878501173055581933725704809e+0004,
 386 };
 387 static const double qs8[6] = {
 388    1.637760268956898345680262381842235272369e+0002,
 389    8.098344946564498460163123708054674227492e+0003,
 390    1.425382914191204905277585267143216379136e+0005,
 391    8.033092571195144136565231198526081387047e+0005,
 392    8.405015798190605130722042369969184811488e+0005,
 393   -3.438992935378666373204500729736454421006e+0005,
 394 };
 395
 396 static const double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
 397    1.840859635945155400568380711372759921179e-0011,
 398    7.324217666126847411304688081129741939255e-0002,
 399    5.835635089620569401157245917610984757296e+0000,
 400    1.351115772864498375785526599119895942361e+0002,
 401    1.027243765961641042977177679021711341529e+0003,
 402    1.989977858646053872589042328678602481924e+0003,
 403 };
 404 static const double qs5[6] = {
 405    8.277661022365377058749454444343415524509e+0001,
 406    2.077814164213929827140178285401017305309e+0003,
 407    1.884728877857180787101956800212453218179e+0004,
 408    5.675111228949473657576693406600265778689e+0004,
 409    3.597675384251145011342454247417399490174e+0004,
 410   -5.354342756019447546671440667961399442388e+0003,
 411 };
 412
 413 static const double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
 414    4.377410140897386263955149197672576223054e-0009,
 415    7.324111800429115152536250525131924283018e-0002,
 416    3.344231375161707158666412987337679317358e+0000,
 417    4.262184407454126175974453269277100206290e+0001,
 418    1.708080913405656078640701512007621675724e+0002,
 419    1.667339486966511691019925923456050558293e+0002,
 420 };
 421 static const double qs3[6] = {
 422    4.875887297245871932865584382810260676713e+0001,
 423    7.096892210566060535416958362640184894280e+0002,
 424    3.704148226201113687434290319905207398682e+0003,
 425    6.460425167525689088321109036469797462086e+0003,
 426    2.516333689203689683999196167394889715078e+0003,
 427   -1.492474518361563818275130131510339371048e+0002,
 428 };
 429
 430 static const double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
 431    1.504444448869832780257436041633206366087e-0007,
 432    7.322342659630792930894554535717104926902e-0002,
 433    1.998191740938159956838594407540292600331e+0000,
 434    1.449560293478857407645853071687125850962e+0001,
 435    3.166623175047815297062638132537957315395e+0001,
 436    1.625270757109292688799540258329430963726e+0001,
 437 };
 438 static const double qs2[6] = {
 439    3.036558483552191922522729838478169383969e+0001,
 440    2.693481186080498724211751445725708524507e+0002,
 441    8.447837575953201460013136756723746023736e+0002,
 442    8.829358451124885811233995083187666981299e+0002,
 443    2.126663885117988324180482985363624996652e+0002,
 444   -5.310954938826669402431816125780738924463e+0000,
 445 };
 446
 447 static double
 448 qzero(double x)
 449 {
 450         const double *p,*q;
 451         double s,r,z;
 452         if (x >= 8.00)                     {p = qr8; q= qs8;}
 453         else if (x >= 4.54545211791992188) {p = qr5; q= qs5;}
 454         else if (x >= 2.85714149475097656) {p = qr3; q= qs3;}
 455         else /* if (x >= 2.00) */          {p = qr2; q= qs2;}
 456         z = one/(x*x);
 457         r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
 458         s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
 459         return (-.125 + r/s)/x;
 460 }