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