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