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           b = (t * __ieee754_j0l (x) / b);
 290         }
 291     }
 292   if (sgn == 1)
 293     return -b;
 294   else
 295     return b;
 296 }
 297
 298 #ifdef __STDC__
 299 long double
 300 __ieee754_ynl (int n, long double x)
 301 #else
 302 long double
 303 __ieee754_ynl (n, x)
 304      int n;
 305      long double x;
 306 #endif
 307 {
 308   u_int32_t se;
 309   int32_t i, ix;
 310   int32_t sign;
 311   long double a, b, temp;
 312   ieee854_long_double_shape_type u;
 313
 314   u.value = x;
 315   se = u.parts32.w0;
 316   ix = se & 0x7fffffff;
 317
 318   /* if Y(n,NaN) is NaN */
 319   if (ix >= 0x7ff00000)
 320     {
 321       if ((u.parts32.w0 & 0xfffff) | u.parts32.w1
 322           | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3)
 323         return x + x;
 324     }
 325   if (x <= 0.0L)
 326     {
 327       if (x == 0.0L)
 328         return -HUGE_VALL + x;
 329       if (se & 0x80000000)
 330         return zero / (zero * x);
 331     }
 332   sign = 1;
 333   if (n < 0)
 334     {
 335       n = -n;
 336       sign = 1 - ((n & 1) << 1);
 337     }
 338   if (n == 0)
 339     return (__ieee754_y0l (x));
 340   if (n == 1)
 341     return (sign * __ieee754_y1l (x));
 342   if (ix >= 0x7ff00000)
 343     return zero;
 344   if (ix >= 0x52D00000)
 345     {                           /* x > 2**302 */
 346
 347       /* ??? See comment above on the possible futility of this.  */
 348
 349       /* (x >> n**2)
 350        *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 351        *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 352        *      Let s=sin(x), c=cos(x),
 353        *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 354        *
 355        *             n    sin(xn)*sqt2    cos(xn)*sqt2
 356        *          ----------------------------------
 357        *             0     s-c             c+s
 358        *             1    -s-c            -c+s
 359        *             2    -s+c            -c-s
 360        *             3     s+c             c-s
 361        */
 362       long double s;
 363       long double c;
 364       __sincosl (x, &s, &c);
 365       switch (n & 3)
 366         {
 367         case 0:
 368           temp = s - c;
 369           break;
 370         case 1:
 371           temp = -s - c;
 372           break;
 373         case 2:
 374           temp = -s + c;
 375           break;
 376         case 3:
 377           temp = s + c;
 378           break;
 379         }
 380       b = invsqrtpi * temp / __ieee754_sqrtl (x);
 381     }
 382   else
 383     {
 384       a = __ieee754_y0l (x);
 385       b = __ieee754_y1l (x);
 386       /* quit if b is -inf */
 387       u.value = b;
 388       se = u.parts32.w0 & 0xfff00000;
 389       for (i = 1; i < n && se != 0xfff00000; i++)
 390         {
 391           temp = b;
 392           b = ((long double) (i + i) / x) * b - a;
 393           u.value = b;
 394           se = u.parts32.w0 & 0xfff00000;
 395           a = temp;
 396         }
 397     }
 398   if (sign > 0)
 399     return b;
 400   else
 401     return -b;
 402 }