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