sysdeps/ieee754/ldbl-96/e_j1l.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_j1(x), __ieee754_y1(x)
  34  * Bessel function of the first and second kinds of order zero.
  35  * Method -- j1(x):
  36  *      1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
  37  *      2. Reduce x to |x| since j1(x)=-j1(-x),  and
  38  *         for x in (0,2)
  39  *              j1(x) = x/2 + x*z*R0/S0,  where z = x*x;
  40  *         for x in (2,inf)
  41  *              j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
  42  *              y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
  43  *         where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
  44  *         as follow:
  45  *              cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
  46  *                      =  1/sqrt(2) * (sin(x) - cos(x))
  47  *              sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
  48  *                      = -1/sqrt(2) * (sin(x) + cos(x))
  49  *         (To avoid cancellation, use
  50  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  51  *          to compute the worse one.)
  52  *
  53  *      3 Special cases
  54  *              j1(nan)= nan
  55  *              j1(0) = 0
  56  *              j1(inf) = 0
  57  *
  58  * Method -- y1(x):
  59  *      1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
  60  *      2. For x<2.
  61  *         Since
  62  *              y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
  63  *         therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
  64  *         We use the following function to approximate y1,
  65  *              y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
  66  *         Note: For tiny x, 1/x dominate y1 and hence
  67  *              y1(tiny) = -2/pi/tiny
  68  *      3. For x>=2.
  69  *              y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
  70  *         where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
  71  *         by method mentioned above.
  72  */
  73
  74 #include "math.h"
  75 #include "math_private.h"
  76
  77 static long double pone (long double), qone (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   /* J1(x) = .5 x + x x^2 R(x^2) / S(x^2)
  86      0 <= x <= 2
  87      Peak relative error 4.5e-21 */
  88 R[5] = {
  89     -9.647406112428107954753770469290757756814E7L,
  90     2.686288565865230690166454005558203955564E6L,
  91     -3.689682683905671185891885948692283776081E4L,
  92     2.195031194229176602851429567792676658146E2L,
  93     -5.124499848728030297902028238597308971319E-1L,
  94 },
  95
  96   S[4] =
  97 {
  98   1.543584977988497274437410333029029035089E9L,
  99   2.133542369567701244002565983150952549520E7L,
 100   1.394077011298227346483732156167414670520E5L,
 101   5.252401789085732428842871556112108446506E2L,
 102   /*  1.000000000000000000000000000000000000000E0L, */
 103 };
 104
 105 static const long double zero = 0.0;
 106
 107
 108 long double
 109 __ieee754_j1l (long double x)
 110 {
 111   long double z, c, r, s, ss, cc, u, v, y;
 112   int32_t ix;
 113   u_int32_t se;
 114
 115   GET_LDOUBLE_EXP (se, x);
 116   ix = se & 0x7fff;
 117   if (__builtin_expect (ix >= 0x7fff, 0))
 118     return one / x;
 119   y = fabsl (x);
 120   if (ix >= 0x4000)
 121     {                           /* |x| >= 2.0 */
 122       __sincosl (y, &s, &c);
 123       ss = -s - c;
 124       cc = s - c;
 125       if (ix < 0x7ffe)
 126         {                       /* make sure y+y not overflow */
 127           z = __cosl (y + y);
 128           if ((s * c) > zero)
 129             cc = z / ss;
 130           else
 131             ss = z / cc;
 132         }
 133       /*
 134        * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
 135        * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
 136        */
 137       if (__builtin_expect (ix > 0x4080, 0))
 138         z = (invsqrtpi * cc) / __ieee754_sqrtl (y);
 139       else
 140         {
 141           u = pone (y);
 142           v = qone (y);
 143           z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (y);
 144         }
 145       if (se & 0x8000)
 146         return -z;
 147       else
 148         return z;
 149     }
 150   if (__builtin_expect (ix < 0x3fde, 0)) /* |x| < 2^-33 */
 151     {
 152       if (huge + x > one)
 153         return 0.5 * x;         /* inexact if x!=0 necessary */
 154     }
 155   z = x * x;
 156   r = z * (R[0] + z * (R[1]+ z * (R[2] + z * (R[3] + z * R[4]))));
 157   s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z)));
 158   r *= x;
 159   return (x * 0.5 + r / s);
 160 }
 161 strong_alias (__ieee754_j1l, __j1l_finite)
 162
 163
 164 /* Y1(x) = 2/pi * (log(x) * j1(x) - 1/x) + x R(x^2)
 165    0 <= x <= 2
 166    Peak relative error 2.3e-23 */
 167 static const long double U0[6] = {
 168   -5.908077186259914699178903164682444848615E10L,
 169   1.546219327181478013495975514375773435962E10L,
 170   -6.438303331169223128870035584107053228235E8L,
 171   9.708540045657182600665968063824819371216E6L,
 172   -6.138043997084355564619377183564196265471E4L,
 173   1.418503228220927321096904291501161800215E2L,
 174 };
 175 static const long double V0[5] = {
 176   3.013447341682896694781964795373783679861E11L,
 177   4.669546565705981649470005402243136124523E9L,
 178   3.595056091631351184676890179233695857260E7L,
 179   1.761554028569108722903944659933744317994E5L,
 180   5.668480419646516568875555062047234534863E2L,
 181   /*  1.000000000000000000000000000000000000000E0L, */
 182 };
 183
 184
 185 long double
 186 __ieee754_y1l (long double x)
 187 {
 188   long double z, s, c, ss, cc, u, v;
 189   int32_t ix;
 190   u_int32_t se, i0, i1;
 191
 192   GET_LDOUBLE_WORDS (se, i0, i1, x);
 193   ix = se & 0x7fff;
 194   /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
 195   if (__builtin_expect (se & 0x8000, 0))
 196     return zero / (zero * x);
 197   if (__builtin_expect (ix >= 0x7fff, 0))
 198     return one / (x + x * x);
 199   if (__builtin_expect ((i0 | i1) == 0, 0))
 200     return -HUGE_VALL + x;  /* -inf and overflow exception.  */
 201   if (ix >= 0x4000)
 202     {                           /* |x| >= 2.0 */
 203       __sincosl (x, &s, &c);
 204       ss = -s - c;
 205       cc = s - c;
 206       if (ix < 0x7fe00000)
 207         {                       /* make sure x+x not overflow */
 208           z = __cosl (x + x);
 209           if ((s * c) > zero)
 210             cc = z / ss;
 211           else
 212             ss = z / cc;
 213         }
 214       /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
 215        * where x0 = x-3pi/4
 216        *      Better formula:
 217        *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
 218        *                      =  1/sqrt(2) * (sin(x) - cos(x))
 219        *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
 220        *                      = -1/sqrt(2) * (cos(x) + sin(x))
 221        * To avoid cancellation, use
 222        *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 223        * to compute the worse one.
 224        */
 225       if (__builtin_expect (ix > 0x4080, 0))
 226         z = (invsqrtpi * ss) / __ieee754_sqrtl (x);
 227       else
 228         {
 229           u = pone (x);
 230           v = qone (x);
 231           z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x);
 232         }
 233       return z;
 234     }
 235   if (__builtin_expect (ix <= 0x3fbe, 0))
 236     {                           /* x < 2**-65 */
 237       return (-tpi / x);
 238     }
 239   z = x * x;
 240  u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * (U0[4] + z * U0[5]))));
 241  v = V0[0] + z * (V0[1] + z * (V0[2] + z * (V0[3] + z * (V0[4] + z))));
 242   return (x * (u / v) +
 243           tpi * (__ieee754_j1l (x) * __ieee754_logl (x) - one / x));
 244 }
 245 strong_alias (__ieee754_y1l, __y1l_finite)
 246
 247
 248 /* For x >= 8, the asymptotic expansions of pone is
 249  *      1 + 15/128 s^2 - 4725/2^15 s^4 - ...,   where s = 1/x.
 250  * We approximate pone by
 251  *      pone(x) = 1 + (R/S)
 252  */
 253
 254 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 255    P1(x) = 1 + z^2 R(z^2), z=1/x
 256    8 <= x <= inf  (0 <= z <= 0.125)
 257    Peak relative error 5.2e-22  */
 258
 259 static const long double pr8[7] = {
 260   8.402048819032978959298664869941375143163E-9L,
 261   1.813743245316438056192649247507255996036E-6L,
 262   1.260704554112906152344932388588243836276E-4L,
 263   3.439294839869103014614229832700986965110E-3L,
 264   3.576910849712074184504430254290179501209E-2L,
 265   1.131111483254318243139953003461511308672E-1L,
 266   4.480715825681029711521286449131671880953E-2L,
 267 };
 268 static const long double ps8[6] = {
 269   7.169748325574809484893888315707824924354E-8L,
 270   1.556549720596672576431813934184403614817E-5L,
 271   1.094540125521337139209062035774174565882E-3L,
 272   3.060978962596642798560894375281428805840E-2L,
 273   3.374146536087205506032643098619414507024E-1L,
 274   1.253830208588979001991901126393231302559E0L,
 275   /* 1.000000000000000000000000000000000000000E0L, */
 276 };
 277
 278 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 279    P1(x) = 1 + z^2 R(z^2), z=1/x
 280    4.54541015625 <= x <= 8
 281    Peak relative error 7.7e-22  */
 282 static const long double pr5[7] = {
 283   4.318486887948814529950980396300969247900E-7L,
 284   4.715341880798817230333360497524173929315E-5L,
 285   1.642719430496086618401091544113220340094E-3L,
 286   2.228688005300803935928733750456396149104E-2L,
 287   1.142773760804150921573259605730018327162E-1L,
 288   1.755576530055079253910829652698703791957E-1L,
 289   3.218803858282095929559165965353784980613E-2L,
 290 };
 291 static const long double ps5[6] = {
 292   3.685108812227721334719884358034713967557E-6L,
 293   4.069102509511177498808856515005792027639E-4L,
 294   1.449728676496155025507893322405597039816E-2L,
 295   2.058869213229520086582695850441194363103E-1L,
 296   1.164890985918737148968424972072751066553E0L,
 297   2.274776933457009446573027260373361586841E0L,
 298   /*  1.000000000000000000000000000000000000000E0L,*/
 299 };
 300
 301 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 302    P1(x) = 1 + z^2 R(z^2), z=1/x
 303    2.85711669921875 <= x <= 4.54541015625
 304    Peak relative error 6.5e-21  */
 305 static const long double pr3[7] = {
 306   1.265251153957366716825382654273326407972E-5L,
 307   8.031057269201324914127680782288352574567E-4L,
 308   1.581648121115028333661412169396282881035E-2L,
 309   1.179534658087796321928362981518645033967E-1L,
 310   3.227936912780465219246440724502790727866E-1L,
 311   2.559223765418386621748404398017602935764E-1L,
 312   2.277136933287817911091370397134882441046E-2L,
 313 };
 314 static const long double ps3[6] = {
 315   1.079681071833391818661952793568345057548E-4L,
 316   6.986017817100477138417481463810841529026E-3L,
 317   1.429403701146942509913198539100230540503E-1L,
 318   1.148392024337075609460312658938700765074E0L,
 319   3.643663015091248720208251490291968840882E0L,
 320   3.990702269032018282145100741746633960737E0L,
 321   /* 1.000000000000000000000000000000000000000E0L, */
 322 };
 323
 324 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 325    P1(x) = 1 + z^2 R(z^2), z=1/x
 326    2 <= x <= 2.85711669921875
 327    Peak relative error 3.5e-21  */
 328 static const long double pr2[7] = {
 329   2.795623248568412225239401141338714516445E-4L,
 330   1.092578168441856711925254839815430061135E-2L,
 331   1.278024620468953761154963591853679640560E-1L,
 332   5.469680473691500673112904286228351988583E-1L,
 333   8.313769490922351300461498619045639016059E-1L,
 334   3.544176317308370086415403567097130611468E-1L,
 335   1.604142674802373041247957048801599740644E-2L,
 336 };
 337 static const long double ps2[6] = {
 338   2.385605161555183386205027000675875235980E-3L,
 339   9.616778294482695283928617708206967248579E-2L,
 340   1.195215570959693572089824415393951258510E0L,
 341   5.718412857897054829999458736064922974662E0L,
 342   1.065626298505499086386584642761602177568E1L,
 343   6.809140730053382188468983548092322151791E0L,
 344  /* 1.000000000000000000000000000000000000000E0L, */
 345 };
 346
 347
 348 static long double
 349 pone (long double x)
 350 {
 351   const long double *p, *q;
 352   long double z, r, s;
 353   int32_t ix;
 354   u_int32_t se, i0, i1;
 355
 356   GET_LDOUBLE_WORDS (se, i0, i1, x);
 357   ix = se & 0x7fff;
 358   if (ix >= 0x4002) /* x >= 8 */
 359     {
 360       p = pr8;
 361       q = ps8;
 362     }
 363   else
 364     {
 365       i1 = (ix << 16) | (i0 >> 16);
 366       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 367         {
 368           p = pr5;
 369           q = ps5;
 370         }
 371       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 372         {
 373           p = pr3;
 374           q = ps3;
 375         }
 376       else if (ix >= 0x4000)    /* x better be >= 2 */
 377         {
 378           p = pr2;
 379           q = ps2;
 380         }
 381     }
 382   z = one / (x * x);
 383  r = p[0] + z * (p[1] +
 384                  z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 385  s = q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
 386   return one + z * r / s;
 387 }
 388
 389
 390 /* For x >= 8, the asymptotic expansions of qone is
 391  *      3/8 s - 105/1024 s^3 - ..., where s = 1/x.
 392  * We approximate pone by
 393  *      qone(x) = s*(0.375 + (R/S))
 394  */
 395
 396 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 397    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 398    8 <= x <= inf
 399    Peak relative error 8.3e-22 */
 400
 401 static const long double qr8[7] = {
 402   -5.691925079044209246015366919809404457380E-10L,
 403   -1.632587664706999307871963065396218379137E-7L,
 404   -1.577424682764651970003637263552027114600E-5L,
 405   -6.377627959241053914770158336842725291713E-4L,
 406   -1.087408516779972735197277149494929568768E-2L,
 407   -6.854943629378084419631926076882330494217E-2L,
 408   -1.055448290469180032312893377152490183203E-1L,
 409 };
 410 static const long double qs8[7] = {
 411   5.550982172325019811119223916998393907513E-9L,
 412   1.607188366646736068460131091130644192244E-6L,
 413   1.580792530091386496626494138334505893599E-4L,
 414   6.617859900815747303032860443855006056595E-3L,
 415   1.212840547336984859952597488863037659161E-1L,
 416   9.017885953937234900458186716154005541075E-1L,
 417   2.201114489712243262000939120146436167178E0L,
 418   /* 1.000000000000000000000000000000000000000E0L, */
 419 };
 420
 421 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 422    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 423    4.54541015625 <= x <= 8
 424    Peak relative error 4.1e-22 */
 425 static const long double qr5[7] = {
 426   -6.719134139179190546324213696633564965983E-8L,
 427   -9.467871458774950479909851595678622044140E-6L,
 428   -4.429341875348286176950914275723051452838E-4L,
 429   -8.539898021757342531563866270278505014487E-3L,
 430   -6.818691805848737010422337101409276287170E-2L,
 431   -1.964432669771684034858848142418228214855E-1L,
 432   -1.333896496989238600119596538299938520726E-1L,
 433 };
 434 static const long double qs5[7] = {
 435   6.552755584474634766937589285426911075101E-7L,
 436   9.410814032118155978663509073200494000589E-5L,
 437   4.561677087286518359461609153655021253238E-3L,
 438   9.397742096177905170800336715661091535805E-2L,
 439   8.518538116671013902180962914473967738771E-1L,
 440   3.177729183645800174212539541058292579009E0L,
 441   4.006745668510308096259753538973038902990E0L,
 442   /* 1.000000000000000000000000000000000000000E0L, */
 443 };
 444
 445 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 446    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 447    2.85711669921875 <= x <= 4.54541015625
 448    Peak relative error 2.2e-21 */
 449 static const long double qr3[7] = {
 450   -3.618746299358445926506719188614570588404E-6L,
 451   -2.951146018465419674063882650970344502798E-4L,
 452   -7.728518171262562194043409753656506795258E-3L,
 453   -8.058010968753999435006488158237984014883E-2L,
 454   -3.356232856677966691703904770937143483472E-1L,
 455   -4.858192581793118040782557808823460276452E-1L,
 456   -1.592399251246473643510898335746432479373E-1L,
 457 };
 458 static const long double qs3[7] = {
 459   3.529139957987837084554591421329876744262E-5L,
 460   2.973602667215766676998703687065066180115E-3L,
 461   8.273534546240864308494062287908662592100E-2L,
 462   9.613359842126507198241321110649974032726E-1L,
 463   4.853923697093974370118387947065402707519E0L,
 464   1.002671608961669247462020977417828796933E1L,
 465   7.028927383922483728931327850683151410267E0L,
 466   /* 1.000000000000000000000000000000000000000E0L, */
 467 };
 468
 469 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 470    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 471    2 <= x <= 2.85711669921875
 472    Peak relative error 6.9e-22 */
 473 static const long double qr2[7] = {
 474   -1.372751603025230017220666013816502528318E-4L,
 475   -6.879190253347766576229143006767218972834E-3L,
 476   -1.061253572090925414598304855316280077828E-1L,
 477   -6.262164224345471241219408329354943337214E-1L,
 478   -1.423149636514768476376254324731437473915E0L,
 479   -1.087955310491078933531734062917489870754E0L,
 480   -1.826821119773182847861406108689273719137E-1L,
 481 };
 482 static const long double qs2[7] = {
 483   1.338768933634451601814048220627185324007E-3L,
 484   7.071099998918497559736318523932241901810E-2L,
 485   1.200511429784048632105295629933382142221E0L,
 486   8.327301713640367079030141077172031825276E0L,
 487   2.468479301872299311658145549931764426840E1L,
 488   2.961179686096262083509383820557051621644E1L,
 489   1.201402313144305153005639494661767354977E1L,
 490  /* 1.000000000000000000000000000000000000000E0L, */
 491 };
 492
 493
 494 static long double
 495 qone (long double x)
 496 {
 497   const long double *p, *q;
 498   static long double s, r, z;
 499   int32_t ix;
 500   u_int32_t se, i0, i1;
 501
 502   GET_LDOUBLE_WORDS (se, i0, i1, x);
 503   ix = se & 0x7fff;
 504   if (ix >= 0x4002)             /* x >= 8 */
 505     {
 506       p = qr8;
 507       q = qs8;
 508     }
 509   else
 510     {
 511       i1 = (ix << 16) | (i0 >> 16);
 512       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 513         {
 514           p = qr5;
 515           q = qs5;
 516         }
 517       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 518         {
 519           p = qr3;
 520           q = qs3;
 521         }
 522       else if (ix >= 0x4000)    /* x better be >= 2 */
 523         {
 524           p = qr2;
 525           q = qs2;
 526         }
 527     }
 528   z = one / (x * x);
 529   r =
 530     p[0] + z * (p[1] +
 531                 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 532   s =
 533     q[0] + z * (q[1] +
 534                 z * (q[2] +
 535                      z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
 536   return (.375 + z * r / s) / x;
 537 }