sysdeps/ieee754/ldbl-128ibm/e_jnl.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 /* Modifications for 128-bit long double 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 /*
  34  * __ieee754_jn(n, x), __ieee754_yn(n, x)
  35  * floating point Bessel's function of the 1st and 2nd kind
  36  * of order n
  37  *
  38  * Special cases:
  39  *      y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
  40  *      y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
  41  * Note 2. About jn(n,x), yn(n,x)
  42  *      For n=0, j0(x) is called,
  43  *      for n=1, j1(x) is called,
  44  *      for n<x, forward recursion us used starting
  45  *      from values of j0(x) and j1(x).
  46  *      for n>x, a continued fraction approximation to
  47  *      j(n,x)/j(n-1,x) is evaluated and then backward
  48  *      recursion is used starting from a supposed value
  49  *      for j(n,x). The resulting value of j(0,x) is
  50  *      compared with the actual value to correct the
  51  *      supposed value of j(n,x).
  52  *
  53  *      yn(n,x) is similar in all respects, except
  54  *      that forward recursion is used for all
  55  *      values of n>1.
  56  *
  57  */
  58
  59 #include <errno.h>
  60 #include <float.h>
  61 #include <math.h>
  62 #include <math_private.h>
  63
  64 static const long double
  65   invsqrtpi = 5.6418958354775628694807945156077258584405E-1L,
  66   two = 2.0e0L,
  67   one = 1.0e0L,
  68   zero = 0.0L;
  69
  70
  71 long double
  72 __ieee754_jnl (int n, long double x)
  73 {
  74   uint32_t se, lx;
  75   int32_t i, ix, sgn;
  76   long double a, b, temp, di, ret;
  77   long double z, w;
  78   double xhi;
  79
  80
  81   /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
  82    * Thus, J(-n,x) = J(n,-x)
  83    */
  84
  85   xhi = ldbl_high (x);
  86   EXTRACT_WORDS (se, lx, xhi);
  87   ix = se & 0x7fffffff;
  88
  89   /* if J(n,NaN) is NaN */
  90   if (ix >= 0x7ff00000)
  91     {
  92       if (((ix - 0x7ff00000) | lx) != 0)
  93         return x + x;
  94     }
  95
  96   if (n < 0)
  97     {
  98       n = -n;
  99       x = -x;
 100       se ^= 0x80000000;
 101     }
 102   if (n == 0)
 103     return (__ieee754_j0l (x));
 104   if (n == 1)
 105     return (__ieee754_j1l (x));
 106   sgn = (n & 1) & (se >> 31);   /* even n -- 0, odd n -- sign(x) */
 107   x = fabsl (x);
 108
 109   {
 110     SET_RESTORE_ROUNDL (FE_TONEAREST);
 111     if (x == 0.0L || ix >= 0x7ff00000)  /* if x is 0 or inf */
 112       return sgn == 1 ? -zero : zero;
 113     else if ((long double) n <= x)
 114       {
 115         /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
 116         if (ix >= 0x52d00000)
 117           {                     /* x > 2**302 */
 118
 119             /* ??? Could use an expansion for large x here.  */
 120
 121             /* (x >> n**2)
 122              *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 123              *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 124              *      Let s=sin(x), c=cos(x),
 125              *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 126              *
 127              *             n    sin(xn)*sqt2    cos(xn)*sqt2
 128              *          ----------------------------------
 129              *             0     s-c             c+s
 130              *             1    -s-c            -c+s
 131              *             2    -s+c            -c-s
 132              *             3     s+c             c-s
 133              */
 134             long double s;
 135             long double c;
 136             __sincosl (x, &s, &c);
 137             switch (n & 3)
 138               {
 139               case 0:
 140                 temp = c + s;
 141                 break;
 142               case 1:
 143                 temp = -c + s;
 144                 break;
 145               case 2:
 146                 temp = -c - s;
 147                 break;
 148               case 3:
 149                 temp = c - s;
 150                 break;
 151               }
 152             b = invsqrtpi * temp / __ieee754_sqrtl (x);
 153           }
 154         else
 155           {
 156             a = __ieee754_j0l (x);
 157             b = __ieee754_j1l (x);
 158             for (i = 1; i < n; i++)
 159               {
 160                 temp = b;
 161                 b = b * ((long double) (i + i) / x) - a;        /* avoid underflow */
 162                 a = temp;
 163               }
 164           }
 165       }
 166     else
 167       {
 168         if (ix < 0x3e100000)
 169           {                     /* x < 2**-29 */
 170             /* x is tiny, return the first Taylor expansion of J(n,x)
 171              * J(n,x) = 1/n!*(x/2)^n  - ...
 172              */
 173             if (n >= 33)                /* underflow, result < 10^-300 */
 174               b = zero;
 175             else
 176               {
 177                 temp = x * 0.5;
 178                 b = temp;
 179                 for (a = one, i = 2; i <= n; i++)
 180                   {
 181                     a *= (long double) i;       /* a = n! */
 182                     b *= temp;  /* b = (x/2)^n */
 183                   }
 184                 b = b / a;
 185               }
 186           }
 187         else
 188           {
 189             /* use backward recurrence */
 190             /*                      x      x^2      x^2
 191              *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
 192              *                      2n  - 2(n+1) - 2(n+2)
 193              *
 194              *                      1      1        1
 195              *  (for large x)   =  ----  ------   ------   .....
 196              *                      2n   2(n+1)   2(n+2)
 197              *                      -- - ------ - ------ -
 198              *                       x     x         x
 199              *
 200              * Let w = 2n/x and h=2/x, then the above quotient
 201              * is equal to the continued fraction:
 202              *                  1
 203              *      = -----------------------
 204              *                     1
 205              *         w - -----------------
 206              *                        1
 207              *              w+h - ---------
 208              *                     w+2h - ...
 209              *
 210              * To determine how many terms needed, let
 211              * Q(0) = w, Q(1) = w(w+h) - 1,
 212              * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
 213              * When Q(k) > 1e4      good for single
 214              * When Q(k) > 1e9      good for double
 215              * When Q(k) > 1e17     good for quadruple
 216              */
 217             /* determine k */
 218             long double t, v;
 219             long double q0, q1, h, tmp;
 220             int32_t k, m;
 221             w = (n + n) / (long double) x;
 222             h = 2.0L / (long double) x;
 223             q0 = w;
 224             z = w + h;
 225             q1 = w * z - 1.0L;
 226             k = 1;
 227             while (q1 < 1.0e17L)
 228               {
 229                 k += 1;
 230                 z += h;
 231                 tmp = z * q1 - q0;
 232                 q0 = q1;
 233                 q1 = tmp;
 234               }
 235             m = n + n;
 236             for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
 237               t = one / (i / x - t);
 238             a = t;
 239             b = one;
 240             /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
 241              *  Hence, if n*(log(2n/x)) > ...
 242              *  single 8.8722839355e+01
 243              *  double 7.09782712893383973096e+02
 244              *  long double 1.1356523406294143949491931077970765006170e+04
 245              *  then recurrent value may overflow and the result is
 246              *  likely underflow to zero
 247              */
 248             tmp = n;
 249             v = two / x;
 250             tmp = tmp * __ieee754_logl (fabsl (v * tmp));
 251
 252             if (tmp < 1.1356523406294143949491931077970765006170e+04L)
 253               {
 254                 for (i = n - 1, di = (long double) (i + i); i > 0; i--)
 255                   {
 256                     temp = b;
 257                     b *= di;
 258                     b = b / x - a;
 259                     a = temp;
 260                     di -= two;
 261                   }
 262               }
 263             else
 264               {
 265                 for (i = n - 1, di = (long double) (i + i); i > 0; i--)
 266                   {
 267                     temp = b;
 268                     b *= di;
 269                     b = b / x - a;
 270                     a = temp;
 271                     di -= two;
 272                     /* scale b to avoid spurious overflow */
 273                     if (b > 1e100L)
 274                       {
 275                         a /= b;
 276                         t /= b;
 277                         b = one;
 278                       }
 279                   }
 280               }
 281             /* j0() and j1() suffer enormous loss of precision at and
 282              * near zero; however, we know that their zero points never
 283              * coincide, so just choose the one further away from zero.
 284              */
 285             z = __ieee754_j0l (x);
 286             w = __ieee754_j1l (x);
 287             if (fabsl (z) >= fabsl (w))
 288               b = (t * z / b);
 289             else
 290               b = (t * w / a);
 291           }
 292       }
 293     if (sgn == 1)
 294       ret = -b;
 295     else
 296       ret = b;
 297   }
 298   if (ret == 0)
 299     ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN;
 300   return ret;
 301 }
 302 strong_alias (__ieee754_jnl, __jnl_finite)
 303
 304 long double
 305 __ieee754_ynl (int n, long double x)
 306 {
 307   uint32_t se, lx;
 308   int32_t i, ix;
 309   int32_t sign;
 310   long double a, b, temp, ret;
 311   double xhi;
 312
 313   xhi = ldbl_high (x);
 314   EXTRACT_WORDS (se, lx, xhi);
 315   ix = se & 0x7fffffff;
 316
 317   /* if Y(n,NaN) is NaN */
 318   if (ix >= 0x7ff00000)
 319     {
 320       if (((ix - 0x7ff00000) | lx) != 0)
 321         return x + x;
 322     }
 323   if (x <= 0.0L)
 324     {
 325       if (x == 0.0L)
 326         return ((n < 0 && (n & 1) != 0) ? 1.0L : -1.0L) / 0.0L;
 327       if (se & 0x80000000)
 328         return zero / (zero * x);
 329     }
 330   sign = 1;
 331   if (n < 0)
 332     {
 333       n = -n;
 334       sign = 1 - ((n & 1) << 1);
 335     }
 336   if (n == 0)
 337     return (__ieee754_y0l (x));
 338   {
 339     SET_RESTORE_ROUNDL (FE_TONEAREST);
 340     if (n == 1)
 341       {
 342         ret = sign * __ieee754_y1l (x);
 343         goto out;
 344       }
 345     if (ix >= 0x7ff00000)
 346       return zero;
 347     if (ix >= 0x52D00000)
 348       {                         /* x > 2**302 */
 349
 350         /* ??? See comment above on the possible futility of this.  */
 351
 352         /* (x >> n**2)
 353          *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 354          *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 355          *      Let s=sin(x), c=cos(x),
 356          *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 357          *
 358          *             n    sin(xn)*sqt2    cos(xn)*sqt2
 359          *          ----------------------------------
 360          *             0     s-c             c+s
 361          *             1    -s-c            -c+s
 362          *             2    -s+c            -c-s
 363          *             3     s+c             c-s
 364          */
 365         long double s;
 366         long double c;
 367         __sincosl (x, &s, &c);
 368         switch (n & 3)
 369           {
 370           case 0:
 371             temp = s - c;
 372             break;
 373           case 1:
 374             temp = -s - c;
 375             break;
 376           case 2:
 377             temp = -s + c;
 378             break;
 379           case 3:
 380             temp = s + c;
 381             break;
 382           }
 383         b = invsqrtpi * temp / __ieee754_sqrtl (x);
 384       }
 385     else
 386       {
 387         a = __ieee754_y0l (x);
 388         b = __ieee754_y1l (x);
 389         /* quit if b is -inf */
 390         xhi = ldbl_high (b);
 391         GET_HIGH_WORD (se, xhi);
 392         se &= 0xfff00000;
 393         for (i = 1; i < n && se != 0xfff00000; i++)
 394           {
 395             temp = b;
 396             b = ((long double) (i + i) / x) * b - a;
 397             xhi = ldbl_high (b);
 398             GET_HIGH_WORD (se, xhi);
 399             se &= 0xfff00000;
 400             a = temp;
 401           }
 402       }
 403     /* If B is +-Inf, set up errno accordingly.  */
 404     if (! isfinite (b))
 405       __set_errno (ERANGE);
 406     if (sign > 0)
 407       ret = b;
 408     else
 409       ret = -b;
 410   }
 411  out:
 412   if (isinf (ret))
 413     ret = __copysignl (LDBL_MAX, ret) * LDBL_MAX;
 414   return ret;
 415 }
 416 strong_alias (__ieee754_ynl, __ynl_finite)