/rust/registry/src/index.crates.io-1949cf8c6b5b557f/pxfm-0.1.27/src/hyperbolic/asinhf.rs

Source
/*
 * // Copyright (c) Radzivon Bartoshyk 7/2025. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
use crate::bits::{get_exponent_f64, set_exponent_f64};
use crate::common::f_fmla;
use crate::polyeval::{f_polyeval4, f_polyeval10};

// Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127.
static ONE_OVER_F: [u64; 128] = [
    0x3ff0000000000000,
    0x3fefc07f01fc07f0,
    0x3fef81f81f81f820,
    0x3fef44659e4a4271,
    0x3fef07c1f07c1f08,
    0x3feecc07b301ecc0,
    0x3fee9131abf0b767,
    0x3fee573ac901e574,
    0x3fee1e1e1e1e1e1e,
    0x3fede5d6e3f8868a,
    0x3fedae6076b981db,
    0x3fed77b654b82c34,
    0x3fed41d41d41d41d,
    0x3fed0cb58f6ec074,
    0x3fecd85689039b0b,
    0x3feca4b3055ee191,
    0x3fec71c71c71c71c,
    0x3fec3f8f01c3f8f0,
    0x3fec0e070381c0e0,
    0x3febdd2b899406f7,
    0x3febacf914c1bad0,
    0x3feb7d6c3dda338b,
    0x3feb4e81b4e81b4f,
    0x3feb2036406c80d9,
    0x3feaf286bca1af28,
    0x3feac5701ac5701b,
    0x3fea98ef606a63be,
    0x3fea6d01a6d01a6d,
    0x3fea41a41a41a41a,
    0x3fea16d3f97a4b02,
    0x3fe9ec8e951033d9,
    0x3fe9c2d14ee4a102,
    0x3fe999999999999a,
    0x3fe970e4f80cb872,
    0x3fe948b0fcd6e9e0,
    0x3fe920fb49d0e229,
    0x3fe8f9c18f9c18fa,
    0x3fe8d3018d3018d3,
    0x3fe8acb90f6bf3aa,
    0x3fe886e5f0abb04a,
    0x3fe8618618618618,
    0x3fe83c977ab2bedd,
    0x3fe8181818181818,
    0x3fe7f405fd017f40,
    0x3fe7d05f417d05f4,
    0x3fe7ad2208e0ecc3,
    0x3fe78a4c8178a4c8,
    0x3fe767dce434a9b1,
    0x3fe745d1745d1746,
    0x3fe724287f46debc,
    0x3fe702e05c0b8170,
    0x3fe6e1f76b4337c7,
    0x3fe6c16c16c16c17,
    0x3fe6a13cd1537290,
    0x3fe6816816816817,
    0x3fe661ec6a5122f9,
    0x3fe642c8590b2164,
    0x3fe623fa77016240,
    0x3fe6058160581606,
    0x3fe5e75bb8d015e7,
    0x3fe5c9882b931057,
    0x3fe5ac056b015ac0,
    0x3fe58ed2308158ed,
    0x3fe571ed3c506b3a,
    0x3fe5555555555555,
    0x3fe5390948f40feb,
    0x3fe51d07eae2f815,
    0x3fe5015015015015,
    0x3fe4e5e0a72f0539,
    0x3fe4cab88725af6e,
    0x3fe4afd6a052bf5b,
    0x3fe49539e3b2d067,
    0x3fe47ae147ae147b,
    0x3fe460cbc7f5cf9a,
    0x3fe446f86562d9fb,
    0x3fe42d6625d51f87,
    0x3fe4141414141414,
    0x3fe3fb013fb013fb,
    0x3fe3e22cbce4a902,
    0x3fe3c995a47babe7,
    0x3fe3b13b13b13b14,
    0x3fe3991c2c187f63,
    0x3fe3813813813814,
    0x3fe3698df3de0748,
    0x3fe3521cfb2b78c1,
    0x3fe33ae45b57bcb2,
    0x3fe323e34a2b10bf,
    0x3fe30d190130d190,
    0x3fe2f684bda12f68,
    0x3fe2e025c04b8097,
    0x3fe2c9fb4d812ca0,
    0x3fe2b404ad012b40,
    0x3fe29e4129e4129e,
    0x3fe288b01288b013,
    0x3fe27350b8812735,
    0x3fe25e22708092f1,
    0x3fe2492492492492,
    0x3fe23456789abcdf,
    0x3fe21fb78121fb78,
    0x3fe20b470c67c0d9,
    0x3fe1f7047dc11f70,
    0x3fe1e2ef3b3fb874,
    0x3fe1cf06ada2811d,
    0x3fe1bb4a4046ed29,
    0x3fe1a7b9611a7b96,
    0x3fe19453808ca29c,
    0x3fe1811811811812,
    0x3fe16e0689427379,
    0x3fe15b1e5f75270d,
    0x3fe1485f0e0acd3b,
    0x3fe135c81135c811,
    0x3fe12358e75d3033,
    0x3fe1111111111111,
    0x3fe0fef010fef011,
    0x3fe0ecf56be69c90,
    0x3fe0db20a88f4696,
    0x3fe0c9714fbcda3b,
    0x3fe0b7e6ec259dc8,
    0x3fe0a6810a6810a7,
    0x3fe0953f39010954,
    0x3fe0842108421084,
    0x3fe073260a47f7c6,
    0x3fe0624dd2f1a9fc,
    0x3fe05197f7d73404,
    0x3fe0410410410410,
    0x3fe03091b51f5e1a,
    0x3fe0204081020408,
    0x3fe0101010101010,
];

// Lookup table for log(f) = log(1 + n*2^(-7)) where n = 0..127.
static LOG_F: [u64; 128] = [
    0x0000000000000000,
    0x3f7fe02a6b106788,
    0x3f8fc0a8b0fc03e3,
    0x3f97b91b07d5b11a,
    0x3f9f829b0e783300,
    0x3fa39e87b9febd5f,
    0x3fa77458f632dcfc,
    0x3fab42dd711971be,
    0x3faf0a30c01162a6,
    0x3fb16536eea37ae0,
    0x3fb341d7961bd1d0,
    0x3fb51b073f06183f,
    0x3fb6f0d28ae56b4b,
    0x3fb8c345d6319b20,
    0x3fba926d3a4ad563,
    0x3fbc5e548f5bc743,
    0x3fbe27076e2af2e5,
    0x3fbfec9131dbeaba,
    0x3fc0d77e7cd08e59,
    0x3fc1b72ad52f67a0,
    0x3fc29552f81ff523,
    0x3fc371fc201e8f74,
    0x3fc44d2b6ccb7d1e,
    0x3fc526e5e3a1b437,
    0x3fc5ff3070a793d3,
    0x3fc6d60fe719d21c,
    0x3fc7ab890210d909,
    0x3fc87fa06520c910,
    0x3fc9525a9cf456b4,
    0x3fca23bc1fe2b563,
    0x3fcaf3c94e80bff2,
    0x3fcbc286742d8cd6,
    0x3fcc8ff7c79a9a21,
    0x3fcd5c216b4fbb91,
    0x3fce27076e2af2e5,
    0x3fcef0adcbdc5936,
    0x3fcfb9186d5e3e2a,
    0x3fd0402594b4d040,
    0x3fd0a324e27390e3,
    0x3fd1058bf9ae4ad5,
    0x3fd1675cababa60e,
    0x3fd1c898c16999fa,
    0x3fd22941fbcf7965,
    0x3fd2895a13de86a3,
    0x3fd2e8e2bae11d30,
    0x3fd347dd9a987d54,
    0x3fd3a64c556945e9,
    0x3fd404308686a7e3,
    0x3fd4618bc21c5ec2,
    0x3fd4be5f957778a0,
    0x3fd51aad872df82d,
    0x3fd5767717455a6c,
    0x3fd5d1bdbf5809ca,
    0x3fd62c82f2b9c795,
    0x3fd686c81e9b14ae,
    0x3fd6e08eaa2ba1e3,
    0x3fd739d7f6bbd006,
    0x3fd792a55fdd47a2,
    0x3fd7eaf83b82afc3,
    0x3fd842d1da1e8b17,
    0x3fd89a3386c1425a,
    0x3fd8f11e873662c7,
    0x3fd947941c2116fa,
    0x3fd99d958117e08a,
    0x3fd9f323ecbf984b,
    0x3fda484090e5bb0a,
    0x3fda9cec9a9a0849,
    0x3fdaf1293247786b,
    0x3fdb44f77bcc8f62,
    0x3fdb9858969310fb,
    0x3fdbeb4d9da71b7b,
    0x3fdc3dd7a7cdad4d,
    0x3fdc8ff7c79a9a21,
    0x3fdce1af0b85f3eb,
    0x3fdd32fe7e00ebd5,
    0x3fdd83e7258a2f3e,
    0x3fddd46a04c1c4a0,
    0x3fde24881a7c6c26,
    0x3fde744261d68787,
    0x3fdec399d2468cc0,
    0x3fdf128f5faf06ec,
    0x3fdf6123fa7028ac,
    0x3fdfaf588f78f31e,
    0x3fdffd2e0857f498,
    0x3fe02552a5a5d0fe,
    0x3fe04bdf9da926d2,
    0x3fe0723e5c1cdf40,
    0x3fe0986f4f573520,
    0x3fe0be72e4252a82,
    0x3fe0e44985d1cc8b,
    0x3fe109f39e2d4c96,
    0x3fe12f719593efbc,
    0x3fe154c3d2f4d5e9,
    0x3fe179eabbd899a0,
    0x3fe19ee6b467c96e,
    0x3fe1c3b81f713c24,
    0x3fe1e85f5e7040d0,
    0x3fe20cdcd192ab6d,
    0x3fe23130d7bebf42,
    0x3fe2555bce98f7cb,
    0x3fe2795e1289b11a,
    0x3fe29d37fec2b08a,
    0x3fe2c0e9ed448e8b,
    0x3fe2e47436e40268,
    0x3fe307d7334f10be,
    0x3fe32b1339121d71,
    0x3fe34e289d9ce1d3,
    0x3fe37117b54747b5,
    0x3fe393e0d3562a19,
    0x3fe3b68449fffc22,
    0x3fe3d9026a7156fa,
    0x3fe3fb5b84d16f42,
    0x3fe41d8fe84672ae,
    0x3fe43f9fe2f9ce67,
    0x3fe4618bc21c5ec2,
    0x3fe48353d1ea88df,
    0x3fe4a4f85db03ebb,
    0x3fe4c679afccee39,
    0x3fe4e7d811b75bb0,
    0x3fe50913cc01686b,
    0x3fe52a2d265bc5aa,
    0x3fe54b2467999497,
    0x3fe56bf9d5b3f399,
    0x3fe58cadb5cd7989,
    0x3fe5ad404c359f2c,
    0x3fe5cdb1dc6c1764,
    0x3fe5ee02a9241675,
    0x3fe60e32f44788d8,
];

#[inline]
pub(crate) fn log_eval(x: f64) -> f64 {
    let ex = get_exponent_f64(x);

    // p1 is the leading 7 bits of mx, i.e.
    // p1 * 2^(-7) <= m_x < (p1 + 1) * 2^(-7).
    let p1 = ((x.to_bits() & ((1u64 << 52) - 1)) >> (52 - 7)) as i32;

    // Set bs to (1 + (mx - p1*2^(-7))
    let mut bs = x.to_bits() & (((1u64 << 52) - 1) >> 7);
    const EXP_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
    bs = set_exponent_f64(bs, EXP_BIAS);
    // dx = (mx - p1*2^(-7)) / (1 + p1*2^(-7)).
    let dx = (f64::from_bits(bs) - 1.0) * f64::from_bits(ONE_OVER_F[p1 as usize]);

    // Minimax polynomial of log(1 + dx) generated by Sollya with:
    // > P = fpminimax(log(1 + x)/x, 6, [|D...|], [0, 2^-7]);
    const COEFFS: [u64; 6] = [
        0xbfdffffffffffffc,
        0x3fd5555555552dde,
        0xbfcffffffefe562d,
        0x3fc9999817d3a50f,
        0xbfc554317b3f67a5,
        0x3fc1dc5c45e09c18,
    ];
    let dx2 = dx * dx;
    let c1 = f_fmla(dx, f64::from_bits(COEFFS[1]), f64::from_bits(COEFFS[0]));
    let c2 = f_fmla(dx, f64::from_bits(COEFFS[3]), f64::from_bits(COEFFS[2]));
    let c3 = f_fmla(dx, f64::from_bits(COEFFS[5]), f64::from_bits(COEFFS[4]));

    let p = f_polyeval4(dx2, dx, c1, c2, c3);
    f_fmla(
        ex as f64,
        /*log(2)*/ f64::from_bits(0x3fe62e42fefa39ef),
        f64::from_bits(LOG_F[p1 as usize]) + p,
    )
}

