Print this page
11175 libm should use signbit() correctly
11188 c99 math macros should return strictly backward compatible values

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/Q/jnl.c
          +++ new/usr/src/lib/libm/common/Q/jnl.c
↓ open down ↓ 60 lines elided ↑ open up ↑
  61   61  
  62   62  #define GENERIC long double
  63   63  
  64   64  static const GENERIC
  65   65  invsqrtpi = 5.641895835477562869480794515607725858441e-0001L,
  66   66  two  = 2.0L,
  67   67  zero = 0.0L,
  68   68  one  = 1.0L;
  69   69  
  70   70  GENERIC
  71      -jnl(n, x) int n; GENERIC x; {
       71 +jnl(int n, GENERIC x)
       72 +{
  72   73          int i, sgn;
  73   74          GENERIC a, b, temp, z, w;
  74   75  
  75   76          /*
  76   77           * J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
  77   78           * Thus, J(-n,x) = J(n,-x)
  78   79           */
  79   80          if (n < 0) {
  80   81                  n = -n;
  81   82                  x = -x;
  82   83          }
  83   84          if (n == 0)
  84   85                  return (j0l(x));
  85   86          if (n == 1)
  86   87                  return (j1l(x));
  87   88          if (x != x)
  88   89                  return (x+x);
  89   90          if ((n&1) == 0)
  90      -                sgn = 0;                        /* even n */
       91 +                sgn = 0;                        /* even n */
  91   92          else
  92   93                  sgn = signbitl(x);      /* old n  */
  93   94          x = fabsl(x);
  94   95          if (x == zero || !finitel(x)) b = zero;
  95   96          else if ((GENERIC)n <= x) {
  96   97                                          /*
  97   98                                           * Safe to use
  98   99                                           * J(n+1,x)=2n/x *J(n,x)-J(n-1,x)
  99  100                                           */
 100      -            if (x > 1.0e91L) {
      101 +                if (x > 1.0e91L) {
 101  102                                  /*
 102  103                                   * x >> n**2
 103  104                                   *  Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 104  105                                   *   Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 105  106                                   *   Let s=sin(x), c=cos(x),
 106  107                                   *      xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 107  108                                   *
 108  109                                   *         n    sin(xn)*sqt2    cos(xn)*sqt2
 109  110                                   *      ----------------------------------
 110  111                                   *         0     s-c             c+s
 111      -                                 *         1    -s-c            -c+s
      112 +                                 *         1    -s-c            -c+s
 112  113                                   *         2    -s+c            -c-s
 113  114                                   *         3     s+c             c-s
 114  115                                   */
 115      -                switch (n&3) {
 116      -                    case 0: temp =  cosl(x)+sinl(x); break;
 117      -                    case 1: temp = -cosl(x)+sinl(x); break;
 118      -                    case 2: temp = -cosl(x)-sinl(x); break;
 119      -                    case 3: temp =  cosl(x)-sinl(x); break;
 120      -                }
 121      -                b = invsqrtpi*temp/sqrtl(x);
 122      -            } else {
      116 +                        switch (n&3) {
      117 +                        case 0:
      118 +                                temp =  cosl(x)+sinl(x);
      119 +                                break;
      120 +                        case 1:
      121 +                                temp = -cosl(x)+sinl(x);
      122 +                                break;
      123 +                        case 2:
      124 +                                temp = -cosl(x)-sinl(x);
      125 +                                break;
      126 +                        case 3:
      127 +                                temp =  cosl(x)-sinl(x);
      128 +                                break;
      129 +                        }
      130 +                        b = invsqrtpi*temp/sqrtl(x);
      131 +                } else {
 123  132                          a = j0l(x);
 124  133                          b = j1l(x);
 125  134                          for (i = 1; i < n; i++) {
 126      -                    temp = b;
 127      -                    b = b*((GENERIC)(i+i)/x) - a; /* avoid underflow */
 128      -                    a = temp;
      135 +                                temp = b;
      136 +                                /* avoid underflow */
      137 +                                b = b*((GENERIC)(i+i)/x) - a;
      138 +                                a = temp;
 129  139                          }
 130      -            }
 131      -        } else {
 132      -            if (x < 1e-17L) {   /* use J(n,x) = 1/n!*(x/2)^n */
 133      -                b = powl(0.5L*x, (GENERIC)n);
 134      -                if (b != zero) {
 135      -                    for (a = one, i = 1; i <= n; i++) a *= (GENERIC)i;
 136      -                    b = b/a;
 137      -                }
 138      -            } else {
 139      -                /* use backward recurrence */
 140      -                /*
 141      -                 *                      x      x^2      x^2
 142      -                 *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
 143      -                 *                      2n  - 2(n+1) - 2(n+2)
 144      -                 *
 145      -                 *                      1      1        1
 146      -                 *  (for large x)   =  ----  ------   ------   .....
 147      -                 *                      2n   2(n+1)   2(n+2)
 148      -                 *                      -- - ------ - ------ -
 149      -                 *                       x     x         x
 150      -                 *
 151      -                 * Let w = 2n/x and h=2/x, then the above quotient
 152      -                 * is equal to the continued fraction:
 153      -                 *                  1
 154      -                 *      = -----------------------
 155      -                 *                     1
 156      -                 *         w - -----------------
 157      -                 *                        1
 158      -                 *              w+h - ---------
 159      -                 *                     w+2h - ...
 160      -                 *
 161      -                 * To determine how many terms needed, let
 162      -                 * Q(0) = w, Q(1) = w(w+h) - 1,
 163      -                 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
 164      -                 * When Q(k) > 1e4      good for single
 165      -                 * When Q(k) > 1e9      good for double
 166      -                 * When Q(k) > 1e17     good for quaduple
 167      -                 */
 168      -            /* determin k */
 169      -                GENERIC t, v;
 170      -                double q0, q1, h, tmp; int k, m;
 171      -                w  = (n+n)/(double)x; h = 2.0/(double)x;
 172      -                q0 = w;  z = w+h; q1 = w*z - 1.0; k = 1;
 173      -                while (q1 < 1.0e17) {
 174      -                        k += 1; z += h;
 175      -                        tmp = z*q1 - q0;
 176      -                        q0 = q1;
 177      -                        q1 = tmp;
 178  140                  }
 179      -                m = n+n;
 180      -                for (t = zero, i = 2*(n+k); i >= m; i -= 2) t = one/(i/x-t);
 181      -                a = t;
 182      -                b = one;
 183      -                /*
 184      -                 * estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
 185      -                 * hence, if n*(log(2n/x)) > ...
 186      -                 *  single 8.8722839355e+01
 187      -                 *  double 7.09782712893383973096e+02
 188      -                 *  long double 1.1356523406294143949491931077970765006170e+04
 189      -                 *  then recurrent value may overflow and the result is
 190      -                 *  likely underflow to zero
 191      -                 */
 192      -                tmp = n;
 193      -                v = two/x;
 194      -                tmp = tmp*logl(fabsl(v*tmp));
 195      -                if (tmp < 1.1356523406294143949491931077970765e+04L) {
 196      -                                for (i = n-1; i > 0; i--) {
 197      -                        temp = b;
 198      -                        b = ((i+i)/x)*b - a;
 199      -                        a = temp;
 200      -                    }
      141 +        } else {
      142 +                if (x < 1e-17L) {       /* use J(n,x) = 1/n!*(x/2)^n */
      143 +                        b = powl(0.5L*x, (GENERIC)n);
      144 +                        if (b != zero) {
      145 +                                for (a = one, i = 1; i <= n; i++)
      146 +                                        a *= (GENERIC)i;
      147 +                                b = b/a;
      148 +                        }
 201  149                  } else {
      150 +                        /* use backward recurrence */
      151 +                        /*
      152 +                         *                      x      x^2      x^2
      153 +                         *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
      154 +                         *                      2n  - 2(n+1) - 2(n+2)
      155 +                         *
      156 +                         *                      1      1        1
      157 +                         *  (for large x)   =  ----  ------   ------   .....
      158 +                         *                      2n   2(n+1)   2(n+2)
      159 +                         *                      -- - ------ - ------ -
      160 +                         *                       x     x         x
      161 +                         *
      162 +                         * Let w = 2n/x and h=2/x, then the above quotient
      163 +                         * is equal to the continued fraction:
      164 +                         *                  1
      165 +                         *      = -----------------------
      166 +                         *                     1
      167 +                         *         w - -----------------
      168 +                         *                        1
      169 +                         *              w+h - ---------
      170 +                         *                     w+2h - ...
      171 +                         *
      172 +                         * To determine how many terms needed, let
      173 +                         * Q(0) = w, Q(1) = w(w+h) - 1,
      174 +                         * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
      175 +                         * When Q(k) > 1e4      good for single
      176 +                         * When Q(k) > 1e9      good for double
      177 +                         * When Q(k) > 1e17     good for quaduple
      178 +                         */
      179 +                        /* determine k */
      180 +                        GENERIC t, v;
      181 +                        double q0, q1, h, tmp;
      182 +                        int k, m;
      183 +                        w  = (n+n)/(double)x;
      184 +                        h = 2.0/(double)x;
      185 +                        q0 = w;
      186 +                        z = w+h;
      187 +                        q1 = w*z - 1.0;
      188 +                        k = 1;
      189 +                        while (q1 < 1.0e17) {
      190 +                                k += 1;
      191 +                                z += h;
      192 +                                tmp = z*q1 - q0;
      193 +                                q0 = q1;
      194 +                                q1 = tmp;
      195 +                        }
      196 +                        m = n+n;
      197 +                        for (t = zero, i = 2*(n+k); i >= m; i -= 2)
      198 +                                t = one/(i/x-t);
      199 +                        a = t;
      200 +                        b = one;
      201 +                        /*
      202 +                         * estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
      203 +                         * hence, if n*(log(2n/x)) > ...
      204 +                         *  single 8.8722839355e+01
      205 +                         *  double 7.09782712893383973096e+02
      206 +                         *  long double 1.1356523406294143949491931077970765006170e+04
      207 +                         *  then recurrent value may overflow and the result is
      208 +                         *  likely underflow to zero
      209 +                         */
      210 +                        tmp = n;
      211 +                        v = two/x;
      212 +                        tmp = tmp*logl(fabsl(v*tmp));
      213 +                        if (tmp < 1.1356523406294143949491931077970765e+04L) {
      214 +                                for (i = n-1; i > 0; i--) {
      215 +                                        temp = b;
      216 +                                        b = ((i+i)/x)*b - a;
      217 +                                        a = temp;
      218 +                                }
      219 +                        } else {
 202  220                                  for (i = n-1; i > 0; i--) {
 203      -                        temp = b;
 204      -                        b = ((i+i)/x)*b - a;
 205      -                        a = temp;
 206      -                        if (b > 1e1000L) {
 207      -                            a /= b;
 208      -                            t /= b;
 209      -                            b  = 1.0;
 210      -                        }
 211      -                    }
      221 +                                        temp = b;
      222 +                                        b = ((i+i)/x)*b - a;
      223 +                                        a = temp;
      224 +                                        if (b > 1e1000L) {
      225 +                                                a /= b;
      226 +                                                t /= b;
      227 +                                                b  = 1.0;
      228 +                                        }
      229 +                                }
      230 +                        }
      231 +                        b = (t*j0l(x)/b);
 212  232                  }
 213      -                b = (t*j0l(x)/b);
 214      -            }
 215  233          }
 216      -        if (sgn == 1)
      234 +        if (sgn != 0)
 217  235                  return (-b);
 218  236          else
 219  237                  return (b);
 220  238  }
 221  239  
 222      -GENERIC ynl(n, x)
 223      -int n; GENERIC x; {
      240 +GENERIC
      241 +ynl(int n, GENERIC x)
      242 +{
 224  243          int i;
 225  244          int sign;
 226  245          GENERIC a, b, temp;
 227  246  
 228  247          if (x != x)
 229  248                  return (x+x);
 230  249          if (x <= zero) {
 231  250                  if (x == zero)
 232  251                          return (-one/zero);
 233  252                  else
↓ open down ↓ 4 lines elided ↑ open up ↑
 238  257                  n = -n;
 239  258                  if ((n&1) == 1) sign = -1;
 240  259          }
 241  260          if (n == 0)
 242  261                  return (y0l(x));
 243  262          if (n == 1)
 244  263                  return (sign*y1l(x));
 245  264          if (!finitel(x))
 246  265                  return (zero);
 247  266  
 248      -        if (x > 1.0e91L) {      /* x >> n**2
 249      -                                    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 250      -                                    Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 251      -                                    Let s = sin(x), c = cos(x),
 252      -                                        xn = x-(2n+1)*pi/4, sqt2 = sqrt(2), then
 253      -
 254      -                                           n    sin(xn)*sqt2    cos(xn)*sqt2
 255      -                                        ----------------------------------
 256      -                                           0     s-c             c+s
 257      -                                           1    -s-c            -c+s
 258      -                                           2    -s+c            -c-s
 259      -                                           3     s+c             c-s
 260      -                                 */
      267 +        if (x > 1.0e91L) {
      268 +                /*
      269 +                 * x >> n**2
      270 +                 *   Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
      271 +                 *   Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
      272 +                 *   Let s = sin(x), c = cos(x),
      273 +                 *      xn = x-(2n+1)*pi/4, sqt2 = sqrt(2), then
      274 +                 *
      275 +                 *         n    sin(xn)*sqt2    cos(xn)*sqt2
      276 +                 *      ----------------------------------
      277 +                 *         0     s-c             c+s
      278 +                 *         1    -s-c            -c+s
      279 +                 *         2    -s+c            -c-s
      280 +                 *         3     s+c             c-s
      281 +                 */
 261  282                  switch (n&3) {
 262      -                    case 0: temp =  sinl(x)-cosl(x); break;
 263      -                    case 1: temp = -sinl(x)-cosl(x); break;
 264      -                    case 2: temp = -sinl(x)+cosl(x); break;
 265      -                    case 3: temp =  sinl(x)+cosl(x); break;
      283 +                case 0:
      284 +                        temp =  sinl(x)-cosl(x);
      285 +                        break;
      286 +                case 1:
      287 +                        temp = -sinl(x)-cosl(x);
      288 +                        break;
      289 +                case 2:
      290 +                        temp = -sinl(x)+cosl(x);
      291 +                        break;
      292 +                case 3:
      293 +                        temp =  sinl(x)+cosl(x);
      294 +                        break;
 266  295                  }
 267  296                  b = invsqrtpi*temp/sqrtl(x);
 268  297          } else {
 269  298                  a = y0l(x);
 270  299                  b = y1l(x);
 271  300                  /*
 272  301                   * fix 1262058 and take care of non-default rounding
 273  302                   */
 274  303                  for (i = 1; i < n; i++) {
 275  304                          temp = b;
↓ open down ↓ 12 lines elided ↑ open up ↑
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX