usr/src/lib/libm/common/complex/catan.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 __catan = catan
  31
  32 /* INDENT OFF */
  33 /*
  34  * dcomplex catan(dcomplex z);
  35  *
  36  * If
  37  *     z = x + iy,
  38  *
  39  * then
  40  *          1       (    2x     )    1                2    2
  41  * Re w  =  - arctan(-----------)  = - ATAN2(2x, 1 - x  - y )
  42  *          2       (     2    2)    2
  43  *                  (1 - x  - y )
  44  *
  45  *               ( 2         2)
  46  *          1    (x  +  (y+1) )      1                  4y
  47  * Im w  =  - log(------------) .=  --- log [ 1 + ------------- ]
  48  *          4    ( 2         2)      4              2         2
  49  *               (x  +  (y-1) )                    x  +  (y-1)
  50  *
  51  *                 2    16  3                         y
  52  *         = t - 2t   + -- t  - ..., where t = -----------------
  53  *                      3                      x*x + (y-1)*(y-1)
  54  *
  55  * Note that: if catan( x, y) = ( u, v), then
  56  *               catan(-x, y) = (-u, v)
  57  *               catan( x,-y) = ( u,-v)
  58  *
  59  * Also,   catan(x,y) = -i*catanh(-y,x), or
  60  *        catanh(x,y) =  i*catan(-y,x)
  61  * So, if catanh(y,x) = (v,u), then catan(x,y) = -i*(-v,u) = (u,v), i.e.,
  62  *        catan(x,y) = (u,v)
  63  *
  64  * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
  65  *    catan( 0  , 0   ) =  (0    ,  0   )
  66  *    catan( NaN, 0   ) =  (NaN  ,  0   )
  67  *    catan( 0  , 1   ) =  (0    ,  +inf) with divide-by-zero
  68  *    catan( inf, y   ) =  (pi/2 ,  0   ) for finite +y
  69  *    catan( NaN, y   ) =  (NaN  ,  NaN ) with invalid for finite y != 0
  70  *    catan( x  , inf ) =  (pi/2 ,  0   ) for finite +x
  71  *    catan( inf, inf ) =  (pi/2 ,  0   )
  72  *    catan( NaN, inf ) =  (NaN  ,  0   )
  73  *    catan( x  , NaN ) =  (NaN  ,  NaN ) with invalid for finite x
  74  *    catan( inf, NaN ) =  (pi/2 ,  +-0 )
  75  */
  76 /* INDENT ON */
  77
  78 #include "libm.h"               /* atan/atan2/fabs/log/log1p */
  79 #include "complex_wrapper.h"
  80
  81 /* INDENT OFF */
  82 static const double
  83         pi_2 = 1.570796326794896558e+00,
  84         zero = 0.0,
  85         half = 0.5,
  86         two = 2.0,
  87         ln2 = 6.931471805599453094172321214581765680755e-0001,
  88         one = 1.0;
  89 /* INDENT ON */
  90
  91 dcomplex
  92 catan(dcomplex z) {
  93         dcomplex ans;
  94         double x, y, ax, ay, t;
  95         int hx, hy, ix, iy;
  96         unsigned lx, ly;
  97
  98         x = D_RE(z);
  99         y = D_IM(z);
 100         ax = fabs(x);
 101         ay = fabs(y);
 102         hx = HI_WORD(x);
 103         lx = LO_WORD(x);
 104         hy = HI_WORD(y);
 105         ly = LO_WORD(y);
 106         ix = hx & 0x7fffffff;
 107         iy = hy & 0x7fffffff;
 108
 109         /* x is inf or NaN */
 110         if (ix >= 0x7ff00000) {
 111                 if (ISINF(ix, lx)) {
 112                         D_RE(ans) = pi_2;
 113                         D_IM(ans) = zero;
 114                 } else {
 115                         D_RE(ans) = x + x;
 116                         if ((iy | ly) == 0 || (ISINF(iy, ly)))
 117                                 D_IM(ans) = zero;
 118                         else
 119                                 D_IM(ans) = (fabs(y) - ay) / (fabs(y) - ay);
 120                 }
 121         } else if (iy >= 0x7ff00000) {
 122                 /* y is inf or NaN */
 123                 if (ISINF(iy, ly)) {
 124                         D_RE(ans) = pi_2;
 125                         D_IM(ans) = zero;
 126                 } else {
 127                         D_RE(ans) = (fabs(x) - ax) / (fabs(x) - ax);
 128                         D_IM(ans) = y;
 129                 }
 130         } else if ((ix | lx) == 0) {
 131                 /* INDENT OFF */
 132                 /*
 133                  * x = 0
 134                  *      1                            1
 135                  * A = --- * atan2(2x, 1-x*x-y*y) = --- atan2(0,1-|y|)
 136                  *      2                            2
 137                  *
 138                  *     1     [  (y+1)*(y+1) ]   1          2      1         2y
 139                  * B = - log [ ------------ ] = - log (1+ ---) or - log(1+ ----)
 140                  *     4     [  (y-1)*(y-1) ]   2         y-1     2         1-y
 141                  */
 142                 /* INDENT ON */
 143                 t = one - ay;
 144                 if (((iy - 0x3ff00000) | ly) == 0) {
 145                         /* y=1: catan(0,1)=(0,+inf) with 1/0 signal */
 146                         D_IM(ans) = ay / ax;
 147                         D_RE(ans) = zero;
 148                 } else if (iy >= 0x3ff00000) {  /* y>1 */
 149                         D_IM(ans) = half * log1p(two / (-t));
 150                         D_RE(ans) = pi_2;
 151                 } else {                /* y<1 */
 152                         D_IM(ans) = half * log1p((ay + ay) / t);
 153                         D_RE(ans) = zero;
 154                 }
 155         } else if (iy < 0x3e200000 || ((ix - iy) >> 20) >= 30) {
 156         /* INDENT OFF */
 157         /*
 158          * Tiny y (relative to 1+|x|)
 159          *     |y| < E*(1+|x|)
 160          * where E=2**-29, -35, -60 for double, double extended, quad precision
 161          *
 162          *      1                           [ x<=1:   atan(x)
 163          * A = --- * atan2(2x, 1-x*x-y*y) ~ [       1                 1+x
 164          *      2                           [ x>=1: - atan2(2,(1-x)*(-----))
 165          *                                          2                  x
 166          *
 167          *                               y/x
 168          * B ~ t*(1-2t), where t = ----------------- is tiny
 169          *                         x + (y-1)*(y-1)/x
 170          */
 171                 /* INDENT ON */
 172                 if (ix < 0x3ff00000)
 173                         D_RE(ans) = atan(ax);
 174                 else
 175                         D_RE(ans) = half * atan2(two, (one - ax) * (one +
 176                                 one / ax));
 177                 if ((iy | ly) == 0) {
 178                         D_IM(ans) = ay;
 179                 } else {
 180                         if (ix < 0x3e200000)
 181                                 t = ay / ((ay - one) * (ay - one));
 182                         else if (ix > 0x41c00000)
 183                                 t = (ay / ax) / ax;
 184                         else
 185                                 t = ay / (ax * ax + (ay - one) * (ay - one));
 186                         D_IM(ans) = t * (one - (t + t));
 187                 }
 188         } else if (iy >= 0x41c00000 && ((iy - ix) >> 20) >= 30) {
 189                 /* INDENT OFF */
 190                 /*
 191                  * Huge y relative to 1+|x|
 192                  *            |y| > Einv*(1+|x|), where Einv~2**(prec/2+3),
 193                  *            1
 194                  *       A ~ --- * atan2(2x, -y*y) ~ pi/2
 195                  *            2
 196                  *                                     y
 197                  *       B ~ t*(1-2t), where t = --------------- is tiny
 198                  *                                (y-1)*(y-1)
 199                  */
 200                 /* INDENT ON */
 201                 D_RE(ans) = pi_2;
 202                 t = (ay / (ay - one)) / (ay - one);
 203                 D_IM(ans) = t * (one - (t + t));
 204         } else if (((iy - 0x3ff00000) | ly) == 0) {
 205                 /* INDENT OFF */
 206                 /*
 207                  * y = 1
 208                  *      1                       1
 209                  * A = --- * atan2(2x, -x*x) = --- atan2(2,-x)
 210                  *      2                       2
 211                  *
 212                  *     1     [x*x + 4]   1          4     [ 0.5(log2-logx) if
 213                  * B = - log [-------] = - log (1+ ---) = [ |x|<E, else 0.25*
 214                  *     4     [  x*x  ]   4         x*x    [ log1p((2/x)*(2/x))
 215                  */
 216                 /* INDENT ON */
 217                 D_RE(ans) = half * atan2(two, -ax);
 218                 if (ix < 0x3e200000)
 219                         D_IM(ans) = half * (ln2 - log(ax));
 220                 else {
 221                         t = two / ax;
 222                         D_IM(ans) = 0.25 * log1p(t * t);
 223                 }
 224         } else if (ix >= 0x43900000) {
 225                 /* INDENT OFF */
 226                 /*
 227                  * Huge x:
 228                  * when |x| > 1/E^2,
 229                  *      1                           pi
 230                  * A ~ --- * atan2(2x, -x*x-y*y) ~ ---
 231                  *      2                           2
 232                  *                               y                 y/x
 233                  * B ~ t*(1-2t), where t = --------------- = (-------------- )/x
 234                  *                         x*x+(y-1)*(y-1)     1+((y-1)/x)^2
 235                  */
 236                 /* INDENT ON */
 237                 D_RE(ans) = pi_2;
 238                 t = ((ay / ax) / (one + ((ay - one) / ax) * ((ay - one) /
 239                         ax))) / ax;
 240                 D_IM(ans) = t * (one - (t + t));
 241         } else if (ix < 0x38b00000) {
 242                 /* INDENT OFF */
 243                 /*
 244                  * Tiny x:
 245                  * when |x| < E^4,  (note that y != 1)
 246                  *      1                            1
 247                  * A = --- * atan2(2x, 1-x*x-y*y) ~ --- * atan2(2x,(1-y)*(1+y))
 248                  *      2                            2
 249                  *
 250                  *     1     [(y+1)*(y+1)]   1          2      1         2y
 251                  * B = - log [-----------] = - log (1+ ---) or - log(1+ ----)
 252                  *     4     [(y-1)*(y-1)]   2         y-1     2         1-y
 253                  */
 254                 /* INDENT ON */
 255                 D_RE(ans) = half * atan2(ax + ax, (one - ay) * (one + ay));
 256                 if (iy >= 0x3ff00000)
 257                         D_IM(ans) = half * log1p(two / (ay - one));
 258                 else
 259                         D_IM(ans) = half * log1p((ay + ay) / (one - ay));
 260         } else {
 261                 /* INDENT OFF */
 262                 /*
 263                  * normal x,y
 264                  *      1
 265                  * A = --- * atan2(2x, 1-x*x-y*y)
 266                  *      2
 267                  *
 268                  *     1     [x*x+(y+1)*(y+1)]   1               4y
 269                  * B = - log [---------------] = - log (1+ -----------------)
 270                  *     4     [x*x+(y-1)*(y-1)]   4         x*x + (y-1)*(y-1)
 271                  */
 272                 /* INDENT ON */
 273                 t = one - ay;
 274                 if (iy >= 0x3fe00000 && iy < 0x40000000) {
 275                         /* y close to 1 */
 276                         D_RE(ans) = half * (atan2((ax + ax), (t * (one + ay) -
 277                                 ax * ax)));
 278                 } else if (ix >= 0x3fe00000 && ix < 0x40000000) {
 279                         /* x close to 1 */
 280                         D_RE(ans) = half * atan2((ax + ax), ((one - ax) *
 281                                 (one + ax) - ay * ay));
 282                 } else
 283                         D_RE(ans) = half * atan2((ax + ax), ((one - ax * ax) -
 284                                 ay * ay));
 285                 D_IM(ans) = 0.25 * log1p((4.0 * ay) / (ax * ax + t * t));
 286         }
 287         if (hx < 0)
 288                 D_RE(ans) = -D_RE(ans);
 289         if (hy < 0)
 290                 D_IM(ans) = -D_IM(ans);
 291         return (ans);
 292 }