blob: 113d5c103c03c733db433f693a6ead597615eb8f [file] [log] [blame]
Charles L. Doriana0690b62010-02-01 22:21:40 -08001// Copyright 2010 The Go Authors. All rights reserved.
2// Use of this source code is governed by a BSD-style
3// license that can be found in the LICENSE file.
4
5package math
6
Charles L. Doriana0690b62010-02-01 22:21:40 -08007// The original C code, the long comment, and the constants
8// below are from FreeBSD's /usr/src/lib/msun/src/e_atanh.c
9// and came with this notice. The go code is a simplified
10// version of the original C.
11//
12// ====================================================
13// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
14//
15// Developed at SunPro, a Sun Microsystems, Inc. business.
16// Permission to use, copy, modify, and distribute this
17// software is freely granted, provided that this notice
18// is preserved.
19// ====================================================
20//
21//
22// __ieee754_atanh(x)
23// Method :
24// 1. Reduce x to positive by atanh(-x) = -atanh(x)
25// 2. For x>=0.5
26// 1 2x x
27// atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
28// 2 1 - x 1 - x
29//
30// For x<0.5
31// atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
32//
33// Special cases:
34// atanh(x) is NaN if |x| > 1 with signal;
35// atanh(NaN) is that NaN with no signal;
36// atanh(+-1) is +-INF with signal.
37//
38
Charles L. Dorianf2734872012-04-06 14:01:12 -040039// Atanh returns the inverse hyperbolic tangent of x.
Charles L. Doriana0690b62010-02-01 22:21:40 -080040//
41// Special cases are:
Charles L. Doriana0690b62010-02-01 22:21:40 -080042// Atanh(1) = +Inf
Charles L. Dorian94b03422011-12-08 17:07:13 -050043// Atanh(±0) = ±0
Charles L. Doriana0690b62010-02-01 22:21:40 -080044// Atanh(-1) = -Inf
Charles L. Dorian94b03422011-12-08 17:07:13 -050045// Atanh(x) = NaN if x < -1 or x > 1
Charles L. Doriana0690b62010-02-01 22:21:40 -080046// Atanh(NaN) = NaN
47func Atanh(x float64) float64 {
Charles L. Dorian3c3e68b2010-04-09 14:37:33 -070048 const NearZero = 1.0 / (1 << 28) // 2**-28
Charles L. Doriana0690b62010-02-01 22:21:40 -080049 // special cases
50 switch {
Luuk van Dijk8dd3de42012-02-01 16:08:31 +010051 case x < -1 || x > 1 || IsNaN(x):
Charles L. Doriana0690b62010-02-01 22:21:40 -080052 return NaN()
53 case x == 1:
54 return Inf(1)
55 case x == -1:
56 return Inf(-1)
57 }
58 sign := false
59 if x < 0 {
60 x = -x
61 sign = true
62 }
63 var temp float64
64 switch {
65 case x < NearZero:
66 temp = x
67 case x < 0.5:
68 temp = x + x
69 temp = 0.5 * Log1p(temp+temp*x/(1-x))
70 default:
71 temp = 0.5 * Log1p((x+x)/(1-x))
72 }
73 if sign {
74 temp = -temp
75 }
76 return temp
77}