pxfm/hyperbolic/
tanhf.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 6/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::common::{f_fmla, f_fmlaf};
30
31/// Hyperbolic tangent
32///
33/// Max found ULP 0.4999994
34#[inline]
35pub fn f_tanhf(x: f32) -> f32 {
36    let z = x as f64;
37    let t = x.to_bits();
38    let ux = t;
39    let e = ux.wrapping_shr(23) & 0xff;
40    if e == 0xff {
41        if ux << 9 != 0 {
42            return x + x;
43        } // x = nan
44        const IR: [f32; 2] = [1.0, -1.0];
45        return IR[ux.wrapping_shr(31) as usize]; // x = +-inf
46    }
47    if e < 115 {
48        // |x| < 2^-13
49        if e < 102 {
50            // |x| < 2^-26
51            if ux.wrapping_shl(1) == 0 {
52                return x;
53            }
54            let res = f_fmlaf(-x, x.abs(), x);
55            return res;
56        }
57        let x2 = x * x;
58        return f_fmlaf(x, -f64::from_bits(0x3fd5555560000000) as f32 * x2, x);
59    }
60    if ux.wrapping_shl(1) > (0x41102cb3u32 << 1) {
61        return f32::copysign(1.0, x) - f32::copysign(f64::from_bits(0x3e60000000000000) as f32, x);
62    }
63    let z2 = z * z;
64    let z4 = z2 * z2;
65    let z8 = z4 * z4;
66    const CN: [u64; 8] = [
67        0x3ff0000000000000,
68        0x3fc30877b8b72d33,
69        0x3f7694aa09ae9e5e,
70        0x3f14101377abb729,
71        0x3e9e0392b1db0018,
72        0x3e12533756e546f7,
73        0x3d6d62e5abe6ae8a,
74        0x3c9b06be534182de,
75    ];
76    const CD: [u64; 8] = [
77        0x3ff0000000000000,
78        0x3fded99131b0ebea,
79        0x3fa0d27ed6c95a69,
80        0x3f47cbdaca0e9fcc,
81        0x3edb4e60b892578e,
82        0x3e5a6f707c5c71ab,
83        0x3dc35a8b6e2cd94c,
84        0x3d0ca8230677aa01,
85    ];
86    let mut n0 = f_fmla(z2, f64::from_bits(CN[1]), f64::from_bits(CN[0]));
87    let n2 = f_fmla(z2, f64::from_bits(CN[3]), f64::from_bits(CN[2]));
88    let mut n4 = f_fmla(z2, f64::from_bits(CN[5]), f64::from_bits(CN[4]));
89    let n6 = f_fmla(z2, f64::from_bits(CN[7]), f64::from_bits(CN[6]));
90    n0 = f_fmla(z4, n2, n0);
91    n4 = f_fmla(z4, n6, n4);
92    n0 = f_fmla(z8, n4, n0);
93    let mut d0 = f_fmla(z2, f64::from_bits(CD[1]), f64::from_bits(CD[0]));
94    let d2 = f_fmla(z2, f64::from_bits(CD[3]), f64::from_bits(CD[2]));
95    let mut d4 = f_fmla(z2, f64::from_bits(CD[5]), f64::from_bits(CD[4]));
96    let d6 = f_fmla(z2, f64::from_bits(CD[7]), f64::from_bits(CD[6]));
97    d0 = f_fmla(z4, d2, d0);
98    d4 = f_fmla(z4, d6, d4);
99    d0 = f_fmla(z8, d4, d0);
100    let r = z * n0 / d0;
101    r as f32
102}
103
104#[cfg(test)]
105mod tests {
106    use super::*;
107
108    #[test]
109    fn test_tanhf() {
110        assert_eq!(f_tanhf(-0.5), -0.46211717);
111        assert_eq!(f_tanhf(0.5), 0.46211717);
112        assert_eq!(f_tanhf(7.), 0.99999833);
113    }
114}