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, FpError, FpResult},
11    fbig::FBig,
12    math::cache::{reborrow_cache, ConstCache},
13    repr::{Context, Repr, Word},
14    round::{Round, Rounded},
15    utils::ceil_usize,
16};
17use core::cmp::Ordering;
18
19impl<const B: Word> EstimatedLog2 for Repr<B> {
20    // currently a Word has at most 64 bits, so log2() < f32::MAX
21    fn log2_bounds(&self) -> (f32, f32) {
22        if self.significand.is_zero() {
23            return (f32::NEG_INFINITY, f32::NEG_INFINITY);
24        }
25
26        // log(s*B^e) = log(s) + e*log(B)
27        let (logs_lb, logs_ub) = self.significand.log2_bounds();
28        let (logb_lb, logb_ub) = if B.is_power_of_two() {
29            let log = B.trailing_zeros() as f32;
30            (log, log)
31        } else {
32            B.log2_bounds()
33        };
34        let e = self.exponent as f32;
35        let (lb, ub) = if self.exponent >= 0 {
36            (logs_lb + e * logb_lb, logs_ub + e * logb_ub)
37        } else {
38            (logs_lb + e * logb_ub, logs_ub + e * logb_lb)
39        };
40        (next_down(lb), next_up(ub))
41    }
42
43    fn log2_est(&self) -> f32 {
44        let logs = self.significand.log2_est();
45        let logb = if B.is_power_of_two() {
46            B.trailing_zeros() as f32
47        } else {
48            B.log2_est()
49        };
50        logs + self.exponent as f32 * logb
51    }
52}
53
54impl<R: Round, const B: Word> EstimatedLog2 for FBig<R, B> {
55    #[inline]
56    fn log2_bounds(&self) -> (f32, f32) {
57        self.repr.log2_bounds()
58    }
59
60    #[inline]
61    fn log2_est(&self) -> f32 {
62        self.repr.log2_est()
63    }
64}
65
66impl<R: Round, const B: Word> FBig<R, B> {
67    /// Calculate the natural logarithm function (`log(x)`) on the float number.
68    ///
69    /// # Examples
70    ///
71    /// ```
72    /// # use core::str::FromStr;
73    /// # use dashu_base::ParseError;
74    /// # use dashu_float::DBig;
75    /// let a = DBig::from_str("1.234")?;
76    /// assert_eq!(a.ln(), DBig::from_str("0.2103")?);
77    /// # Ok::<(), ParseError>(())
78    /// ```
79    #[inline]
80    pub fn ln(&self) -> Self {
81        self.context.unwrap_fp(self.context.ln(&self.repr, None))
82    }
83
84    /// Calculate the natural logarithm function (`log(x+1)`) on the float number
85    ///
86    /// # Examples
87    ///
88    /// ```
89    /// # use core::str::FromStr;
90    /// # use dashu_base::ParseError;
91    /// # use dashu_float::DBig;
92    /// let a = DBig::from_str("0.1234")?;
93    /// assert_eq!(a.ln_1p(), DBig::from_str("0.11636")?);
94    /// # Ok::<(), ParseError>(())
95    /// ```
96    #[inline]
97    pub fn ln_1p(&self) -> Self {
98        self.context.unwrap_fp(self.context.ln_1p(&self.repr, None))
99    }
100}
101
102impl<R: Round> Context<R> {
103    /// Calculate log(2)
104    ///
105    /// The precision of the output will be larger than self.precision
106    #[inline]
107    fn ln2<const B: Word>(&self, cache: Option<&mut ConstCache>) -> FBig<R, B> {
108        if let Some(c) = cache {
109            return c.ln2::<B, R>(self.precision);
110        }
111        // log(2) = 4L(6) + 2L(99)
112        // see formula (24) from Gourdon, Xavier, and Pascal Sebah.
113        // "The Logarithmic Constant: Log 2." (2004)
114        4 * self.iacoth(6.into()) + 2 * self.iacoth(99.into())
115    }
116
117    /// Calculate log(10)
118    ///
119    /// The precision of the output will be larger than self.precision
120    #[inline]
121    fn ln10<const B: Word>(&self, cache: Option<&mut ConstCache>) -> FBig<R, B> {
122        if let Some(c) = cache {
123            return c.ln10::<B, R>(self.precision);
124        }
125        // log(10) = log(2) + log(5) = 3log(2) + 2L(9)
126        3 * self.ln2(None) + 2 * self.iacoth(9.into())
127    }
128
129    /// Calculate log(B), for internal use only
130    ///
131    /// The precision of the output will be larger than self.precision
132    #[inline]
133    pub(crate) fn ln_base<const B: Word>(&self, cache: Option<&mut ConstCache>) -> FBig<R, B> {
134        if let Some(c) = cache {
135            return c.ln_base::<B, R>(self.precision);
136        }
137        match B {
138            2 => self.ln2(None),
139            10 => self.ln10(None),
140            i if i.is_power_of_two() => self.ln2(None) * i.trailing_zeros(),
141            _ => self.unwrap_fp(self.ln(&Repr::new(Repr::<B>::BASE.into(), 0), None)),
142        }
143    }
144
145    /// Calculate L(n) = acoth(n) = atanh(1/n) = 1/2 log((n+1)/(n-1)), given by the
146    /// series
147    ///
148    /// ```text
149    ///                1     n + 1              1
150    ///   atanh(1/n) = — log(—————) = Σ   ——————————————————
151    ///                2     n - 1   i≥0 n^(2i+1) · (2i+1)
152    /// ```
153    ///
154    /// This method is intended to be used in logarithm calculation,
155    /// so the precision of the output will be larger than desired precision.
156    ///
157    /// Evaluated by binary splitting (see [`iacoth_bs`][crate::math::cache::iacoth_bs]):
158    /// the exact integer tree state `(P, Q, T)` over `[1, N)` satisfies
159    /// `L(n) = (Q + T)/(n·Q)`, with `Q` kept at O(p) digits by the ratio-form
160    /// term recurrence.
161    fn iacoth<const B: Word>(&self, n: IBig) -> FBig<R, B> {
162        let n: u32 = (&n).try_into().expect("iacoth argument must fit in u32");
163
164        // number of series terms until r_k < B^{-p}:  (2k+1)·log_B(n) > p.
165        // The count is generously over-provisioned, so a truncating cast stands in
166        // for a ceiling.
167        let log_b_n = n.log2_est() / B.log2_est();
168        let num_terms = (self.precision as f32 / (2.0 * log_b_n)) as usize + 10;
169
170        let (_p, q, t) = crate::math::cache::iacoth_bs(n, 1, num_terms + 1);
171
172        // L(n) = (Q + T) / (n·Q). Extra guard digits absorb the division's rounding
173        // (the binary-splitting state is exact, so only this single round loses anything).
174        let guard_digits = ceil_usize(self.precision.log2_est() / B.log2_est());
175        let work_context = Self::new(self.precision + guard_digits + 2);
176
177        let num = work_context.convert_int::<B>(q.as_ibig() + &t).value();
178        let denom = work_context.convert_int::<B>(IBig::from(n) * &q).value();
179        num / denom
180    }
181
182    /// Calculate the natural logarithm function (`log(x)`) on the float number under this context.
183    ///
184    /// # Examples
185    ///
186    /// ```
187    /// # use core::str::FromStr;
188    /// # use dashu_base::ParseError;
189    /// # use dashu_float::DBig;
190    /// use dashu_base::Approximation::*;
191    /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
192    ///
193    /// let context = Context::<HalfAway>::new(2);
194    /// let a = DBig::from_str("1.234")?;
195    /// assert_eq!(context.ln(&a.repr(), None), Ok(Inexact(DBig::from_str("0.21")?, NoOp)));
196    /// # Ok::<(), ParseError>(())
197    /// ```
198    #[inline]
199    pub fn ln<const B: Word>(
200        &self,
201        x: &Repr<B>,
202        cache: Option<&mut ConstCache>,
203    ) -> FpResult<FBig<R, B>> {
204        if x.is_infinite() {
205            return Err(FpError::InfiniteInput);
206        }
207        if x.significand.is_zero() {
208            // ln(±0) = -inf (a value, not an error)
209            return Ok(Exact(FBig::new(Repr::neg_infinity(), *self)));
210        }
211        if x.sign() == Sign::Negative {
212            return Err(FpError::OutOfDomain);
213        }
214        Ok(self.ln_internal(x, false, cache))
215    }
216
217    /// Calculate the natural logarithm function (`log(x+1)`) on the float number under this context.
218    ///
219    /// # Examples
220    ///
221    /// ```
222    /// # use core::str::FromStr;
223    /// # use dashu_base::ParseError;
224    /// # use dashu_float::DBig;
225    /// use dashu_base::Approximation::*;
226    /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
227    ///
228    /// let context = Context::<HalfAway>::new(2);
229    /// let a = DBig::from_str("0.1234")?;
230    /// assert_eq!(context.ln_1p(&a.repr(), None), Ok(Inexact(DBig::from_str("0.12")?, AddOne)));
231    /// # Ok::<(), ParseError>(())
232    /// ```
233    #[inline]
234    pub fn ln_1p<const B: Word>(
235        &self,
236        x: &Repr<B>,
237        cache: Option<&mut ConstCache>,
238    ) -> FpResult<FBig<R, B>> {
239        if x.is_infinite() {
240            return Err(FpError::InfiniteInput);
241        }
242        // Domain of ln_1p is x > -1. x == -1 gives -inf; x < -1 is out of domain.
243        if x.sign() == Sign::Negative && !x.significand.is_zero() {
244            match FBig::<R, B>::new(x.clone(), *self).abs_cmp(&FBig::ONE) {
245                Ordering::Greater => return Err(FpError::OutOfDomain), // x < -1
246                Ordering::Equal => return Ok(Exact(FBig::new(Repr::neg_infinity(), *self))),
247                _ => {}
248            }
249        }
250        Ok(self.ln_internal(x, true, cache))
251    }
252
253    fn ln_internal<const B: Word>(
254        &self,
255        x: &Repr<B>,
256        one_plus: bool,
257        mut cache: Option<&mut ConstCache>,
258    ) -> Rounded<FBig<R, B>> {
259        assert_finite(x);
260        assert_limited_precision(self.precision);
261
262        if !one_plus && x.is_one() {
263            return Exact(FBig::ZERO); // ln(1) = +0
264        }
265        if one_plus && x.significand.is_zero() {
266            // ln_1p(±0) = ±0
267            let zero = if x.is_neg_zero() {
268                FBig::new(Repr::neg_zero(), *self)
269            } else {
270                FBig::ZERO
271            };
272            return Exact(zero);
273        }
274
275        // A simple algorithm:
276        // - let log(x) = log(x/2^s) + slog2 where s = floor(log2(x))
277        // - such that x*2^s is close to but larger than 1 (and x*2^s < 2)
278        let guard_digits = ceil_usize(self.precision.log2_est() / B.log2_est()) + 2;
279        let mut work_precision = self.precision + guard_digits + one_plus as usize;
280        let context = Context::<R>::new(work_precision);
281        let x = FBig::new(context.repr_round_ref(x).value(), context);
282
283        // When one_plus is true and |x| < 1/B, the input is fed into the Maclaurin without scaling
284        let no_scaling = one_plus && x.log2_est() < -B.log2_est();
285
286        let (s, mut x_scaled) = if no_scaling {
287            (0, x)
288        } else {
289            let x = if one_plus { x + FBig::ONE } else { x };
290
291            let log2 = x.log2_bounds().0;
292            let s = log2 as isize - (log2 < 0.) as isize; // floor(log2(x))
293
294            let x_scaled = if B == 2 {
295                x >> s
296            } else if s > 0 {
297                x / (IBig::ONE << s as usize)
298            } else {
299                x * (IBig::ONE << (-s) as usize)
300            };
301            debug_assert!(x_scaled >= FBig::<R, B>::ONE);
302            (s, x_scaled)
303        };
304
305        if s < 0 || x_scaled.repr.sign() == Sign::Negative {
306            // when s or x_scaled is negative, the final addition is actually a subtraction,
307            // therefore we need to double the precision to get the correct result
308            work_precision += self.precision;
309            x_scaled.context.precision = work_precision;
310        };
311        let work_context = Context::new(work_precision);
312
313        // after the number is scaled to nearly one, use Maclaurin series on log(x) = 2atanh(z):
314        // let z = (x-1)/(x+1) < 1, log(x) = 2atanh(z) = 2Σ(z²ⁱ⁺¹/(2i+1)) for i = 1,3,5,...
315        // similar to iacoth, the required iterations stop at i = -p/2log_B(z), and we need log_B(i) guard bits
316        let z = if no_scaling {
317            let d = &x_scaled + (FBig::ONE + FBig::ONE);
318            x_scaled / d
319        } else {
320            (&x_scaled - FBig::ONE) / (x_scaled + FBig::ONE)
321        };
322        let z2 = z.sqr();
323        let mut pow = z.clone();
324        let mut sum = z;
325
326        let mut k: usize = 3;
327        loop {
328            pow *= &z2;
329
330            let increase = &pow / work_context.convert_int::<B>(k.into()).value();
331            if increase.abs_cmp(&sum.sub_ulp()).is_le() {
332                break;
333            }
334
335            sum += increase;
336            k += 2;
337        }
338
339        // compose the logarithm of the original number
340        let result: FBig<R, B> = if no_scaling {
341            2 * sum
342        } else {
343            2 * sum + s * work_context.ln2::<B>(reborrow_cache(&mut cache))
344        };
345        result.with_precision(self.precision)
346    }
347}
348
349#[cfg(test)]
350mod tests {
351    use super::*;
352    use crate::round::mode;
353
354    #[test]
355    fn test_ln_zero_is_neg_infinity() {
356        let ctx = Context::<mode::HalfEven>::new(53);
357        let r = ctx.ln::<2>(&Repr::<2>::zero(), None).unwrap().value();
358        assert!(r.repr().is_infinite());
359        assert_eq!(r.repr().sign(), Sign::Negative);
360    }
361
362    #[test]
363    fn test_iacoth() {
364        let context = Context::<mode::Zero>::new(10);
365        let binary_6 = context.iacoth::<2>(6.into()).with_precision(10).value();
366        assert_eq!(binary_6.repr.significand, IBig::from(689));
367        let decimal_6 = context.iacoth::<10>(6.into()).with_precision(10).value();
368        assert_eq!(decimal_6.repr.significand, IBig::from(1682361183));
369
370        let context = Context::<mode::Zero>::new(40);
371        let decimal_6 = context.iacoth::<10>(6.into()).with_precision(40).value();
372        assert_eq!(
373            decimal_6.repr.significand,
374            IBig::from_str_radix("1682361183106064652522967051084960450557", 10).unwrap()
375        );
376
377        let context = Context::<mode::Zero>::new(201);
378        let binary_6 = context.iacoth::<2>(6.into()).with_precision(201).value();
379        assert_eq!(
380            binary_6.repr.significand,
381            IBig::from_str_radix(
382                "2162760151454160450909229890833066944953539957685348083415205",
383                10
384            )
385            .unwrap()
386        );
387    }
388
389    #[test]
390    fn test_ln2_ln10() {
391        let context = Context::<mode::Zero>::new(45);
392        let decimal_ln2 = context.ln2::<10>(None).with_precision(45).value();
393        assert_eq!(
394            decimal_ln2.repr.significand,
395            IBig::from_str_radix("693147180559945309417232121458176568075500134", 10).unwrap()
396        );
397        let decimal_ln10 = context.ln10::<10>(None).with_precision(45).value();
398        assert_eq!(
399            decimal_ln10.repr.significand,
400            IBig::from_str_radix("230258509299404568401799145468436420760110148", 10).unwrap()
401        );
402
403        let context = Context::<mode::Zero>::new(180);
404        let binary_ln2 = context.ln2::<2>(None).with_precision(180).value();
405        assert_eq!(
406            binary_ln2.repr.significand,
407            IBig::from_str_radix("1062244963371879310175186301324412638028404515790072203", 10)
408                .unwrap()
409        );
410        let binary_ln10 = context.ln10::<2>(None).with_precision(180).value();
411        assert_eq!(
412            binary_ln10.repr.significand,
413            IBig::from_str_radix("882175346869410758689845931257775553286341791676474847", 10)
414                .unwrap()
415        );
416    }
417}
dashu_float/log.rs

dashu_float/
log.rs
