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 __ctanh = ctanh
  32 
  33 
  34 /*
  35  * dcomplex ctanh(dcomplex z);
  36  *
  37  *            tanh x  + i tan y      sinh 2x  +  i sin 2y
  38  * ctanh z = --------------------- = --------------------
  39  *           1 + i tanh(x)tan(y)       cosh 2x + cos 2y
  40  *
  41  * For |x| >= prec/2 (14,28,34,60 for single, double, double extended, quad),
  42  * we use
  43  *
  44  *                         1   2x                              2 sin 2y
  45  *    cosh 2x = sinh 2x = --- e    and hence  ctanh z = 1 + i -----------;
  46  *                         2                                       2x
  47  *                                                                e
  48  *
  49  * otherwise, to avoid cancellation, for |x| < prec/2,
  50  *                              2x     2
  51  *                            (e   - 1)        2       2
  52  *    cosh 2x + cos 2y = 1 + ------------ + cos y - sin y
  53  *                                 2x
  54  *                              2 e
  55  *
  56  *                        1    2x     2  -2x         2
  57  *                     = --- (e   - 1)  e     + 2 cos y
  58  *                        2
  59  * and
  60  *
  61  *                  [            2x      ]
  62  *               1  [  2x       e   - 1  ]
  63  *    sinh 2x = --- [ e  - 1 + --------- ]
  64  *               2  [               2x   ]
  65  *                  [              e     ]
  66  *                                             2x
  67  * Implementation notes:  let t = expm1(2x) = e   - 1, then
  68  *
  69  *                     1    [  t*t         2  ]              1    [      t  ]
  70  * cosh 2x + cos 2y = --- * [ ----- + 4 cos y ];  sinh 2x = --- * [ t + --- ]
  71  *                     2    [  t+1            ]              2    [     t+1 ]
  72  *
  73  * Hence,
  74  *
  75  *
  76  *                        t*t+2t                  [4(t+1)(cos y)]*(sin y)
  77  *    ctanh z = --------------------------- + i --------------------------
  78  *               t*t+[4(t+1)(cos y)](cos y)     t*t+[4(t+1)(cos y)](cos y)
  79  *
  80  * EXCEPTION (conform to ISO/IEC 9899:1999(E)):
  81  *      ctanh(0,0)=(0,0)
  82  *      ctanh(x,inf) = (NaN,NaN) for finite x
  83  *      ctanh(x,NaN) = (NaN,NaN) for finite x
  84  *      ctanh(inf,y) = 1+ i*0*sin(2y) for positive-signed finite y
  85  *      ctanh(inf,inf) = (1, +-0)
  86  *      ctanh(inf,NaN) = (1, +-0)
  87  *      ctanh(NaN,0) = (NaN,0)
  88  *      ctanh(NaN,y) = (NaN,NaN) for non-zero y
  89  *      ctanh(NaN,NaN) = (NaN,NaN)
  90  */
  91 
  92 #include "libm.h"                       /* exp/expm1/fabs/sin/tanh/sincos */
  93 #include "complex_wrapper.h"
  94 
  95 static const double four = 4.0, two = 2.0, one = 1.0, zero = 0.0;
  96 
  97 dcomplex
  98 ctanh(dcomplex z)
  99 {
 100         double t, r, v, u, x, y, S, C;
 101         int hx, ix, lx, hy, iy, ly;
 102         dcomplex ans;
 103 
 104         x = D_RE(z);
 105         y = D_IM(z);
 106         hx = HI_WORD(x);
 107         lx = LO_WORD(x);
 108         ix = hx & 0x7fffffff;
 109         hy = HI_WORD(y);
 110         ly = LO_WORD(y);
 111         iy = hy & 0x7fffffff;
 112         x = fabs(x);
 113         y = fabs(y);
 114 
 115         if ((iy | ly) == 0) {   /* ctanh(x,0) = (x,0) for x = 0 or NaN */
 116                 D_RE(ans) = tanh(x);
 117                 D_IM(ans) = zero;
 118         } else if (iy >= 0x7ff00000) { /* y is inf or NaN */
 119                 if (ix < 0x7ff00000) { /* catanh(finite x,inf/nan) is nan */
 120                         D_RE(ans) = D_IM(ans) = y - y;
 121                 } else if (((ix - 0x7ff00000) | lx) == 0) {     /* x is inf */
 122                         D_RE(ans) = one;
 123                         D_IM(ans) = zero;
 124                 } else {
 125                         D_RE(ans) = x + y;
 126                         D_IM(ans) = y - y;
 127                 }
 128         } else if (ix >= 0x403c0000) {
 129                 /*
 130                  * |x| > 28 = prec/2 (14,28,34,60)
 131                  * ctanh z ~ 1 + i (sin2y)/(exp(2x))
 132                  */
 133                 D_RE(ans) = one;
 134 
 135                 if (iy < 0x7fe00000) {       /* t = sin(2y) */
 136                         S = sin(y + y);
 137                 } else {
 138                         (void) sincos(y, &S, &C);
 139                         S = (S + S) * C;
 140                 }
 141 
 142                 if (ix >= 0x7fe00000) {              /* |x| > max/2 */
 143                         if (ix >= 0x7ff00000) {      /* |x| is inf or NaN */
 144                                 if (((ix - 0x7ff00000) | lx) != 0)
 145                                         D_RE(ans) = D_IM(ans) = x + y;
 146                                 /* x is NaN */
 147                                 else
 148                                         D_IM(ans) = zero * S;   /* x is inf */
 149                         } else {
 150                                 D_IM(ans) = S * exp(-x);        /* underflow */
 151                         }
 152                 } else {
 153                         D_IM(ans) = (S + S) * exp(-(x + x));
 154                 }
 155 
 156                 /* 2 sin 2y / exp(2x) */
 157         } else {
 158                 /* BEGIN CSTYLED */
 159                 /*
 160                  *                        t*t+2t
 161                  *    ctanh z = --------------------------- +
 162                  *               t*t+[4(t+1)(cos y)](cos y)
 163                  *
 164                  *                  [4(t+1)(cos y)]*(sin y)
 165                  *              i --------------------------
 166                  *                t*t+[4(t+1)(cos y)](cos y)
 167                  */
 168                 /* END CSTYLED */
 169                 (void) sincos(y, &S, &C);
 170                 t = expm1(x + x);
 171                 r = (four * C) * (t + one);
 172                 u = t * t;
 173                 v = one / (u + r * C);
 174                 D_RE(ans) = (u + two * t) * v;
 175                 D_IM(ans) = (r * S) * v;
 176         }
 177 
 178         if (hx < 0)
 179                 D_RE(ans) = -D_RE(ans);
 180 
 181         if (hy < 0)
 182                 D_IM(ans) = -D_IM(ans);
 183 
 184         return (ans);
 185 }