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