dashu_float/
log.rs

1use dashu_base::{
2    utils::{next_down, next_up},
3    AbsOrd,
4    Approximation::*,
5    EstimatedLog2, Sign,
6};
7use dashu_int::IBig;
8
9use crate::{
10    error::{assert_finite, assert_limited_precision},
11    fbig::FBig,
12    repr::{Context, Repr, Word},
13    round::{Round, Rounded},
14};
15
16impl<const B: Word> EstimatedLog2 for Repr<B> {
17    // currently a Word has at most 64 bits, so log2() < f32::MAX
18    fn log2_bounds(&self) -> (f32, f32) {
19        if self.significand.is_zero() {
20            return (f32::NEG_INFINITY, f32::NEG_INFINITY);
21        }
22
23        // log(s*B^e) = log(s) + e*log(B)
24        let (logs_lb, logs_ub) = self.significand.log2_bounds();
25        let (logb_lb, logb_ub) = if B.is_power_of_two() {
26            let log = B.trailing_zeros() as f32;
27            (log, log)
28        } else {
29            B.log2_bounds()
30        };
31        let e = self.exponent as f32;
32        let (lb, ub) = if self.exponent >= 0 {
33            (logs_lb + e * logb_lb, logs_ub + e * logb_ub)
34        } else {
35            (logs_lb + e * logb_ub, logs_ub + e * logb_lb)
36        };
37        (next_down(lb), next_up(ub))
38    }
39
40    fn log2_est(&self) -> f32 {
41        let logs = self.significand.log2_est();
42        let logb = if B.is_power_of_two() {
43            B.trailing_zeros() as f32
44        } else {
45            B.log2_est()
46        };
47        logs + self.exponent as f32 * logb
48    }
49}
50
51impl<R: Round, const B: Word> EstimatedLog2 for FBig<R, B> {
52    #[inline]
53    fn log2_bounds(&self) -> (f32, f32) {
54        self.repr.log2_bounds()
55    }
56
57    #[inline]
58    fn log2_est(&self) -> f32 {
59        self.repr.log2_est()
60    }
61}
62
63impl<R: Round, const B: Word> FBig<R, B> {
64    /// Calculate the natural logarithm function (`log(x)`) on the float number.
65    ///
66    /// # Examples
67    ///
68    /// ```
69    /// # use core::str::FromStr;
70    /// # use dashu_base::ParseError;
71    /// # use dashu_float::DBig;
72    /// let a = DBig::from_str("1.234")?;
73    /// assert_eq!(a.ln(), DBig::from_str("0.2103")?);
74    /// # Ok::<(), ParseError>(())
75    /// ```
76    #[inline]
77    pub fn ln(&self) -> Self {
78        self.context.ln(&self.repr).value()
79    }
80
81    /// Calculate the natural logarithm function (`log(x+1)`) on the float number
82    ///
83    /// # Examples
84    ///
85    /// ```
86    /// # use core::str::FromStr;
87    /// # use dashu_base::ParseError;
88    /// # use dashu_float::DBig;
89    /// let a = DBig::from_str("0.1234")?;
90    /// assert_eq!(a.ln_1p(), DBig::from_str("0.11636")?);
91    /// # Ok::<(), ParseError>(())
92    /// ```
93    #[inline]
94    pub fn ln_1p(&self) -> Self {
95        self.context.ln_1p(&self.repr).value()
96    }
97}
98
99impl<R: Round> Context<R> {
100    /// Calculate log(2)
101    ///
102    /// The precision of the output will be larger than self.precision
103    #[inline]
104    fn ln2<const B: Word>(&self) -> FBig<R, B> {
105        // log(2) = 4L(6) + 2L(99)
106        // see formula (24) from Gourdon, Xavier, and Pascal Sebah.
107        // "The Logarithmic Constant: Log 2." (2004)
108        4 * self.iacoth(6.into()) + 2 * self.iacoth(99.into())
109    }
110
111    /// Calculate log(2)
112    ///
113    /// The precision of the output will be larger than self.precision
114    #[inline]
115    fn ln10<const B: Word>(&self) -> FBig<R, B> {
116        // log(10) = log(2) + log(5) = 3log(2) + 2L(9)
117        3 * self.ln2() + 2 * self.iacoth(9.into())
118    }
119
120    /// Calculate log(B), for internal use only
121    ///
122    /// The precision of the output will be larger than self.precision
123    #[inline]
124    pub(crate) fn ln_base<const B: Word>(&self) -> FBig<R, B> {
125        match B {
126            2 => self.ln2(),
127            10 => self.ln10(),
128            i if i.is_power_of_two() => self.ln2() * i.trailing_zeros(),
129            _ => self.ln(&Repr::new(Repr::<B>::BASE.into(), 0)).value(),
130        }
131    }
132
133    /// Calculate L(n) = acoth(n) = atanh(1/n) = 1/2 log((n+1)/(n-1))
134    ///
135    /// This method is intended to be used in logarithm calculation,
136    /// so the precision of the output will be larger than desired precision.
137    fn iacoth<const B: Word>(&self, n: IBig) -> FBig<R, B> {
138        /*
139         * use Maclaurin series:
140         *       1    1     n+1             1
141         * atanh(—) = — log(———) =  Σ  ———————————
142         *       n    2     n-1    i≥0 n²ⁱ⁺¹(2i+1)
143         *
144         * Therefore to achieve precision B^p, the series should be stopped at
145         *    n²ⁱ⁺¹(2i+1) / n = B^p
146         * => 2i*ln(n) + ln(2i+1) = p ln(B)
147         * ~> 2i*ln(n) = p ln(B)
148         * => 2i = p/log_B(n)
149         *
150         * There will be i summations when calculating the series, to prevent
151         * loss of significant, we needs about log_B(i) guard digits.
152         *    log_B(i)
153         * <= log_B(p/2log_B(n))
154         *  = log_B(p/2) - log_B(log_B(n))
155         * <= log_B(p/2)
156         */
157
158        // extras digits are added to ensure precise result
159        // TODO: test if we can use log_B(p/2log_B(n)) directly
160        let guard_digits = (self.precision.log2_est() / B.log2_est()) as usize;
161        let work_context = Self::new(self.precision + guard_digits + 2);
162
163        let n = work_context.convert_int(n).value();
164        let inv = FBig::ONE / n;
165        let inv2 = inv.sqr();
166        let mut sum = inv.clone();
167        let mut pow = inv;
168
169        let mut k: usize = 3;
170        loop {
171            pow *= &inv2;
172
173            let increase = &pow / work_context.convert_int::<B>(k.into()).value();
174            if increase < sum.sub_ulp() {
175                return sum;
176            }
177
178            sum += increase;
179            k += 2;
180        }
181    }
182
183    /// Calculate the natural logarithm function (`log(x)`) on the float number under this context.
184    ///
185    /// # Examples
186    ///
187    /// ```
188    /// # use core::str::FromStr;
189    /// # use dashu_base::ParseError;
190    /// # use dashu_float::DBig;
191    /// use dashu_base::Approximation::*;
192    /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
193    ///
194    /// let context = Context::<HalfAway>::new(2);
195    /// let a = DBig::from_str("1.234")?;
196    /// assert_eq!(context.ln(&a.repr()), Inexact(DBig::from_str("0.21")?, NoOp));
197    /// # Ok::<(), ParseError>(())
198    /// ```
199    #[inline]
200    pub fn ln<const B: Word>(&self, x: &Repr<B>) -> Rounded<FBig<R, B>> {
201        self.ln_internal(x, false)
202    }
203
204    /// Calculate the natural logarithm function (`log(x+1)`) on the float number under this context.
205    ///
206    /// # Examples
207    ///
208    /// ```
209    /// # use core::str::FromStr;
210    /// # use dashu_base::ParseError;
211    /// # use dashu_float::DBig;
212    /// use dashu_base::Approximation::*;
213    /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
214    ///
215    /// let context = Context::<HalfAway>::new(2);
216    /// let a = DBig::from_str("0.1234")?;
217    /// assert_eq!(context.ln_1p(&a.repr()), Inexact(DBig::from_str("0.12")?, AddOne));
218    /// # Ok::<(), ParseError>(())
219    /// ```
220    #[inline]
221    pub fn ln_1p<const B: Word>(&self, x: &Repr<B>) -> Rounded<FBig<R, B>> {
222        self.ln_internal(x, true)
223    }
224
225    fn ln_internal<const B: Word>(&self, x: &Repr<B>, one_plus: bool) -> Rounded<FBig<R, B>> {
226        assert_finite(x);
227        assert_limited_precision(self.precision);
228
229        if (one_plus && x.is_zero()) || (!one_plus && x.is_one()) {
230            return Exact(FBig::ZERO);
231        }
232
233        // A simple algorithm:
234        // - let log(x) = log(x/2^s) + slog2 where s = floor(log2(x))
235        // - such that x*2^s is close to but larger than 1 (and x*2^s < 2)
236        let guard_digits = (self.precision.log2_est() / B.log2_est()) as usize + 2;
237        let mut work_precision = self.precision + guard_digits + one_plus as usize;
238        let context = Context::<R>::new(work_precision);
239        let x = FBig::new(context.repr_round_ref(x).value(), context);
240
241        // When one_plus is true and |x| < 1/B, the input is fed into the Maclaurin without scaling
242        let no_scaling = one_plus && x.log2_est() < -B.log2_est();
243
244        let (s, mut x_scaled) = if no_scaling {
245            (0, x)
246        } else {
247            let x = if one_plus { x + FBig::ONE } else { x };
248
249            let log2 = x.log2_bounds().0;
250            let s = log2 as isize - (log2 < 0.) as isize; // floor(log2(x))
251
252            let x_scaled = if B == 2 {
253                x >> s
254            } else if s > 0 {
255                x / (IBig::ONE << s as usize)
256            } else {
257                x * (IBig::ONE << (-s) as usize)
258            };
259            debug_assert!(x_scaled >= FBig::<R, B>::ONE);
260            (s, x_scaled)
261        };
262
263        if s < 0 || x_scaled.repr.sign() == Sign::Negative {
264            // when s or x_scaled is negative, the final addition is actually a subtraction,
265            // therefore we need to double the precision to get the correct result
266            work_precision += self.precision;
267            x_scaled.context.precision = work_precision;
268        };
269        let work_context = Context::new(work_precision);
270
271        // after the number is scaled to nearly one, use Maclaurin series on log(x) = 2atanh(z):
272        // let z = (x-1)/(x+1) < 1, log(x) = 2atanh(z) = 2Σ(z²ⁱ⁺¹/(2i+1)) for i = 1,3,5,...
273        // similar to iacoth, the required iterations stop at i = -p/2log_B(z), and we need log_B(i) guard bits
274        let z = if no_scaling {
275            let d = &x_scaled + (FBig::ONE + FBig::ONE);
276            x_scaled / d
277        } else {
278            (&x_scaled - FBig::ONE) / (x_scaled + FBig::ONE)
279        };
280        let z2 = z.sqr();
281        let mut pow = z.clone();
282        let mut sum = z;
283
284        let mut k: usize = 3;
285        loop {
286            pow *= &z2;
287
288            let increase = &pow / work_context.convert_int::<B>(k.into()).value();
289            if increase.abs_cmp(&sum.sub_ulp()).is_le() {
290                break;
291            }
292
293            sum += increase;
294            k += 2;
295        }
296
297        // compose the logarithm of the original number
298        let result: FBig<R, B> = if no_scaling {
299            2 * sum
300        } else {
301            2 * sum + s * work_context.ln2()
302        };
303        result.with_precision(self.precision)
304    }
305}
306
307#[cfg(test)]
308mod tests {
309    use super::*;
310    use crate::round::mode;
311
312    #[test]
313    fn test_iacoth() {
314        let context = Context::<mode::Zero>::new(10);
315        let binary_6 = context.iacoth::<2>(6.into()).with_precision(10).value();
316        assert_eq!(binary_6.repr.significand, IBig::from(689));
317        let decimal_6 = context.iacoth::<10>(6.into()).with_precision(10).value();
318        assert_eq!(decimal_6.repr.significand, IBig::from(1682361183));
319
320        let context = Context::<mode::Zero>::new(40);
321        let decimal_6 = context.iacoth::<10>(6.into()).with_precision(40).value();
322        assert_eq!(
323            decimal_6.repr.significand,
324            IBig::from_str_radix("1682361183106064652522967051084960450557", 10).unwrap()
325        );
326
327        let context = Context::<mode::Zero>::new(201);
328        let binary_6 = context.iacoth::<2>(6.into()).with_precision(201).value();
329        assert_eq!(
330            binary_6.repr.significand,
331            IBig::from_str_radix(
332                "2162760151454160450909229890833066944953539957685348083415205",
333                10
334            )
335            .unwrap()
336        );
337    }
338
339    #[test]
340    fn test_ln2_ln10() {
341        let context = Context::<mode::Zero>::new(45);
342        let decimal_ln2 = context.ln2::<10>().with_precision(45).value();
343        assert_eq!(
344            decimal_ln2.repr.significand,
345            IBig::from_str_radix("693147180559945309417232121458176568075500134", 10).unwrap()
346        );
347        let decimal_ln10 = context.ln10::<10>().with_precision(45).value();
348        assert_eq!(
349            decimal_ln10.repr.significand,
350            IBig::from_str_radix("230258509299404568401799145468436420760110148", 10).unwrap()
351        );
352
353        let context = Context::<mode::Zero>::new(180);
354        let binary_ln2 = context.ln2::<2>().with_precision(180).value();
355        assert_eq!(
356            binary_ln2.repr.significand,
357            IBig::from_str_radix("1062244963371879310175186301324412638028404515790072203", 10)
358                .unwrap()
359        );
360        let binary_ln10 = context.ln10::<2>().with_precision(180).value();
361        assert_eq!(
362            binary_ln10.repr.significand,
363            IBig::from_str_radix("882175346869410758689845931257775553286341791676474847", 10)
364                .unwrap()
365        );
366    }
367}
dashu_float/log.rs

dashu_float/
log.rs
