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