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