sysdeps/ieee754/ldbl-128/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 <math.h>
  61 #include <math_private.h>
  62
  63 static const long double
  64   invsqrtpi = 5.6418958354775628694807945156077258584405E-1L,
  65   two = 2.0e0L,
  66   one = 1.0e0L,
  67   zero = 0.0L;
  68
  69
  70 long double
  71 __ieee754_jnl (int n, long double x)
  72 {
  73   u_int32_t se;
  74   int32_t i, ix, sgn;
  75   long double a, b, temp, di;
  76   long double z, w;
  77   ieee854_long_double_shape_type u;
  78
  79
  80   /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
  81    * Thus, J(-n,x) = J(n,-x)
  82    */
  83
  84   u.value = x;
  85   se = u.parts32.w0;
  86   ix = se & 0x7fffffff;
  87
  88   /* if J(n,NaN) is NaN */
  89   if (ix >= 0x7fff0000)
  90     {
  91       if ((u.parts32.w0 & 0xffff) | u.parts32.w1 | u.parts32.w2 | 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 >= 0x7fff0000)    /* 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 >= 0x412D0000)
 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 < 0x3fc60000)
 166         {                       /* x < 2**-57 */
 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 >= 400)         /* underflow, result < 10^-4952 */
 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 >= 0x7fff0000)
 312     {
 313       if ((u.parts32.w0 & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3)
 314         return x + x;
 315     }
 316   if (x <= 0.0L)
 317     {
 318       if (x == 0.0L)
 319         return ((n < 0 && (n & 1) != 0) ? 1.0L : -1.0L) / 0.0L;
 320       if (se & 0x80000000)
 321         return zero / (zero * x);
 322     }
 323   sign = 1;
 324   if (n < 0)
 325     {
 326       n = -n;
 327       sign = 1 - ((n & 1) << 1);
 328     }
 329   if (n == 0)
 330     return (__ieee754_y0l (x));
 331   if (n == 1)
 332     return (sign * __ieee754_y1l (x));
 333   if (ix >= 0x7fff0000)
 334     return zero;
 335   if (ix >= 0x412D0000)
 336     {                           /* x > 2**302 */
 337
 338       /* ??? See comment above on the possible futility of this.  */
 339
 340       /* (x >> n**2)
 341        *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 342        *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 343        *      Let s=sin(x), c=cos(x),
 344        *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 345        *
 346        *             n    sin(xn)*sqt2    cos(xn)*sqt2
 347        *          ----------------------------------
 348        *             0     s-c             c+s
 349        *             1    -s-c            -c+s
 350        *             2    -s+c            -c-s
 351        *             3     s+c             c-s
 352        */
 353       long double s;
 354       long double c;
 355       __sincosl (x, &s, &c);
 356       switch (n & 3)
 357         {
 358         case 0:
 359           temp = s - c;
 360           break;
 361         case 1:
 362           temp = -s - c;
 363           break;
 364         case 2:
 365           temp = -s + c;
 366           break;
 367         case 3:
 368           temp = s + c;
 369           break;
 370         }
 371       b = invsqrtpi * temp / __ieee754_sqrtl (x);
 372     }
 373   else
 374     {
 375       a = __ieee754_y0l (x);
 376       b = __ieee754_y1l (x);
 377       /* quit if b is -inf */
 378       u.value = b;
 379       se = u.parts32.w0 & 0xffff0000;
 380       for (i = 1; i < n && se != 0xffff0000; i++)
 381         {
 382           temp = b;
 383           b = ((long double) (i + i) / x) * b - a;
 384           u.value = b;
 385           se = u.parts32.w0 & 0xffff0000;
 386           a = temp;
 387         }
 388     }
 389   /* If B is +-Inf, set up errno accordingly.  */
 390   if (! __finitel (b))
 391     __set_errno (ERANGE);
 392   if (sign > 0)
 393     return b;
 394   else
 395     return -b;
 396 }
 397 strong_alias (__ieee754_ynl, __ynl_finite)