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