nsz Git - musl/blob - src/math/jnf.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/e_jnf.c */
   2 /*
   3  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   4  */
   5 /*
   6  * ====================================================
   7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   8  *
   9  * Developed at SunPro, a Sun Microsystems, Inc. business.
  10  * Permission to use, copy, modify, and distribute this
  11  * software is freely granted, provided that this notice
  12  * is preserved.
  13  * ====================================================
  14  */
  15
  16 #define _GNU_SOURCE
  17 #include "libm.h"
  18
  19 float jnf(int n, float x)
  20 {
  21         int32_t i,hx,ix, sgn;
  22         float a, b, temp, di;
  23         float z, w;
  24
  25         /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
  26          * Thus, J(-n,x) = J(n,-x)
  27          */
  28         GET_FLOAT_WORD(hx, x);
  29         ix = 0x7fffffff & hx;
  30         /* if J(n,NaN) is NaN */
  31         if (ix > 0x7f800000)
  32                 return x+x;
  33         if (n < 0) {
  34                 n = -n;
  35                 x = -x;
  36                 hx ^= 0x80000000;
  37         }
  38         if (n == 0) return j0f(x);
  39         if (n == 1) return j1f(x);
  40
  41         sgn = (n&1)&(hx>>31);  /* even n -- 0, odd n -- sign(x) */
  42         x = fabsf(x);
  43         if (ix == 0 || ix >= 0x7f800000)  /* if x is 0 or inf */
  44                 b = 0.0f;
  45         else if((float)n <= x) {
  46                 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
  47                 a = j0f(x);
  48                 b = j1f(x);
  49                 for (i=1; i<n; i++){
  50                         temp = b;
  51                         b = b*((float)(i+i)/x) - a; /* avoid underflow */
  52                         a = temp;
  53                 }
  54         } else {
  55                 if (ix < 0x30800000) { /* x < 2**-29 */
  56                         /* x is tiny, return the first Taylor expansion of J(n,x)
  57                          * J(n,x) = 1/n!*(x/2)^n  - ...
  58                          */
  59                         if (n > 33)  /* underflow */
  60                                 b = 0.0f;
  61                         else {
  62                                 temp = 0.5f * x;
  63                                 b = temp;
  64                                 for (a=1.0f,i=2; i<=n; i++) {
  65                                         a *= (float)i;    /* a = n! */
  66                                         b *= temp;        /* b = (x/2)^n */
  67                                 }
  68                                 b = b/a;
  69                         }
  70                 } else {
  71                         /* use backward recurrence */
  72                         /*                      x      x^2      x^2
  73                          *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
  74                          *                      2n  - 2(n+1) - 2(n+2)
  75                          *
  76                          *                      1      1        1
  77                          *  (for large x)   =  ----  ------   ------   .....
  78                          *                      2n   2(n+1)   2(n+2)
  79                          *                      -- - ------ - ------ -
  80                          *                       x     x         x
  81                          *
  82                          * Let w = 2n/x and h=2/x, then the above quotient
  83                          * is equal to the continued fraction:
  84                          *                  1
  85                          *      = -----------------------
  86                          *                     1
  87                          *         w - -----------------
  88                          *                        1
  89                          *              w+h - ---------
  90                          *                     w+2h - ...
  91                          *
  92                          * To determine how many terms needed, let
  93                          * Q(0) = w, Q(1) = w(w+h) - 1,
  94                          * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
  95                          * When Q(k) > 1e4      good for single
  96                          * When Q(k) > 1e9      good for double
  97                          * When Q(k) > 1e17     good for quadruple
  98                          */
  99                         /* determine k */
 100                         float t,v;
 101                         float q0,q1,h,tmp;
 102                         int32_t k,m;
 103
 104                         w = (n+n)/x;
 105                         h = 2.0f/x;
 106                         z = w+h;
 107                         q0 = w;
 108                         q1 = w*z - 1.0f;
 109                         k = 1;
 110                         while (q1 < 1.0e9f) {
 111                                 k += 1;
 112                                 z += h;
 113                                 tmp = z*q1 - q0;
 114                                 q0 = q1;
 115                                 q1 = tmp;
 116                         }
 117                         m = n+n;
 118                         for (t=0.0f, i = 2*(n+k); i>=m; i -= 2)
 119                                 t = 1.0f/(i/x-t);
 120                         a = t;
 121                         b = 1.0f;
 122                         /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
 123                          *  Hence, if n*(log(2n/x)) > ...
 124                          *  single 8.8722839355e+01
 125                          *  double 7.09782712893383973096e+02
 126                          *  long double 1.1356523406294143949491931077970765006170e+04
 127                          *  then recurrent value may overflow and the result is
 128                          *  likely underflow to zero
 129                          */
 130                         tmp = n;
 131                         v = 2.0f/x;
 132                         tmp = tmp*logf(fabsf(v*tmp));
 133                         if (tmp < 88.721679688f) {
 134                                 for (i=n-1,di=(float)(i+i); i>0; i--) {
 135                                         temp = b;
 136                                         b *= di;
 137                                         b = b/x - a;
 138                                         a = temp;
 139                                         di -= 2.0f;
 140                                 }
 141                         } else {
 142                                 for (i=n-1,di=(float)(i+i); i>0; i--){
 143                                         temp = b;
 144                                         b *= di;
 145                                         b = b/x - a;
 146                                         a = temp;
 147                                         di -= 2.0f;
 148                                         /* scale b to avoid spurious overflow */
 149                                         if (b > 1e10f) {
 150                                                 a /= b;
 151                                                 t /= b;
 152                                                 b = 1.0f;
 153                                         }
 154                                 }
 155                         }
 156                         z = j0f(x);
 157                         w = j1f(x);
 158                         if (fabsf(z) >= fabsf(w))
 159                                 b = t*z/b;
 160                         else
 161                                 b = t*w/a;
 162                 }
 163         }
 164         if (sgn == 1) return -b;
 165         return b;
 166 }
 167
 168 float ynf(int n, float x)
 169 {
 170         int32_t i,hx,ix,ib;
 171         int32_t sign;
 172         float a, b, temp;
 173
 174         GET_FLOAT_WORD(hx, x);
 175         ix = 0x7fffffff & hx;
 176         /* if Y(n,NaN) is NaN */
 177         if (ix > 0x7f800000)
 178                 return x+x;
 179         if (ix == 0)
 180                 return -1.0f/0.0f;
 181         if (hx < 0)
 182                 return 0.0f/0.0f;
 183         sign = 1;
 184         if (n < 0) {
 185                 n = -n;
 186                 sign = 1 - ((n&1)<<1);
 187         }
 188         if (n == 0)
 189                 return y0f(x);
 190         if (n == 1)
 191                 return sign*y1f(x);
 192         if (ix == 0x7f800000)
 193                 return 0.0f;
 194
 195         a = y0f(x);
 196         b = y1f(x);
 197         /* quit if b is -inf */
 198         GET_FLOAT_WORD(ib,b);
 199         for (i = 1; i < n && ib != 0xff800000; i++){
 200                 temp = b;
 201                 b = ((float)(i+i)/x)*b - a;
 202                 GET_FLOAT_WORD(ib, b);
 203                 a = temp;
 204         }
 205         if (sign > 0)
 206                 return b;
 207         return -b;
 208 }