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 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __cacos = cacos
  32 
  33 
  34 /*
  35  * dcomplex cacos(dcomplex z);
  36  *
  37  * Alogrithm
  38  * (based on T.E.Hull, Thomas F. Fairgrieve and Ping Tak Peter Tang's
  39  * paper "Implementing the Complex Arcsine and Arccosine Functins Using
  40  * Exception Handling", ACM TOMS, Vol 23, pp 299-335)
  41  *
  42  * The principal value of complex inverse cosine function cacos(z),
  43  * where z = x+iy, can be defined by
  44  *
  45  *      cacos(z) = acos(B) - i sign(y) log (A + sqrt(A*A-1)),
  46  *
  47  * where the log function is the natural log, and
  48  *             ____________           ____________
  49  *       1    /     2    2      1    /     2    2
  50  *  A = ---  / (x+1)  + y   +  ---  / (x-1)  + y
  51  *       2 \/                   2 \/
  52  *             ____________           ____________
  53  *       1    /     2    2      1    /     2    2
  54  *  B = ---  / (x+1)  + y   -  ---  / (x-1)  + y   .
  55  *       2 \/                   2 \/
  56  *
  57  * The Branch cuts are on the real line from -inf to -1 and from 1 to inf.
  58  * The real and imaginary parts are based on Abramowitz and Stegun
  59  * [Handbook of Mathematic Functions, 1972].  The sign of the imaginary
  60  * part is chosen to be the generally considered the principal value of
  61  * this function.
  62  *
  63  * Notes:1. A is the average of the distances from z to the points (1,0)
  64  *          and (-1,0) in the complex z-plane, and in particular A>=1.
  65  *       2. B is in [-1,1], and A*B = x
  66  *
  67  * Basic relations
  68  *    cacos(conj(z)) = conj(cacos(z))
  69  *    cacos(-z)      = pi   - cacos(z)
  70  *    cacos( z)      = pi/2 - casin(z)
  71  *
  72  * Special cases (conform to ISO/IEC 9899:1999(E)):
  73  *    cacos(+-0  + i y  ) = pi/2 - i y for y is +-0, +-inf, NaN
  74  *    cacos( x   + i inf) = pi/2 - i inf for all x
  75  *    cacos( x   + i NaN) = NaN  + i NaN with invalid for non-zero finite x
  76  *    cacos(-inf + i y  ) = pi   - i inf for finite +y
  77  *    cacos( inf + i y  ) = 0    - i inf for finite +y
  78  *    cacos(-inf + i inf) = 3pi/4- i inf
  79  *    cacos( inf + i inf) = pi/4 - i inf
  80  *    cacos(+-inf+ i NaN) = NaN  - i inf (sign of imaginary is unspecified)
  81  *    cacos(NaN  + i y  ) = NaN  + i NaN with invalid for finite y
  82  *    cacos(NaN  + i inf) = NaN  - i inf
  83  *    cacos(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 + ------------ ~ x,  if x > 1
 121  *                    [      2(x+1)(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 - -----------) ~ x,  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(x+1)(x-1)
 143  *      Thus
 144  *                            [ acos(x) - i y/sqrt((x-1)*(x+1)),      if x < 1,
 145  *                            [
 146  *              cacos(x+i*y)~ [ 0 - i 0,                              if x = 1,
 147  *                            [
 148  *                            [ y/sqrt(x*x-1) - i log(x+sqrt(x*x-1)), if x > 1.
 149  *
 150  *      Note: y/sqrt(x*x-1) ~ y/x when x >= 2**26.
 151  *  case 3. y < 4 sqrt(u), where u = minimum normal x.
 152  *      After case 1 and 2, this will only occurs when x=1. When x=1, we have
 153  *         A = (sqrt(4+y*y)+y)/2 ~ 1 + y/2 + y^2/8 + ...
 154  *      and
 155  *         B = 1/A = 1 - y/2 + y^2/8 + ...
 156  *      Since
 157  *         cos(sqrt(y)) ~ 1 - y/2 + ...
 158  *      we have, for the real part,
 159  *         acos(B) ~ acos(1 - y/2) ~ sqrt(y)
 160  *      For the imaginary part,
 161  *         log(A+sqrt(A*A-1)) ~ log(1+y/2+sqrt(2*y/2))
 162  *                            = log(1+y/2+sqrt(y))
 163  *                            = (y/2+sqrt(y)) - (y/2+sqrt(y))^2/2 + ...
 164  *                            ~ sqrt(y) - y*(sqrt(y)+y/2)/2
 165  *                            ~ sqrt(y)
 166  *
 167  *  case 4. y >= (x+1)/ulp(0.5). In this case, A ~ y and B ~ x/y. Thus
 168  *         real part = acos(B) ~ pi/2
 169  *      and
 170  *         imag part = log(y+sqrt(y*y-one))
 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 = acos(B) = atan(y/x),
 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 = acos(B) ~ pi/2
 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  *      cacos(z) = acos(B) - i sign(y) log (A + sqrt(A*A-1)),
 193  */
 194 
 195 #include "libm.h"
 196 #include "complex_wrapper.h"
 197 
 198 static const double zero = 0.0,
 199         one = 1.0,
 200         E = 1.11022302462515654042e-16, /* 2**-53 */
 201         ln2 = 6.93147180559945286227e-01,
 202         pi = 3.1415926535897931159979634685,
 203         pi_l = 1.224646799147353177e-16,
 204         pi_2 = 1.570796326794896558e+00,
 205         pi_2_l = 6.123233995736765886e-17,
 206         pi_4 = 0.78539816339744827899949,
 207         pi_4_l = 3.061616997868382943e-17,
 208         pi3_4 = 2.356194490192344836998,
 209         pi3_4_l = 9.184850993605148829195e-17,
 210         Foursqrtu = 5.96667258496016539463e-154,        /* 2**(-509) */
 211         Acrossover = 1.5,
 212         Bcrossover = 0.6417,
 213         half = 0.5;
 214 
 215 
 216 dcomplex
 217 cacos(dcomplex z)
 218 {
 219         double x, y, t, R, S, A, Am1, B, y2, xm1, xp1, Apx;
 220         int ix, iy, hx, hy;
 221         unsigned lx, ly;
 222         dcomplex ans;
 223 
 224         x = D_RE(z);
 225         y = D_IM(z);
 226         hx = HI_WORD(x);
 227         lx = LO_WORD(x);
 228         hy = HI_WORD(y);
 229         ly = LO_WORD(y);
 230         ix = hx & 0x7fffffff;
 231         iy = hy & 0x7fffffff;
 232 
 233         /* x is 0 */
 234         if ((ix | lx) == 0) {
 235                 if (((iy | ly) == 0) || (iy >= 0x7ff00000)) {
 236                         D_RE(ans) = pi_2;
 237                         D_IM(ans) = -y;
 238                         return (ans);
 239                 }
 240         }
 241 
 242         /* |y| is inf or NaN */
 243         if (iy >= 0x7ff00000) {
 244                 if (ISINF(iy, ly)) {    /* cacos(x + i inf) = pi/2  - i inf */
 245                         D_IM(ans) = -y;
 246 
 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 
 260                         if (ISINF(ix, lx))
 261                                 D_IM(ans) = -fabs(x);
 262                         else
 263                                 D_IM(ans) = y;
 264                 }
 265 
 266                 return (ans);
 267         }
 268 
 269         x = fabs(x);
 270         y = fabs(y);
 271 
 272         /* x is inf or NaN */
 273         if (ix >= 0x7ff00000) {              /* x is inf or NaN */
 274                 if (ISINF(ix, lx)) {    /* x is INF */
 275                         D_IM(ans) = -x;
 276 
 277                         if (iy >= 0x7ff00000) {
 278                                 if (ISINF(iy, ly)) {
 279                                         /*
 280                                          * cacos(inf + i inf) = pi/4 - i inf
 281                                          * cacos(-inf+ i inf) =3pi/4 - i inf
 282                                          */
 283                                         if (hx >= 0)
 284                                                 D_RE(ans) = pi_4 + pi_4_l;
 285                                         else
 286                                                 D_RE(ans) = pi3_4 + pi3_4_l;
 287                                 } else {
 288                                         /*
 289                                          * cacos(inf + i NaN) = NaN  - i inf
 290                                          */
 291                                         D_RE(ans) = y + y;
 292                                 }
 293                         } else
 294                         /*
 295                          * cacos(inf + iy ) = 0  - i inf
 296                          * cacos(-inf+ iy  ) = pi - i inf
 297                          */
 298                         if (hx >= 0) {
 299                                 D_RE(ans) = zero;
 300                         } else {
 301                                 D_RE(ans) = pi + pi_l;
 302                         }
 303                 } else {                /* x is NaN */
 304 
 305                         /*
 306                          * cacos(NaN + i inf) = NaN  - i inf
 307                          * cacos(NaN + i y  ) = NaN  + i NaN
 308                          * cacos(NaN + i NaN) = NaN  + i NaN
 309                          */
 310                         D_RE(ans) = x + y;
 311 
 312                         if (iy >= 0x7ff00000)
 313                                 D_IM(ans) = -y;
 314                         else
 315                                 D_IM(ans) = x;
 316                 }
 317 
 318                 if (hy < 0)
 319                         D_IM(ans) = -D_IM(ans);
 320 
 321                 return (ans);
 322         }
 323 
 324         if ((iy | ly) == 0) {           /* region 1: y=0 */
 325                 if (ix < 0x3ff00000) {       /* |x| < 1 */
 326                         D_RE(ans) = acos(x);
 327                         D_IM(ans) = zero;
 328                 } else {
 329                         D_RE(ans) = zero;
 330 
 331                         if (ix >= 0x43500000) {              /* |x| >= 2**54 */
 332                                 D_IM(ans) = ln2 + log(x);
 333                         } else if (ix >= 0x3ff80000) {       /* x > Acrossover */
 334                                 D_IM(ans) = log(x + sqrt((x - one) * (x +
 335                                     one)));
 336                         } else {
 337                                 xm1 = x - one;
 338                                 D_IM(ans) = log1p(xm1 + sqrt(xm1 * (x + one)));
 339                         }
 340                 }
 341         } else if (y <= E * fabs(x - one)) { /* region 2: y < tiny*|x-1| */
 342                 if (ix < 0x3ff00000) {               /* x < 1 */
 343                         D_RE(ans) = acos(x);
 344                         D_IM(ans) = y / sqrt((one + x) * (one - x));
 345                 } else if (ix >= 0x43500000) {       /* |x| >= 2**54 */
 346                         D_RE(ans) = y / x;
 347                         D_IM(ans) = ln2 + log(x);
 348                 } else {
 349                         t = sqrt((x - one) * (x + one));
 350                         D_RE(ans) = y / t;
 351 
 352                         if (ix >= 0x3ff80000)        /* x > Acrossover */
 353                                 D_IM(ans) = log(x + t);
 354                         else
 355                                 D_IM(ans) = log1p((x - one) + t);
 356                 }
 357         } else if (y < Foursqrtu) {  /* region 3 */
 358                 t = sqrt(y);
 359                 D_RE(ans) = t;
 360                 D_IM(ans) = t;
 361         } else if (E * y - one >= x) {                               /* region 4 */
 362                 D_RE(ans) = pi_2;
 363                 D_IM(ans) = ln2 + log(y);
 364         } else if (ix >= 0x5fc00000 || iy >= 0x5fc00000) {        /* x,y>2**509 */
 365                 /* region 5: x+1 or y is very large (>= sqrt(max)/8) */
 366                 t = x / y;
 367                 D_RE(ans) = atan(y / x);
 368                 D_IM(ans) = ln2 + log(y) + half * log1p(t * t);
 369         } else if (x < Foursqrtu) {
 370                 /* region 6: x is very small, < 4sqrt(min) */
 371                 D_RE(ans) = pi_2;
 372                 A = sqrt(one + y * y);
 373 
 374                 if (iy >= 0x3ff80000)        /* if y > Acrossover */
 375                         D_IM(ans) = log(y + A);
 376                 else
 377                         D_IM(ans) = half * log1p((y + y) * (y + A));
 378         } else {                        /* safe region */
 379                 y2 = y * y;
 380                 xp1 = x + one;
 381                 xm1 = x - one;
 382                 R = sqrt(xp1 * xp1 + y2);
 383                 S = sqrt(xm1 * xm1 + y2);
 384                 A = half * (R + S);
 385                 B = x / A;
 386 
 387                 if (B <= Bcrossover) {
 388                         D_RE(ans) = acos(B);
 389                 } else {        /* use atan and an accurate approx to a-x */
 390                         Apx = A + x;
 391 
 392                         if (x <= one)
 393                                 D_RE(ans) = atan(sqrt(half * Apx * (y2 / (R +
 394                                     xp1) + (S - xm1))) / x);
 395                         else
 396                                 D_RE(ans) = atan((y * sqrt(half * (Apx / (R +
 397                                     xp1) + Apx / (S + xm1)))) / x);
 398                 }
 399 
 400                 if (A <= Acrossover) {
 401                         /* use log1p and an accurate approx to A-1 */
 402                         if (x < one)
 403                                 Am1 = half * (y2 / (R + xp1) + y2 / (S - xm1));
 404                         else
 405                                 Am1 = half * (y2 / (R + xp1) + (S + xm1));
 406 
 407                         D_IM(ans) = log1p(Am1 + sqrt(Am1 * (A + one)));
 408                 } else {
 409                         D_IM(ans) = log(A + sqrt(A * A - one));
 410                 }
 411         }
 412 
 413         if (hx < 0)
 414                 D_RE(ans) = pi - D_RE(ans);
 415 
 416         if (hy >= 0)
 417                 D_IM(ans) = -D_IM(ans);
 418 
 419         return (ans);
 420 }