Print this page
11210 libm should be cstyle(1ONBLD) clean

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/complex/casin.c
          +++ new/usr/src/lib/libm/common/complex/casin.c
↓ open down ↓ 14 lines elided ↑ open up ↑
  15   15   * If applicable, add the following below this CDDL HEADER, with the
  16   16   * fields enclosed by brackets "[]" replaced with your own identifying
  17   17   * information: Portions Copyright [yyyy] [name of copyright owner]
  18   18   *
  19   19   * CDDL HEADER END
  20   20   */
  21   21  
  22   22  /*
  23   23   * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24   24   */
       25 +
  25   26  /*
  26   27   * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27   28   * Use is subject to license terms.
  28   29   */
  29   30  
  30   31  #pragma weak __casin = casin
  31   32  
  32      -/* INDENT OFF */
       33 +
  33   34  /*
  34   35   * dcomplex casin(dcomplex z);
  35   36   *
  36   37   * Alogrithm
  37   38   * (based on T.E.Hull, Thomas F. Fairgrieve and Ping Tak Peter Tang's
  38   39   * paper "Implementing the Complex Arcsine and Arccosine Functins Using
  39   40   * Exception Handling", ACM TOMS, Vol 23, pp 299-335)
  40   41   *
  41   42   * The principal value of complex inverse sine function casin(z),
  42   43   * where z = x+iy, can be defined by
  43   44   *
  44      - *      casin(z) = asin(B) + i sign(y) log (A + sqrt(A*A-1)),
       45 + *      casin(z) = asin(B) + i sign(y) log (A + sqrt(A*A-1)),
  45   46   *
  46   47   * where the log function is the natural log, and
  47   48   *             ____________           ____________
  48   49   *       1    /     2    2      1    /     2    2
  49   50   *  A = ---  / (x+1)  + y   +  ---  / (x-1)  + y
  50   51   *       2 \/                   2 \/
  51   52   *             ____________           ____________
  52   53   *       1    /     2    2      1    /     2    2
  53   54   *  B = ---  / (x+1)  + y   -  ---  / (x-1)  + y   .
  54   55   *       2 \/                   2 \/
↓ open down ↓ 88 lines elided ↑ open up ↑
 143  144   *      Thus
 144  145   *                            [ asin(x) + i y/sqrt((x-1)*(x+1)), if x <  1
 145  146   *              casin(x+i*y)=[
 146  147   *                            [ pi/2    + i log(x+sqrt(x*x-1)),  if x >= 1
 147  148   *
 148  149   *  case 3. y < 4 sqrt(u), where u = minimum normal x.
 149  150   *      After case 1 and 2, this will only occurs when x=1. When x=1, we have
 150  151   *         A = (sqrt(4+y*y)+y)/2 ~ 1 + y/2 + y^2/8 + ...
 151  152   *      and
 152  153   *         B = 1/A = 1 - y/2 + y^2/8 + ...
 153      - *      Since
      154 + *      Since
 154  155   *         asin(x) = pi/2-2*asin(sqrt((1-x)/2))
 155  156   *         asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
 156  157   *      we have, for the real part asin(B),
 157  158   *         asin(1-y/2) ~ pi/2 - 2 asin(sqrt(y/4))
 158  159   *                     ~ pi/2 - sqrt(y)
 159  160   *      For the imaginary part,
 160  161   *         log(A+sqrt(A*A-1)) ~ log(1+y/2+sqrt(2*y/2))
 161  162   *                            = log(1+y/2+sqrt(y))
 162  163   *                            = (y/2+sqrt(y)) - (y/2+sqrt(y))^2/2 + ...
 163  164   *                            ~ sqrt(y) - y*(sqrt(y)+y/2)/2
 164  165   *                            ~ sqrt(y)
 165  166   *
 166  167   *  case 4. y >= (x+1)ulp(0.5). In this case, A ~ y and B ~ x/y. Thus
 167  168   *         real part = asin(B) ~ x/y (be careful, x/y may underflow)
 168      - *      and
      169 + *      and
 169  170   *         imag part = log(y+sqrt(y*y-one))
 170  171   *
 171  172   *
 172  173   *  case 5. Both x and y are large: x and y > sqrt(M)/8, where M = maximum x
 173  174   *      In this case,
 174  175   *         A ~ sqrt(x*x+y*y)
 175  176   *         B ~ x/sqrt(x*x+y*y).
 176  177   *      Thus
 177  178   *         real part = asin(B) = atan(x/y),
 178  179   *         imag part = log(A+sqrt(A*A-1)) ~ log(2A)
↓ open down ↓ 3 lines elided ↑ open up ↑
 182  183   *  case 6. x < 4 sqrt(u). In this case, we have
 183  184   *          A ~ sqrt(1+y*y), B = x/sqrt(1+y*y).
 184  185   *      Since B is tiny, we have
 185  186   *          real part = asin(B) ~ B = x/sqrt(1+y*y)
 186  187   *          imag part = log(A+sqrt(A*A-1)) = log (A+sqrt(y*y))
 187  188   *                    = log(y+sqrt(1+y*y))
 188  189   *                    = 0.5*log(y^2+2ysqrt(1+y^2)+1+y^2)
 189  190   *                    = 0.5*log(1+2y(y+sqrt(1+y^2)));
 190  191   *                    = 0.5*log1p(2y(y+A));
 191  192   *
 192      - *      casin(z) = asin(B) + i sign(y) log (A + sqrt(A*A-1)),
      193 + *      casin(z) = asin(B) + i sign(y) log (A + sqrt(A*A-1)),
 193  194   */
 194      -/* INDENT ON */
 195  195  
 196      -#include "libm.h"               /* asin/atan/fabs/log/log1p/sqrt */
      196 +#include "libm.h"                       /* asin/atan/fabs/log/log1p/sqrt */
 197  197  #include "complex_wrapper.h"
 198  198  
 199      -/* INDENT OFF */
 200      -static const double
 201      -        zero = 0.0,
      199 +static const double zero = 0.0,
 202  200          one = 1.0,
 203      -        E = 1.11022302462515654042e-16,                 /* 2**-53 */
      201 +        E = 1.11022302462515654042e-16, /* 2**-53 */
 204  202          ln2 = 6.93147180559945286227e-01,
 205  203          pi_2 = 1.570796326794896558e+00,
 206  204          pi_2_l = 6.123233995736765886e-17,
 207  205          pi_4 = 7.85398163397448278999e-01,
 208      -        Foursqrtu = 5.96667258496016539463e-154,        /* 2**(-509) */
      206 +        Foursqrtu = 5.96667258496016539463e-154, /* 2**(-509) */
 209  207          Acrossover = 1.5,
 210  208          Bcrossover = 0.6417,
 211  209          half = 0.5;
 212      -/* INDENT ON */
      210 +
 213  211  
 214  212  dcomplex
 215      -casin(dcomplex z) {
      213 +casin(dcomplex z)
      214 +{
 216  215          double x, y, t, R, S, A, Am1, B, y2, xm1, xp1, Apx;
 217  216          int ix, iy, hx, hy;
 218  217          unsigned lx, ly;
 219  218          dcomplex ans;
 220  219  
 221  220          x = D_RE(z);
 222  221          y = D_IM(z);
 223  222          hx = HI_WORD(x);
 224  223          lx = LO_WORD(x);
 225  224          hy = HI_WORD(y);
 226  225          ly = LO_WORD(y);
 227  226          ix = hx & 0x7fffffff;
 228  227          iy = hy & 0x7fffffff;
 229  228          x = fabs(x);
 230  229          y = fabs(y);
 231  230  
 232  231          /* special cases */
 233  232  
 234  233          /* x is inf or NaN */
 235      -        if (ix >= 0x7ff00000) { /* x is inf or NaN */
      234 +        if (ix >= 0x7ff00000) {         /* x is inf or NaN */
 236  235                  if (ISINF(ix, lx)) {    /* x is INF */
 237  236                          D_IM(ans) = x;
      237 +
 238  238                          if (iy >= 0x7ff00000) {
 239  239                                  if (ISINF(iy, ly))
 240  240                                          /* casin(inf + i inf) = pi/4 + i inf */
 241  241                                          D_RE(ans) = pi_4;
 242  242                                  else    /* casin(inf + i NaN) = NaN  + i inf  */
 243  243                                          D_RE(ans) = y + y;
 244      -                        } else  /* casin(inf + iy) = pi/2 + i inf */
      244 +                        } else { /* casin(inf + iy) = pi/2 + i inf */
 245  245                                  D_RE(ans) = pi_2;
      246 +                        }
 246  247                  } else {                /* x is NaN */
 247  248                          if (iy >= 0x7ff00000) {
 248      -                                /* INDENT OFF */
      249 +
 249  250                                  /*
 250  251                                   * casin(NaN + i inf) = NaN + i inf
 251  252                                   * casin(NaN + i NaN) = NaN + i NaN
 252  253                                   */
 253      -                                /* INDENT ON */
 254  254                                  D_IM(ans) = y + y;
 255  255                                  D_RE(ans) = x + x;
 256  256                          } else {
 257  257                                  /* casin(NaN + i y ) = NaN  + i NaN */
 258  258                                  D_IM(ans) = D_RE(ans) = x + y;
 259  259                          }
 260  260                  }
      261 +
 261  262                  if (hx < 0)
 262  263                          D_RE(ans) = -D_RE(ans);
      264 +
 263  265                  if (hy < 0)
 264  266                          D_IM(ans) = -D_IM(ans);
      267 +
 265  268                  return (ans);
 266  269          }
 267  270  
 268  271          /* casin(+0 + i 0  ) =  0   + i 0. */
 269  272          if ((ix | lx | iy | ly) == 0)
 270  273                  return (z);
 271  274  
 272      -        if (iy >= 0x7ff00000) { /* y is inf or NaN */
      275 +        if (iy >= 0x7ff00000) {         /* y is inf or NaN */
 273  276                  if (ISINF(iy, ly)) {    /* casin(x + i inf) =  0   + i inf */
 274  277                          D_IM(ans) = y;
 275  278                          D_RE(ans) = zero;
 276  279                  } else {                /* casin(x + i NaN) = NaN  + i NaN */
 277  280                          D_IM(ans) = x + y;
      281 +
 278  282                          if ((ix | lx) == 0)
 279  283                                  D_RE(ans) = x;
 280  284                          else
 281  285                                  D_RE(ans) = y;
 282  286                  }
      287 +
 283  288                  if (hx < 0)
 284  289                          D_RE(ans) = -D_RE(ans);
      290 +
 285  291                  if (hy < 0)
 286  292                          D_IM(ans) = -D_IM(ans);
      293 +
 287  294                  return (ans);
 288  295          }
 289  296  
 290      -        if ((iy | ly) == 0) {   /* region 1: y=0 */
      297 +        if ((iy | ly) == 0) {           /* region 1: y=0 */
 291  298                  if (ix < 0x3ff00000) {  /* |x| < 1 */
 292  299                          D_RE(ans) = asin(x);
 293  300                          D_IM(ans) = zero;
 294  301                  } else {
 295  302                          D_RE(ans) = pi_2;
 296      -                        if (ix >= 0x43500000)   /* |x| >= 2**54 */
      303 +
      304 +                        if (ix >= 0x43500000) {         /* |x| >= 2**54 */
 297  305                                  D_IM(ans) = ln2 + log(x);
 298      -                        else if (ix >= 0x3ff80000)      /* x > Acrossover */
      306 +                        } else if (ix >= 0x3ff80000) {  /* x > Acrossover */
 299  307                                  D_IM(ans) = log(x + sqrt((x - one) * (x +
 300      -                                        one)));
 301      -                        else {
      308 +                                    one)));
      309 +                        } else {
 302  310                                  xm1 = x - one;
 303  311                                  D_IM(ans) = log1p(xm1 + sqrt(xm1 * (x + one)));
 304  312                          }
 305  313                  }
 306  314          } else if (y <= E * fabs(x - one)) {    /* region 2: y < tiny*|x-1| */
 307      -                if (ix < 0x3ff00000) {  /* x < 1 */
      315 +                if (ix < 0x3ff00000) {          /* x < 1 */
 308  316                          D_RE(ans) = asin(x);
 309  317                          D_IM(ans) = y / sqrt((one + x) * (one - x));
 310  318                  } else {
 311  319                          D_RE(ans) = pi_2;
 312      -                        if (ix >= 0x43500000) { /* |x| >= 2**54 */
      320 +
      321 +                        if (ix >= 0x43500000)           /* |x| >= 2**54 */
 313  322                                  D_IM(ans) = ln2 + log(x);
 314      -                        } else if (ix >= 0x3ff80000)    /* x > Acrossover */
      323 +                        else if (ix >= 0x3ff80000)      /* x > Acrossover */
 315  324                                  D_IM(ans) = log(x + sqrt((x - one) * (x +
 316      -                                        one)));
      325 +                                    one)));
 317  326                          else
 318  327                                  D_IM(ans) = log1p((x - one) + sqrt((x - one) *
 319      -                                        (x + one)));
      328 +                                    (x + one)));
 320  329                  }
 321  330          } else if (y < Foursqrtu) {     /* region 3 */
 322  331                  t = sqrt(y);
 323  332                  D_RE(ans) = pi_2 - (t - pi_2_l);
 324  333                  D_IM(ans) = t;
 325      -        } else if (E * y - one >= x) {  /* region 4 */
 326      -                D_RE(ans) = x / y;      /* need to fix underflow cases */
      334 +        } else if (E * y - one >= x) { /* region 4 */
      335 +                D_RE(ans) = x / y; /* need to fix underflow cases */
 327  336                  D_IM(ans) = ln2 + log(y);
 328  337          } else if (ix >= 0x5fc00000 || iy >= 0x5fc00000) {      /* x,y>2**509 */
 329  338                  /* region 5: x+1 or y is very large (>= sqrt(max)/8) */
 330  339                  t = x / y;
 331  340                  D_RE(ans) = atan(t);
 332  341                  D_IM(ans) = ln2 + log(y) + half * log1p(t * t);
 333  342          } else if (x < Foursqrtu) {
 334  343                  /* region 6: x is very small, < 4sqrt(min) */
 335  344                  A = sqrt(one + y * y);
 336  345                  D_RE(ans) = x / A;      /* may underflow */
      346 +
 337  347                  if (iy >= 0x3ff80000)   /* if y > Acrossover */
 338  348                          D_IM(ans) = log(y + A);
 339  349                  else
 340  350                          D_IM(ans) = half * log1p((y + y) * (y + A));
 341      -        } else {        /* safe region */
      351 +        } else {                        /* safe region */
 342  352                  y2 = y * y;
 343  353                  xp1 = x + one;
 344  354                  xm1 = x - one;
 345  355                  R = sqrt(xp1 * xp1 + y2);
 346  356                  S = sqrt(xm1 * xm1 + y2);
 347  357                  A = half * (R + S);
 348  358                  B = x / A;
 349  359  
 350      -                if (B <= Bcrossover)
      360 +                if (B <= Bcrossover) {
 351  361                          D_RE(ans) = asin(B);
 352      -                else {          /* use atan and an accurate approx to a-x */
      362 +                } else {        /* use atan and an accurate approx to a-x */
 353  363                          Apx = A + x;
      364 +
 354  365                          if (x <= one)
 355  366                                  D_RE(ans) = atan(x / sqrt(half * Apx * (y2 /
 356      -                                        (R + xp1) + (S - xm1))));
      367 +                                    (R + xp1) + (S - xm1))));
 357  368                          else
 358  369                                  D_RE(ans) = atan(x / (y * sqrt(half * (Apx /
 359      -                                        (R + xp1) + Apx / (S + xm1)))));
      370 +                                    (R + xp1) + Apx / (S + xm1)))));
 360  371                  }
      372 +
 361  373                  if (A <= Acrossover) {
 362  374                          /* use log1p and an accurate approx to A-1 */
 363  375                          if (x < one)
 364  376                                  Am1 = half * (y2 / (R + xp1) + y2 / (S - xm1));
 365  377                          else
 366  378                                  Am1 = half * (y2 / (R + xp1) + (S + xm1));
      379 +
 367  380                          D_IM(ans) = log1p(Am1 + sqrt(Am1 * (A + one)));
 368  381                  } else {
 369  382                          D_IM(ans) = log(A + sqrt(A * A - one));
 370  383                  }
 371  384          }
 372  385  
 373  386          if (hx < 0)
 374  387                  D_RE(ans) = -D_RE(ans);
      388 +
 375  389          if (hy < 0)
 376  390                  D_IM(ans) = -D_IM(ans);
 377  391  
 378  392          return (ans);
 379  393  }
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX