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;
  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   if (x == 0.0L || ix >= 0x7ff00000)    /* if x is 0 or inf */
 110     b = zero;
 111   else if ((long double) n <= x)
 112     {
 113       /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
 114       if (ix >= 0x52d00000)
 115         {                       /* x > 2**302 */
 116
 117           /* ??? Could use an expansion for large x here.  */
 118
 119           /* (x >> n**2)
 120            *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 121            *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 122            *      Let s=sin(x), c=cos(x),
 123            *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 124            *
 125            *             n    sin(xn)*sqt2    cos(xn)*sqt2
 126            *          ----------------------------------
 127            *             0     s-c             c+s
 128            *             1    -s-c            -c+s
 129            *             2    -s+c            -c-s
 130            *             3     s+c             c-s
 131            */
 132           long double s;
 133           long double c;
 134           __sincosl (x, &s, &c);
 135           switch (n & 3)
 136             {
 137             case 0:
 138               temp = c + s;
 139               break;
 140             case 1:
 141               temp = -c + s;
 142               break;
 143             case 2:
 144               temp = -c - s;
 145               break;
 146             case 3:
 147               temp = c - s;
 148               break;
 149             }
 150           b = invsqrtpi * temp / __ieee754_sqrtl (x);
 151         }
 152       else
 153         {
 154           a = __ieee754_j0l (x);
 155           b = __ieee754_j1l (x);
 156           for (i = 1; i < n; i++)
 157             {
 158               temp = b;
 159               b = b * ((long double) (i + i) / x) - a;  /* avoid underflow */
 160               a = temp;
 161             }
 162         }
 163     }
 164   else
 165     {
 166       if (ix < 0x3e100000)
 167         {                       /* x < 2**-29 */
 168           /* x is tiny, return the first Taylor expansion of J(n,x)
 169            * J(n,x) = 1/n!*(x/2)^n  - ...
 170            */
 171           if (n >= 33)          /* underflow, result < 10^-300 */
 172             b = zero;
 173           else
 174             {
 175               temp = x * 0.5;
 176               b = temp;
 177               for (a = one, i = 2; i <= n; i++)
 178                 {
 179                   a *= (long double) i; /* a = n! */
 180                   b *= temp;    /* b = (x/2)^n */
 181                 }
 182               b = b / a;
 183             }
 184         }
 185       else
 186         {
 187           /* use backward recurrence */
 188           /*                      x      x^2      x^2
 189            *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
 190            *                      2n  - 2(n+1) - 2(n+2)
 191            *
 192            *                      1      1        1
 193            *  (for large x)   =  ----  ------   ------   .....
 194            *                      2n   2(n+1)   2(n+2)
 195            *                      -- - ------ - ------ -
 196            *                       x     x         x
 197            *
 198            * Let w = 2n/x and h=2/x, then the above quotient
 199            * is equal to the continued fraction:
 200            *                  1
 201            *      = -----------------------
 202            *                     1
 203            *         w - -----------------
 204            *                        1
 205            *              w+h - ---------
 206            *                     w+2h - ...
 207            *
 208            * To determine how many terms needed, let
 209            * Q(0) = w, Q(1) = w(w+h) - 1,
 210            * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
 211            * When Q(k) > 1e4      good for single
 212            * When Q(k) > 1e9      good for double
 213            * When Q(k) > 1e17     good for quadruple
 214            */
 215           /* determine k */
 216           long double t, v;
 217           long double q0, q1, h, tmp;
 218           int32_t k, m;
 219           w = (n + n) / (long double) x;
 220           h = 2.0L / (long double) x;
 221           q0 = w;
 222           z = w + h;
 223           q1 = w * z - 1.0L;
 224           k = 1;
 225           while (q1 < 1.0e17L)
 226             {
 227               k += 1;
 228               z += h;
 229               tmp = z * q1 - q0;
 230               q0 = q1;
 231               q1 = tmp;
 232             }
 233           m = n + n;
 234           for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
 235             t = one / (i / x - t);
 236           a = t;
 237           b = one;
 238           /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
 239            *  Hence, if n*(log(2n/x)) > ...
 240            *  single 8.8722839355e+01
 241            *  double 7.09782712893383973096e+02
 242            *  long double 1.1356523406294143949491931077970765006170e+04
 243            *  then recurrent value may overflow and the result is
 244            *  likely underflow to zero
 245            */
 246           tmp = n;
 247           v = two / x;
 248           tmp = tmp * __ieee754_logl (fabsl (v * tmp));
 249
 250           if (tmp < 1.1356523406294143949491931077970765006170e+04L)
 251             {
 252               for (i = n - 1, di = (long double) (i + i); i > 0; i--)
 253                 {
 254                   temp = b;
 255                   b *= di;
 256                   b = b / x - a;
 257                   a = temp;
 258                   di -= two;
 259                 }
 260             }
 261           else
 262             {
 263               for (i = n - 1, di = (long double) (i + i); i > 0; i--)
 264                 {
 265                   temp = b;
 266                   b *= di;
 267                   b = b / x - a;
 268                   a = temp;
 269                   di -= two;
 270                   /* scale b to avoid spurious overflow */
 271                   if (b > 1e100L)
 272                     {
 273                       a /= b;
 274                       t /= b;
 275                       b = one;
 276                     }
 277                 }
 278             }
 279           /* j0() and j1() suffer enormous loss of precision at and
 280            * near zero; however, we know that their zero points never
 281            * coincide, so just choose the one further away from zero.
 282            */
 283           z = __ieee754_j0l (x);
 284           w = __ieee754_j1l (x);
 285           if (fabsl (z) >= fabsl (w))
 286             b = (t * z / b);
 287           else
 288             b = (t * w / a);
 289         }
 290     }
 291   if (sgn == 1)
 292     return -b;
 293   else
 294     return b;
 295 }
 296 strong_alias (__ieee754_jnl, __jnl_finite)
 297
 298 long double
 299 __ieee754_ynl (int n, long double x)
 300 {
 301   uint32_t se, lx;
 302   int32_t i, ix;
 303   int32_t sign;
 304   long double a, b, temp, ret;
 305   double xhi;
 306
 307   xhi = ldbl_high (x);
 308   EXTRACT_WORDS (se, lx, xhi);
 309   ix = se & 0x7fffffff;
 310
 311   /* if Y(n,NaN) is NaN */
 312   if (ix >= 0x7ff00000)
 313     {
 314       if (((ix - 0x7ff00000) | lx) != 0)
 315         return x + x;
 316     }
 317   if (x <= 0.0L)
 318     {
 319       if (x == 0.0L)
 320         return ((n < 0 && (n & 1) != 0) ? 1.0L : -1.0L) / 0.0L;
 321       if (se & 0x80000000)
 322         return zero / (zero * x);
 323     }
 324   sign = 1;
 325   if (n < 0)
 326     {
 327       n = -n;
 328       sign = 1 - ((n & 1) << 1);
 329     }
 330   if (n == 0)
 331     return (__ieee754_y0l (x));
 332   {
 333     SET_RESTORE_ROUNDL (FE_TONEAREST);
 334     if (n == 1)
 335       {
 336         ret = sign * __ieee754_y1l (x);
 337         goto out;
 338       }
 339     if (ix >= 0x7ff00000)
 340       return zero;
 341     if (ix >= 0x52D00000)
 342       {                         /* x > 2**302 */
 343
 344         /* ??? See comment above on the possible futility of this.  */
 345
 346         /* (x >> n**2)
 347          *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 348          *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 349          *      Let s=sin(x), c=cos(x),
 350          *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 351          *
 352          *             n    sin(xn)*sqt2    cos(xn)*sqt2
 353          *          ----------------------------------
 354          *             0     s-c             c+s
 355          *             1    -s-c            -c+s
 356          *             2    -s+c            -c-s
 357          *             3     s+c             c-s
 358          */
 359         long double s;
 360         long double c;
 361         __sincosl (x, &s, &c);
 362         switch (n & 3)
 363           {
 364           case 0:
 365             temp = s - c;
 366             break;
 367           case 1:
 368             temp = -s - c;
 369             break;
 370           case 2:
 371             temp = -s + c;
 372             break;
 373           case 3:
 374             temp = s + c;
 375             break;
 376           }
 377         b = invsqrtpi * temp / __ieee754_sqrtl (x);
 378       }
 379     else
 380       {
 381         a = __ieee754_y0l (x);
 382         b = __ieee754_y1l (x);
 383         /* quit if b is -inf */
 384         xhi = ldbl_high (b);
 385         GET_HIGH_WORD (se, xhi);
 386         se &= 0xfff00000;
 387         for (i = 1; i < n && se != 0xfff00000; i++)
 388           {
 389             temp = b;
 390             b = ((long double) (i + i) / x) * b - a;
 391             xhi = ldbl_high (b);
 392             GET_HIGH_WORD (se, xhi);
 393             se &= 0xfff00000;
 394             a = temp;
 395           }
 396       }
 397     /* If B is +-Inf, set up errno accordingly.  */
 398     if (! isfinite (b))
 399       __set_errno (ERANGE);
 400     if (sign > 0)
 401       ret = b;
 402     else
 403       ret = -b;
 404   }
 405  out:
 406   if (isinf (ret))
 407     ret = __copysignl (LDBL_MAX, ret) * LDBL_MAX;
 408   return ret;
 409 }
 410 strong_alias (__ieee754_ynl, __ynl_finite)