amath  1.8.5
Simple command line calculator
atanh.c
Go to the documentation of this file.
1 /*-
2  * Copyright (c) 2014-2018 Carsten Sonne Larsen <cs@innolan.net>
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions
7  * are met:
8  * 1. Redistributions of source code must retain the above copyright
9  * notice, this list of conditions and the following disclaimer.
10  * 2. Redistributions in binary form must reproduce the above copyright
11  * notice, this list of conditions and the following disclaimer in the
12  * documentation and/or other materials provided with the distribution.
13  *
14  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
15  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
16  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
17  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
18  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
19  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
20  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
21  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
22  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
23  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24  *
25  * Project homepage:
26  * https://amath.innolan.net
27  *
28  * The original source code can be obtained from:
29  * http://www.netlib.org/fdlibm/e_atanh.c
30  *
31  * =================================================================
32  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
33  *
34  * Developed at SunSoft, a Sun Microsystems, Inc. business.
35  * Permission to use, copy, modify, and distribute this
36  * software is freely granted, provided that this notice
37  * is preserved.
38  * =================================================================
39  */
40 
41 /**
42  * @file atanh.c
43  * @brief Inverse hyperbolic tangent function
44  */
45 
46 #include "prim.h"
47 
48 static const double one = 1.0, huge = 1e300;
49 static double zero = 0.0;
50 
51 /**
52  * @brief Inverse hyperbolic tangent function
53  * @details
54  * <pre>
55  * Method
56  *
57  * 1.Reduced x to positive by atanh(-x) = -atanh(x)
58  * 2.For x>=0.5
59  * 1 2x x
60  * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
61  * 2 1 - x 1 - x
62  *
63  * For x<0.5
64  * atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
65  *
66  * Special cases
67  * atanh(x) is NaN if |x| > 1
68  * atanh(NaN) is that NaN
69  * atanh(+-1) is +-INF
70  * </pre>
71  */
72 double atanh(double x)
73 {
74  double t;
75  int32_t hx, ix;
76  uint32_t lx;
77  GET_HIGH_WORD(hx, x);
78  GET_LOW_WORD(lx, x);
79  ix = hx & 0x7FFFFFFF;
80 
81  // |x| > 1
82  if ((ix | ((lx | (-lx)) >> 31)) > 0x3FF00000)
83  {
84  return NAN;
85  }
86 
87  if (ix == 0x3FF00000)
88  return x / zero;
89  if (ix < 0x3E300000 && (huge + x) > zero)
90  return x; /* x<2**-28 */
91  SET_HIGH_WORD(x, ix); /* x <- |x| */
92  if (ix < 0x3FE00000)
93  { /* x < 0.5 */
94  t = x + x;
95  t = 0.5 * log1p(t + t * x / (one - x));
96  }
97  else
98  t = 0.5 * log1p((x + x) / (one - x));
99 
100  if (hx >= 0)
101  {
102  return t;
103  }
104 
105  return -t;
106 }
#define GET_HIGH_WORD(i, d)
Get the more significant 32 bit int from a double.
Definition: prim.h:167
double log1p(double x)
Definition: log1p.c:122
#define GET_LOW_WORD(i, d)
Get the less significant 32 bit int from a double.
Definition: prim.h:177
static double zero
Definition: atanh.c:49
double atanh(double x)
Inverse hyperbolic tangent function.
Definition: atanh.c:72
#define NAN
Definition: mathr.h:53
static const double huge
Definition: atanh.c:48
static const double one
Definition: atanh.c:48
#define SET_HIGH_WORD(d, v)
Set the more significant 32 bits of a double from an int.
Definition: prim.h:198