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 __cacos = cacos
  31 
  32 /* INDENT OFF */
  33 /*
  34  * dcomplex cacos(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 cosine function cacos(z),
  42  * where z = x+iy, can be defined by
  43  *
  44  *      cacos(z) = acos(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  * Basic relations
  67  *    cacos(conj(z)) = conj(cacos(z))
  68  *    cacos(-z)      = pi   - cacos(z)
  69  *    cacos( z)      = pi/2 - casin(z)
  70  *
  71  * Special cases (conform to ISO/IEC 9899:1999(E)):
  72  *    cacos(+-0  + i y  ) = pi/2 - i y for y is +-0, +-inf, NaN
  73  *    cacos( x   + i inf) = pi/2 - i inf for all x
  74  *    cacos( x   + i NaN) = NaN  + i NaN with invalid for non-zero finite x
  75  *    cacos(-inf + i y  ) = pi   - i inf for finite +y
  76  *    cacos( inf + i y  ) = 0    - i inf for finite +y
  77  *    cacos(-inf + i inf) = 3pi/4- i inf
  78  *    cacos( inf + i inf) = pi/4 - i inf
  79  *    cacos(+-inf+ i NaN) = NaN  - i inf (sign of imaginary is unspecified)
  80  *    cacos(NaN  + i y  ) = NaN  + i NaN with invalid for finite y
  81  *    cacos(NaN  + i inf) = NaN  - i inf
  82  *    cacos(NaN  + i NaN) = NaN  + i NaN
  83  *
  84  * Special Regions (better formula for accuracy and for avoiding spurious
  85  * overflow or underflow) (all x and y are assumed nonnegative):
  86  *  case 1: y = 0
  87  *  case 2: tiny y relative to x-1: y <= ulp(0.5)*|x-1|
  88  *  case 3: tiny y: y < 4 sqrt(u), where u = minimum normal number
  89  *  case 4: huge y relative to x+1: y >= (1+x)/ulp(0.5)
  90  *  case 5: huge x and y: x and y >= sqrt(M)/8, where M = maximum normal number
  91  *  case 6: tiny x: x < 4 sqrt(u)
  92  *  --------
  93  *  case        1 & 2. y=0 or y/|x-1| is tiny. We have
  94  *             ____________              _____________
  95  *            /      2    2             /       y    2
  96  *           / (x+-1)  + y   =  |x+-1| / 1 + (------)
  97  *         \/                        \/       |x+-1|
  98  *
  99  *                                            1     y    2
 100  *                           ~  |x+-1| ( 1 + --- (------)  )
 101  *                                            2   |x+-1|
 102  *
 103  *                                          2
 104  *                                         y
 105  *                           = |x+-1| + --------.
 106  *                                      2|x+-1|
 107  *
 108  *      Consequently, it is not difficult to see that
 109  *                                 2
 110  *                                y
 111  *                    [ 1 + ------------ ,     if x < 1,
 112  *                    [      2(1+x)(1-x)
 113  *                    [
 114  *                    [
 115  *                    [ x,                     if x = 1 (y = 0),
 116  *                    [
 117  *              A ~=  [             2
 118  *                    [        x * y
 119  *                    [ x + ------------ ~ x,  if x > 1
 120  *                    [      2(x+1)(x-1)
 121  *
 122  *      and hence
 123  *                      ______                                 2
 124  *                     / 2                    y               y
 125  *               A + \/ A  - 1  ~  1 + ---------------- + -----------, if x < 1,
 126  *                                     sqrt((x+1)(1-x))   2(x+1)(1-x)
 127  *
 128  *
 129  *                              ~  x + sqrt((x-1)*(x+1)),             if x >= 1.
 130  *
 131  *                                         2
 132  *                                        y
 133  *                          [ x(1 - -----------) ~ x,  if x < 1,
 134  *                          [       2(1+x)(1-x)
 135  *              B = x/A  ~  [
 136  *                          [ 1,                       if x = 1,
 137  *                          [
 138  *                          [           2
 139  *                          [          y
 140  *                          [ 1 - ------------ ,       if x > 1,
 141  *                          [      2(x+1)(x-1)
 142  *      Thus
 143  *                            [ acos(x) - i y/sqrt((x-1)*(x+1)),      if x < 1,
 144  *                            [
 145  *              cacos(x+i*y)~ [ 0 - i 0,                              if x = 1,
 146  *                            [
 147  *                            [ y/sqrt(x*x-1) - i log(x+sqrt(x*x-1)), if x > 1.
 148  *
 149  *      Note: y/sqrt(x*x-1) ~ y/x when x >= 2**26.
 150  *  case 3. y < 4 sqrt(u), where u = minimum normal x.
 151  *      After case 1 and 2, this will only occurs when x=1. When x=1, we have
 152  *         A = (sqrt(4+y*y)+y)/2 ~ 1 + y/2 + y^2/8 + ...
 153  *      and
 154  *         B = 1/A = 1 - y/2 + y^2/8 + ...
 155  *      Since
 156  *         cos(sqrt(y)) ~ 1 - y/2 + ...
 157  *      we have, for the real part,
 158  *         acos(B) ~ acos(1 - y/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 = acos(B) ~ pi/2
 168  *      and
 169  *         imag part = log(y+sqrt(y*y-one))
 170  *
 171  *  case 5. Both x and y are large: x and y > sqrt(M)/8, where M = maximum x
 172  *      In this case,
 173  *         A ~ sqrt(x*x+y*y)
 174  *         B ~ x/sqrt(x*x+y*y).
 175  *      Thus
 176  *         real part = acos(B) = atan(y/x),
 177  *         imag part = log(A+sqrt(A*A-1)) ~ log(2A)
 178  *                   = log(2) + 0.5*log(x*x+y*y)
 179  *                   = log(2) + log(y) + 0.5*log(1+(x/y)^2)
 180  *
 181  *  case 6. x < 4 sqrt(u). In this case, we have
 182  *          A ~ sqrt(1+y*y), B = x/sqrt(1+y*y).
 183  *      Since B is tiny, we have
 184  *          real part = acos(B) ~ pi/2
 185  *          imag part = log(A+sqrt(A*A-1)) = log (A+sqrt(y*y))
 186  *                    = log(y+sqrt(1+y*y))
 187  *                    = 0.5*log(y^2+2ysqrt(1+y^2)+1+y^2)
 188  *                    = 0.5*log(1+2y(y+sqrt(1+y^2)));
 189  *                    = 0.5*log1p(2y(y+A));
 190  *
 191  *      cacos(z) = acos(B) - i sign(y) log (A + sqrt(A*A-1)),
 192  */
 193 /* INDENT ON */
 194 
 195 #include "libm.h"
 196 #include "complex_wrapper.h"
 197 
 198 /* INDENT OFF */
 199 static const double
 200         zero = 0.0,
 201         one = 1.0,
 202         E = 1.11022302462515654042e-16,                 /* 2**-53 */
 203         ln2 = 6.93147180559945286227e-01,
 204         pi = 3.1415926535897931159979634685,
 205         pi_l = 1.224646799147353177e-16,
 206         pi_2 = 1.570796326794896558e+00,
 207         pi_2_l = 6.123233995736765886e-17,
 208         pi_4 = 0.78539816339744827899949,
 209         pi_4_l = 3.061616997868382943e-17,
 210         pi3_4 = 2.356194490192344836998,
 211         pi3_4_l = 9.184850993605148829195e-17,
 212         Foursqrtu = 5.96667258496016539463e-154,        /* 2**(-509) */
 213         Acrossover = 1.5,
 214         Bcrossover = 0.6417,
 215         half = 0.5;
 216 /* INDENT ON */
 217 
 218 dcomplex
 219 cacos(dcomplex z) {
 220         double x, y, t, R, S, A, Am1, B, y2, xm1, xp1, Apx;
 221         int ix, iy, hx, hy;
 222         unsigned lx, ly;
 223         dcomplex ans;
 224 
 225         x = D_RE(z);
 226         y = D_IM(z);
 227         hx = HI_WORD(x);
 228         lx = LO_WORD(x);
 229         hy = HI_WORD(y);
 230         ly = LO_WORD(y);
 231         ix = hx & 0x7fffffff;
 232         iy = hy & 0x7fffffff;
 233 
 234         /* x is 0 */
 235         if ((ix | lx) == 0) {
 236                 if (((iy | ly) == 0) || (iy >= 0x7ff00000)) {
 237                         D_RE(ans) = pi_2;
 238                         D_IM(ans) = -y;
 239                         return (ans);
 240                 }
 241         }
 242 
 243         /* |y| is inf or NaN */
 244         if (iy >= 0x7ff00000) {
 245                 if (ISINF(iy, ly)) {    /* cacos(x + i inf) = pi/2  - i inf */
 246                         D_IM(ans) = -y;
 247                         if (ix < 0x7ff00000) {
 248                                 D_RE(ans) = pi_2 + pi_2_l;
 249                         } else if (ISINF(ix, lx)) {
 250                                 if (hx >= 0)
 251                                         D_RE(ans) = pi_4 + pi_4_l;
 252                                 else
 253                                         D_RE(ans) = pi3_4 + pi3_4_l;
 254                         } else {
 255                                 D_RE(ans) = x;
 256                         }
 257                 } else {                /* cacos(x + i NaN) = NaN  + i NaN */
 258                         D_RE(ans) = y + x;
 259                         if (ISINF(ix, lx))
 260                                 D_IM(ans) = -fabs(x);
 261                         else
 262                                 D_IM(ans) = y;
 263                 }
 264                 return (ans);
 265         }
 266 
 267         x = fabs(x);
 268         y = fabs(y);
 269 
 270         /* x is inf or NaN */
 271         if (ix >= 0x7ff00000) {      /* x is inf or NaN */
 272                 if (ISINF(ix, lx)) {    /* x is INF */
 273                         D_IM(ans) = -x;
 274                         if (iy >= 0x7ff00000) {
 275                                 if (ISINF(iy, ly)) {
 276                                         /* INDENT OFF */
 277                                         /* cacos(inf + i inf) = pi/4 - i inf */
 278                                         /* cacos(-inf+ i inf) =3pi/4 - i inf */
 279                                         /* INDENT ON */
 280                                         if (hx >= 0)
 281                                                 D_RE(ans) = pi_4 + pi_4_l;
 282                                         else
 283                                                 D_RE(ans) = pi3_4 + pi3_4_l;
 284                                 } else
 285                                         /* INDENT OFF */
 286                                         /* cacos(inf + i NaN) = NaN  - i inf  */
 287                                         /* INDENT ON */
 288                                         D_RE(ans) = y + y;
 289                         } else
 290                                 /* INDENT OFF */
 291                                 /* cacos(inf + iy ) = 0  - i inf */
 292                                 /* cacos(-inf+ iy  ) = pi - i inf */
 293                                 /* INDENT ON */
 294                         if (hx >= 0)
 295                                 D_RE(ans) = zero;
 296                         else
 297                                 D_RE(ans) = pi + pi_l;
 298                 } else {                /* x is NaN */
 299                         /* INDENT OFF */
 300                         /*
 301                          * cacos(NaN + i inf) = NaN  - i inf
 302                          * cacos(NaN + i y  ) = NaN  + i NaN
 303                          * cacos(NaN + i NaN) = NaN  + i NaN
 304                          */
 305                         /* INDENT ON */
 306                         D_RE(ans) = x + y;
 307                         if (iy >= 0x7ff00000) {
 308                                 D_IM(ans) = -y;
 309                         } else {
 310                                 D_IM(ans) = x;
 311                         }
 312                 }
 313                 if (hy < 0)
 314                         D_IM(ans) = -D_IM(ans);
 315                 return (ans);
 316         }
 317 
 318         if ((iy | ly) == 0) {   /* region 1: y=0 */
 319                 if (ix < 0x3ff00000) {       /* |x| < 1 */
 320                         D_RE(ans) = acos(x);
 321                         D_IM(ans) = zero;
 322                 } else {
 323                         D_RE(ans) = zero;
 324                         if (ix >= 0x43500000)        /* |x| >= 2**54 */
 325                                 D_IM(ans) = ln2 + log(x);
 326                         else if (ix >= 0x3ff80000)   /* x > Acrossover */
 327                                 D_IM(ans) = log(x + sqrt((x - one) * (x +
 328                                         one)));
 329                         else {
 330                                 xm1 = x - one;
 331                                 D_IM(ans) = log1p(xm1 + sqrt(xm1 * (x + one)));
 332                         }
 333                 }
 334         } else if (y <= E * fabs(x - one)) { /* region 2: y < tiny*|x-1| */
 335                 if (ix < 0x3ff00000) {       /* x < 1 */
 336                         D_RE(ans) = acos(x);
 337                         D_IM(ans) = y / sqrt((one + x) * (one - x));
 338                 } else if (ix >= 0x43500000) {       /* |x| >= 2**54 */
 339                         D_RE(ans) = y / x;
 340                         D_IM(ans) = ln2 + log(x);
 341                 } else {
 342                         t = sqrt((x - one) * (x + one));
 343                         D_RE(ans) = y / t;
 344                         if (ix >= 0x3ff80000)        /* x > Acrossover */
 345                                 D_IM(ans) = log(x + t);
 346                         else
 347                                 D_IM(ans) = log1p((x - one) + t);
 348                 }
 349         } else if (y < Foursqrtu) {  /* region 3 */
 350                 t = sqrt(y);
 351                 D_RE(ans) = t;
 352                 D_IM(ans) = t;
 353         } else if (E * y - one >= x) {       /* region 4 */
 354                 D_RE(ans) = pi_2;
 355                 D_IM(ans) = ln2 + log(y);
 356         } else if (ix >= 0x5fc00000 || iy >= 0x5fc00000) {        /* x,y>2**509 */
 357                 /* region 5: x+1 or y is very large (>= sqrt(max)/8) */
 358                 t = x / y;
 359                 D_RE(ans) = atan(y / x);
 360                 D_IM(ans) = ln2 + log(y) + half * log1p(t * t);
 361         } else if (x < Foursqrtu) {
 362                 /* region 6: x is very small, < 4sqrt(min) */
 363                 D_RE(ans) = pi_2;
 364                 A = sqrt(one + y * y);
 365                 if (iy >= 0x3ff80000)        /* if y > Acrossover */
 366                         D_IM(ans) = log(y + A);
 367                 else
 368                         D_IM(ans) = half * log1p((y + y) * (y + A));
 369         } else {        /* safe region */
 370                 y2 = y * y;
 371                 xp1 = x + one;
 372                 xm1 = x - one;
 373                 R = sqrt(xp1 * xp1 + y2);
 374                 S = sqrt(xm1 * xm1 + y2);
 375                 A = half * (R + S);
 376                 B = x / A;
 377                 if (B <= Bcrossover)
 378                         D_RE(ans) = acos(B);
 379                 else {          /* use atan and an accurate approx to a-x */
 380                         Apx = A + x;
 381                         if (x <= one)
 382                                 D_RE(ans) = atan(sqrt(half * Apx * (y2 / (R +
 383                                         xp1) + (S - xm1))) / x);
 384                         else
 385                                 D_RE(ans) = atan((y * sqrt(half * (Apx / (R +
 386                                         xp1) + Apx / (S + xm1)))) / x);
 387                 }
 388                 if (A <= Acrossover) {
 389                         /* use log1p and an accurate approx to A-1 */
 390                         if (x < one)
 391                                 Am1 = half * (y2 / (R + xp1) + y2 / (S - xm1));
 392                         else
 393                                 Am1 = half * (y2 / (R + xp1) + (S + xm1));
 394                         D_IM(ans) = log1p(Am1 + sqrt(Am1 * (A + one)));
 395                 } else {
 396                         D_IM(ans) = log(A + sqrt(A * A - one));
 397                 }
 398         }
 399         if (hx < 0)
 400                 D_RE(ans) = pi - D_RE(ans);
 401         if (hy >= 0)
 402                 D_IM(ans) = -D_IM(ans);
 403         return (ans);
 404 }