/// Hyperbolic arcsine function
///
/// Max ULP 0.5
#[inline]
pub fn f_asinhf(x: f32) -> f32 {
    let x_u = x.to_bits();
    let x_abs = x_u & 0x7fff_ffff;

    // |x| <= 2^-3
    if x_abs <= 0x3e80_0000u32 {
        // |x| <= 2^-26
        if x_abs <= 0x3280_0000u32 {
            return if x_abs == 0 {
                x
            } else {
                (x as f64 - f64::from_bits(0x3fc5555555555555) * x as f64 * x as f64 * x as f64)
                    as f32
            };
        }

        let x_d = x as f64;
        let x_sq = x_d * x_d;
        // Generated by Sollya with:
        // R = remez(asinh(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18|], [0, 2^-3]);
        // P = fpminimax(asinh(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18|], [|D...|], [0, 2^-3], absolute, floating, R);
        // See `./notes/asinhf.sollya
        let p = f_polyeval10(
            x_sq,
            0.0,
            f64::from_bits(0xbfc5555555555555),
            f64::from_bits(0x3fb3333333333333),
            f64::from_bits(0xbfa6db6db6db6d8e),
            f64::from_bits(0x3f9f1c71c71beb52),
            f64::from_bits(0xbf96e8ba2e0dde02),
            f64::from_bits(0x3f91c4ec071a2f97),
            f64::from_bits(0xbf8c9966fc6b6fda),
            f64::from_bits(0x3f879da45ad06ce8),
            f64::from_bits(0xbf82b3657f620c14),
        );
        return f_fmla(x_d, p, x_d) as f32;
    }

    static SIGN: [f64; 2] = [1.0, -1.0];
    let x_sign = SIGN[(x_u >> 31) as usize];
    let x_d = x as f64;

    if !x.is_finite() {
        return x;
    }

    // asinh(x) = log(x + sqrt(x^2 + 1))
    (x_sign * log_eval(f_fmla(x_d, x_sign, f_fmla(x_d, x_d, 1.0).sqrt()))) as f32
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_asinhf() {
        assert_eq!(f_asinhf(-0.24319594), -0.2408603);
        assert_eq!(f_asinhf(2.0), 1.4436355);
        assert_eq!(f_asinhf(-2.0), -1.4436355);
        assert_eq!(f_asinhf(1.054234), 0.9192077);
    }
}

