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 __catanl = catanl
  32 
  33 
  34 /*
  35  * ldcomplex catanl(ldcomplex z);
  36  *
  37  * Atan(z) return A + Bi where,
  38  *            1
  39  *      A = --- * atan2(2x, 1-x*x-y*y)
  40  *            2
  41  *
  42  *            1      [ x*x + (y+1)*(y+1) ]   1               4y
  43  *       B = --- log [ ----------------- ] = - log (1+ -----------------)
  44  *            4      [ x*x + (y-1)*(y-1) ]   4         x*x + (y-1)*(y-1)
  45  *
  46  *                 2    16  3                         y
  47  *         = t - 2t   + -- t  - ..., where t = -----------------
  48  *                      3                      x*x + (y-1)*(y-1)
  49  * Proof:
  50  * Let w = atan(z=x+yi) = A + B i. Then tan(w) = z.
  51  * Since sin(w) = (exp(iw)-exp(-iw))/(2i), cos(w)=(exp(iw)+exp(-iw))/(2),
  52  * Let p = exp(iw), then z = tan(w) = ((p-1/p)/(p+1/p))/i, or
  53  * iz = (p*p-1)/(p*p+1), or, after simplification,
  54  *      p*p = (1+iz)/(1-iz)                                 ... (1)
  55  * LHS of (1) = exp(2iw) = exp(2i(A+Bi)) = exp(-2B)*exp(2iA)
  56  *            = exp(-2B)*(cos(2A)+i*sin(2A))                ... (2)
  57  *              1-y+ix   (1-y+ix)*(1+y+ix)   1-x*x-y*y + 2xi
  58  * RHS of (1) = ------ = ----------------- = --------------- ... (3)
  59  *              1+y-ix    (1+y)**2 + x**2    (1+y)**2 + x**2
  60  *
  61  * Comparing the real and imaginary parts of (2) and (3), we have:
  62  *      cos(2A) : 1-x*x-y*y = sin(2A) : 2x
  63  * and hence
  64  *      tan(2A) = 2x/(1-x*x-y*y), or
  65  *      A = 0.5 * atan2(2x, 1-x*x-y*y)                      ... (4)
  66  *
  67  * For the imaginary part B, Note that |p*p| = exp(-2B), and
  68  *      |1+iz|   |i-z|   hypot(x,(y-1))
  69  *       |----| = |---| = --------------
  70  *      |1-iz|   |i+z|   hypot(x,(y+1))
  71  * Thus
  72  *                 x*x + (y+1)*(y+1)
  73  *      exp(4B) = -----------------, or
  74  *                 x*x + (y-1)*(y-1)
  75  *
  76  *            1     [x^2+(y+1)^2]   1             4y
  77  *       B =  - log [-----------] = - log(1+ -------------)  ... (5)
  78  *            4     [x^2+(y-1)^2]   4         x^2+(y-1)^2
  79  *
  80  * QED.
  81  *
  82  * Note that: if catan( x, y) = ( u, v), then
  83  *               catan(-x, y) = (-u, v)
  84  *               catan( x,-y) = ( u,-v)
  85  *
  86  * Also,   catan(x,y) = -i*catanh(-y,x), or
  87  *        catanh(x,y) =  i*catan(-y,x)
  88  * So, if catanh(y,x) = (v,u), then catan(x,y) = -i*(-v,u) = (u,v), i.e.,
  89  *         catan(x,y) = (u,v)
  90  *
  91  * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
  92  *    catan( 0  , 0   ) =  (0    ,  0   )
  93  *    catan( NaN, 0   ) =  (NaN  ,  0   )
  94  *    catan( 0  , 1   ) =  (0    ,  +inf) with divide-by-zero
  95  *    catan( inf, y   ) =  (pi/2 ,  0   ) for finite +y
  96  *    catan( NaN, y   ) =  (NaN  ,  NaN ) with invalid for finite y != 0
  97  *    catan( x  , inf ) =  (pi/2 ,  0   ) for finite +x
  98  *    catan( inf, inf ) =  (pi/2 ,  0   )
  99  *    catan( NaN, inf ) =  (NaN  ,  0   )
 100  *    catan( x  , NaN ) =  (NaN  ,  NaN ) with invalid for finite x
 101  *    catan( inf, NaN ) =  (pi/2 ,  +-0 )
 102  */
 103 
 104 #include "libm.h"       /* atan2l/atanl/fabsl/isinfl/iszerol/log1pl/logl */
 105 #include "complex_wrapper.h"
 106 #include "longdouble.h"
 107 
 108 /* BEGIN CSTYLED */
 109 static const long double zero = 0.0L,
 110         one = 1.0L,
 111         two = 2.0L,
 112         half = 0.5L,
 113         ln2 = 6.931471805599453094172321214581765680755e-0001L,
 114         pi_2 = 1.570796326794896619231321691639751442098584699687552910487472L,
 115 #if defined(__x86)
 116         E = 2.910383045673370361328125000000000000000e-11L,             /* 2**-35 */
 117         Einv = 3.435973836800000000000000000000000000000e+10L;  /* 2**+35 */
 118 #else
 119         E = 8.673617379884035472059622406959533691406e-19L,             /* 2**-60 */
 120         Einv = 1.152921504606846976000000000000000000000e18L;           /* 2**+60 */
 121 #endif
 122 /* END CSTYLED */
 123 
 124 ldcomplex
 125 catanl(ldcomplex z)
 126 {
 127         ldcomplex ans;
 128         long double x, y, t1, ax, ay, t;
 129         int hx, hy, ix, iy;
 130 
 131         x = LD_RE(z);
 132         y = LD_IM(z);
 133         ax = fabsl(x);
 134         ay = fabsl(y);
 135         hx = HI_XWORD(x);
 136         hy = HI_XWORD(y);
 137         ix = hx & 0x7fffffff;
 138         iy = hy & 0x7fffffff;
 139 
 140         /* x is inf or NaN */
 141         if (ix >= 0x7fff0000) {
 142                 if (isinfl(x)) {
 143                         LD_RE(ans) = pi_2;
 144                         LD_IM(ans) = zero;
 145                 } else {
 146                         LD_RE(ans) = x + x;
 147 
 148                         if (iszerol(y) || (isinfl(y)))
 149                                 LD_IM(ans) = zero;
 150                         else
 151                                 LD_IM(ans) = (fabsl(y) - ay) / (fabsl(y) - ay);
 152                 }
 153         } else if (iy >= 0x7fff0000) {
 154                 /* y is inf or NaN */
 155                 if (isinfl(y)) {
 156                         LD_RE(ans) = pi_2;
 157                         LD_IM(ans) = zero;
 158                 } else {
 159                         LD_RE(ans) = (fabsl(x) - ax) / (fabsl(x) - ax);
 160                         LD_IM(ans) = y;
 161                 }
 162         } else if (iszerol(x)) {
 163                 /* BEGIN CSTYLED */
 164                 /*
 165                  * x = 0
 166                  *      1                            1
 167                  * A = --- * atan2(2x, 1-x*x-y*y) = --- atan2(0,1-|y|)
 168                  *      2                            2
 169                  *
 170                  *     1     [ (y+1)*(y+1) ]   1          2      1         2y
 171                  * B = - log [ ----------- ] = - log (1+ ---) or - log(1+ ----)
 172                  *     4     [ (y-1)*(y-1) ]   2         y-1     2         1-y
 173                  */
 174                 /* END CSTYLED */
 175                 t = one - ay;
 176 
 177                 if (ay == one) {
 178                         /* y=1: catan(0,1)=(0,+inf) with 1/0 signal */
 179                         LD_IM(ans) = ay / ax;
 180                         LD_RE(ans) = zero;
 181                 } else if (ay > one) {       /* y>1 */
 182                         LD_IM(ans) = half * log1pl(two / (-t));
 183                         LD_RE(ans) = pi_2;
 184                 } else {                /* y<1 */
 185                         LD_IM(ans) = half * log1pl((ay + ay) / t);
 186                         LD_RE(ans) = zero;
 187                 }
 188         } else if (ay < E * (one + ax)) {
 189                 /* BEGIN CSTYLED */
 190                 /*
 191                  * Tiny y (relative to 1+|x|)
 192                  *     |y| < E*(1+|x|)
 193                  * where E=2**-29, -35, -60 for double, extended, quad precision
 194                  *
 195                  *     1                         [x<=1:   atan(x)
 196                  * A = - * atan2(2x,1-x*x-y*y) ~ [      1                 1+x
 197                  *     2                         [x>=1: - atan2(2,(1-x)*(-----))
 198                  *                                      2                  x
 199                  *
 200                  *                               y/x
 201                  * B ~ t*(1-2t), where t = ----------------- is tiny
 202                  *                         x + (y-1)*(y-1)/x
 203                  *
 204                  *                           y
 205                  * (when x < 2**-60, t = ----------- )
 206                  *                       (y-1)*(y-1)
 207                  */
 208                 /* END CSTYLED */
 209                 if (ay == zero) {
 210                         LD_IM(ans) = ay;
 211                 } else {
 212                         t1 = ay - one;
 213 
 214                         if (ix < 0x3fc30000)
 215                                 t = ay / (t1 * t1);
 216                         else if (ix > 0x403b0000)
 217                                 t = (ay / ax) / ax;
 218                         else
 219                                 t = ay / (ax * ax + t1 * t1);
 220 
 221                         LD_IM(ans) = t * (one - two * t);
 222                 }
 223 
 224                 if (ix < 0x3fff0000)
 225                         LD_RE(ans) = atanl(ax);
 226                 else
 227                         LD_RE(ans) = half * atan2l(two, (one - ax) * (one +
 228                             one / ax));
 229         } else if (ay > Einv * (one + ax)) {
 230                 /* BEGIN CSTYLED */
 231                 /*
 232                  * Huge y relative to 1+|x|
 233                  *     |y| > Einv*(1+|x|), where Einv~2**(prec/2+3),
 234                  *      1
 235                  * A ~ --- * atan2(2x, -y*y) ~ pi/2
 236                  *      2
 237                  *                               y
 238                  * B ~ t*(1-2t), where t = --------------- is tiny
 239                  *                          (y-1)*(y-1)
 240                  */
 241                 /* END CSTYLED */
 242                 LD_RE(ans) = pi_2;
 243                 t = (ay / (ay - one)) / (ay - one);
 244                 LD_IM(ans) = t * (one - (t + t));
 245         } else if (ay == one) {
 246                 /* BEGIN CSTYLED */
 247                 /*
 248                  * y=1
 249                  *     1                      1
 250                  * A = - * atan2(2x, -x*x) = --- atan2(2,-x)
 251                  *     2                      2
 252                  *
 253                  *     1     [ x*x+4]   1          4     [ 0.5(log2-logx) if
 254                  * B = - log [ -----] = - log (1+ ---) = [ |x|<E, else 0.25*
 255                  *     4     [  x*x ]   4         x*x    [ log1p((2/x)*(2/x))
 256                  */
 257                 /* END CSTYLED */
 258                 LD_RE(ans) = half * atan2l(two, -ax);
 259 
 260                 if (ax < E) {
 261                         LD_IM(ans) = half * (ln2 - logl(ax));
 262                 } else {
 263                         t = two / ax;
 264                         LD_IM(ans) = 0.25L * log1pl(t * t);
 265                 }
 266         } else if (ax > Einv * Einv) {
 267                 /* BEGIN CSTYLED */
 268                 /*
 269                  * Huge x:
 270                  * when |x| > 1/E^2,
 271                  *      1                           pi
 272                  * A ~ --- * atan2(2x, -x*x-y*y) ~ ---
 273                  *      2                           2
 274                  *                               y                 y/x
 275                  * B ~ t*(1-2t), where t = --------------- = (-------------- )/x
 276                  *                         x*x+(y-1)*(y-1)     1+((y-1)/x)^2
 277                  */
 278                 /* END CSTYLED */
 279                 LD_RE(ans) = pi_2;
 280                 t = ((ay / ax) / (one + ((ay - one) / ax) * ((ay - one) /
 281                     ax))) / ax;
 282                 LD_IM(ans) = t * (one - (t + t));
 283         } else if (ax < E * E * E * E) {
 284                 /* BEGIN CSTYLED */
 285                 /*
 286                  * Tiny x:
 287                  * when |x| < E^4,  (note that y != 1)
 288                  *      1                            1
 289                  * A = --- * atan2(2x, 1-x*x-y*y) ~ --- * atan2(2x,1-y*y)
 290                  *      2                            2
 291                  *
 292                  *     1     [ (y+1)*(y+1) ]   1          2      1         2y
 293                  * B = - log [ ----------- ] = - log (1+ ---) or - log(1+ ----)
 294                  *     4     [ (y-1)*(y-1) ]   2         y-1     2         1-y
 295                  */
 296                 /* END CSTYLED */
 297                 LD_RE(ans) = half * atan2l(ax + ax, (one - ay) * (one + ay));
 298 
 299                 if (ay > one)                /* y>1 */
 300                         LD_IM(ans) = half * log1pl(two / (ay - one));
 301                 else                    /* y<1 */
 302                         LD_IM(ans) = half * log1pl((ay + ay) / (one - ay));
 303         } else {
 304                 /* BEGIN CSTYLED */
 305                 /*
 306                  * normal x,y
 307                  *      1
 308                  * A = --- * atan2(2x, 1-x*x-y*y)
 309                  *      2
 310                  *
 311                  *     1     [ x*x+(y+1)*(y+1) ]   1               4y
 312                  * B = - log [ --------------- ] = - log (1+ -----------------)
 313                  *     4     [ x*x+(y-1)*(y-1) ]   4         x*x + (y-1)*(y-1)
 314                  */
 315                 /* END CSTYLED */
 316                 t = one - ay;
 317 
 318                 if (iy >= 0x3ffe0000 && iy < 0x40000000) {
 319                         /* y close to 1 */
 320                         LD_RE(ans) = half * (atan2l((ax + ax), (t * (one + ay) -
 321                             ax * ax)));
 322                 } else if (ix >= 0x3ffe0000 && ix < 0x40000000) {
 323                         /* x close to 1 */
 324                         LD_RE(ans) = half * atan2l((ax + ax), ((one - ax) *
 325                             (one + ax) - ay * ay));
 326                 } else {
 327                         LD_RE(ans) = half * atan2l((ax + ax), ((one - ax * ax) -
 328                             ay * ay));
 329                 }
 330 
 331                 LD_IM(ans) = 0.25L * log1pl((4.0L * ay) / (ax * ax + t * t));
 332         }
 333 
 334         if (hx < 0)
 335                 LD_RE(ans) = -LD_RE(ans);
 336 
 337         if (hy < 0)
 338                 LD_IM(ans) = -LD_IM(ans);
 339 
 340         return (ans);
 341 }