Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 1 | // 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 | |
| 5 | package math |
| 6 | |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 7 | // 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. Dorian | f273487 | 2012-04-06 14:01:12 -0400 | [diff] [blame] | 39 | // Atanh returns the inverse hyperbolic tangent of x. |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 40 | // |
| 41 | // Special cases are: |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 42 | // Atanh(1) = +Inf |
Charles L. Dorian | 94b0342 | 2011-12-08 17:07:13 -0500 | [diff] [blame] | 43 | // Atanh(±0) = ±0 |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 44 | // Atanh(-1) = -Inf |
Charles L. Dorian | 94b0342 | 2011-12-08 17:07:13 -0500 | [diff] [blame] | 45 | // Atanh(x) = NaN if x < -1 or x > 1 |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 46 | // Atanh(NaN) = NaN |
| 47 | func Atanh(x float64) float64 { |
Charles L. Dorian | 3c3e68b | 2010-04-09 14:37:33 -0700 | [diff] [blame] | 48 | const NearZero = 1.0 / (1 << 28) // 2**-28 |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 49 | // special cases |
| 50 | switch { |
Luuk van Dijk | 8dd3de4 | 2012-02-01 16:08:31 +0100 | [diff] [blame] | 51 | case x < -1 || x > 1 || IsNaN(x): |
Charles L. Dorian | a0690b6 | 2010-02-01 22:21:40 -0800 | [diff] [blame] | 52 | 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 | } |