pxfm/logs/
log_dd.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::double_double::DoubleDouble;
30use crate::logs::log_dd_coeffs::LOG_NEG_DD;
31use crate::polyeval::f_polyeval7;
32use crate::pow_tables::POW_INVERSE;
33
34#[inline(always)]
35pub(crate) fn log_poly(z: f64) -> DoubleDouble {
36    /*
37      See ./notes/dd_log.sollya
38    */
39    const P: [(u64, u64); 9] = [
40        (0x3c755555556a4311, 0x3fd5555555555555),
41        (0x3b8ffa82859b4000, 0xbfd0000000000000),
42        (0xbc699b6b44796cd4, 0x3fc999999999999a),
43        (0x3c489642b3424250, 0xbfc5555555555563),
44        (0x3c60e3ef2c7fc443, 0x3fc24924924924ad),
45        (0x3c535269fe0ce5df, 0xbfbfffffffc0b25c),
46        (0x3c24f14f55e95858, 0x3fbc71c71c187ea4),
47        (0x3c5b15951c5e1a17, 0xbfb999d71e5042d5),
48        (0x3c5ecb0133e43410, 0x3fb74615b842a94d),
49    ];
50    let x2 = DoubleDouble::from_exact_mult(z, z);
51    let mut t = DoubleDouble::quick_mul_f64_add(
52        DoubleDouble::from_bit_pair(P[8]),
53        z,
54        DoubleDouble::from_bit_pair(P[7]),
55    );
56    t = DoubleDouble::quick_mul_f64_add(t, z, DoubleDouble::from_bit_pair(P[6]));
57    t = DoubleDouble::quick_mul_f64_add(t, z, DoubleDouble::from_bit_pair(P[5]));
58    t = DoubleDouble::quick_mul_f64_add(t, z, DoubleDouble::from_bit_pair(P[4]));
59    t = DoubleDouble::quick_mul_f64_add(t, z, DoubleDouble::from_bit_pair(P[3]));
60    t = DoubleDouble::quick_mul_f64_add(t, z, DoubleDouble::from_bit_pair(P[2]));
61    t = DoubleDouble::quick_mul_f64_add(t, z, DoubleDouble::from_bit_pair(P[1]));
62    t = DoubleDouble::quick_mul_f64_add(t, z, DoubleDouble::from_bit_pair(P[0]));
63    t = DoubleDouble::quick_mult(t, x2);
64    t = DoubleDouble::quick_mult_f64(t, z);
65    DoubleDouble::mul_f64_add(x2, -0.5, t)
66}
67
68#[inline(always)]
69pub(crate) fn log_poly_fast(z: f64) -> DoubleDouble {
70    /*
71      See ./notes/dd_log_fast.sollya
72    */
73    let x2 = DoubleDouble::from_exact_mult(z, z);
74    let r = f_polyeval7(
75        z,
76        f64::from_bits(0xbfc5555555555555),
77        f64::from_bits(0x3fc24924924924aa),
78        f64::from_bits(0xbfc000000000bc96),
79        f64::from_bits(0x3fbc71c71c202bf0),
80        f64::from_bits(0xbfb9999839f1aa36),
81        f64::from_bits(0x3fb74612adef67e0),
82        f64::from_bits(0xbfb5cb6ab20b8efa),
83    );
84    let mut p = DoubleDouble::f64_mul_f64_add(
85        z,
86        r,
87        DoubleDouble::from_bit_pair((0xbc699b293fa3344b, 0x3fc999999999999a)),
88    );
89    p = DoubleDouble::mul_f64_add(
90        p,
91        z,
92        DoubleDouble::from_bit_pair((0xbb3b08905e500000, 0xbfd0000000000000)),
93    );
94    p = DoubleDouble::mul_f64_add(
95        p,
96        z,
97        DoubleDouble::from_bit_pair((0x3c7555555567a1b0, 0x3fd5555555555555)),
98    );
99    let mut t = DoubleDouble::quick_mult(x2, p);
100    t = DoubleDouble::quick_mult_f64(t, z);
101    DoubleDouble::mul_f64_add(x2, -0.5, t)
102}
103
104#[inline]
105pub(crate) fn log_dd(x: f64) -> DoubleDouble {
106    let x_u = x.to_bits();
107    let mut m = x_u & 0xfffffffffffff;
108    let mut e: i64 = ((x_u >> 52) & 0x7ff) as i64;
109
110    let t;
111    if e != 0 {
112        t = m | (0x3ffu64 << 52);
113        m = m.wrapping_add(1u64 << 52);
114        e -= 0x3ff;
115    } else {
116        /* x is a subnormal double  */
117        let k = m.leading_zeros() - 11;
118
119        e = -0x3fei64 - k as i64;
120        m = m.wrapping_shl(k);
121        t = m | (0x3ffu64 << 52);
122    }
123
124    /* now |x| = 2^_e*_t = 2^(_e-52)*m with 1 <= _t < 2,
125    and 2^52 <= _m < 2^53 */
126
127    //   log(x) = log(t) + E · log(2)
128    let mut t = f64::from_bits(t);
129
130    // If m > sqrt(2) we divide it by 2 so ensure 1/sqrt(2) < t < sqrt(2)
131    let c: usize = (m >= 0x16a09e667f3bcd) as usize;
132    static CY: [f64; 2] = [1.0, 0.5];
133    static CM: [u64; 2] = [44, 45];
134
135    e = e.wrapping_add(c as i64);
136    let be = e;
137    let i = m >> CM[c];
138    t *= CY[c];
139
140    let r = f64::from_bits(POW_INVERSE[(i - 181) as usize]);
141    let log_r = DoubleDouble::from_bit_pair(LOG_NEG_DD[(i - 181) as usize]);
142
143    let z = f64::mul_add(r, t, -1.0);
144
145    const LOG2_DD: DoubleDouble = DoubleDouble::new(
146        f64::from_bits(0x3c7abc9e3b39803f),
147        f64::from_bits(0x3fe62e42fefa39ef),
148    );
149
150    let tt = DoubleDouble::mul_f64_add(LOG2_DD, be as f64, log_r);
151
152    let v = DoubleDouble::full_add_f64(tt, z);
153    let p = log_poly(z);
154    DoubleDouble::f64_add(v.hi, DoubleDouble::new(v.lo + p.lo, p.hi))
155}
156
157#[inline]
158pub(crate) fn log_dd_fast(x: f64) -> DoubleDouble {
159    let x_u = x.to_bits();
160    let mut m = x_u & 0xfffffffffffff;
161    let mut e: i64 = ((x_u >> 52) & 0x7ff) as i64;
162
163    let t;
164    if e != 0 {
165        t = m | (0x3ffu64 << 52);
166        m = m.wrapping_add(1u64 << 52);
167        e -= 0x3ff;
168    } else {
169        /* x is a subnormal double  */
170        let k = m.leading_zeros() - 11;
171
172        e = -0x3fei64 - k as i64;
173        m = m.wrapping_shl(k);
174        t = m | (0x3ffu64 << 52);
175    }
176
177    /* now |x| = 2^_e*_t = 2^(_e-52)*m with 1 <= _t < 2,
178    and 2^52 <= _m < 2^53 */
179
180    //   log(x) = log(t) + E · log(2)
181    let mut t = f64::from_bits(t);
182
183    // If m > sqrt(2) we divide it by 2 so ensure 1/sqrt(2) < t < sqrt(2)
184    let c: usize = (m >= 0x16a09e667f3bcd) as usize;
185    static CY: [f64; 2] = [1.0, 0.5];
186    static CM: [u64; 2] = [44, 45];
187
188    e = e.wrapping_add(c as i64);
189    let be = e;
190    let i = m >> CM[c];
191    t *= CY[c];
192
193    let r = f64::from_bits(POW_INVERSE[(i - 181) as usize]);
194    let log_r = DoubleDouble::from_bit_pair(LOG_NEG_DD[(i - 181) as usize]);
195
196    let z = f64::mul_add(r, t, -1.0);
197
198    const LOG2_DD: DoubleDouble = DoubleDouble::new(
199        f64::from_bits(0x3c7abc9e3b39803f),
200        f64::from_bits(0x3fe62e42fefa39ef),
201    );
202
203    let tt = DoubleDouble::mul_f64_add(LOG2_DD, be as f64, log_r);
204
205    let v = DoubleDouble::full_add_f64(tt, z);
206    let p = log_poly_fast(z);
207    DoubleDouble::f64_add(v.hi, DoubleDouble::new(v.lo + p.lo, p.hi))
208}
209
210#[cfg(test)]
211mod tests {
212    use super::*;
213
214    #[test]
215    fn test_log_dd() {
216        assert_eq!(log_dd(std::f64::consts::E).to_f64(), 1.);
217    }
218}
pxfm/logs/log_dd.rs

pxfm/logs/
log_dd.rs
