pxfm/hyperbolic/
asinhf.rs

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 */
29use crate::bits::{get_exponent_f64, set_exponent_f64};
30use crate::common::f_fmla;
31use crate::polyeval::{f_polyeval4, f_polyeval10};
32
33// Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127.
34static 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.
166static 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]
298pub(crate) fn log_eval(x: f64) -> f64 {
299    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    let p1 = ((x.to_bits() & ((1u64 << 52) - 1)) >> (52 - 7)) as i32;
304
305    // Set bs to (1 + (mx - p1*2^(-7))
306    let mut bs = x.to_bits() & (((1u64 << 52) - 1) >> 7);
307    const EXP_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
308    bs = set_exponent_f64(bs, EXP_BIAS);
309    // dx = (mx - p1*2^(-7)) / (1 + p1*2^(-7)).
310    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    let dx2 = dx * dx;
323    let c1 = f_fmla(dx, f64::from_bits(COEFFS[1]), f64::from_bits(COEFFS[0]));
324    let c2 = f_fmla(dx, f64::from_bits(COEFFS[3]), f64::from_bits(COEFFS[2]));
325    let c3 = f_fmla(dx, f64::from_bits(COEFFS[5]), f64::from_bits(COEFFS[4]));
326
327    let p = f_polyeval4(dx2, dx, c1, c2, c3);
328    f_fmla(
329        ex as f64,
330        /*log(2)*/ f64::from_bits(0x3fe62e42fefa39ef),
331        f64::from_bits(LOG_F[p1 as usize]) + p,
332    )
333}
334
335/// Hyperbolic arcsine function
336///
337/// Max ULP 0.5
338#[inline]
339pub fn f_asinhf(x: f32) -> f32 {
340    let x_u = x.to_bits();
341    let x_abs = x_u & 0x7fff_ffff;
342
343    // |x| <= 2^-3
344    if x_abs <= 0x3e80_0000u32 {
345        // |x| <= 2^-26
346        if x_abs <= 0x3280_0000u32 {
347            return if x_abs == 0 {
348                x
349            } else {
350                (x as f64 - f64::from_bits(0x3fc5555555555555) * x as f64 * x as f64 * x as f64)
351                    as f32
352            };
353        }
354
355        let x_d = x as f64;
356        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        let p = f_polyeval10(
362            x_sq,
363            0.0,
364            f64::from_bits(0xbfc5555555555555),
365            f64::from_bits(0x3fb3333333333333),
366            f64::from_bits(0xbfa6db6db6db6d8e),
367            f64::from_bits(0x3f9f1c71c71beb52),
368            f64::from_bits(0xbf96e8ba2e0dde02),
369            f64::from_bits(0x3f91c4ec071a2f97),
370            f64::from_bits(0xbf8c9966fc6b6fda),
371            f64::from_bits(0x3f879da45ad06ce8),
372            f64::from_bits(0xbf82b3657f620c14),
373        );
374        return f_fmla(x_d, p, x_d) as f32;
375    }
376
377    static SIGN: [f64; 2] = [1.0, -1.0];
378    let x_sign = SIGN[(x_u >> 31) as usize];
379    let x_d = x as f64;
380
381    if !x.is_finite() {
382        return x;
383    }
384
385    // asinh(x) = log(x + sqrt(x^2 + 1))
386    (x_sign * log_eval(f_fmla(x_d, x_sign, f_fmla(x_d, x_d, 1.0).sqrt()))) as f32
387}
388
389#[cfg(test)]
390mod 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}
pxfm/hyperbolic/asinhf.rs

pxfm/hyperbolic/
asinhf.rs
