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