Source file src/pkg/math/atanh.go
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(x) calculates the inverse hyperbolic tangent of x. 40 // 41 // Special cases are: 42 // Atanh(x) = NaN if x < -1 or x > 1 43 // Atanh(1) = +Inf 44 // Atanh(-1) = -Inf 45 // Atanh(NaN) = NaN 46 func Atanh(x float64) float64 { 47 const NearZero = 1.0 / (1 << 28) // 2**-28 48 // TODO(rsc): Remove manual inlining of IsNaN 49 // when compiler does it for us 50 // special cases 51 switch { 52 case x < -1 || x > 1 || x != x: // x < -1 || x > 1 || IsNaN(x): 53 return NaN() 54 case x == 1: 55 return Inf(1) 56 case x == -1: 57 return Inf(-1) 58 } 59 sign := false 60 if x < 0 { 61 x = -x 62 sign = true 63 } 64 var temp float64 65 switch { 66 case x < NearZero: 67 temp = x 68 case x < 0.5: 69 temp = x + x 70 temp = 0.5 * Log1p(temp+temp*x/(1-x)) 71 default: 72 temp = 0.5 * Log1p((x+x)/(1-x)) 73 } 74 if sign { 75 temp = -temp 76 } 77 return temp 78 }