usr/src/lib/libm/common/complex/casin.c

   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29
  30 #pragma weak __casin = casin
  31
  32 /* INDENT OFF */
  33 /*
  34  * dcomplex casin(dcomplex z);
  35  *
  36  * Alogrithm
  37  * (based on T.E.Hull, Thomas F. Fairgrieve and Ping Tak Peter Tang's
  38  * paper "Implementing the Complex Arcsine and Arccosine Functins Using
  39  * Exception Handling", ACM TOMS, Vol 23, pp 299-335)
  40  *
  41  * The principal value of complex inverse sine function casin(z),
  42  * where z = x+iy, can be defined by
  43  *
  44  *      casin(z) = asin(B) + i sign(y) log (A + sqrt(A*A-1)),
  45  *
  46  * where the log function is the natural log, and
  47  *             ____________           ____________
  48  *       1    /     2    2      1    /     2    2
  49  *  A = ---  / (x+1)  + y   +  ---  / (x-1)  + y
  50  *       2 \/                   2 \/
  51  *             ____________           ____________
  52  *       1    /     2    2      1    /     2    2
  53  *  B = ---  / (x+1)  + y   -  ---  / (x-1)  + y   .
  54  *       2 \/                   2 \/
  55  *
  56  * The Branch cuts are on the real line from -inf to -1 and from 1 to inf.
  57  * The real and imaginary parts are based on Abramowitz and Stegun
  58  * [Handbook of Mathematic Functions, 1972].  The sign of the imaginary
  59  * part is chosen to be the generally considered the principal value of
  60  * this function.
  61  *
  62  * Notes:1. A is the average of the distances from z to the points (1,0)
  63  *          and (-1,0) in the complex z-plane, and in particular A>=1.
  64  *       2. B is in [-1,1], and A*B = x.
  65  *
  66  * Special notes: if casin( x, y) = ( u, v), then
  67  *                  casin(-x, y) = (-u, v),
  68  *                  casin( x,-y) = ( u,-v),
  69  *    in general, we have casin(conj(z))     =  conj(casin(z))
  70  *                       casin(-z)          = -casin(z)
  71  *                       casin(z)           =  pi/2 - cacos(z)
  72  *
  73  * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
  74  *    casin( 0 + i 0   ) =  0    + i 0
  75  *    casin( 0 + i NaN ) =  0    + i NaN
  76  *    casin( x + i inf ) =  0    + i inf for finite x
  77  *    casin( x + i NaN ) =  NaN  + i NaN with invalid for finite x != 0
  78  *    casin(inf + iy   ) =  pi/2 + i inf finite y
  79  *    casin(inf + i inf) =  pi/4 + i inf
  80  *    casin(inf + i NaN) =  NaN  + i inf
  81  *    casin(NaN + i y  ) =  NaN  + i NaN for finite y
  82  *    casin(NaN + i inf) =  NaN  + i inf
  83  *    casin(NaN + i NaN) =  NaN  + i NaN
  84  *
  85  * Special Regions (better formula for accuracy and for avoiding spurious
  86  * overflow or underflow) (all x and y are assumed nonnegative):
  87  *  case 1: y = 0
  88  *  case 2: tiny y relative to x-1: y <= ulp(0.5)*|x-1|
  89  *  case 3: tiny y: y < 4 sqrt(u), where u = minimum normal number
  90  *  case 4: huge y relative to x+1: y >= (1+x)/ulp(0.5)
  91  *  case 5: huge x and y: x and y >= sqrt(M)/8, where M = maximum normal number
  92  *  case 6: tiny x: x < 4 sqrt(u)
  93  *  --------
  94  *  case        1 & 2. y=0 or y/|x-1| is tiny. We have
  95  *             ____________              _____________
  96  *            /      2    2             /       y    2
  97  *           / (x+-1)  + y   =  |x+-1| / 1 + (------)
  98  *         \/                        \/       |x+-1|
  99  *
 100  *                                            1      y   2
 101  *                           ~  |x+-1| ( 1 + --- (------)  )
 102  *                                            2   |x+-1|
 103  *
 104  *                                           2
 105  *                                          y
 106  *                           =  |x+-1| + --------.
 107  *                                       2|x+-1|
 108  *
 109  *      Consequently, it is not difficult to see that
 110  *                                 2
 111  *                                y
 112  *                    [ 1 + ------------ ,  if x < 1,
 113  *                    [      2(1+x)(1-x)
 114  *                    [
 115  *                    [
 116  *                    [ x,                 if x = 1 (y = 0),
 117  *                    [
 118  *              A ~=  [             2
 119  *                    [        x * y
 120  *                    [ x + ------------ ,  if x > 1
 121  *                    [      2(1+x)(x-1)
 122  *
 123  *      and hence
 124  *                      ______                                 2
 125  *                     / 2                    y               y
 126  *               A + \/ A  - 1  ~  1 + ---------------- + -----------, if x < 1,
 127  *                                     sqrt((x+1)(1-x))   2(x+1)(1-x)
 128  *
 129  *
 130  *                             ~  x + sqrt((x-1)*(x+1)),              if x >= 1.
 131  *
 132  *                                         2
 133  *                                        y
 134  *                          [ x(1 - ------------), if x < 1,
 135  *                          [       2(1+x)(1-x)
 136  *              B = x/A  ~  [
 137  *                          [ 1,                  if x = 1,
 138  *                          [
 139  *                          [           2
 140  *                          [          y
 141  *                          [ 1 - ------------ ,   if x > 1,
 142  *                          [      2(1+x)(1-x)
 143  *      Thus
 144  *                            [ asin(x) + i y/sqrt((x-1)*(x+1)), if x <  1
 145  *              casin(x+i*y)=[
 146  *                            [ pi/2    + i log(x+sqrt(x*x-1)),  if x >= 1
 147  *
 148  *  case 3. y < 4 sqrt(u), where u = minimum normal x.
 149  *      After case 1 and 2, this will only occurs when x=1. When x=1, we have
 150  *         A = (sqrt(4+y*y)+y)/2 ~ 1 + y/2 + y^2/8 + ...
 151  *      and
 152  *         B = 1/A = 1 - y/2 + y^2/8 + ...
 153  *      Since
 154  *         asin(x) = pi/2-2*asin(sqrt((1-x)/2))
 155  *         asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
 156  *      we have, for the real part asin(B),
 157  *         asin(1-y/2) ~ pi/2 - 2 asin(sqrt(y/4))
 158  *                     ~ pi/2 - sqrt(y)
 159  *      For the imaginary part,
 160  *         log(A+sqrt(A*A-1)) ~ log(1+y/2+sqrt(2*y/2))
 161  *                            = log(1+y/2+sqrt(y))
 162  *                            = (y/2+sqrt(y)) - (y/2+sqrt(y))^2/2 + ...
 163  *                            ~ sqrt(y) - y*(sqrt(y)+y/2)/2
 164  *                            ~ sqrt(y)
 165  *
 166  *  case 4. y >= (x+1)ulp(0.5). In this case, A ~ y and B ~ x/y. Thus
 167  *         real part = asin(B) ~ x/y (be careful, x/y may underflow)
 168  *      and
 169  *         imag part = log(y+sqrt(y*y-one))
 170  *
 171  *
 172  *  case 5. Both x and y are large: x and y > sqrt(M)/8, where M = maximum x
 173  *      In this case,
 174  *         A ~ sqrt(x*x+y*y)
 175  *         B ~ x/sqrt(x*x+y*y).
 176  *      Thus
 177  *         real part = asin(B) = atan(x/y),
 178  *         imag part = log(A+sqrt(A*A-1)) ~ log(2A)
 179  *                   = log(2) + 0.5*log(x*x+y*y)
 180  *                   = log(2) + log(y) + 0.5*log(1+(x/y)^2)
 181  *
 182  *  case 6. x < 4 sqrt(u). In this case, we have
 183  *          A ~ sqrt(1+y*y), B = x/sqrt(1+y*y).
 184  *      Since B is tiny, we have
 185  *          real part = asin(B) ~ B = x/sqrt(1+y*y)
 186  *          imag part = log(A+sqrt(A*A-1)) = log (A+sqrt(y*y))
 187  *                    = log(y+sqrt(1+y*y))
 188  *                    = 0.5*log(y^2+2ysqrt(1+y^2)+1+y^2)
 189  *                    = 0.5*log(1+2y(y+sqrt(1+y^2)));
 190  *                    = 0.5*log1p(2y(y+A));
 191  *
 192  *      casin(z) = asin(B) + i sign(y) log (A + sqrt(A*A-1)),
 193  */
 194 /* INDENT ON */
 195
 196 #include "libm.h"               /* asin/atan/fabs/log/log1p/sqrt */
 197 #include "complex_wrapper.h"
 198
 199 /* INDENT OFF */
 200 static const double
 201         zero = 0.0,
 202         one = 1.0,
 203         E = 1.11022302462515654042e-16,                 /* 2**-53 */
 204         ln2 = 6.93147180559945286227e-01,
 205         pi_2 = 1.570796326794896558e+00,
 206         pi_2_l = 6.123233995736765886e-17,
 207         pi_4 = 7.85398163397448278999e-01,
 208         Foursqrtu = 5.96667258496016539463e-154,        /* 2**(-509) */
 209         Acrossover = 1.5,
 210         Bcrossover = 0.6417,
 211         half = 0.5;
 212 /* INDENT ON */
 213
 214 dcomplex
 215 casin(dcomplex z) {
 216         double x, y, t, R, S, A, Am1, B, y2, xm1, xp1, Apx;
 217         int ix, iy, hx, hy;
 218         unsigned lx, ly;
 219         dcomplex ans;
 220
 221         x = D_RE(z);
 222         y = D_IM(z);
 223         hx = HI_WORD(x);
 224         lx = LO_WORD(x);
 225         hy = HI_WORD(y);
 226         ly = LO_WORD(y);
 227         ix = hx & 0x7fffffff;
 228         iy = hy & 0x7fffffff;
 229         x = fabs(x);
 230         y = fabs(y);
 231
 232         /* special cases */
 233
 234         /* x is inf or NaN */
 235         if (ix >= 0x7ff00000) { /* x is inf or NaN */
 236                 if (ISINF(ix, lx)) {    /* x is INF */
 237                         D_IM(ans) = x;
 238                         if (iy >= 0x7ff00000) {
 239                                 if (ISINF(iy, ly))
 240                                         /* casin(inf + i inf) = pi/4 + i inf */
 241                                         D_RE(ans) = pi_4;
 242                                 else    /* casin(inf + i NaN) = NaN  + i inf  */
 243                                         D_RE(ans) = y + y;
 244                         } else  /* casin(inf + iy) = pi/2 + i inf */
 245                                 D_RE(ans) = pi_2;
 246                 } else {                /* x is NaN */
 247                         if (iy >= 0x7ff00000) {
 248                                 /* INDENT OFF */
 249                                 /*
 250                                  * casin(NaN + i inf) = NaN + i inf
 251                                  * casin(NaN + i NaN) = NaN + i NaN
 252                                  */
 253                                 /* INDENT ON */
 254                                 D_IM(ans) = y + y;
 255                                 D_RE(ans) = x + x;
 256                         } else {
 257                                 /* casin(NaN + i y ) = NaN  + i NaN */
 258                                 D_IM(ans) = D_RE(ans) = x + y;
 259                         }
 260                 }
 261                 if (hx < 0)
 262                         D_RE(ans) = -D_RE(ans);
 263                 if (hy < 0)
 264                         D_IM(ans) = -D_IM(ans);
 265                 return (ans);
 266         }
 267
 268         /* casin(+0 + i 0  ) =  0   + i 0. */
 269         if ((ix | lx | iy | ly) == 0)
 270                 return (z);
 271
 272         if (iy >= 0x7ff00000) { /* y is inf or NaN */
 273                 if (ISINF(iy, ly)) {    /* casin(x + i inf) =  0   + i inf */
 274                         D_IM(ans) = y;
 275                         D_RE(ans) = zero;
 276                 } else {                /* casin(x + i NaN) = NaN  + i NaN */
 277                         D_IM(ans) = x + y;
 278                         if ((ix | lx) == 0)
 279                                 D_RE(ans) = x;
 280                         else
 281                                 D_RE(ans) = y;
 282                 }
 283                 if (hx < 0)
 284                         D_RE(ans) = -D_RE(ans);
 285                 if (hy < 0)
 286                         D_IM(ans) = -D_IM(ans);
 287                 return (ans);
 288         }
 289
 290         if ((iy | ly) == 0) {   /* region 1: y=0 */
 291                 if (ix < 0x3ff00000) {  /* |x| < 1 */
 292                         D_RE(ans) = asin(x);
 293                         D_IM(ans) = zero;
 294                 } else {
 295                         D_RE(ans) = pi_2;
 296                         if (ix >= 0x43500000)   /* |x| >= 2**54 */
 297                                 D_IM(ans) = ln2 + log(x);
 298                         else if (ix >= 0x3ff80000)      /* x > Acrossover */
 299                                 D_IM(ans) = log(x + sqrt((x - one) * (x +
 300                                         one)));
 301                         else {
 302                                 xm1 = x - one;
 303                                 D_IM(ans) = log1p(xm1 + sqrt(xm1 * (x + one)));
 304                         }
 305                 }
 306         } else if (y <= E * fabs(x - one)) {    /* region 2: y < tiny*|x-1| */
 307                 if (ix < 0x3ff00000) {  /* x < 1 */
 308                         D_RE(ans) = asin(x);
 309                         D_IM(ans) = y / sqrt((one + x) * (one - x));
 310                 } else {
 311                         D_RE(ans) = pi_2;
 312                         if (ix >= 0x43500000) { /* |x| >= 2**54 */
 313                                 D_IM(ans) = ln2 + log(x);
 314                         } else if (ix >= 0x3ff80000)    /* x > Acrossover */
 315                                 D_IM(ans) = log(x + sqrt((x - one) * (x +
 316                                         one)));
 317                         else
 318                                 D_IM(ans) = log1p((x - one) + sqrt((x - one) *
 319                                         (x + one)));
 320                 }
 321         } else if (y < Foursqrtu) {     /* region 3 */
 322                 t = sqrt(y);
 323                 D_RE(ans) = pi_2 - (t - pi_2_l);
 324                 D_IM(ans) = t;
 325         } else if (E * y - one >= x) {  /* region 4 */
 326                 D_RE(ans) = x / y;      /* need to fix underflow cases */
 327                 D_IM(ans) = ln2 + log(y);
 328         } else if (ix >= 0x5fc00000 || iy >= 0x5fc00000) {      /* x,y>2**509 */
 329                 /* region 5: x+1 or y is very large (>= sqrt(max)/8) */
 330                 t = x / y;
 331                 D_RE(ans) = atan(t);
 332                 D_IM(ans) = ln2 + log(y) + half * log1p(t * t);
 333         } else if (x < Foursqrtu) {
 334                 /* region 6: x is very small, < 4sqrt(min) */
 335                 A = sqrt(one + y * y);
 336                 D_RE(ans) = x / A;      /* may underflow */
 337                 if (iy >= 0x3ff80000)   /* if y > Acrossover */
 338                         D_IM(ans) = log(y + A);
 339                 else
 340                         D_IM(ans) = half * log1p((y + y) * (y + A));
 341         } else {        /* safe region */
 342                 y2 = y * y;
 343                 xp1 = x + one;
 344                 xm1 = x - one;
 345                 R = sqrt(xp1 * xp1 + y2);
 346                 S = sqrt(xm1 * xm1 + y2);
 347                 A = half * (R + S);
 348                 B = x / A;
 349
 350                 if (B <= Bcrossover)
 351                         D_RE(ans) = asin(B);
 352                 else {          /* use atan and an accurate approx to a-x */
 353                         Apx = A + x;
 354                         if (x <= one)
 355                                 D_RE(ans) = atan(x / sqrt(half * Apx * (y2 /
 356                                         (R + xp1) + (S - xm1))));
 357                         else
 358                                 D_RE(ans) = atan(x / (y * sqrt(half * (Apx /
 359                                         (R + xp1) + Apx / (S + xm1)))));
 360                 }
 361                 if (A <= Acrossover) {
 362                         /* use log1p and an accurate approx to A-1 */
 363                         if (x < one)
 364                                 Am1 = half * (y2 / (R + xp1) + y2 / (S - xm1));
 365                         else
 366                                 Am1 = half * (y2 / (R + xp1) + (S + xm1));
 367                         D_IM(ans) = log1p(Am1 + sqrt(Am1 * (A + one)));
 368                 } else {
 369                         D_IM(ans) = log(A + sqrt(A * A - one));
 370                 }
 371         }
 372
 373         if (hx < 0)
 374                 D_RE(ans) = -D_RE(ans);
 375         if (hy < 0)
 376                 D_IM(ans) = -D_IM(ans);
 377
 378         return (ans);
 379 }