dashu_float/exp.rs
1use core::convert::TryInto;
2
3use crate::{
4 error::{assert_finite, assert_limited_precision, panic_power_negative_base},
5 fbig::FBig,
6 repr::{Context, Repr, Word},
7 round::{Round, Rounded},
8};
9use dashu_base::{AbsOrd, Approximation::*, BitTest, DivRemEuclid, EstimatedLog2, Sign};
10use dashu_int::IBig;
11
12impl<R: Round, const B: Word> FBig<R, B> {
13 /// Raise the floating point number to an integer power.
14 ///
15 /// # Examples
16 ///
17 /// ```
18 /// # use dashu_base::ParseError;
19 /// # use dashu_float::DBig;
20 /// let a = DBig::from_str_native("-1.234")?;
21 /// assert_eq!(a.powi(10.into()), DBig::from_str_native("8.188")?);
22 /// # Ok::<(), ParseError>(())
23 /// ```
24 #[inline]
25 pub fn powi(&self, exp: IBig) -> FBig<R, B> {
26 self.context.powi(&self.repr, exp).value()
27 }
28
29 /// Raise the floating point number to an floating point power.
30 ///
31 /// # Examples
32 ///
33 /// ```
34 /// # use dashu_base::ParseError;
35 /// # use dashu_float::DBig;
36 /// let x = DBig::from_str_native("1.23")?;
37 /// let y = DBig::from_str_native("-4.56")?;
38 /// assert_eq!(x.powf(&y), DBig::from_str_native("0.389")?);
39 /// # Ok::<(), ParseError>(())
40 /// ```
41 #[inline]
42 pub fn powf(&self, exp: &Self) -> Self {
43 let context = Context::max(self.context, exp.context);
44 context.powf(&self.repr, &exp.repr).value()
45 }
46
47 /// Calculate the exponential function (`eˣ`) on the floating point number.
48 ///
49 /// # Examples
50 ///
51 /// ```
52 /// # use dashu_base::ParseError;
53 /// # use dashu_float::DBig;
54 /// let a = DBig::from_str_native("-1.234")?;
55 /// assert_eq!(a.exp(), DBig::from_str_native("0.2911")?);
56 /// # Ok::<(), ParseError>(())
57 /// ```
58 #[inline]
59 pub fn exp(&self) -> FBig<R, B> {
60 self.context.exp(&self.repr).value()
61 }
62
63 /// Calculate the exponential minus one function (`eˣ-1`) on the floating point number.
64 ///
65 /// # Examples
66 ///
67 /// ```
68 /// # use dashu_base::ParseError;
69 /// # use dashu_float::DBig;
70 /// let a = DBig::from_str_native("-0.1234")?;
71 /// assert_eq!(a.exp_m1(), DBig::from_str_native("-0.11609")?);
72 /// # Ok::<(), ParseError>(())
73 /// ```
74 #[inline]
75 pub fn exp_m1(&self) -> FBig<R, B> {
76 self.context.exp_m1(&self.repr).value()
77 }
78}
79
80// TODO: give the exact formulation of required guard bits
81
82impl<R: Round> Context<R> {
83 /// Raise the floating point number to an integer power under this context.
84 ///
85 /// # Examples
86 ///
87 /// ```
88 /// # use dashu_base::ParseError;
89 /// # use dashu_float::DBig;
90 /// use dashu_base::Approximation::*;
91 /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
92 ///
93 /// let context = Context::<HalfAway>::new(2);
94 /// let a = DBig::from_str_native("-1.234")?;
95 /// assert_eq!(context.powi(&a.repr(), 10.into()), Inexact(DBig::from_str_native("8.2")?, AddOne));
96 /// # Ok::<(), ParseError>(())
97 /// ```
98 ///
99 /// # Panics
100 ///
101 /// Panics if the precision is unlimited and the exponent is negative. In this case, the exact
102 /// result is likely to have infinite digits.
103 pub fn powi<const B: Word>(&self, base: &Repr<B>, exp: IBig) -> Rounded<FBig<R, B>> {
104 assert_finite(base);
105
106 let (exp_sign, exp) = exp.into_parts();
107 if exp_sign == Sign::Negative {
108 // if the exponent is negative, then negate the exponent
109 // note that do the inverse at last requires less guard bits
110 assert_limited_precision(self.precision); // TODO: we can allow this if the inverse is exact (only when significand is one?)
111
112 let guard_bits = self.precision.bit_len() * 2; // heuristic
113 let rev_context = Context::<R::Reverse>::new(self.precision + guard_bits);
114 let pow = rev_context.powi(base, exp.into()).value();
115 let inv = rev_context.repr_div(Repr::one(), pow.repr);
116 let repr = inv.and_then(|v| self.repr_round(v));
117 return repr.map(|v| FBig::new(v, *self));
118 }
119 if exp.is_zero() {
120 return Exact(FBig::ONE);
121 } else if exp.is_one() {
122 let repr = self.repr_round_ref(base);
123 return repr.map(|v| FBig::new(v, *self));
124 }
125
126 let work_context = if self.is_limited() {
127 // increase working precision when the exponent is large
128 let guard_digits = exp.bit_len() + self.precision.bit_len(); // heuristic
129 Context::<R>::new(self.precision + guard_digits)
130 } else {
131 Context::<R>::new(0)
132 };
133
134 // binary exponentiation from left to right
135 let mut p = exp.bit_len() - 2;
136 let mut res = work_context.sqr(base);
137 loop {
138 if exp.bit(p) {
139 res = res.and_then(|v| work_context.mul(v.repr(), base));
140 }
141 if p == 0 {
142 break;
143 }
144 p -= 1;
145 res = res.and_then(|v| work_context.sqr(v.repr()));
146 }
147
148 res.and_then(|v| v.with_precision(self.precision))
149 }
150
151 /// Raise the floating point number to an floating point power under this context.
152 ///
153 /// Note that this method will not rely on [FBig::powi] even if the `exp` is actually an integer.
154 ///
155 /// # Examples
156 ///
157 /// ```
158 /// # use dashu_base::ParseError;
159 /// # use dashu_float::DBig;
160 /// use dashu_base::Approximation::*;
161 /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
162 ///
163 /// let context = Context::<HalfAway>::new(2);
164 /// let x = DBig::from_str_native("1.23")?;
165 /// let y = DBig::from_str_native("-4.56")?;
166 /// assert_eq!(context.powf(&x.repr(), &y.repr()), Inexact(DBig::from_str_native("0.39")?, AddOne));
167 /// # Ok::<(), ParseError>(())
168 /// ```
169 ///
170 /// # Panics
171 ///
172 /// Panics if the precision is unlimited.
173 pub fn powf<const B: Word>(&self, base: &Repr<B>, exp: &Repr<B>) -> Rounded<FBig<R, B>> {
174 assert_finite(base);
175 assert_limited_precision(self.precision); // TODO: we can allow it if exp is integer
176
177 // shortcuts
178 if exp.is_zero() {
179 return Exact(FBig::ONE);
180 } else if exp.is_one() {
181 let repr = self.repr_round_ref(base);
182 return repr.map(|v| FBig::new(v, *self));
183 } else if base.is_zero() {
184 return Exact(FBig::ZERO);
185 }
186 if base.sign() == Sign::Negative {
187 // TODO: we should allow negative base when exp is an integer
188 panic_power_negative_base()
189 }
190
191 // x^y = exp(y*ln(x)), use a simple rule for guard bits
192 let guard_digits = 10 + self.precision.log2_est() as usize;
193 let work_context = Context::<R>::new(self.precision + guard_digits);
194
195 let res = work_context
196 .ln(base)
197 .and_then(|v| work_context.mul(&v.repr, exp))
198 .and_then(|v| work_context.exp(&v.repr));
199 res.and_then(|v| v.with_precision(self.precision))
200 }
201
202 /// Calculate the exponential function (`eˣ`) on the floating point number under this context.
203 ///
204 /// # Examples
205 ///
206 /// ```
207 /// # use dashu_base::ParseError;
208 /// # use dashu_float::DBig;
209 /// use dashu_base::Approximation::*;
210 /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
211 ///
212 /// let context = Context::<HalfAway>::new(2);
213 /// let a = DBig::from_str_native("-1.234")?;
214 /// assert_eq!(context.exp(&a.repr()), Inexact(DBig::from_str_native("0.29")?, NoOp));
215 /// # Ok::<(), ParseError>(())
216 /// ```
217 #[inline]
218 pub fn exp<const B: Word>(&self, x: &Repr<B>) -> Rounded<FBig<R, B>> {
219 self.exp_internal(x, false)
220 }
221
222 /// Calculate the exponential minus one function (`eˣ-1`) on the floating point number under this context.
223 ///
224 /// # Examples
225 ///
226 /// ```
227 /// # use dashu_base::ParseError;
228 /// # use dashu_float::DBig;
229 /// use dashu_base::Approximation::*;
230 /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
231 ///
232 /// let context = Context::<HalfAway>::new(2);
233 /// let a = DBig::from_str_native("-0.1234")?;
234 /// assert_eq!(context.exp_m1(&a.repr()), Inexact(DBig::from_str_native("-0.12")?, SubOne));
235 /// # Ok::<(), ParseError>(())
236 /// ```
237 #[inline]
238 pub fn exp_m1<const B: Word>(&self, x: &Repr<B>) -> Rounded<FBig<R, B>> {
239 self.exp_internal(x, true)
240 }
241
242 // TODO: change reduction to (x - s log2) / 2ⁿ, so that the final powering is always base 2, and doesn't depends on powi.
243 // the powering exp(r)^(2ⁿ) could be optimized by noticing (1+x)^2 - 1 = x^2 + 2x
244 // consider this change after having a benchmark
245
246 fn exp_internal<const B: Word>(&self, x: &Repr<B>, minus_one: bool) -> Rounded<FBig<R, B>> {
247 assert_finite(x);
248 assert_limited_precision(self.precision);
249
250 if x.is_zero() {
251 return match minus_one {
252 false => Exact(FBig::ONE),
253 true => Exact(FBig::ZERO),
254 };
255 }
256
257 // A simple algorithm:
258 // - let r = (x - s logB) / Bⁿ, where s = floor(x / logB), such that r < B⁻ⁿ.
259 // - if the target precision is p digits, then there're only about p/m terms in Tyler series
260 // - finally, exp(x) = Bˢ * exp(r)^(Bⁿ)
261 // - the optimal n is √p as given by MPFR
262
263 // Maclaurin series: exp(r) = 1 + Σ(rⁱ/i!)
264 // There will be about p/log_B(r) summations when calculating the series, to prevent
265 // loss of significant, we needs about log_B(p) guard digits.
266 let series_guard_digits = (self.precision.log2_est() / B.log2_est()) as usize + 2;
267 let pow_guard_digits = (self.precision.bit_len() as f32 * B.log2_est() * 2.) as usize; // heuristic
268 let work_precision;
269
270 // When minus_one is true and |x| < 1/B, the input is fed into the Maclaurin series without scaling
271 let no_scaling = minus_one && x.log2_est() < -B.log2_est();
272 let (s, n, r) = if no_scaling {
273 // if minus_one is true and x is already small (x < 1/B),
274 // then directly evaluate the Maclaurin series without scaling
275 if x.sign() == Sign::Negative {
276 // extra digits are required to prevent cancellation during the summation
277 work_precision = self.precision + 2 * series_guard_digits;
278 } else {
279 work_precision = self.precision + series_guard_digits;
280 }
281 let context = Context::<R>::new(work_precision);
282 (0, 0, FBig::new(context.repr_round_ref(x).value(), context))
283 } else {
284 work_precision = self.precision + series_guard_digits + pow_guard_digits;
285 let context = Context::<R>::new(work_precision);
286 let x = FBig::new(context.repr_round_ref(x).value(), context);
287 let logb = context.ln_base::<B>();
288 let (s, r) = x.div_rem_euclid(logb);
289
290 // here m is roughly equal to sqrt(self.precision)
291 let n = 1usize << (self.precision.bit_len() / 2);
292 let s: isize = s.try_into().expect("exponent is too large");
293 (s, n, r)
294 };
295
296 let r = r >> n as isize;
297 let mut factorial = IBig::ONE;
298 let mut pow = r.clone();
299 let mut sum = if no_scaling {
300 r.clone()
301 } else {
302 FBig::ONE + &r
303 };
304
305 let mut k = 2;
306 loop {
307 factorial *= k;
308 pow *= &r;
309
310 let increase = &pow / &factorial;
311 if increase.abs_cmp(&sum.sub_ulp()).is_le() {
312 break;
313 }
314 sum += increase;
315 k += 1;
316 }
317
318 if no_scaling {
319 sum.with_precision(self.precision)
320 } else if minus_one {
321 // add extra digits to compensate for the subtraction
322 Context::<R>::new(self.precision + self.precision / 8 + 1) // heuristic
323 .powi(sum.repr(), Repr::<B>::BASE.pow(n).into())
324 .map(|v| (v << s) - FBig::ONE)
325 .and_then(|v| v.with_precision(self.precision))
326 } else {
327 self.powi(sum.repr(), Repr::<B>::BASE.pow(n).into())
328 .map(|v| v << s)
329 }
330 }
331}