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 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 39 // Atanh returns the inverse hyperbolic tangent of x. 40 // 41 // Special cases are: 42 // Atanh(1) = +Inf 43 // Atanh(±0) = ±0 44 // Atanh(-1) = -Inf 45 // Atanh(x) = NaN if x < -1 or x > 1 46 // Atanh(NaN) = NaN 47 func Atanh(x float64) float64 { 48 if haveArchAtanh { 49 return archAtanh(x) 50 } 51 return atanh(x) 52 } 53 54 func atanh(x float64) float64 { 55 const NearZero = 1.0 / (1 << 28) // 2**-28 56 // special cases 57 switch { 58 case x < -1 || x > 1 || IsNaN(x): 59 return NaN() 60 case x == 1: 61 return Inf(1) 62 case x == -1: 63 return Inf(-1) 64 } 65 sign := false 66 if x < 0 { 67 x = -x 68 sign = true 69 } 70 var temp float64 71 switch { 72 case x < NearZero: 73 temp = x 74 case x < 0.5: 75 temp = x + x 76 temp = 0.5 * Log1p(temp+temp*x/(1-x)) 77 default: 78 temp = 0.5 * Log1p((x+x)/(1-x)) 79 } 80 if sign { 81 temp = -temp 82 } 83 return temp 84 } 85