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