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