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