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 <errno.h>
  60 #include "quadmath-imp.h"
  61
  62 static const __float128
  63   invsqrtpi = 5.6418958354775628694807945156077258584405E-1Q,
  64   two = 2.0e0Q,
  65   one = 1.0e0Q,
  66   zero = 0.0Q;
  67
  68
  69 __float128
  70 jnq (int n, __float128 x)
  71 {
  72   uint32_t se;
  73   int32_t i, ix, sgn;
  74   __float128 a, b, temp, di;
  75   __float128 z, w;
  76   ieee854_float128 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.words32.w0;
  85   ix = se & 0x7fffffff;
  86
  87   /* if J(n,NaN) is NaN */
  88   if (ix >= 0x7fff0000)
  89     {
  90       if ((u.words32.w0 & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.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 (j0q (x));
 102   if (n == 1)
 103     return (j1q (x));
 104   sgn = (n & 1) & (se >> 31);   /* even n -- 0, odd n -- sign(x) */
 105   x = fabsq (x);
 106
 107   if (x == 0.0Q || ix >= 0x7fff0000)    /* if x is 0 or inf */
 108     b = zero;
 109   else if ((__float128) 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           __float128 s;
 131           __float128 c;
 132           sincosq (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 / sqrtq (x);
 149         }
 150       else
 151         {
 152           a = j0q (x);
 153           b = j1q (x);
 154           for (i = 1; i < n; i++)
 155             {
 156               temp = b;
 157               b = b * ((__float128) (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 *= (__float128) 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           __float128 t, v;
 215           __float128 q0, q1, h, tmp;
 216           int32_t k, m;
 217           w = (n + n) / (__float128) x;
 218           h = 2.0Q / (__float128) x;
 219           q0 = w;
 220           z = w + h;
 221           q1 = w * z - 1.0Q;
 222           k = 1;
 223           while (q1 < 1.0e17Q)
 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            *  __float128 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 * logq (fabsq (v * tmp));
 247
 248           if (tmp < 1.1356523406294143949491931077970765006170e+04Q)
 249             {
 250               for (i = n - 1, di = (__float128) (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 = (__float128) (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 > 1e100Q)
 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 = j0q (x);
 282           w = j1q (x);
 283           if (fabsq (z) >= fabsq (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
 295 __float128
 296 ynq (int n, __float128 x)
 297 {
 298   uint32_t se;
 299   int32_t i, ix;
 300   int32_t sign;
 301   __float128 a, b, temp;
 302   ieee854_float128 u;
 303
 304   u.value = x;
 305   se = u.words32.w0;
 306   ix = se & 0x7fffffff;
 307
 308   /* if Y(n,NaN) is NaN */
 309   if (ix >= 0x7fff0000)
 310     {
 311       if ((u.words32.w0 & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3)
 312         return x + x;
 313     }
 314   if (x <= 0.0Q)
 315     {
 316       if (x == 0.0Q)
 317         return -HUGE_VALQ + x;
 318       if (se & 0x80000000)
 319         return zero / (zero * x);
 320     }
 321   sign = 1;
 322   if (n < 0)
 323     {
 324       n = -n;
 325       sign = 1 - ((n & 1) << 1);
 326     }
 327   if (n == 0)
 328     return (y0q (x));
 329   if (n == 1)
 330     return (sign * y1q (x));
 331   if (ix >= 0x7fff0000)
 332     return zero;
 333   if (ix >= 0x412D0000)
 334     {                           /* x > 2**302 */
 335
 336       /* ??? See comment above on the possible futility of this.  */
 337
 338       /* (x >> n**2)
 339        *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 340        *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 341        *      Let s=sin(x), c=cos(x),
 342        *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 343        *
 344        *             n    sin(xn)*sqt2    cos(xn)*sqt2
 345        *          ----------------------------------
 346        *             0     s-c             c+s
 347        *             1    -s-c            -c+s
 348        *             2    -s+c            -c-s
 349        *             3     s+c             c-s
 350        */
 351       __float128 s;
 352       __float128 c;
 353       sincosq (x, &s, &c);
 354       switch (n & 3)
 355         {
 356         case 0:
 357           temp = s - c;
 358           break;
 359         case 1:
 360           temp = -s - c;
 361           break;
 362         case 2:
 363           temp = -s + c;
 364           break;
 365         case 3:
 366           temp = s + c;
 367           break;
 368         }
 369       b = invsqrtpi * temp / sqrtq (x);
 370     }
 371   else
 372     {
 373       a = y0q (x);
 374       b = y1q (x);
 375       /* quit if b is -inf */
 376       u.value = b;
 377       se = u.words32.w0 & 0xffff0000;
 378       for (i = 1; i < n && se != 0xffff0000; i++)
 379         {
 380           temp = b;
 381           b = ((__float128) (i + i) / x) * b - a;
 382           u.value = b;
 383           se = u.words32.w0 & 0xffff0000;
 384           a = temp;
 385         }
 386     }
 387   /* If B is +-Inf, set up errno accordingly.  */
 388   if (! finiteq (b))
 389     errno = ERANGE;
 390   if (sign > 0)
 391     return b;
 392   else
 393     return -b;
 394 }