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 "math.h"
  60 #include "math_private.h"
  61
  62 static const long double
  63   invsqrtpi = 5.6418958354775628694807945156077258584405E-1L,
  64   two = 2.0e0L,
  65   one = 1.0e0L,
  66   zero = 0.0L;
  67
  68
  69 long double
  70 __ieee754_jnl (int n, long double x)
  71 {
  72   u_int32_t se;
  73   int32_t i, ix, sgn;
  74   long double a, b, temp, di;
  75   long double z, w;
  76   ieee854_long_double_shape_type u;
  77
  78
  79   /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
  80    * Thus, J(-n,x) = J(n,-x)
  81    */
  82
  83   u.value = x;
  84   se = u.parts32.w0;
  85   ix = se & 0x7fffffff;
  86
  87   /* if J(n,NaN) is NaN */
  88   if (ix >= 0x7ff00000)
  89     {
  90       if ((u.parts32.w0 & 0xfffff) | u.parts32.w1
  91           | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3)
  92         return x + x;
  93     }
  94
  95   if (n < 0)
  96     {
  97       n = -n;
  98       x = -x;
  99       se ^= 0x80000000;
 100     }
 101   if (n == 0)
 102     return (__ieee754_j0l (x));
 103   if (n == 1)
 104     return (__ieee754_j1l (x));
 105   sgn = (n & 1) & (se >> 31);   /* even n -- 0, odd n -- sign(x) */
 106   x = fabsl (x);
 107
 108   if (x == 0.0L || ix >= 0x7ff00000)    /* if x is 0 or inf */
 109     b = zero;
 110   else if ((long double) n <= x)
 111     {
 112       /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
 113       if (ix >= 0x52d00000)
 114         {                       /* x > 2**302 */
 115
 116           /* ??? Could use an expansion for large x here.  */
 117
 118           /* (x >> n**2)
 119            *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 120            *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 121            *      Let s=sin(x), c=cos(x),
 122            *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 123            *
 124            *             n    sin(xn)*sqt2    cos(xn)*sqt2
 125            *          ----------------------------------
 126            *             0     s-c             c+s
 127            *             1    -s-c            -c+s
 128            *             2    -s+c            -c-s
 129            *             3     s+c             c-s
 130            */
 131           long double s;
 132           long double c;
 133           __sincosl (x, &s, &c);
 134           switch (n & 3)
 135             {
 136             case 0:
 137               temp = c + s;
 138               break;
 139             case 1:
 140               temp = -c + s;
 141               break;
 142             case 2:
 143               temp = -c - s;
 144               break;
 145             case 3:
 146               temp = c - s;
 147               break;
 148             }
 149           b = invsqrtpi * temp / __ieee754_sqrtl (x);
 150         }
 151       else
 152         {
 153           a = __ieee754_j0l (x);
 154           b = __ieee754_j1l (x);
 155           for (i = 1; i < n; i++)
 156             {
 157               temp = b;
 158               b = b * ((long double) (i + i) / x) - a;  /* avoid underflow */
 159               a = temp;
 160             }
 161         }
 162     }
 163   else
 164     {
 165       if (ix < 0x3e100000)
 166         {                       /* x < 2**-29 */
 167           /* x is tiny, return the first Taylor expansion of J(n,x)
 168            * J(n,x) = 1/n!*(x/2)^n  - ...
 169            */
 170           if (n >= 33)          /* underflow, result < 10^-300 */
 171             b = zero;
 172           else
 173             {
 174               temp = x * 0.5;
 175               b = temp;
 176               for (a = one, i = 2; i <= n; i++)
 177                 {
 178                   a *= (long double) i; /* a = n! */
 179                   b *= temp;    /* b = (x/2)^n */
 180                 }
 181               b = b / a;
 182             }
 183         }
 184       else
 185         {
 186           /* use backward recurrence */
 187           /*                      x      x^2      x^2
 188            *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
 189            *                      2n  - 2(n+1) - 2(n+2)
 190            *
 191            *                      1      1        1
 192            *  (for large x)   =  ----  ------   ------   .....
 193            *                      2n   2(n+1)   2(n+2)
 194            *                      -- - ------ - ------ -
 195            *                       x     x         x
 196            *
 197            * Let w = 2n/x and h=2/x, then the above quotient
 198            * is equal to the continued fraction:
 199            *                  1
 200            *      = -----------------------
 201            *                     1
 202            *         w - -----------------
 203            *                        1
 204            *              w+h - ---------
 205            *                     w+2h - ...
 206            *
 207            * To determine how many terms needed, let
 208            * Q(0) = w, Q(1) = w(w+h) - 1,
 209            * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
 210            * When Q(k) > 1e4      good for single
 211            * When Q(k) > 1e9      good for double
 212            * When Q(k) > 1e17     good for quadruple
 213            */
 214           /* determine k */
 215           long double t, v;
 216           long double q0, q1, h, tmp;
 217           int32_t k, m;
 218           w = (n + n) / (long double) x;
 219           h = 2.0L / (long double) x;
 220           q0 = w;
 221           z = w + h;
 222           q1 = w * z - 1.0L;
 223           k = 1;
 224           while (q1 < 1.0e17L)
 225             {
 226               k += 1;
 227               z += h;
 228               tmp = z * q1 - q0;
 229               q0 = q1;
 230               q1 = tmp;
 231             }
 232           m = n + n;
 233           for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
 234             t = one / (i / x - t);
 235           a = t;
 236           b = one;
 237           /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
 238            *  Hence, if n*(log(2n/x)) > ...
 239            *  single 8.8722839355e+01
 240            *  double 7.09782712893383973096e+02
 241            *  long double 1.1356523406294143949491931077970765006170e+04
 242            *  then recurrent value may overflow and the result is
 243            *  likely underflow to zero
 244            */
 245           tmp = n;
 246           v = two / x;
 247           tmp = tmp * __ieee754_logl (fabsl (v * tmp));
 248
 249           if (tmp < 1.1356523406294143949491931077970765006170e+04L)
 250             {
 251               for (i = n - 1, di = (long double) (i + i); i > 0; i--)
 252                 {
 253                   temp = b;
 254                   b *= di;
 255                   b = b / x - a;
 256                   a = temp;
 257                   di -= two;
 258                 }
 259             }
 260           else
 261             {
 262               for (i = n - 1, di = (long double) (i + i); i > 0; i--)
 263                 {
 264                   temp = b;
 265                   b *= di;
 266                   b = b / x - a;
 267                   a = temp;
 268                   di -= two;
 269                   /* scale b to avoid spurious overflow */
 270                   if (b > 1e100L)
 271                     {
 272                       a /= b;
 273                       t /= b;
 274                       b = one;
 275                     }
 276                 }
 277             }
 278           /* j0() and j1() suffer enormous loss of precision at and
 279            * near zero; however, we know that their zero points never
 280            * coincide, so just choose the one further away from zero.
 281            */
 282           z = __ieee754_j0l (x);
 283           w = __ieee754_j1l (x);
 284           if (fabsl (z) >= fabsl (w))
 285             b = (t * z / b);
 286           else
 287             b = (t * w / a);
 288         }
 289     }
 290   if (sgn == 1)
 291     return -b;
 292   else
 293     return b;
 294 }
 295 strong_alias (__ieee754_jnl, __jnl_finite)
 296
 297 long double
 298 __ieee754_ynl (int n, long double x)
 299 {
 300   u_int32_t se;
 301   int32_t i, ix;
 302   int32_t sign;
 303   long double a, b, temp;
 304   ieee854_long_double_shape_type u;
 305
 306   u.value = x;
 307   se = u.parts32.w0;
 308   ix = se & 0x7fffffff;
 309
 310   /* if Y(n,NaN) is NaN */
 311   if (ix >= 0x7ff00000)
 312     {
 313       if ((u.parts32.w0 & 0xfffff) | u.parts32.w1
 314           | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3)
 315         return x + x;
 316     }
 317   if (x <= 0.0L)
 318     {
 319       if (x == 0.0L)
 320         return -HUGE_VALL + x;
 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   if (n == 1)
 333     return (sign * __ieee754_y1l (x));
 334   if (ix >= 0x7ff00000)
 335     return zero;
 336   if (ix >= 0x52D00000)
 337     {                           /* x > 2**302 */
 338
 339       /* ??? See comment above on the possible futility of this.  */
 340
 341       /* (x >> n**2)
 342        *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 343        *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 344        *      Let s=sin(x), c=cos(x),
 345        *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 346        *
 347        *             n    sin(xn)*sqt2    cos(xn)*sqt2
 348        *          ----------------------------------
 349        *             0     s-c             c+s
 350        *             1    -s-c            -c+s
 351        *             2    -s+c            -c-s
 352        *             3     s+c             c-s
 353        */
 354       long double s;
 355       long double c;
 356       __sincosl (x, &s, &c);
 357       switch (n & 3)
 358         {
 359         case 0:
 360           temp = s - c;
 361           break;
 362         case 1:
 363           temp = -s - c;
 364           break;
 365         case 2:
 366           temp = -s + c;
 367           break;
 368         case 3:
 369           temp = s + c;
 370           break;
 371         }
 372       b = invsqrtpi * temp / __ieee754_sqrtl (x);
 373     }
 374   else
 375     {
 376       a = __ieee754_y0l (x);
 377       b = __ieee754_y1l (x);
 378       /* quit if b is -inf */
 379       u.value = b;
 380       se = u.parts32.w0 & 0xfff00000;
 381       for (i = 1; i < n && se != 0xfff00000; i++)
 382         {
 383           temp = b;
 384           b = ((long double) (i + i) / x) * b - a;
 385           u.value = b;
 386           se = u.parts32.w0 & 0xfff00000;
 387           a = temp;
 388         }
 389     }
 390   if (sign > 0)
 391     return b;
 392   else
 393     return -b;
 394 }
 395 strong_alias (__ieee754_ynl, __ynl_finite)