1 use super::log1p; 2 3 /* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */ 4 /// Inverse hyperbolic tangent (f64) 5 /// 6 /// Calculates the inverse hyperbolic tangent of `x`. 7 /// Is defined as `log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2`. 8 pub fn atanh(x: f64) -> f64 { 9 let u = x.to_bits(); 10 let e = ((u >> 52) as usize) & 0x7ff; 11 let sign = (u >> 63) != 0; 12 13 /* |x| */ 14 let mut y = f64::from_bits(u & 0x7fff_ffff_ffff_ffff); 15 16 if e < 0x3ff - 1 { 17 if e < 0x3ff - 32 { 18 /* handle underflow */ 19 if e == 0 { 20 force_eval!(y as f32); 21 } 22 } else { 23 /* |x| < 0.5, up to 1.7ulp error */ 24 y = 0.5 * log1p(2.0 * y + 2.0 * y * y / (1.0 - y)); 25 } 26 } else { 27 /* avoid overflow */ 28 y = 0.5 * log1p(2.0 * (y / (1.0 - y))); atan2f(y: f32, x: f32) -> f3229 } 30 31 if sign { 32 -y 33 } else { 34 y 35 } 36 } 37