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   /* ix >= 0x4000 for all calls to this function.  */
 360   if (ix >= 0x4002)
 361     {
 362       p = pR8;
 363       q = pS8;
 364     }                           /* x >= 8 */
 365   else
 366     {
 367       i1 = (ix << 16) | (i0 >> 16);
 368       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 369         {
 370           p = pR5;
 371           q = pS5;
 372         }
 373       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 374         {
 375           p = pR3;
 376           q = pS3;
 377         }
 378       else      /* x >= 2 */
 379         {
 380           p = pR2;
 381           q = pS2;
 382         }
 383     }
 384   z = one / (x * x);
 385   r =
 386     p[0] + z * (p[1] +
 387                 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 388   s =
 389     q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
 390   return (one + z * r / s);
 391 }
 392
 393
 394 /* For x >= 8, the asymptotic expansions of qzero is
 395  *      -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
 396  * We approximate qzero by
 397  *      qzero(x) = s*(-.125 + R(s^2) / S(s^2))
 398  */
 399 static const long double qR8[7] = {
 400   /* 8 <= x <= inf
 401      peak relative error 2.23e-21 */
 402   3.001267180483191397885272640777189348008E-10L,
 403   8.693186311430836495238494289942413810121E-8L,
 404   8.496875536711266039522937037850596580686E-6L,
 405   3.482702869915288984296602449543513958409E-4L,
 406   6.036378380706107692863811938221290851352E-3L,
 407   3.881970028476167836382607922840452192636E-2L,
 408   6.132191514516237371140841765561219149638E-2L,
 409 };
 410 static const long double qS8[7] = {
 411   4.097730123753051126914971174076227600212E-9L,
 412   1.199615869122646109596153392152131139306E-6L,
 413   1.196337580514532207793107149088168946451E-4L,
 414   5.099074440112045094341500497767181211104E-3L,
 415   9.577420799632372483249761659674764460583E-2L,
 416   7.385243015344292267061953461563695918646E-1L,
 417   1.917266424391428937962682301561699055943E0L,
 418   /* 1.000000000000000000000000000000000000000E0L, */
 419 };
 420
 421 static const long double qR5[7] = {
 422   /* 4.54541015625 <= x <= 8
 423      peak relative error 1.03e-21 */
 424   3.406256556438974327309660241748106352137E-8L,
 425   4.855492710552705436943630087976121021980E-6L,
 426   2.301011739663737780613356017352912281980E-4L,
 427   4.500470249273129953870234803596619899226E-3L,
 428   3.651376459725695502726921248173637054828E-2L,
 429   1.071578819056574524416060138514508609805E-1L,
 430   7.458950172851611673015774675225656063757E-2L,
 431 };
 432 static const long double qS5[7] = {
 433   4.650675622764245276538207123618745150785E-7L,
 434   6.773573292521412265840260065635377164455E-5L,
 435   3.340711249876192721980146877577806687714E-3L,
 436   7.036218046856839214741678375536970613501E-2L,
 437   6.569599559163872573895171876511377891143E-1L,
 438   2.557525022583599204591036677199171155186E0L,
 439   3.457237396120935674982927714210361269133E0L,
 440   /* 1.000000000000000000000000000000000000000E0L,*/
 441 };
 442
 443 static const long double qR3[7] = {
 444   /* 2.85711669921875 <= x <= 4.54541015625
 445      peak relative error 5.24e-21 */
 446   1.749459596550816915639829017724249805242E-6L,
 447   1.446252487543383683621692672078376929437E-4L,
 448   3.842084087362410664036704812125005761859E-3L,
 449   4.066369994699462547896426554180954233581E-2L,
 450   1.721093619117980251295234795188992722447E-1L,
 451   2.538595333972857367655146949093055405072E-1L,
 452   8.560591367256769038905328596020118877936E-2L,
 453 };
 454 static const long double qS3[7] = {
 455   2.388596091707517488372313710647510488042E-5L,
 456   2.048679968058758616370095132104333998147E-3L,
 457   5.824663198201417760864458765259945181513E-2L,
 458   6.953906394693328750931617748038994763958E-1L,
 459   3.638186936390881159685868764832961092476E0L,
 460   7.900169524705757837298990558459547842607E0L,
 461   5.992718532451026507552820701127504582907E0L,
 462   /* 1.000000000000000000000000000000000000000E0L, */
 463 };
 464
 465 static const long double qR2[7] = {
 466   /* 2 <= x <= 2.85711669921875
 467      peak relative error 1.58e-21  */
 468   6.306524405520048545426928892276696949540E-5L,
 469   3.209606155709930950935893996591576624054E-3L,
 470   5.027828775702022732912321378866797059604E-2L,
 471   3.012705561838718956481911477587757845163E-1L,
 472   6.960544893905752937420734884995688523815E-1L,
 473   5.431871999743531634887107835372232030655E-1L,
 474   9.447736151202905471899259026430157211949E-2L,
 475 };
 476 static const long double qS2[7] = {
 477   8.610579901936193494609755345106129102676E-4L,
 478   4.649054352710496997203474853066665869047E-2L,
 479   8.104282924459837407218042945106320388339E-1L,
 480   5.807730930825886427048038146088828206852E0L,
 481   1.795310145936848873627710102199881642939E1L,
 482   2.281313316875375733663657188888110605044E1L,
 483   1.011242067883822301487154844458322200143E1L,
 484   /* 1.000000000000000000000000000000000000000E0L, */
 485 };
 486
 487 static long double
 488 qzero (long double x)
 489 {
 490   const long double *p, *q;
 491   long double s, r, z;
 492   int32_t ix;
 493   u_int32_t se, i0, i1;
 494
 495   GET_LDOUBLE_WORDS (se, i0, i1, x);
 496   ix = se & 0x7fff;
 497   /* ix >= 0x4000 for all calls to this function.  */
 498   if (ix >= 0x4002)             /* x >= 8 */
 499     {
 500       p = qR8;
 501       q = qS8;
 502     }
 503   else
 504     {
 505       i1 = (ix << 16) | (i0 >> 16);
 506       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 507         {
 508           p = qR5;
 509           q = qS5;
 510         }
 511       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 512         {
 513           p = qR3;
 514           q = qS3;
 515         }
 516       else      /* x >= 2 */
 517         {
 518           p = qR2;
 519           q = qS2;
 520         }
 521     }
 522   z = one / (x * x);
 523   r =
 524     p[0] + z * (p[1] +
 525                 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 526   s =
 527     q[0] + z * (q[1] +
 528                 z * (q[2] +
 529                      z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
 530   return (-.125 + z * r / s) / x;
 531 }