pxfm/compound/
powm1.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 8/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::{is_integer, is_odd_integer};
30use crate::double_double::DoubleDouble;
31use crate::exponents::{EXPM1_T0, EXPM1_T1, ldexp};
32use crate::pow_exec::pow_log_1;
33use crate::round_ties_even::RoundTiesEven;
34
35/// Computes x^y - 1
36pub fn f_powm1(x: f64, y: f64) -> f64 {
37    let ax: u64 = x.to_bits().wrapping_shl(1);
38    let ay: u64 = y.to_bits().wrapping_shl(1);
39
40    // filter out exceptional cases
41    if ax == 0 || ax >= 0x7ffu64 << 53 || ay == 0 || ay >= 0x7ff64 << 53 {
42        if x.is_nan() || y.is_nan() {
43            return f64::NAN;
44        }
45
46        // Handle infinities
47        if x.is_infinite() {
48            return if x.is_sign_positive() {
49                if y.is_infinite() {
50                    return f64::INFINITY;
51                } else if y > 0.0 {
52                    f64::INFINITY // inf^positive -> inf
53                } else if y < 0.0 {
54                    -1.0 // inf^negative -> 0, so powm1 = -1
55                } else {
56                    f64::NAN // inf^0 -> NaN (0^0 conventionally 1, inf^0 = NaN)
57                }
58            } else {
59                // x = -inf
60                if y.is_infinite() {
61                    return -1.0;
62                }
63                if is_integer(y) {
64                    // Negative base: (-inf)^even = +inf, (-inf)^odd = -inf
65                    let pow = if y as i32 % 2 == 0 {
66                        f64::INFINITY
67                    } else {
68                        f64::NEG_INFINITY
69                    };
70                    pow - 1.0
71                } else {
72                    f64::NAN // Negative base with non-integer exponent
73                }
74            };
75        }
76
77        // Handle y infinite
78        if y.is_infinite() {
79            return if x.abs() > 1.0 {
80                if y.is_sign_positive() {
81                    f64::INFINITY
82                } else {
83                    -1.0
84                }
85            } else if x.abs() < 1.0 {
86                if y.is_sign_positive() {
87                    -1.0
88                } else {
89                    f64::INFINITY
90                }
91            } else {
92                // |x| == 1
93                f64::NAN // 1^inf or -1^inf is undefined
94            };
95        }
96
97        // Handle zero base
98        if x == 0.0 {
99            return if y > 0.0 {
100                -1.0 // 0^positive -> 0, powm1 = -1
101            } else if y < 0.0 {
102                f64::INFINITY // 0^negative -> inf
103            } else {
104                0.0 // 0^0 -> conventionally 1, powm1 = 0
105            };
106        }
107    }
108
109    let y_integer = is_integer(y);
110
111    let mut negative_parity: bool = false;
112
113    let mut x = x;
114
115    // Handle negative base with non-integer exponent
116    if x < 0.0 {
117        if !y_integer {
118            return f64::NAN; // x < 0 and non-integer y
119        }
120        x = x.abs();
121        if is_odd_integer(y) {
122            negative_parity = true;
123        }
124    }
125
126    let (mut l, _) = pow_log_1(x);
127    l = DoubleDouble::from_exact_add(l.hi, l.lo);
128
129    let r = DoubleDouble::quick_mult_f64(l, y);
130    if r.hi < -37.42994775023705 {
131        // underflow
132        return -1.;
133    }
134    let res = powm1_expm1_1(r);
135    // For x < 0 and integer y = n:
136    // if n is even: x^n = |x|^n → powm1 = |x|^n - 1 (same sign as res).
137    // if n is odd: x^n = -|x|^n → powm1 = -|x|^n - 1 = - (|x|^n + 1).
138    if negative_parity {
139        DoubleDouble::full_add_f64(-res, -2.).to_f64()
140    } else {
141        res.to_f64()
142    }
143}
144
145#[inline]
146pub(crate) fn powm1_expm1_1(r: DoubleDouble) -> DoubleDouble {
147    let ax = r.hi.to_bits() & 0x7fffffffffffffffu64;
148
149    const LOG2H: f64 = f64::from_bits(0x3f262e42fefa39ef);
150    const LOG2L: f64 = f64::from_bits(0x3bbabc9e3b39803f);
151
152    if ax <= 0x3f80000000000000 {
153        // |x| < 2^-7
154        if ax < 0x3970000000000000 {
155            // |x| < 2^-104
156            return r;
157        }
158        let d = crate::pow_exec::expm1_poly_dd_tiny(r);
159        return d;
160    }
161
162    const INVLOG2: f64 = f64::from_bits(0x40b71547652b82fe);
163
164    let k = (r.hi * INVLOG2).round_ties_even_finite();
165
166    let z = DoubleDouble::mul_f64_add(DoubleDouble::new(LOG2L, LOG2H), -k, r);
167
168    let bk = unsafe { k.to_int_unchecked::<i64>() }; /* Note: k is an integer, this is just a conversion. */
169    let mk = (bk >> 12) + 0x3ff;
170    let i2 = (bk >> 6) & 0x3f;
171    let i1 = bk & 0x3f;
172
173    let t0 = DoubleDouble::from_bit_pair(EXPM1_T0[i2 as usize]);
174    let t1 = DoubleDouble::from_bit_pair(EXPM1_T1[i1 as usize]);
175    let tbh = DoubleDouble::quick_mult(t1, t0);
176    let mut de = tbh;
177    // exp(k)=2^k*exp(r) + (2^k - 1)
178    let q = crate::pow_exec::expm1_poly_fast(z);
179    de = DoubleDouble::quick_mult(de, q);
180    de = DoubleDouble::add(tbh, de);
181
182    let ie = mk - 0x3ff;
183
184    let off: f64 = f64::from_bits((2048i64 + 1023i64).wrapping_sub(ie).wrapping_shl(52) as u64);
185
186    let e: f64;
187    if ie < 53 {
188        let fhz = DoubleDouble::from_exact_add(off, de.hi);
189        de.hi = fhz.hi;
190        e = fhz.lo;
191    } else if ie < 104 {
192        let fhz = DoubleDouble::from_exact_add(de.hi, off);
193        de.hi = fhz.hi;
194        e = fhz.lo;
195    } else {
196        e = 0.;
197    }
198    de.lo += e;
199    de.hi = ldexp(de.to_f64(), ie as i32);
200    de.lo = 0.;
201    de
202}
203
204#[cfg(test)]
205mod tests {
206    use super::*;
207    #[test]
208    fn test_powm1() {
209        assert_eq!(f_powm1(f64::INFINITY, f64::INFINITY), f64::INFINITY);
210        assert_eq!(f_powm1(50850368932909610000000000., 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000023201985303960773), 1.3733470789307166e-303);
211        assert_eq!(f_powm1(-3.375, -9671689000000000000000000.), -1.);
212        assert_eq!(f_powm1(1.83329e-40, 2.4645883e-32), -2.255031542428047e-30);
213        assert_eq!(f_powm1(3., 2.), 8.);
214        assert_eq!(f_powm1(3., 3.), 26.);
215        assert_eq!(f_powm1(5., 2.), 24.);
216        assert_eq!(f_powm1(5., -2.), 1. / 25. - 1.);
217        assert_eq!(f_powm1(-5., 2.), 24.);
218        assert_eq!(f_powm1(-5., 3.), -126.);
219        assert_eq!(
220            f_powm1(196560., 0.000000000000000000000000000000000000001193773),
221            1.4550568430468268e-38
222        );
223    }
224}
pxfm/compound/powm1.rs

pxfm/compound/
powm1.rs
