sysdeps/ieee754/ldbl-96/e_j0l.c

   1 /*
   2  * ====================================================
   3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   4  *
   5  * Developed at SunPro, a Sun Microsystems, Inc. business.
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10  */
  11
  12 /* Long double expansions are
  13   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
  14   and are incorporated herein by permission of the author.  The author
  15   reserves the right to distribute this material elsewhere under different
  16   copying permissions.  These modifications are distributed here under
  17   the following terms:
  18
  19     This library is free software; you can redistribute it and/or
  20     modify it under the terms of the GNU Lesser General Public
  21     License as published by the Free Software Foundation; either
  22     version 2.1 of the License, or (at your option) any later version.
  23
  24     This library is distributed in the hope that it will be useful,
  25     but WITHOUT ANY WARRANTY; without even the implied warranty of
  26     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  27     Lesser General Public License for more details.
  28
  29     You should have received a copy of the GNU Lesser General Public
  30     License along with this library; if not, see
  31     <http://www.gnu.org/licenses/>.  */
  32
  33 /* __ieee754_j0(x), __ieee754_y0(x)
  34  * Bessel function of the first and second kinds of order zero.
  35  * Method -- j0(x):
  36  *      1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
  37  *      2. Reduce x to |x| since j0(x)=j0(-x),  and
  38  *         for x in (0,2)
  39  *              j0(x) = 1 - z/4 + z^2*R0/S0,  where z = x*x;
  40  *         for x in (2,inf)
  41  *              j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
  42  *         where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
  43  *         as follow:
  44  *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
  45  *                      = 1/sqrt(2) * (cos(x) + sin(x))
  46  *              sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
  47  *                      = 1/sqrt(2) * (sin(x) - cos(x))
  48  *         (To avoid cancellation, use
  49  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  50  *          to compute the worse one.)
  51  *
  52  *      3 Special cases
  53  *              j0(nan)= nan
  54  *              j0(0) = 1
  55  *              j0(inf) = 0
  56  *
  57  * Method -- y0(x):
  58  *      1. For x<2.
  59  *         Since
  60  *              y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
  61  *         therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
  62  *         We use the following function to approximate y0,
  63  *              y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
  64  *
  65  *         Note: For tiny x, U/V = u0 and j0(x)~1, hence
  66  *              y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
  67  *      2. For x>=2.
  68  *              y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
  69  *         where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
  70  *         by the method mentioned above.
  71  *      3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
  72  */
  73
  74 #include <math.h>
  75 #include <math_private.h>
  76
  77 static long double pzero (long double), qzero (long double);
  78
  79 static const long double
  80   huge = 1e4930L,
  81   one = 1.0L,
  82   invsqrtpi = 5.6418958354775628694807945156077258584405e-1L,
  83   tpi = 6.3661977236758134307553505349005744813784e-1L,
  84
  85   /* J0(x) = 1 - x^2 / 4 + x^4 R0(x^2) / S0(x^2)
  86      0 <= x <= 2
  87      peak relative error 1.41e-22 */
  88   R[5] = {
  89   4.287176872744686992880841716723478740566E7L,
  90   -6.652058897474241627570911531740907185772E5L,
  91   7.011848381719789863458364584613651091175E3L,
  92   -3.168040850193372408702135490809516253693E1L,
  93   6.030778552661102450545394348845599300939E-2L,
  94 },
  95
  96  S[4] = {
  97    2.743793198556599677955266341699130654342E9L,
  98    3.364330079384816249840086842058954076201E7L,
  99    1.924119649412510777584684927494642526573E5L,
 100    6.239282256012734914211715620088714856494E2L,
 101    /*   1.000000000000000000000000000000000000000E0L,*/
 102 };
 103
 104 static const long double zero = 0.0;
 105
 106 long double
 107 __ieee754_j0l (long double x)
 108 {
 109   long double z, s, c, ss, cc, r, u, v;
 110   int32_t ix;
 111   u_int32_t se;
 112
 113   GET_LDOUBLE_EXP (se, x);
 114   ix = se & 0x7fff;
 115   if (__glibc_unlikely (ix >= 0x7fff))
 116     return one / (x * x);
 117   x = fabsl (x);
 118   if (ix >= 0x4000)             /* |x| >= 2.0 */
 119     {
 120       __sincosl (x, &s, &c);
 121       ss = s - c;
 122       cc = s + c;
 123       if (ix < 0x7ffe)
 124         {                       /* make sure x+x not overflow */
 125           z = -__cosl (x + x);
 126           if ((s * c) < zero)
 127             cc = z / ss;
 128           else
 129             ss = z / cc;
 130         }
 131       /*
 132        * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 133        * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 134        */
 135       if (__glibc_unlikely (ix > 0x4080))       /* 2^129 */
 136         z = (invsqrtpi * cc) / __ieee754_sqrtl (x);
 137       else
 138         {
 139           u = pzero (x);
 140           v = qzero (x);
 141           z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (x);
 142         }
 143       return z;
 144     }
 145   if (__glibc_unlikely (ix < 0x3fef))       /* |x| < 2**-16 */
 146     {
 147       /* raise inexact if x != 0 */
 148       math_force_eval (huge + x);
 149       if (ix < 0x3fde) /* |x| < 2^-33 */
 150         return one;
 151       else
 152         return one - 0.25 * x * x;
 153     }
 154   z = x * x;
 155   r = z * (R[0] + z * (R[1] + z * (R[2] + z * (R[3] + z * R[4]))));
 156   s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z)));
 157   if (ix < 0x3fff)
 158     {                           /* |x| < 1.00 */
 159       return (one - 0.25 * z + z * (r / s));
 160     }
 161   else
 162     {
 163       u = 0.5 * x;
 164       return ((one + u) * (one - u) + z * (r / s));
 165     }
 166 }
 167 strong_alias (__ieee754_j0l, __j0l_finite)
 168
 169
 170 /* y0(x) = 2/pi ln(x) J0(x) + U(x^2)/V(x^2)
 171    0 < x <= 2
 172    peak relative error 1.7e-21 */
 173 static const long double
 174 U[6] = {
 175   -1.054912306975785573710813351985351350861E10L,
 176   2.520192609749295139432773849576523636127E10L,
 177   -1.856426071075602001239955451329519093395E9L,
 178   4.079209129698891442683267466276785956784E7L,
 179   -3.440684087134286610316661166492641011539E5L,
 180   1.005524356159130626192144663414848383774E3L,
 181 };
 182 static const long double
 183 V[5] = {
 184   1.429337283720789610137291929228082613676E11L,
 185   2.492593075325119157558811370165695013002E9L,
 186   2.186077620785925464237324417623665138376E7L,
 187   1.238407896366385175196515057064384929222E5L,
 188   4.693924035211032457494368947123233101664E2L,
 189   /*  1.000000000000000000000000000000000000000E0L */
 190 };
 191
 192 long double
 193 __ieee754_y0l (long double x)
 194 {
 195   long double z, s, c, ss, cc, u, v;
 196   int32_t ix;
 197   u_int32_t se, i0, i1;
 198
 199   GET_LDOUBLE_WORDS (se, i0, i1, x);
 200   ix = se & 0x7fff;
 201   /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
 202   if (__glibc_unlikely (se & 0x8000))
 203     return zero / (zero * x);
 204   if (__glibc_unlikely (ix >= 0x7fff))
 205     return one / (x + x * x);
 206   if (__glibc_unlikely ((i0 | i1) == 0))
 207     return -HUGE_VALL + x;  /* -inf and overflow exception.  */
 208   if (ix >= 0x4000)
 209     {                           /* |x| >= 2.0 */
 210
 211       /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
 212        * where x0 = x-pi/4
 213        *      Better formula:
 214        *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
 215        *                      =  1/sqrt(2) * (sin(x) + cos(x))
 216        *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
 217        *                      =  1/sqrt(2) * (sin(x) - cos(x))
 218        * To avoid cancellation, use
 219        *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 220        * to compute the worse one.
 221        */
 222       __sincosl (x, &s, &c);
 223       ss = s - c;
 224       cc = s + c;
 225       /*
 226        * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 227        * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 228        */
 229       if (ix < 0x7ffe)
 230         {                       /* make sure x+x not overflow */
 231           z = -__cosl (x + x);
 232           if ((s * c) < zero)
 233             cc = z / ss;
 234           else
 235             ss = z / cc;
 236         }
 237       if (__glibc_unlikely (ix > 0x4080))       /* 1e39 */
 238         z = (invsqrtpi * ss) / __ieee754_sqrtl (x);
 239       else
 240         {
 241           u = pzero (x);
 242           v = qzero (x);
 243           z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x);
 244         }
 245       return z;
 246     }
 247   if (__glibc_unlikely (ix <= 0x3fde))       /* x < 2^-33 */
 248     {
 249       z = -7.380429510868722527629822444004602747322E-2L
 250         + tpi * __ieee754_logl (x);
 251       return z;
 252     }
 253   z = x * x;
 254   u = U[0] + z * (U[1] + z * (U[2] + z * (U[3] + z * (U[4] + z * U[5]))));
 255   v = V[0] + z * (V[1] + z * (V[2] + z * (V[3] + z * (V[4] + z))));
 256   return (u / v + tpi * (__ieee754_j0l (x) * __ieee754_logl (x)));
 257 }
 258 strong_alias (__ieee754_y0l, __y0l_finite)
 259
 260 /* The asymptotic expansions of pzero is
 261  *      1 - 9/128 s^2 + 11025/98304 s^4 - ...,  where s = 1/x.
 262  * For x >= 2, We approximate pzero by
 263  *      pzero(x) = 1 + s^2 R(s^2) / S(s^2)
 264  */
 265 static const long double pR8[7] = {
 266   /* 8 <= x <= inf
 267      Peak relative error 4.62 */
 268   -4.094398895124198016684337960227780260127E-9L,
 269   -8.929643669432412640061946338524096893089E-7L,
 270   -6.281267456906136703868258380673108109256E-5L,
 271   -1.736902783620362966354814353559382399665E-3L,
 272   -1.831506216290984960532230842266070146847E-2L,
 273   -5.827178869301452892963280214772398135283E-2L,
 274   -2.087563267939546435460286895807046616992E-2L,
 275 };
 276 static const long double pS8[6] = {
 277   5.823145095287749230197031108839653988393E-8L,
 278   1.279281986035060320477759999428992730280E-5L,
 279   9.132668954726626677174825517150228961304E-4L,
 280   2.606019379433060585351880541545146252534E-2L,
 281   2.956262215119520464228467583516287175244E-1L,
 282   1.149498145388256448535563278632697465675E0L,
 283   /* 1.000000000000000000000000000000000000000E0L, */
 284 };
 285
 286 static const long double pR5[7] = {
 287   /* 4.54541015625 <= x <= 8
 288      Peak relative error 6.51E-22 */
 289   -2.041226787870240954326915847282179737987E-7L,
 290   -2.255373879859413325570636768224534428156E-5L,
 291   -7.957485746440825353553537274569102059990E-4L,
 292   -1.093205102486816696940149222095559439425E-2L,
 293   -5.657957849316537477657603125260701114646E-2L,
 294   -8.641175552716402616180994954177818461588E-2L,
 295   -1.354654710097134007437166939230619726157E-2L,
 296 };
 297 static const long double pS5[6] = {
 298   2.903078099681108697057258628212823545290E-6L,
 299   3.253948449946735405975737677123673867321E-4L,
 300   1.181269751723085006534147920481582279979E-2L,
 301   1.719212057790143888884745200257619469363E-1L,
 302   1.006306498779212467670654535430694221924E0L,
 303   2.069568808688074324555596301126375951502E0L,
 304   /* 1.000000000000000000000000000000000000000E0L, */
 305 };
 306
 307 static const long double pR3[7] = {
 308   /* 2.85711669921875 <= x <= 4.54541015625
 309      peak relative error 5.25e-21 */
 310   -5.755732156848468345557663552240816066802E-6L,
 311   -3.703675625855715998827966962258113034767E-4L,
 312   -7.390893350679637611641350096842846433236E-3L,
 313   -5.571922144490038765024591058478043873253E-2L,
 314   -1.531290690378157869291151002472627396088E-1L,
 315   -1.193350853469302941921647487062620011042E-1L,
 316   -8.567802507331578894302991505331963782905E-3L,
 317 };
 318 static const long double pS3[6] = {
 319   8.185931139070086158103309281525036712419E-5L,
 320   5.398016943778891093520574483111255476787E-3L,
 321   1.130589193590489566669164765853409621081E-1L,
 322   9.358652328786413274673192987670237145071E-1L,
 323   3.091711512598349056276917907005098085273E0L,
 324   3.594602474737921977972586821673124231111E0L,
 325   /* 1.000000000000000000000000000000000000000E0L, */
 326 };
 327
 328 static const long double pR2[7] = {
 329   /* 2 <= x <= 2.85711669921875
 330      peak relative error 2.64e-21 */
 331   -1.219525235804532014243621104365384992623E-4L,
 332   -4.838597135805578919601088680065298763049E-3L,
 333   -5.732223181683569266223306197751407418301E-2L,
 334   -2.472947430526425064982909699406646503758E-1L,
 335   -3.753373645974077960207588073975976327695E-1L,
 336   -1.556241316844728872406672349347137975495E-1L,
 337   -5.355423239526452209595316733635519506958E-3L,
 338 };
 339 static const long double pS2[6] = {
 340   1.734442793664291412489066256138894953823E-3L,
 341   7.158111826468626405416300895617986926008E-2L,
 342   9.153839713992138340197264669867993552641E-1L,
 343   4.539209519433011393525841956702487797582E0L,
 344   8.868932430625331650266067101752626253644E0L,
 345   6.067161890196324146320763844772857713502E0L,
 346   /* 1.000000000000000000000000000000000000000E0L, */
 347 };
 348
 349 static long double
 350 pzero (long double x)
 351 {
 352   const long double *p, *q;
 353   long double z, r, s;
 354   int32_t ix;
 355   u_int32_t se, i0, i1;
 356
 357   GET_LDOUBLE_WORDS (se, i0, i1, x);
 358   ix = se & 0x7fff;
 359   if (ix >= 0x4002)
 360     {
 361       p = pR8;
 362       q = pS8;
 363     }                           /* x >= 8 */
 364   else
 365     {
 366       i1 = (ix << 16) | (i0 >> 16);
 367       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 368         {
 369           p = pR5;
 370           q = pS5;
 371         }
 372       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 373         {
 374           p = pR3;
 375           q = pS3;
 376         }
 377       else if (ix >= 0x4000)    /* x better be >= 2 */
 378         {
 379           p = pR2;
 380           q = pS2;
 381         }
 382     }
 383   z = one / (x * x);
 384   r =
 385     p[0] + z * (p[1] +
 386                 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 387   s =
 388     q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
 389   return (one + z * r / s);
 390 }
 391
 392
 393 /* For x >= 8, the asymptotic expansions of qzero is
 394  *      -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
 395  * We approximate qzero by
 396  *      qzero(x) = s*(-.125 + R(s^2) / S(s^2))
 397  */
 398 static const long double qR8[7] = {
 399   /* 8 <= x <= inf
 400      peak relative error 2.23e-21 */
 401   3.001267180483191397885272640777189348008E-10L,
 402   8.693186311430836495238494289942413810121E-8L,
 403   8.496875536711266039522937037850596580686E-6L,
 404   3.482702869915288984296602449543513958409E-4L,
 405   6.036378380706107692863811938221290851352E-3L,
 406   3.881970028476167836382607922840452192636E-2L,
 407   6.132191514516237371140841765561219149638E-2L,
 408 };
 409 static const long double qS8[7] = {
 410   4.097730123753051126914971174076227600212E-9L,
 411   1.199615869122646109596153392152131139306E-6L,
 412   1.196337580514532207793107149088168946451E-4L,
 413   5.099074440112045094341500497767181211104E-3L,
 414   9.577420799632372483249761659674764460583E-2L,
 415   7.385243015344292267061953461563695918646E-1L,
 416   1.917266424391428937962682301561699055943E0L,
 417   /* 1.000000000000000000000000000000000000000E0L, */
 418 };
 419
 420 static const long double qR5[7] = {
 421   /* 4.54541015625 <= x <= 8
 422      peak relative error 1.03e-21 */
 423   3.406256556438974327309660241748106352137E-8L,
 424   4.855492710552705436943630087976121021980E-6L,
 425   2.301011739663737780613356017352912281980E-4L,
 426   4.500470249273129953870234803596619899226E-3L,
 427   3.651376459725695502726921248173637054828E-2L,
 428   1.071578819056574524416060138514508609805E-1L,
 429   7.458950172851611673015774675225656063757E-2L,
 430 };
 431 static const long double qS5[7] = {
 432   4.650675622764245276538207123618745150785E-7L,
 433   6.773573292521412265840260065635377164455E-5L,
 434   3.340711249876192721980146877577806687714E-3L,
 435   7.036218046856839214741678375536970613501E-2L,
 436   6.569599559163872573895171876511377891143E-1L,
 437   2.557525022583599204591036677199171155186E0L,
 438   3.457237396120935674982927714210361269133E0L,
 439   /* 1.000000000000000000000000000000000000000E0L,*/
 440 };
 441
 442 static const long double qR3[7] = {
 443   /* 2.85711669921875 <= x <= 4.54541015625
 444      peak relative error 5.24e-21 */
 445   1.749459596550816915639829017724249805242E-6L,
 446   1.446252487543383683621692672078376929437E-4L,
 447   3.842084087362410664036704812125005761859E-3L,
 448   4.066369994699462547896426554180954233581E-2L,
 449   1.721093619117980251295234795188992722447E-1L,
 450   2.538595333972857367655146949093055405072E-1L,
 451   8.560591367256769038905328596020118877936E-2L,
 452 };
 453 static const long double qS3[7] = {
 454   2.388596091707517488372313710647510488042E-5L,
 455   2.048679968058758616370095132104333998147E-3L,
 456   5.824663198201417760864458765259945181513E-2L,
 457   6.953906394693328750931617748038994763958E-1L,
 458   3.638186936390881159685868764832961092476E0L,
 459   7.900169524705757837298990558459547842607E0L,
 460   5.992718532451026507552820701127504582907E0L,
 461   /* 1.000000000000000000000000000000000000000E0L, */
 462 };
 463
 464 static const long double qR2[7] = {
 465   /* 2 <= x <= 2.85711669921875
 466      peak relative error 1.58e-21  */
 467   6.306524405520048545426928892276696949540E-5L,
 468   3.209606155709930950935893996591576624054E-3L,
 469   5.027828775702022732912321378866797059604E-2L,
 470   3.012705561838718956481911477587757845163E-1L,
 471   6.960544893905752937420734884995688523815E-1L,
 472   5.431871999743531634887107835372232030655E-1L,
 473   9.447736151202905471899259026430157211949E-2L,
 474 };
 475 static const long double qS2[7] = {
 476   8.610579901936193494609755345106129102676E-4L,
 477   4.649054352710496997203474853066665869047E-2L,
 478   8.104282924459837407218042945106320388339E-1L,
 479   5.807730930825886427048038146088828206852E0L,
 480   1.795310145936848873627710102199881642939E1L,
 481   2.281313316875375733663657188888110605044E1L,
 482   1.011242067883822301487154844458322200143E1L,
 483   /* 1.000000000000000000000000000000000000000E0L, */
 484 };
 485
 486 static long double
 487 qzero (long double x)
 488 {
 489   const long double *p, *q;
 490   long double s, r, z;
 491   int32_t ix;
 492   u_int32_t se, i0, i1;
 493
 494   GET_LDOUBLE_WORDS (se, i0, i1, x);
 495   ix = se & 0x7fff;
 496   if (ix >= 0x4002)             /* x >= 8 */
 497     {
 498       p = qR8;
 499       q = qS8;
 500     }
 501   else
 502     {
 503       i1 = (ix << 16) | (i0 >> 16);
 504       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 505         {
 506           p = qR5;
 507           q = qS5;
 508         }
 509       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 510         {
 511           p = qR3;
 512           q = qS3;
 513         }
 514       else if (ix >= 0x4000)    /* x better be >= 2 */
 515         {
 516           p = qR2;
 517           q = qS2;
 518         }
 519     }
 520   z = one / (x * x);
 521   r =
 522     p[0] + z * (p[1] +
 523                 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 524   s =
 525     q[0] + z * (q[1] +
 526                 z * (q[2] +
 527                      z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
 528   return (-.125 + z * r / s) / x;
 529 }