github.com/twelsh-aw/go/src@v0.0.0-20230516233729-a56fe86a7c81/math/atanh.go (about) 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 // 43 // Atanh(1) = +Inf 44 // Atanh(±0) = ±0 45 // Atanh(-1) = -Inf 46 // Atanh(x) = NaN if x < -1 or x > 1 47 // Atanh(NaN) = NaN 48 func Atanh(x float64) float64 { 49 if haveArchAtanh { 50 return archAtanh(x) 51 } 52 return atanh(x) 53 } 54 55 func atanh(x float64) float64 { 56 const NearZero = 1.0 / (1 << 28) // 2**-28 57 // special cases 58 switch { 59 case x < -1 || x > 1 || IsNaN(x): 60 return NaN() 61 case x == 1: 62 return Inf(1) 63 case x == -1: 64 return Inf(-1) 65 } 66 sign := false 67 if x < 0 { 68 x = -x 69 sign = true 70 } 71 var temp float64 72 switch { 73 case x < NearZero: 74 temp = x 75 case x < 0.5: 76 temp = x + x 77 temp = 0.5 * Log1p(temp+temp*x/(1-x)) 78 default: 79 temp = 0.5 * Log1p((x+x)/(1-x)) 80 } 81 if sign { 82 temp = -temp 83 } 84 return temp 85 }