1 use crate::prelude::*; 2 use num_traits::pow::Pow; 3 4 // Tolerance for inaccuracies when calculating exp 5 const EXP_TOLERANCE: Decimal = Decimal::from_parts(2, 0, 0, false, 7); 6 // Approximation of 1/ln(10) = 0.4342944819032518276511289189 7 const LN10_INVERSE: Decimal = Decimal::from_parts_raw(1763037029, 1670682625, 235431510, 1835008); 8 // Total iterations of taylor series for Trig. 9 const TRIG_SERIES_UPPER_BOUND: usize = 6; 10 // PI / 8 11 const EIGHTH_PI: Decimal = Decimal::from_parts_raw(2822163429, 3244459792, 212882598, 1835008); 12 13 // Table representing {index}! 14 const FACTORIAL: [Decimal; 28] = [ 15 Decimal::from_parts(1, 0, 0, false, 0), 16 Decimal::from_parts(1, 0, 0, false, 0), 17 Decimal::from_parts(2, 0, 0, false, 0), 18 Decimal::from_parts(6, 0, 0, false, 0), 19 Decimal::from_parts(24, 0, 0, false, 0), 20 // 5! 21 Decimal::from_parts(120, 0, 0, false, 0), 22 Decimal::from_parts(720, 0, 0, false, 0), 23 Decimal::from_parts(5040, 0, 0, false, 0), 24 Decimal::from_parts(40320, 0, 0, false, 0), 25 Decimal::from_parts(362880, 0, 0, false, 0), 26 // 10! 27 Decimal::from_parts(3628800, 0, 0, false, 0), 28 Decimal::from_parts(39916800, 0, 0, false, 0), 29 Decimal::from_parts(479001600, 0, 0, false, 0), 30 Decimal::from_parts(1932053504, 1, 0, false, 0), 31 Decimal::from_parts(1278945280, 20, 0, false, 0), 32 // 15! 33 Decimal::from_parts(2004310016, 304, 0, false, 0), 34 Decimal::from_parts(2004189184, 4871, 0, false, 0), 35 Decimal::from_parts(4006445056, 82814, 0, false, 0), 36 Decimal::from_parts(3396534272, 1490668, 0, false, 0), 37 Decimal::from_parts(109641728, 28322707, 0, false, 0), 38 // 20! 39 Decimal::from_parts(2192834560, 566454140, 0, false, 0), 40 Decimal::from_parts(3099852800, 3305602358, 2, false, 0), 41 Decimal::from_parts(3772252160, 4003775155, 60, false, 0), 42 Decimal::from_parts(862453760, 1892515369, 1401, false, 0), 43 Decimal::from_parts(3519021056, 2470695900, 33634, false, 0), 44 // 25! 45 Decimal::from_parts(2076180480, 1637855376, 840864, false, 0), 46 Decimal::from_parts(2441084928, 3929534124, 21862473, false, 0), 47 Decimal::from_parts(1484783616, 3018206259, 590286795, false, 0), 48 ]; 49 50 /// Trait exposing various mathematical operations that can be applied using a Decimal. This is only 51 /// present when the `maths` feature has been enabled. 52 pub trait MathematicalOps { 53 /// The estimated exponential function, e<sup>x</sup>. Stops calculating when it is within 54 /// tolerance of roughly `0.0000002`. exp(&self) -> Decimal55 fn exp(&self) -> Decimal; 56 57 /// The estimated exponential function, e<sup>x</sup>. Stops calculating when it is within 58 /// tolerance of roughly `0.0000002`. Returns `None` on overflow. checked_exp(&self) -> Option<Decimal>59 fn checked_exp(&self) -> Option<Decimal>; 60 61 /// The estimated exponential function, e<sup>x</sup> using the `tolerance` provided as a hint 62 /// as to when to stop calculating. A larger tolerance will cause the number to stop calculating 63 /// sooner at the potential cost of a slightly less accurate result. exp_with_tolerance(&self, tolerance: Decimal) -> Decimal64 fn exp_with_tolerance(&self, tolerance: Decimal) -> Decimal; 65 66 /// The estimated exponential function, e<sup>x</sup> using the `tolerance` provided as a hint 67 /// as to when to stop calculating. A larger tolerance will cause the number to stop calculating 68 /// sooner at the potential cost of a slightly less accurate result. 69 /// Returns `None` on overflow. checked_exp_with_tolerance(&self, tolerance: Decimal) -> Option<Decimal>70 fn checked_exp_with_tolerance(&self, tolerance: Decimal) -> Option<Decimal>; 71 72 /// Raise self to the given integer exponent: x<sup>y</sup> powi(&self, exp: i64) -> Decimal73 fn powi(&self, exp: i64) -> Decimal; 74 75 /// Raise self to the given integer exponent x<sup>y</sup> returning `None` on overflow. checked_powi(&self, exp: i64) -> Option<Decimal>76 fn checked_powi(&self, exp: i64) -> Option<Decimal>; 77 78 /// Raise self to the given unsigned integer exponent: x<sup>y</sup> powu(&self, exp: u64) -> Decimal79 fn powu(&self, exp: u64) -> Decimal; 80 81 /// Raise self to the given unsigned integer exponent x<sup>y</sup> returning `None` on overflow. checked_powu(&self, exp: u64) -> Option<Decimal>82 fn checked_powu(&self, exp: u64) -> Option<Decimal>; 83 84 /// Raise self to the given floating point exponent: x<sup>y</sup> powf(&self, exp: f64) -> Decimal85 fn powf(&self, exp: f64) -> Decimal; 86 87 /// Raise self to the given floating point exponent x<sup>y</sup> returning `None` on overflow. checked_powf(&self, exp: f64) -> Option<Decimal>88 fn checked_powf(&self, exp: f64) -> Option<Decimal>; 89 90 /// Raise self to the given Decimal exponent: x<sup>y</sup>. If `exp` is not whole then the approximation 91 /// e<sup>y*ln(x)</sup> is used. powd(&self, exp: Decimal) -> Decimal92 fn powd(&self, exp: Decimal) -> Decimal; 93 94 /// Raise self to the given Decimal exponent x<sup>y</sup> returning `None` on overflow. 95 /// If `exp` is not whole then the approximation e<sup>y*ln(x)</sup> is used. checked_powd(&self, exp: Decimal) -> Option<Decimal>96 fn checked_powd(&self, exp: Decimal) -> Option<Decimal>; 97 98 /// The square root of a Decimal. Uses a standard Babylonian method. sqrt(&self) -> Option<Decimal>99 fn sqrt(&self) -> Option<Decimal>; 100 101 /// Calculates the natural logarithm for a Decimal calculated using Taylor's series. ln(&self) -> Decimal102 fn ln(&self) -> Decimal; 103 104 /// Calculates the checked natural logarithm for a Decimal calculated using Taylor's series. 105 /// Returns `None` for negative numbers or zero. checked_ln(&self) -> Option<Decimal>106 fn checked_ln(&self) -> Option<Decimal>; 107 108 /// Calculates the base 10 logarithm of a specified Decimal number. log10(&self) -> Decimal109 fn log10(&self) -> Decimal; 110 111 /// Calculates the checked base 10 logarithm of a specified Decimal number. 112 /// Returns `None` for negative numbers or zero. checked_log10(&self) -> Option<Decimal>113 fn checked_log10(&self) -> Option<Decimal>; 114 115 /// Abramowitz Approximation of Error Function from [wikipedia](https://en.wikipedia.org/wiki/Error_function#Numerical_approximations) erf(&self) -> Decimal116 fn erf(&self) -> Decimal; 117 118 /// The Cumulative distribution function for a Normal distribution norm_cdf(&self) -> Decimal119 fn norm_cdf(&self) -> Decimal; 120 121 /// The Probability density function for a Normal distribution. norm_pdf(&self) -> Decimal122 fn norm_pdf(&self) -> Decimal; 123 124 /// The Probability density function for a Normal distribution returning `None` on overflow. checked_norm_pdf(&self) -> Option<Decimal>125 fn checked_norm_pdf(&self) -> Option<Decimal>; 126 127 /// Computes the sine of a number (in radians). 128 /// Panics upon overflow. sin(&self) -> Decimal129 fn sin(&self) -> Decimal; 130 131 /// Computes the checked sine of a number (in radians). checked_sin(&self) -> Option<Decimal>132 fn checked_sin(&self) -> Option<Decimal>; 133 134 /// Computes the cosine of a number (in radians). 135 /// Panics upon overflow. cos(&self) -> Decimal136 fn cos(&self) -> Decimal; 137 138 /// Computes the checked cosine of a number (in radians). checked_cos(&self) -> Option<Decimal>139 fn checked_cos(&self) -> Option<Decimal>; 140 141 /// Computes the tangent of a number (in radians). 142 /// Panics upon overflow or upon approaching a limit. tan(&self) -> Decimal143 fn tan(&self) -> Decimal; 144 145 /// Computes the checked tangent of a number (in radians). 146 /// Returns None on limit. checked_tan(&self) -> Option<Decimal>147 fn checked_tan(&self) -> Option<Decimal>; 148 } 149 150 impl MathematicalOps for Decimal { exp(&self) -> Decimal151 fn exp(&self) -> Decimal { 152 self.exp_with_tolerance(EXP_TOLERANCE) 153 } 154 checked_exp(&self) -> Option<Decimal>155 fn checked_exp(&self) -> Option<Decimal> { 156 self.checked_exp_with_tolerance(EXP_TOLERANCE) 157 } 158 exp_with_tolerance(&self, tolerance: Decimal) -> Decimal159 fn exp_with_tolerance(&self, tolerance: Decimal) -> Decimal { 160 match self.checked_exp_with_tolerance(tolerance) { 161 Some(d) => d, 162 None => { 163 if self.is_sign_negative() { 164 panic!("Exp underflowed") 165 } else { 166 panic!("Exp overflowed") 167 } 168 } 169 } 170 } 171 checked_exp_with_tolerance(&self, tolerance: Decimal) -> Option<Decimal>172 fn checked_exp_with_tolerance(&self, tolerance: Decimal) -> Option<Decimal> { 173 if self.is_zero() { 174 return Some(Decimal::ONE); 175 } 176 if self.is_sign_negative() { 177 let mut flipped = *self; 178 flipped.set_sign_positive(true); 179 let exp = flipped.checked_exp_with_tolerance(tolerance)?; 180 return Decimal::ONE.checked_div(exp); 181 } 182 183 let mut term = *self; 184 let mut result = self.checked_add(Decimal::ONE)?; 185 186 for factorial in FACTORIAL.iter().skip(2) { 187 term = self.checked_mul(term)?; 188 let next = result + (term / factorial); 189 let diff = (next - result).abs(); 190 result = next; 191 if diff <= tolerance { 192 break; 193 } 194 } 195 196 Some(result) 197 } 198 powi(&self, exp: i64) -> Decimal199 fn powi(&self, exp: i64) -> Decimal { 200 match self.checked_powi(exp) { 201 Some(result) => result, 202 None => panic!("Pow overflowed"), 203 } 204 } 205 checked_powi(&self, exp: i64) -> Option<Decimal>206 fn checked_powi(&self, exp: i64) -> Option<Decimal> { 207 // For negative exponents we change x^-y into 1 / x^y. 208 // Otherwise, we calculate a standard unsigned exponent 209 if exp >= 0 { 210 return self.checked_powu(exp as u64); 211 } 212 213 // Get the unsigned exponent 214 let exp = exp.unsigned_abs(); 215 let pow = match self.checked_powu(exp) { 216 Some(v) => v, 217 None => return None, 218 }; 219 Decimal::ONE.checked_div(pow) 220 } 221 powu(&self, exp: u64) -> Decimal222 fn powu(&self, exp: u64) -> Decimal { 223 match self.checked_powu(exp) { 224 Some(result) => result, 225 None => panic!("Pow overflowed"), 226 } 227 } 228 checked_powu(&self, exp: u64) -> Option<Decimal>229 fn checked_powu(&self, exp: u64) -> Option<Decimal> { 230 match exp { 231 0 => Some(Decimal::ONE), 232 1 => Some(*self), 233 2 => self.checked_mul(*self), 234 _ => { 235 // Get the squared value 236 let squared = match self.checked_mul(*self) { 237 Some(s) => s, 238 None => return None, 239 }; 240 // Square self once and make an infinite sized iterator of the square. 241 let iter = core::iter::repeat(squared); 242 243 // We then take half of the exponent to create a finite iterator and then multiply those together. 244 let mut product = Decimal::ONE; 245 for x in iter.take((exp >> 1) as usize) { 246 match product.checked_mul(x) { 247 Some(r) => product = r, 248 None => return None, 249 }; 250 } 251 252 // If the exponent is odd we still need to multiply once more 253 if exp & 0x1 > 0 { 254 match self.checked_mul(product) { 255 Some(p) => product = p, 256 None => return None, 257 } 258 } 259 product.normalize_assign(); 260 Some(product) 261 } 262 } 263 } 264 powf(&self, exp: f64) -> Decimal265 fn powf(&self, exp: f64) -> Decimal { 266 match self.checked_powf(exp) { 267 Some(result) => result, 268 None => panic!("Pow overflowed"), 269 } 270 } 271 checked_powf(&self, exp: f64) -> Option<Decimal>272 fn checked_powf(&self, exp: f64) -> Option<Decimal> { 273 let exp = match Decimal::from_f64(exp) { 274 Some(f) => f, 275 None => return None, 276 }; 277 self.checked_powd(exp) 278 } 279 powd(&self, exp: Decimal) -> Decimal280 fn powd(&self, exp: Decimal) -> Decimal { 281 match self.checked_powd(exp) { 282 Some(result) => result, 283 None => panic!("Pow overflowed"), 284 } 285 } 286 checked_powd(&self, exp: Decimal) -> Option<Decimal>287 fn checked_powd(&self, exp: Decimal) -> Option<Decimal> { 288 if exp.is_zero() { 289 return Some(Decimal::ONE); 290 } 291 if self.is_zero() { 292 return Some(Decimal::ZERO); 293 } 294 if self.is_one() { 295 return Some(Decimal::ONE); 296 } 297 if exp.is_one() { 298 return Some(*self); 299 } 300 301 // If the scale is 0 then it's a trivial calculation 302 let exp = exp.normalize(); 303 if exp.scale() == 0 { 304 if exp.mid() != 0 || exp.hi() != 0 { 305 // Exponent way too big 306 return None; 307 } 308 309 return if exp.is_sign_negative() { 310 self.checked_powi(-(exp.lo() as i64)) 311 } else { 312 self.checked_powu(exp.lo() as u64) 313 }; 314 } 315 316 // We do some approximations since we've got a decimal exponent. 317 // For positive bases: a^b = exp(b*ln(a)) 318 let negative = self.is_sign_negative(); 319 let e = match self.abs().ln().checked_mul(exp) { 320 Some(e) => e, 321 None => return None, 322 }; 323 let mut result = e.checked_exp()?; 324 result.set_sign_negative(negative); 325 Some(result) 326 } 327 sqrt(&self) -> Option<Decimal>328 fn sqrt(&self) -> Option<Decimal> { 329 if self.is_sign_negative() { 330 return None; 331 } 332 333 if self.is_zero() { 334 return Some(Decimal::ZERO); 335 } 336 337 // Start with an arbitrary number as the first guess 338 let mut result = self / Decimal::TWO; 339 // Too small to represent, so we start with self 340 // Future iterations could actually avoid using a decimal altogether and use a buffered 341 // vector, only combining back into a decimal on return 342 if result.is_zero() { 343 result = *self; 344 } 345 let mut last = result + Decimal::ONE; 346 347 // Keep going while the difference is larger than the tolerance 348 let mut circuit_breaker = 0; 349 while last != result { 350 circuit_breaker += 1; 351 assert!(circuit_breaker < 1000, "geo mean circuit breaker"); 352 353 last = result; 354 result = (result + self / result) / Decimal::TWO; 355 } 356 357 Some(result) 358 } 359 360 #[cfg(feature = "maths-nopanic")] ln(&self) -> Decimal361 fn ln(&self) -> Decimal { 362 match self.checked_ln() { 363 Some(result) => result, 364 None => Decimal::ZERO, 365 } 366 } 367 368 #[cfg(not(feature = "maths-nopanic"))] ln(&self) -> Decimal369 fn ln(&self) -> Decimal { 370 match self.checked_ln() { 371 Some(result) => result, 372 None => { 373 if self.is_sign_negative() { 374 panic!("Unable to calculate ln for negative numbers") 375 } else if self.is_zero() { 376 panic!("Unable to calculate ln for zero") 377 } else { 378 panic!("Calculation of ln failed for unknown reasons") 379 } 380 } 381 } 382 } 383 checked_ln(&self) -> Option<Decimal>384 fn checked_ln(&self) -> Option<Decimal> { 385 if self.is_sign_negative() || self.is_zero() { 386 return None; 387 } 388 if self.is_one() { 389 return Some(Decimal::ZERO); 390 } 391 392 // Approximate using Taylor Series 393 let mut x = *self; 394 let mut count = 0; 395 while x >= Decimal::ONE { 396 x *= Decimal::E_INVERSE; 397 count += 1; 398 } 399 while x <= Decimal::E_INVERSE { 400 x *= Decimal::E; 401 count -= 1; 402 } 403 x -= Decimal::ONE; 404 if x.is_zero() { 405 return Some(Decimal::new(count, 0)); 406 } 407 let mut result = Decimal::ZERO; 408 let mut iteration = 0; 409 let mut y = Decimal::ONE; 410 let mut last = Decimal::ONE; 411 while last != result && iteration < 100 { 412 iteration += 1; 413 last = result; 414 y *= -x; 415 result += y / Decimal::new(iteration, 0); 416 } 417 Some(Decimal::new(count, 0) - result) 418 } 419 420 #[cfg(feature = "maths-nopanic")] log10(&self) -> Decimal421 fn log10(&self) -> Decimal { 422 match self.checked_log10() { 423 Some(result) => result, 424 None => Decimal::ZERO, 425 } 426 } 427 428 #[cfg(not(feature = "maths-nopanic"))] log10(&self) -> Decimal429 fn log10(&self) -> Decimal { 430 match self.checked_log10() { 431 Some(result) => result, 432 None => { 433 if self.is_sign_negative() { 434 panic!("Unable to calculate log10 for negative numbers") 435 } else if self.is_zero() { 436 panic!("Unable to calculate log10 for zero") 437 } else { 438 panic!("Calculation of log10 failed for unknown reasons") 439 } 440 } 441 } 442 } 443 checked_log10(&self) -> Option<Decimal>444 fn checked_log10(&self) -> Option<Decimal> { 445 use crate::ops::array::{div_by_u32, is_all_zero}; 446 // Early exits 447 if self.is_sign_negative() || self.is_zero() { 448 return None; 449 } 450 if self.is_one() { 451 return Some(Decimal::ZERO); 452 } 453 // This uses a very basic method for calculating log10. We know the following is true: 454 // log10(n) = ln(n) / ln(10) 455 // From this we can perform some small optimizations: 456 // 1. ln(10) is a constant 457 // 2. Multiplication is faster than division, so we can pre-calculate the constant 1/ln(10) 458 // This allows us to then simplify log10(n) to: 459 // log10(n) = C * ln(n) 460 // Before doing all of this however, we see if there are simple calculations to be made. 461 let mut working = self.mantissa_array3(); 462 let mut result = 0; 463 let mut invalid_early_exit = false; 464 while !is_all_zero(&working) { 465 let remainder = div_by_u32(&mut working, 10u32); 466 if remainder != 0 { 467 invalid_early_exit = true; 468 break; 469 } 470 result += 1; 471 if working[2] == 0 && working[1] == 0 && working[0] == 1 { 472 break; 473 } 474 } 475 if !invalid_early_exit { 476 return Some((result - self.scale()).into()); 477 } 478 479 self.checked_ln().map(|result| LN10_INVERSE * result) 480 } 481 erf(&self) -> Decimal482 fn erf(&self) -> Decimal { 483 if self.is_sign_positive() { 484 let one = &Decimal::ONE; 485 486 let xa1 = self * Decimal::from_parts(705230784, 0, 0, false, 10); 487 let xa2 = self.powi(2) * Decimal::from_parts(422820123, 0, 0, false, 10); 488 let xa3 = self.powi(3) * Decimal::from_parts(92705272, 0, 0, false, 10); 489 let xa4 = self.powi(4) * Decimal::from_parts(1520143, 0, 0, false, 10); 490 let xa5 = self.powi(5) * Decimal::from_parts(2765672, 0, 0, false, 10); 491 let xa6 = self.powi(6) * Decimal::from_parts(430638, 0, 0, false, 10); 492 493 let sum = one + xa1 + xa2 + xa3 + xa4 + xa5 + xa6; 494 one - (one / sum.powi(16)) 495 } else { 496 -self.abs().erf() 497 } 498 } 499 norm_cdf(&self) -> Decimal500 fn norm_cdf(&self) -> Decimal { 501 (Decimal::ONE + (self / Decimal::from_parts(2318911239, 3292722, 0, false, 16)).erf()) / Decimal::TWO 502 } 503 norm_pdf(&self) -> Decimal504 fn norm_pdf(&self) -> Decimal { 505 match self.checked_norm_pdf() { 506 Some(d) => d, 507 None => panic!("Norm Pdf overflowed"), 508 } 509 } 510 checked_norm_pdf(&self) -> Option<Decimal>511 fn checked_norm_pdf(&self) -> Option<Decimal> { 512 let sqrt2pi = Decimal::from_parts_raw(2133383024, 2079885984, 1358845910, 1835008); 513 let factor = -self.checked_powi(2)?; 514 let factor = factor.checked_div(Decimal::TWO)?; 515 factor.checked_exp()?.checked_div(sqrt2pi) 516 } 517 sin(&self) -> Decimal518 fn sin(&self) -> Decimal { 519 match self.checked_sin() { 520 Some(x) => x, 521 None => panic!("Sin overflowed"), 522 } 523 } 524 checked_sin(&self) -> Option<Decimal>525 fn checked_sin(&self) -> Option<Decimal> { 526 if self.is_zero() { 527 return Some(Decimal::ZERO); 528 } 529 if self.is_sign_negative() { 530 // -Sin(-x) 531 return (-self).checked_sin().map(|x| -x); 532 } 533 if self >= &Decimal::TWO_PI { 534 // Reduce large numbers early - we can do this using rem to constrain to a range 535 let adjusted = self.checked_rem(Decimal::TWO_PI)?; 536 return adjusted.checked_sin(); 537 } 538 if self >= &Decimal::PI { 539 // -Sin(x-π) 540 return (self - Decimal::PI).checked_sin().map(|x| -x); 541 } 542 if self >= &Decimal::QUARTER_PI { 543 // Cos(π2-x) 544 return (Decimal::HALF_PI - self).checked_cos(); 545 } 546 547 // Taylor series: 548 // ∑(n=0 to ∞) : ((−1)^n / (2n + 1)!) * x^(2n + 1) , x∈R 549 // First few expansions: 550 // x^1/1! - x^3/3! + x^5/5! - x^7/7! + x^9/9! 551 let mut result = Decimal::ZERO; 552 for n in 0..TRIG_SERIES_UPPER_BOUND { 553 let x = 2 * n + 1; 554 let element = self.checked_powi(x as i64)?.checked_div(FACTORIAL[x])?; 555 if n & 0x1 == 0 { 556 result += element; 557 } else { 558 result -= element; 559 } 560 } 561 Some(result) 562 } 563 cos(&self) -> Decimal564 fn cos(&self) -> Decimal { 565 match self.checked_cos() { 566 Some(x) => x, 567 None => panic!("Cos overflowed"), 568 } 569 } 570 checked_cos(&self) -> Option<Decimal>571 fn checked_cos(&self) -> Option<Decimal> { 572 if self.is_zero() { 573 return Some(Decimal::ONE); 574 } 575 if self.is_sign_negative() { 576 // Cos(-x) 577 return (-self).checked_cos(); 578 } 579 if self >= &Decimal::TWO_PI { 580 // Reduce large numbers early - we can do this using rem to constrain to a range 581 let adjusted = self.checked_rem(Decimal::TWO_PI)?; 582 return adjusted.checked_cos(); 583 } 584 if self >= &Decimal::PI { 585 // -Cos(x-π) 586 return (self - Decimal::PI).checked_cos().map(|x| -x); 587 } 588 if self >= &Decimal::QUARTER_PI { 589 // Sin(π2-x) 590 return (Decimal::HALF_PI - self).checked_sin(); 591 } 592 593 // Taylor series: 594 // ∑(n=0 to ∞) : ((−1)^n / (2n)!) * x^(2n) , x∈R 595 // First few expansions: 596 // x^0/0! - x^2/2! + x^4/4! - x^6/6! + x^8/8! 597 let mut result = Decimal::ZERO; 598 for n in 0..TRIG_SERIES_UPPER_BOUND { 599 let x = 2 * n; 600 let element = self.checked_powi(x as i64)?.checked_div(FACTORIAL[x])?; 601 if n & 0x1 == 0 { 602 result += element; 603 } else { 604 result -= element; 605 } 606 } 607 Some(result) 608 } 609 tan(&self) -> Decimal610 fn tan(&self) -> Decimal { 611 match self.checked_tan() { 612 Some(x) => x, 613 None => panic!("Tan overflowed"), 614 } 615 } 616 checked_tan(&self) -> Option<Decimal>617 fn checked_tan(&self) -> Option<Decimal> { 618 if self.is_zero() { 619 return Some(Decimal::ZERO); 620 } 621 if self.is_sign_negative() { 622 // -Tan(-x) 623 return (-self).checked_tan().map(|x| -x); 624 } 625 if self >= &Decimal::TWO_PI { 626 // Reduce large numbers early - we can do this using rem to constrain to a range 627 let adjusted = self.checked_rem(Decimal::TWO_PI)?; 628 return adjusted.checked_tan(); 629 } 630 // Reduce to 0 <= x <= PI 631 if self >= &Decimal::PI { 632 // Tan(x-π) 633 return (self - Decimal::PI).checked_tan(); 634 } 635 // Reduce to 0 <= x <= PI/2 636 if self > &Decimal::HALF_PI { 637 // We can use the symmetrical function inside the first quadrant 638 // e.g. tan(x) = -tan((PI/2 - x) + PI/2) 639 return ((Decimal::HALF_PI - self) + Decimal::HALF_PI).checked_tan().map(|x| -x); 640 } 641 642 // It has now been reduced to 0 <= x <= PI/2. If it is >= PI/4 we can make it even smaller 643 // by calculating tan(PI/2 - x) and taking the reciprocal 644 if self > &Decimal::QUARTER_PI { 645 return match (Decimal::HALF_PI - self).checked_tan() { 646 Some(x) => Decimal::ONE.checked_div(x), 647 None => None, 648 }; 649 } 650 651 // Due the way that tan(x) sharply tends towards infinity, we try to optimize 652 // the resulting accuracy by using Trigonometric identity when > PI/8. We do this by 653 // replacing the angle with one that is half as big. 654 if self > &EIGHTH_PI { 655 // Work out tan(x/2) 656 let tan_half = (self / Decimal::TWO).checked_tan()?; 657 // Work out the dividend i.e. 2tan(x/2) 658 let dividend = Decimal::TWO.checked_mul(tan_half)?; 659 660 // Work out the divisor i.e. 1 - tan^2(x/2) 661 let squared = tan_half.checked_mul(tan_half)?; 662 let divisor = Decimal::ONE - squared; 663 // Treat this as infinity 664 if divisor.is_zero() { 665 return None; 666 } 667 return dividend.checked_div(divisor); 668 } 669 670 // Do a polynomial approximation based upon the Maclaurin series. 671 // This can be simplified to something like: 672 // 673 // ∑(n=1,3,5,7,9)(f(n)(0)/n!)x^n 674 // 675 // First few expansions (which we leverage): 676 // (f'(0)/1!)x^1 + (f'''(0)/3!)x^3 + (f'''''(0)/5!)x^5 + (f'''''''/7!)x^7 677 // 678 // x + (1/3)x^3 + (2/15)x^5 + (17/315)x^7 + (62/2835)x^9 + (1382/155925)x^11 679 // 680 // (Generated by https://www.wolframalpha.com/widgets/view.jsp?id=fe1ad8d4f5dbb3cb866d0c89beb527a6) 681 // The more terms, the better the accuracy. This generates accuracy within approx 10^-8 for angles 682 // less than PI/8. 683 const SERIES: [(Decimal, u64); 6] = [ 684 // 1 / 3 685 (Decimal::from_parts_raw(89478485, 347537611, 180700362, 1835008), 3), 686 // 2 / 15 687 (Decimal::from_parts_raw(894784853, 3574988881, 72280144, 1835008), 5), 688 // 17 / 315 689 (Decimal::from_parts_raw(905437054, 3907911371, 2925624, 1769472), 7), 690 // 62 / 2835 691 (Decimal::from_parts_raw(3191872741, 2108928381, 11855473, 1835008), 9), 692 // 1382 / 155925 693 (Decimal::from_parts_raw(3482645539, 2612995122, 4804769, 1835008), 11), 694 // 21844 / 6081075 695 (Decimal::from_parts_raw(4189029078, 2192791200, 1947296, 1835008), 13), 696 ]; 697 let mut result = *self; 698 for (fraction, pow) in SERIES { 699 result += fraction * self.powu(pow); 700 } 701 Some(result) 702 } 703 } 704 705 impl Pow<Decimal> for Decimal { 706 type Output = Decimal; 707 pow(self, rhs: Decimal) -> Self::Output708 fn pow(self, rhs: Decimal) -> Self::Output { 709 MathematicalOps::powd(&self, rhs) 710 } 711 } 712 713 impl Pow<u64> for Decimal { 714 type Output = Decimal; 715 pow(self, rhs: u64) -> Self::Output716 fn pow(self, rhs: u64) -> Self::Output { 717 MathematicalOps::powu(&self, rhs) 718 } 719 } 720 721 impl Pow<i64> for Decimal { 722 type Output = Decimal; 723 pow(self, rhs: i64) -> Self::Output724 fn pow(self, rhs: i64) -> Self::Output { 725 MathematicalOps::powi(&self, rhs) 726 } 727 } 728 729 impl Pow<f64> for Decimal { 730 type Output = Decimal; 731 pow(self, rhs: f64) -> Self::Output732 fn pow(self, rhs: f64) -> Self::Output { 733 MathematicalOps::powf(&self, rhs) 734 } 735 } 736 737 #[cfg(test)] 738 mod test { 739 use super::*; 740 741 #[test] test_factorials()742 fn test_factorials() { 743 assert_eq!("1", FACTORIAL[0].to_string(), "0!"); 744 assert_eq!("1", FACTORIAL[1].to_string(), "1!"); 745 assert_eq!("2", FACTORIAL[2].to_string(), "2!"); 746 assert_eq!("6", FACTORIAL[3].to_string(), "3!"); 747 assert_eq!("24", FACTORIAL[4].to_string(), "4!"); 748 assert_eq!("120", FACTORIAL[5].to_string(), "5!"); 749 assert_eq!("720", FACTORIAL[6].to_string(), "6!"); 750 assert_eq!("5040", FACTORIAL[7].to_string(), "7!"); 751 assert_eq!("40320", FACTORIAL[8].to_string(), "8!"); 752 assert_eq!("362880", FACTORIAL[9].to_string(), "9!"); 753 assert_eq!("3628800", FACTORIAL[10].to_string(), "10!"); 754 assert_eq!("39916800", FACTORIAL[11].to_string(), "11!"); 755 assert_eq!("479001600", FACTORIAL[12].to_string(), "12!"); 756 assert_eq!("6227020800", FACTORIAL[13].to_string(), "13!"); 757 assert_eq!("87178291200", FACTORIAL[14].to_string(), "14!"); 758 assert_eq!("1307674368000", FACTORIAL[15].to_string(), "15!"); 759 assert_eq!("20922789888000", FACTORIAL[16].to_string(), "16!"); 760 assert_eq!("355687428096000", FACTORIAL[17].to_string(), "17!"); 761 assert_eq!("6402373705728000", FACTORIAL[18].to_string(), "18!"); 762 assert_eq!("121645100408832000", FACTORIAL[19].to_string(), "19!"); 763 assert_eq!("2432902008176640000", FACTORIAL[20].to_string(), "20!"); 764 assert_eq!("51090942171709440000", FACTORIAL[21].to_string(), "21!"); 765 assert_eq!("1124000727777607680000", FACTORIAL[22].to_string(), "22!"); 766 assert_eq!("25852016738884976640000", FACTORIAL[23].to_string(), "23!"); 767 assert_eq!("620448401733239439360000", FACTORIAL[24].to_string(), "24!"); 768 assert_eq!("15511210043330985984000000", FACTORIAL[25].to_string(), "25!"); 769 assert_eq!("403291461126605635584000000", FACTORIAL[26].to_string(), "26!"); 770 assert_eq!("10888869450418352160768000000", FACTORIAL[27].to_string(), "27!"); 771 } 772 } 773