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, 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 "quadmath-imp.h"
  60
  61 static const __float128
  62   invsqrtpi = 5.6418958354775628694807945156077258584405E-1Q,
  63   two = 2,
  64   one = 1,
  65   zero = 0;
  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, ret;
  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   {
 107     SET_RESTORE_ROUNDF128 (FE_TONEAREST);
 108     if (x == 0 || ix >= 0x7fff0000)     /* if x is 0 or inf */
 109       return sgn == 1 ? -zero : zero;
 110     else if ((__float128) n <= x)
 111       {
 112         /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
 113         if (ix >= 0x412D0000)
 114           {                     /* x > 2**302 */
 115
 116             /* ??? Could use an expansion for large x here.  */
 117
 118             /* (x >> n**2)
 119              *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 120              *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 121              *      Let s=sin(x), c=cos(x),
 122              *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 123              *
 124              *             n    sin(xn)*sqt2    cos(xn)*sqt2
 125              *          ----------------------------------
 126              *             0     s-c             c+s
 127              *             1    -s-c            -c+s
 128              *             2    -s+c            -c-s
 129              *             3     s+c             c-s
 130              */
 131             __float128 s;
 132             __float128 c;
 133             sincosq (x, &s, &c);
 134             switch (n & 3)
 135               {
 136               case 0:
 137                 temp = c + s;
 138                 break;
 139               case 1:
 140                 temp = -c + s;
 141                 break;
 142               case 2:
 143                 temp = -c - s;
 144                 break;
 145               case 3:
 146                 temp = c - s;
 147                 break;
 148               }
 149             b = invsqrtpi * temp / sqrtq (x);
 150           }
 151         else
 152           {
 153             a = j0q (x);
 154             b = j1q (x);
 155             for (i = 1; i < n; i++)
 156               {
 157                 temp = b;
 158                 b = b * ((__float128) (i + i) / x) - a; /* avoid underflow */
 159                 a = temp;
 160               }
 161           }
 162       }
 163     else
 164       {
 165         if (ix < 0x3fc60000)
 166           {                     /* x < 2**-57 */
 167             /* x is tiny, return the first Taylor expansion of J(n,x)
 168              * J(n,x) = 1/n!*(x/2)^n  - ...
 169              */
 170             if (n >= 400)               /* underflow, result < 10^-4952 */
 171               b = zero;
 172             else
 173               {
 174                 temp = x * 0.5;
 175                 b = temp;
 176                 for (a = one, i = 2; i <= n; i++)
 177                   {
 178                     a *= (__float128) i;        /* a = n! */
 179                     b *= temp;  /* b = (x/2)^n */
 180                   }
 181                 b = b / a;
 182               }
 183           }
 184         else
 185           {
 186             /* use backward recurrence */
 187             /*                      x      x^2      x^2
 188              *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
 189              *                      2n  - 2(n+1) - 2(n+2)
 190              *
 191              *                      1      1        1
 192              *  (for large x)   =  ----  ------   ------   .....
 193              *                      2n   2(n+1)   2(n+2)
 194              *                      -- - ------ - ------ -
 195              *                       x     x         x
 196              *
 197              * Let w = 2n/x and h=2/x, then the above quotient
 198              * is equal to the continued fraction:
 199              *                  1
 200              *      = -----------------------
 201              *                     1
 202              *         w - -----------------
 203              *                        1
 204              *              w+h - ---------
 205              *                     w+2h - ...
 206              *
 207              * To determine how many terms needed, let
 208              * Q(0) = w, Q(1) = w(w+h) - 1,
 209              * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
 210              * When Q(k) > 1e4      good for single
 211              * When Q(k) > 1e9      good for double
 212              * When Q(k) > 1e17     good for quadruple
 213              */
 214             /* determine k */
 215             __float128 t, v;
 216             __float128 q0, q1, h, tmp;
 217             int32_t k, m;
 218             w = (n + n) / (__float128) x;
 219             h = 2 / (__float128) x;
 220             q0 = w;
 221             z = w + h;
 222             q1 = w * z - 1;
 223             k = 1;
 224             while (q1 < 1.0e17Q)
 225               {
 226                 k += 1;
 227                 z += h;
 228                 tmp = z * q1 - q0;
 229                 q0 = q1;
 230                 q1 = tmp;
 231               }
 232             m = n + n;
 233             for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
 234               t = one / (i / x - t);
 235             a = t;
 236             b = one;
 237             /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
 238              *  Hence, if n*(log(2n/x)) > ...
 239              *  single 8.8722839355e+01
 240              *  double 7.09782712893383973096e+02
 241              *  long double 1.1356523406294143949491931077970765006170e+04
 242              *  then recurrent value may overflow and the result is
 243              *  likely underflow to zero
 244              */
 245             tmp = n;
 246             v = two / x;
 247             tmp = tmp * logq (fabsq (v * tmp));
 248
 249             if (tmp < 1.1356523406294143949491931077970765006170e+04Q)
 250               {
 251                 for (i = n - 1, di = (__float128) (i + i); i > 0; i--)
 252                   {
 253                     temp = b;
 254                     b *= di;
 255                     b = b / x - a;
 256                     a = temp;
 257                     di -= two;
 258                   }
 259               }
 260             else
 261               {
 262                 for (i = n - 1, di = (__float128) (i + i); i > 0; i--)
 263                   {
 264                     temp = b;
 265                     b *= di;
 266                     b = b / x - a;
 267                     a = temp;
 268                     di -= two;
 269                     /* scale b to avoid spurious overflow */
 270                     if (b > 1e100Q)
 271                       {
 272                         a /= b;
 273                         t /= b;
 274                         b = one;
 275                       }
 276                   }
 277               }
 278             /* j0() and j1() suffer enormous loss of precision at and
 279              * near zero; however, we know that their zero points never
 280              * coincide, so just choose the one further away from zero.
 281              */
 282             z = j0q (x);
 283             w = j1q (x);
 284             if (fabsq (z) >= fabsq (w))
 285               b = (t * z / b);
 286             else
 287               b = (t * w / a);
 288           }
 289       }
 290     if (sgn == 1)
 291       ret = -b;
 292     else
 293       ret = b;
 294   }
 295   if (ret == 0)
 296     {
 297       ret = copysignq (FLT128_MIN, ret) * FLT128_MIN;
 298       errno = ERANGE;
 299     }
 300   else
 301     math_check_force_underflow (ret);
 302   return ret;
 303 }
 304
 305
 306 __float128
 307 ynq (int n, __float128 x)
 308 {
 309   uint32_t se;
 310   int32_t i, ix;
 311   int32_t sign;
 312   __float128 a, b, temp, ret;
 313   ieee854_float128 u;
 314
 315   u.value = x;
 316   se = u.words32.w0;
 317   ix = se & 0x7fffffff;
 318
 319   /* if Y(n,NaN) is NaN */
 320   if (ix >= 0x7fff0000)
 321     {
 322       if ((u.words32.w0 & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3)
 323         return x + x;
 324     }
 325   if (x <= 0)
 326     {
 327       if (x == 0)
 328         return ((n < 0 && (n & 1) != 0) ? 1 : -1) / 0.0Q;
 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 (y0q (x));
 340   {
 341     SET_RESTORE_ROUNDF128 (FE_TONEAREST);
 342     if (n == 1)
 343       {
 344         ret = sign * y1q (x);
 345         goto out;
 346       }
 347     if (ix >= 0x7fff0000)
 348       return zero;
 349     if (ix >= 0x412D0000)
 350       {                         /* x > 2**302 */
 351
 352         /* ??? See comment above on the possible futility of this.  */
 353
 354         /* (x >> n**2)
 355          *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 356          *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 357          *      Let s=sin(x), c=cos(x),
 358          *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 359          *
 360          *             n    sin(xn)*sqt2    cos(xn)*sqt2
 361          *          ----------------------------------
 362          *             0     s-c             c+s
 363          *             1    -s-c            -c+s
 364          *             2    -s+c            -c-s
 365          *             3     s+c             c-s
 366          */
 367         __float128 s;
 368         __float128 c;
 369         sincosq (x, &s, &c);
 370         switch (n & 3)
 371           {
 372           case 0:
 373             temp = s - c;
 374             break;
 375           case 1:
 376             temp = -s - c;
 377             break;
 378           case 2:
 379             temp = -s + c;
 380             break;
 381           case 3:
 382             temp = s + c;
 383             break;
 384           }
 385         b = invsqrtpi * temp / sqrtq (x);
 386       }
 387     else
 388       {
 389         a = y0q (x);
 390         b = y1q (x);
 391         /* quit if b is -inf */
 392         u.value = b;
 393         se = u.words32.w0 & 0xffff0000;
 394         for (i = 1; i < n && se != 0xffff0000; i++)
 395           {
 396             temp = b;
 397             b = ((__float128) (i + i) / x) * b - a;
 398             u.value = b;
 399             se = u.words32.w0 & 0xffff0000;
 400             a = temp;
 401           }
 402       }
 403     /* If B is +-Inf, set up errno accordingly.  */
 404     if (! finiteq (b))
 405       errno = ERANGE;
 406     if (sign > 0)
 407       ret = b;
 408     else
 409       ret = -b;
 410   }
 411  out:
 412   if (isinfq (ret))
 413     ret = copysignq (FLT128_MAX, ret) * FLT128_MAX;
 414   return ret;
 415 }