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