Coverage Report

Created: 2025-12-20 06:48

Line	Count	Source
1		/*
2		* // Copyright (c) Radzivon Bartoshyk 7/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		*/
29		use crate::bits::{get_exponent_f64, set_exponent_f64};
30		use crate::common::f_fmla;
31		use crate::polyeval::{f_polyeval4, f_polyeval10};
32
33		// Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127.
34		static ONE_OVER_F: [u64; 128] = [
35		0x3ff0000000000000,
36		0x3fefc07f01fc07f0,
37		0x3fef81f81f81f820,
38		0x3fef44659e4a4271,
39		0x3fef07c1f07c1f08,
40		0x3feecc07b301ecc0,
41		0x3fee9131abf0b767,
42		0x3fee573ac901e574,
43		0x3fee1e1e1e1e1e1e,
44		0x3fede5d6e3f8868a,
45		0x3fedae6076b981db,
46		0x3fed77b654b82c34,
47		0x3fed41d41d41d41d,
48		0x3fed0cb58f6ec074,
49		0x3fecd85689039b0b,
50		0x3feca4b3055ee191,
51		0x3fec71c71c71c71c,
52		0x3fec3f8f01c3f8f0,
53		0x3fec0e070381c0e0,
54		0x3febdd2b899406f7,
55		0x3febacf914c1bad0,
56		0x3feb7d6c3dda338b,
57		0x3feb4e81b4e81b4f,
58		0x3feb2036406c80d9,
59		0x3feaf286bca1af28,
60		0x3feac5701ac5701b,
61		0x3fea98ef606a63be,
62		0x3fea6d01a6d01a6d,
63		0x3fea41a41a41a41a,
64		0x3fea16d3f97a4b02,
65		0x3fe9ec8e951033d9,
66		0x3fe9c2d14ee4a102,
67		0x3fe999999999999a,
68		0x3fe970e4f80cb872,
69		0x3fe948b0fcd6e9e0,
70		0x3fe920fb49d0e229,
71		0x3fe8f9c18f9c18fa,
72		0x3fe8d3018d3018d3,
73		0x3fe8acb90f6bf3aa,
74		0x3fe886e5f0abb04a,
75		0x3fe8618618618618,
76		0x3fe83c977ab2bedd,
77		0x3fe8181818181818,
78		0x3fe7f405fd017f40,
79		0x3fe7d05f417d05f4,
80		0x3fe7ad2208e0ecc3,
81		0x3fe78a4c8178a4c8,
82		0x3fe767dce434a9b1,
83		0x3fe745d1745d1746,
84		0x3fe724287f46debc,
85		0x3fe702e05c0b8170,
86		0x3fe6e1f76b4337c7,
87		0x3fe6c16c16c16c17,
88		0x3fe6a13cd1537290,
89		0x3fe6816816816817,
90		0x3fe661ec6a5122f9,
91		0x3fe642c8590b2164,
92		0x3fe623fa77016240,
93		0x3fe6058160581606,
94		0x3fe5e75bb8d015e7,
95		0x3fe5c9882b931057,
96		0x3fe5ac056b015ac0,
97		0x3fe58ed2308158ed,
98		0x3fe571ed3c506b3a,
99		0x3fe5555555555555,
100		0x3fe5390948f40feb,
101		0x3fe51d07eae2f815,
102		0x3fe5015015015015,
103		0x3fe4e5e0a72f0539,
104		0x3fe4cab88725af6e,
105		0x3fe4afd6a052bf5b,
106		0x3fe49539e3b2d067,
107		0x3fe47ae147ae147b,
108		0x3fe460cbc7f5cf9a,
109		0x3fe446f86562d9fb,
110		0x3fe42d6625d51f87,
111		0x3fe4141414141414,
112		0x3fe3fb013fb013fb,
113		0x3fe3e22cbce4a902,
114		0x3fe3c995a47babe7,
115		0x3fe3b13b13b13b14,
116		0x3fe3991c2c187f63,
117		0x3fe3813813813814,
118		0x3fe3698df3de0748,
119		0x3fe3521cfb2b78c1,
120		0x3fe33ae45b57bcb2,
121		0x3fe323e34a2b10bf,
122		0x3fe30d190130d190,
123		0x3fe2f684bda12f68,
124		0x3fe2e025c04b8097,
125		0x3fe2c9fb4d812ca0,
126		0x3fe2b404ad012b40,
127		0x3fe29e4129e4129e,
128		0x3fe288b01288b013,
129		0x3fe27350b8812735,
130		0x3fe25e22708092f1,
131		0x3fe2492492492492,
132		0x3fe23456789abcdf,
133		0x3fe21fb78121fb78,
134		0x3fe20b470c67c0d9,
135		0x3fe1f7047dc11f70,
136		0x3fe1e2ef3b3fb874,
137		0x3fe1cf06ada2811d,
138		0x3fe1bb4a4046ed29,
139		0x3fe1a7b9611a7b96,
140		0x3fe19453808ca29c,
141		0x3fe1811811811812,
142		0x3fe16e0689427379,
143		0x3fe15b1e5f75270d,
144		0x3fe1485f0e0acd3b,
145		0x3fe135c81135c811,
146		0x3fe12358e75d3033,
147		0x3fe1111111111111,
148		0x3fe0fef010fef011,
149		0x3fe0ecf56be69c90,
150		0x3fe0db20a88f4696,
151		0x3fe0c9714fbcda3b,
152		0x3fe0b7e6ec259dc8,
153		0x3fe0a6810a6810a7,
154		0x3fe0953f39010954,
155		0x3fe0842108421084,
156		0x3fe073260a47f7c6,
157		0x3fe0624dd2f1a9fc,
158		0x3fe05197f7d73404,
159		0x3fe0410410410410,
160		0x3fe03091b51f5e1a,
161		0x3fe0204081020408,
162		0x3fe0101010101010,
163		];
164
165		// Lookup table for log(f) = log(1 + n*2^(-7)) where n = 0..127.
166		static LOG_F: [u64; 128] = [
167		0x0000000000000000,
168		0x3f7fe02a6b106788,
169		0x3f8fc0a8b0fc03e3,
170		0x3f97b91b07d5b11a,
171		0x3f9f829b0e783300,
172		0x3fa39e87b9febd5f,
173		0x3fa77458f632dcfc,
174		0x3fab42dd711971be,
175		0x3faf0a30c01162a6,
176		0x3fb16536eea37ae0,
177		0x3fb341d7961bd1d0,
178		0x3fb51b073f06183f,
179		0x3fb6f0d28ae56b4b,
180		0x3fb8c345d6319b20,
181		0x3fba926d3a4ad563,
182		0x3fbc5e548f5bc743,
183		0x3fbe27076e2af2e5,
184		0x3fbfec9131dbeaba,
185		0x3fc0d77e7cd08e59,
186		0x3fc1b72ad52f67a0,
187		0x3fc29552f81ff523,
188		0x3fc371fc201e8f74,
189		0x3fc44d2b6ccb7d1e,
190		0x3fc526e5e3a1b437,
191		0x3fc5ff3070a793d3,
192		0x3fc6d60fe719d21c,
193		0x3fc7ab890210d909,
194		0x3fc87fa06520c910,
195		0x3fc9525a9cf456b4,
196		0x3fca23bc1fe2b563,
197		0x3fcaf3c94e80bff2,
198		0x3fcbc286742d8cd6,
199		0x3fcc8ff7c79a9a21,
200		0x3fcd5c216b4fbb91,
201		0x3fce27076e2af2e5,
202		0x3fcef0adcbdc5936,
203		0x3fcfb9186d5e3e2a,
204		0x3fd0402594b4d040,
205		0x3fd0a324e27390e3,
206		0x3fd1058bf9ae4ad5,
207		0x3fd1675cababa60e,
208		0x3fd1c898c16999fa,
209		0x3fd22941fbcf7965,
210		0x3fd2895a13de86a3,
211		0x3fd2e8e2bae11d30,
212		0x3fd347dd9a987d54,
213		0x3fd3a64c556945e9,
214		0x3fd404308686a7e3,
215		0x3fd4618bc21c5ec2,
216		0x3fd4be5f957778a0,
217		0x3fd51aad872df82d,
218		0x3fd5767717455a6c,
219		0x3fd5d1bdbf5809ca,
220		0x3fd62c82f2b9c795,
221		0x3fd686c81e9b14ae,
222		0x3fd6e08eaa2ba1e3,
223		0x3fd739d7f6bbd006,
224		0x3fd792a55fdd47a2,
225		0x3fd7eaf83b82afc3,
226		0x3fd842d1da1e8b17,
227		0x3fd89a3386c1425a,
228		0x3fd8f11e873662c7,
229		0x3fd947941c2116fa,
230		0x3fd99d958117e08a,
231		0x3fd9f323ecbf984b,
232		0x3fda484090e5bb0a,
233		0x3fda9cec9a9a0849,
234		0x3fdaf1293247786b,
235		0x3fdb44f77bcc8f62,
236		0x3fdb9858969310fb,
237		0x3fdbeb4d9da71b7b,
238		0x3fdc3dd7a7cdad4d,
239		0x3fdc8ff7c79a9a21,
240		0x3fdce1af0b85f3eb,
241		0x3fdd32fe7e00ebd5,
242		0x3fdd83e7258a2f3e,
243		0x3fddd46a04c1c4a0,
244		0x3fde24881a7c6c26,
245		0x3fde744261d68787,
246		0x3fdec399d2468cc0,
247		0x3fdf128f5faf06ec,
248		0x3fdf6123fa7028ac,
249		0x3fdfaf588f78f31e,
250		0x3fdffd2e0857f498,
251		0x3fe02552a5a5d0fe,
252		0x3fe04bdf9da926d2,
253		0x3fe0723e5c1cdf40,
254		0x3fe0986f4f573520,
255		0x3fe0be72e4252a82,
256		0x3fe0e44985d1cc8b,
257		0x3fe109f39e2d4c96,
258		0x3fe12f719593efbc,
259		0x3fe154c3d2f4d5e9,
260		0x3fe179eabbd899a0,
261		0x3fe19ee6b467c96e,
262		0x3fe1c3b81f713c24,
263		0x3fe1e85f5e7040d0,
264		0x3fe20cdcd192ab6d,
265		0x3fe23130d7bebf42,
266		0x3fe2555bce98f7cb,
267		0x3fe2795e1289b11a,
268		0x3fe29d37fec2b08a,
269		0x3fe2c0e9ed448e8b,
270		0x3fe2e47436e40268,
271		0x3fe307d7334f10be,
272		0x3fe32b1339121d71,
273		0x3fe34e289d9ce1d3,
274		0x3fe37117b54747b5,
275		0x3fe393e0d3562a19,
276		0x3fe3b68449fffc22,
277		0x3fe3d9026a7156fa,
278		0x3fe3fb5b84d16f42,
279		0x3fe41d8fe84672ae,
280		0x3fe43f9fe2f9ce67,
281		0x3fe4618bc21c5ec2,
282		0x3fe48353d1ea88df,
283		0x3fe4a4f85db03ebb,
284		0x3fe4c679afccee39,
285		0x3fe4e7d811b75bb0,
286		0x3fe50913cc01686b,
287		0x3fe52a2d265bc5aa,
288		0x3fe54b2467999497,
289		0x3fe56bf9d5b3f399,
290		0x3fe58cadb5cd7989,
291		0x3fe5ad404c359f2c,
292		0x3fe5cdb1dc6c1764,
293		0x3fe5ee02a9241675,
294		0x3fe60e32f44788d8,
295		];
296
297		#[inline]
298	0	pub(crate) fn log_eval(x: f64) -> f64 {
299	0	let ex = get_exponent_f64(x);
300
301		// p1 is the leading 7 bits of mx, i.e.
302		// p1 * 2^(-7) <= m_x < (p1 + 1) * 2^(-7).
303	0	let p1 = ((x.to_bits() & ((1u64 << 52) - 1)) >> (52 - 7)) as i32;
304
305		// Set bs to (1 + (mx - p1*2^(-7))
306	0	let mut bs = x.to_bits() & (((1u64 << 52) - 1) >> 7);
307		const EXP_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
308	0	bs = set_exponent_f64(bs, EXP_BIAS);
309		// dx = (mx - p12^(-7)) / (1 + p12^(-7)).
310	0	let dx = (f64::from_bits(bs) - 1.0) * f64::from_bits(ONE_OVER_F[p1 as usize]);
311
312		// Minimax polynomial of log(1 + dx) generated by Sollya with:
313		// > P = fpminimax(log(1 + x)/x, 6, [\|D...\|], [0, 2^-7]);
314		const COEFFS: [u64; 6] = [
315		0xbfdffffffffffffc,
316		0x3fd5555555552dde,
317		0xbfcffffffefe562d,
318		0x3fc9999817d3a50f,
319		0xbfc554317b3f67a5,
320		0x3fc1dc5c45e09c18,
321		];
322	0	let dx2 = dx * dx;
323	0	let c1 = f_fmla(dx, f64::from_bits(COEFFS[1]), f64::from_bits(COEFFS[0]));
324	0	let c2 = f_fmla(dx, f64::from_bits(COEFFS[3]), f64::from_bits(COEFFS[2]));
325	0	let c3 = f_fmla(dx, f64::from_bits(COEFFS[5]), f64::from_bits(COEFFS[4]));
326
327	0	let p = f_polyeval4(dx2, dx, c1, c2, c3);
328	0	f_fmla(
329	0	ex as f64,
330	0	/log(2)/ f64::from_bits(0x3fe62e42fefa39ef),
331	0	f64::from_bits(LOG_F[p1 as usize]) + p,
332		)
333	0	}
334
335		/// Hyperbolic arcsine function
336		///
337		/// Max ULP 0.5
338		#[inline]
339	0	pub fn f_asinhf(x: f32) -> f32 {
340	0	let x_u = x.to_bits();
341	0	let x_abs = x_u & 0x7fff_ffff;
342
343		// \|x\| <= 2^-3
344	0	if x_abs <= 0x3e80_0000u32 {
345		// \|x\| <= 2^-26
346	0	if x_abs <= 0x3280_0000u32 {
347	0	return if x_abs == 0 {
348	0	x
349		} else {
350	0	(x as f64 - f64::from_bits(0x3fc5555555555555) * x as f64 * x as f64 * x as f64)
351	0	as f32
352		};
353	0	}
354
355	0	let x_d = x as f64;
356	0	let x_sq = x_d * x_d;
357		// Generated by Sollya with:
358		// R = remez(asinh(x)/x, [\|0, 2, 4, 6, 8, 10, 12, 14, 16, 18\|], [0, 2^-3]);
359		// P = fpminimax(asinh(x)/x, [\|0, 2, 4, 6, 8, 10, 12, 14, 16, 18\|], [\|D...\|], [0, 2^-3], absolute, floating, R);
360		// See `./notes/asinhf.sollya
361	0	let p = f_polyeval10(
362	0	x_sq,
363		0.0,
364	0	f64::from_bits(0xbfc5555555555555),
365	0	f64::from_bits(0x3fb3333333333333),
366	0	f64::from_bits(0xbfa6db6db6db6d8e),
367	0	f64::from_bits(0x3f9f1c71c71beb52),
368	0	f64::from_bits(0xbf96e8ba2e0dde02),
369	0	f64::from_bits(0x3f91c4ec071a2f97),
370	0	f64::from_bits(0xbf8c9966fc6b6fda),
371	0	f64::from_bits(0x3f879da45ad06ce8),
372	0	f64::from_bits(0xbf82b3657f620c14),
373		);
374	0	return f_fmla(x_d, p, x_d) as f32;
375	0	}
376
377		static SIGN: [f64; 2] = [1.0, -1.0];
378	0	let x_sign = SIGN[(x_u >> 31) as usize];
379	0	let x_d = x as f64;
380
381	0	if !x.is_finite() {
382	0	return x;
383	0	}
384
385		// asinh(x) = log(x + sqrt(x^2 + 1))
386	0	(x_sign * log_eval(f_fmla(x_d, x_sign, f_fmla(x_d, x_d, 1.0).sqrt()))) as f32
387	0	}
388
389		#[cfg(test)]
390		mod tests {
391		use super::*;
392
393		#[test]
394		fn test_asinhf() {
395		assert_eq!(f_asinhf(-0.24319594), -0.2408603);
396		assert_eq!(f_asinhf(2.0), 1.4436355);
397		assert_eq!(f_asinhf(-2.0), -1.4436355);
398		assert_eq!(f_asinhf(1.054234), 0.9192077);
399		}
400		}