usr/src/lib/libm/common/complex/csinh.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 __csinh = csinh
  31
  32 /* INDENT OFF */
  33 /*
  34  * dcomplex csinh(dcomplex z);
  35  *
  36  *             z      -z       x                      -x
  37  *            e   -  e        e  (cos(y)+i*sin(y)) - e  (cos(-y)+i*sin(-y))
  38  * sinh z = -------------- =  ---------------------------------------------
  39  *                2                                2
  40  *                     x    -x                x    -x
  41  *           cos(y) ( e  - e  )  + i*sin(y) (e  + e   )
  42  *        = --------------------------------------------
  43  *                               2
  44  *
  45  *        =  cos(y) sinh(x)  + i sin(y) cosh(x)
  46  *
  47  * Implementation Note
  48  * -------------------
  49  *
  50  *             |x|    -|x|   |x|        -2|x|       -2|x|    -P-4
  51  * Note that  e   +- e    = e   ( 1 +- e     ). If e      < 2     , where
  52  *
  53  * P stands for the number of significant bits of the machine precision,
  54  *                                     |x|
  55  * then the result will be rounded to e   . Therefore, we have
  56  *
  57  *                 z
  58  *                e
  59  *     sinh z = -----  if |x| >= (P/2 + 2)*ln2
  60  *                2
  61  *
  62  * EXCEPTION (conform to ISO/IEC 9899:1999(E)):
  63  *      csinh(0,0)=(0,0)
  64  *      csinh(0,inf)=(+-0,NaN)
  65  *      csinh(0,NaN)=(+-0,NaN)
  66  *      csinh(x,inf) = (NaN,NaN) for finite positive x
  67  *      csinh(x,NaN) = (NaN,NaN) for finite non-zero x
  68  *      csinh(inf,0) = (inf, 0)
  69  *      csinh(inf,y) = (inf*cos(y),inf*sin(y)) for positive finite y
  70  *      csinh(inf,inf) = (+-inf,NaN)
  71  *      csinh(inf,NaN) = (+-inf,NaN)
  72  *      csinh(NaN,0) = (NaN,0)
  73  *      csinh(NaN,y) = (NaN,NaN) for non-zero y
  74  *      csinh(NaN,NaN) = (NaN,NaN)
  75  */
  76 /* INDENT ON */
  77
  78 #include "libm.h"               /* cosh/exp/fabs/scalbn/sinh/sincos/__k_cexp */
  79 #include "complex_wrapper.h"
  80
  81 dcomplex
  82 csinh(dcomplex z) {
  83         double t, x, y, S, C;
  84         int hx, ix, lx, hy, iy, ly, n;
  85         dcomplex ans;
  86
  87         x = D_RE(z);
  88         y = D_IM(z);
  89         hx = HI_WORD(x);
  90         lx = LO_WORD(x);
  91         ix = hx & 0x7fffffff;
  92         hy = HI_WORD(y);
  93         ly = LO_WORD(y);
  94         iy = hy & 0x7fffffff;
  95         x = fabs(x);
  96         y = fabs(y);
  97
  98         (void) sincos(y, &S, &C);
  99         if (ix >= 0x403c0000) { /* |x| > 28 = prec/2 (14,28,34,60) */
 100                 if (ix >= 0x40862E42) { /* |x| > 709.78... ~ log(2**1024) */
 101                         if (ix >= 0x7ff00000) { /* |x| is inf or NaN */
 102                                 if ((iy | ly) == 0) {
 103                                         D_RE(ans) = x;
 104                                         D_IM(ans) = y;
 105                                 } else if (iy >= 0x7ff00000) {
 106                                         D_RE(ans) = x;
 107                                         D_IM(ans) = x - y;
 108                                 } else {
 109                                         D_RE(ans) = C * x;
 110                                         D_IM(ans) = S * x;
 111                                 }
 112                         } else {
 113                                 /* return exp(x)=t*2**n */
 114                                 t = __k_cexp(x, &n);
 115                                 D_RE(ans) = scalbn(C * t, n - 1);
 116                                 D_IM(ans) = scalbn(S * t, n - 1);
 117                         }
 118                 } else {
 119                         t = exp(x) * 0.5;
 120                         D_RE(ans) = C * t;
 121                         D_IM(ans) = S * t;
 122                 }
 123         } else {
 124                 if ((ix | lx) == 0) {   /* x = 0, return (0,S) */
 125                         D_RE(ans) = 0.0;
 126                         D_IM(ans) = S;
 127                 } else {
 128                         D_RE(ans) = C * sinh(x);
 129                         D_IM(ans) = S * cosh(x);
 130                 }
 131         }
 132         if (hx < 0)
 133                 D_RE(ans) = -D_RE(ans);
 134         if (hy < 0)
 135                 D_IM(ans) = -D_IM(ans);
 136         return (ans);
 137 }