From 408ef733915f750104a07a4e1e103fb8e4e9a843 Mon Sep 17 00:00:00 2001 From: Sami Farin Date: Thu, 1 Jan 2026 18:30:50 +0200 Subject: [PATCH] pow optimization for big values Configurable values for the algorithm change can be modified here: bigrat.rs // Use exp/ln approximation if: // 1. The denominator is large (e.g. ^0.0001), indicating a complex decimal root. // 2. The numerator is large AND it's a root (e.g. ^210.1 = ^2101/10). let large_root = rhs.den > BigUint::from(1000u64); let large_power = rhs.den > BigUint::from(1u64) && rhs.num > BigUint::from(100u64); Also Dynamic precision for {sin,cos,atan}_series, inv trigs accurate up to 10e369 --- core/Cargo.toml | 3 + core/src/num/bigrat.rs | 958 ++++++++++++++++++++++++++++---- core/tests/integration_tests.rs | 32 +- 3 files changed, 883 insertions(+), 110 deletions(-) diff --git a/core/Cargo.toml b/core/Cargo.toml index 74c67dfd..b5b36743 100644 --- a/core/Cargo.toml +++ b/core/Cargo.toml @@ -9,3 +9,6 @@ keywords.workspace = true categories.workspace = true license.workspace = true readme = "README.md" + +[dependencies] +num-traits = "0.2.19" diff --git a/core/src/num/bigrat.rs b/core/src/num/bigrat.rs index 8949d36a..5d98cadc 100644 --- a/core/src/num/bigrat.rs +++ b/core/src/num/bigrat.rs @@ -7,6 +7,8 @@ use crate::num::{Base, Exact, FormattingStyle, Range, RangeBound}; use crate::result::FResult; use crate::serialize::CborValue; use core::f64; +use num_traits::ToPrimitive; +use std::sync::OnceLock; use std::{cmp, fmt, hash, ops}; pub(crate) mod sign { @@ -120,6 +122,9 @@ impl hash::Hash for BigRat { } impl BigRat { + const INTERNAL_PRECISION: usize = 200; + const FIXED_POINT_GUARD_DIGITS: usize = 20; + pub(crate) fn serialize(&self) -> CborValue { let num = self.num.serialize(self.sign); if self.den == 1.into() { @@ -219,94 +224,294 @@ impl BigRat { } } - #[allow( - clippy::cast_possible_truncation, - clippy::cast_sign_loss, - clippy::cast_precision_loss - )] - pub(crate) fn from_f64(mut f: f64, int: &I) -> FResult { - let negative = f < 0.0; - if negative { - f = -f; + fn high_precision_ln_2() -> Self { + static LN_2: OnceLock = OnceLock::new(); + LN_2.get_or_init(|| { + // ln(2) to ~280 digits + let s = "0.69314718055994530941723212145817656807550013436025525412068000949339362196969471560586332699641868754200148102057068573368552023575813055703267075163507596193072757082837143519030703862389167347112335011536449795523912047517268157493206515552473413952588295045300709532636664265410423915781495204374"; + let split_idx = std::cmp::min(2 + Self::INTERNAL_PRECISION, s.len()); + let slice = &s[..split_idx]; + let int = &crate::interrupt::Never; + let fraction_str = &slice[2..]; + let mut fraction = BigUint::from(0u64); + let ten = BigUint::from(10u64); + for c in fraction_str.chars() { + let digit = u64::from(c.to_digit(10).unwrap()); + fraction = fraction.mul(&ten, int).unwrap().add(&BigUint::from(digit)); + } + let den = BigUint::pow(&ten, &BigUint::from(fraction_str.len() as u64), int).unwrap(); + Self { + sign: Sign::Positive, + num: fraction, + den, + } + }).clone() + } + + fn high_precision_ln_10() -> Self { + static LN_10: OnceLock = OnceLock::new(); + LN_10.get_or_init(|| { + // ln(10) to ~280 digits + let s = "2.30258509299404568401799145468436420760110148862877297603332790096757260967735248023599720508959829834196778404228624863340952546508280675666628736909878168948290720832555468084379989482623319852839350530896537773262884616336622228769821988674654366747440424327436515504893431493939147961940440022210"; + let split_idx = std::cmp::min(2 + Self::INTERNAL_PRECISION, s.len()); + let slice = &s[..split_idx]; + let int = &crate::interrupt::Never; + let integer = BigUint::from(2u64); + let fraction_str = &slice[2..]; + let mut fraction = BigUint::from(0u64); + let ten = BigUint::from(10u64); + for c in fraction_str.chars() { + let digit = u64::from(c.to_digit(10).unwrap()); + fraction = fraction.mul(&ten, int).unwrap().add(&BigUint::from(digit)); + } + let den = BigUint::pow(&ten, &BigUint::from(fraction_str.len() as u64), int).unwrap(); + Self { + sign: Sign::Positive, + num: integer.mul(&den, int).unwrap().add(&fraction), + den, + } + }).clone() + } + + // exp(x) = 1 + x + x^2/2! + ... + fn exp_series(x: &Self, int: &I) -> FResult { + if x.is_definitely_zero() { + return Ok(Self::from(1)); + } + + let scale = Self::fixed_point_scale(); + // We expect x to be small (reduced), so fixed point is fine + let x_fixed = x.num.clone().mul(&scale.clone(), int)?.div(&x.den, int)?; + + if x_fixed == 0.into() { + return Ok(Self::from(1)); } - let i = (f * u64::MAX as f64) as u128; - let part1 = i as u64; - let part2 = (i >> 64) as u64; + + let mut result = scale.clone(); // 1.0 + let mut term = scale.clone(); // 1.0 + + for i in 1..Self::INTERNAL_PRECISION { + term = term.mul(&x_fixed, int)?; + term = term.div(scale, int)?; + term = term.div(&BigUint::from(i as u64), int)?; + + if term == 0.into() { + break; + } + + result = result.add(&term); + } + Ok(Self { - sign: if negative { - Sign::Negative - } else { - Sign::Positive - }, - num: BigUint::from(part1) - .add(&BigUint::from(part2).mul(&BigUint::from(u64::MAX), int)?), - den: BigUint::from(u64::MAX), + sign: Sign::Positive, + num: result, + den: scale.clone(), }) } - // sin works for all real numbers - pub(crate) fn sin(self, int: &I) -> FResult> { - Ok(if self == 0.into() { - Exact::new(Self::from(0), true) - } else { - Exact::new(Self::from_f64(f64::sin(self.into_f64(int)?), int)?, false) - }) - } + // ln((1+z)/(1-z)) = 2(z + z^3/3 + ...) + fn ln_series_inv_tanh(z: Self, int: &I) -> FResult { + if z.is_definitely_zero() { + return Ok(z); + } - // cos works for all real numbers - pub(crate) fn cos(self, int: &I) -> FResult> { - Ok(if self == 0.into() { - Exact::new(Self::from(1), true) - } else { - Exact::new(Self::from_f64(f64::cos(self.into_f64(int)?), int)?, false) + let scale = Self::fixed_point_scale(); + let z_fixed = z.num.clone().mul(&scale.clone(), int)?.div(&z.den, int)?; + + if z_fixed == 0.into() { + return Ok(z); + } + + let mut result = z_fixed.clone(); + let mut term = z_fixed.clone(); + let z_sq = z_fixed.clone().mul(&z_fixed, int)?; + let scale_sq = scale.clone().mul(scale, int)?; + + for i in 1..Self::INTERNAL_PRECISION { + term = term.mul(&z_sq, int)?; + term = term.div(&scale_sq, int)?; + + if term == 0.into() { + break; + } + + let n = 2 * i + 1; + let val = term.clone().div(&BigUint::from(n as u64), int)?; + + result = result.add(&val); + } + + // Multiply by 2 + let result = result.mul(&BigUint::from(2u64), int)?; + + Ok(Self { + sign: Sign::Positive, + num: result, + den: scale.clone(), }) } - // asin, acos and atan only work for values between -1 and 1 - pub(crate) fn asin(self, int: &I) -> FResult { + pub(crate) fn exp(self, int: &I) -> FResult> { + if self.is_definitely_zero() { + return Ok(Exact::new(Self::from(1), true)); + } let one = Self::from(1); - if self > one || self < -one { - return Err(out_of_range(self.fm(int)?, Range::open(-1, 1))); + + let mut x = self.clone(); + let mut invert = false; + if x.sign == Sign::Negative { + x = -x; + invert = true; + } + + // Range reduction: Reduce x until x < 0.5 + let mut k = 0; + // 1/2 + let half = Self { + sign: Sign::Positive, + num: 1.into(), + den: 2.into(), + }; + + while x > half { + x.den = x.den.mul(&2.into(), int)?; + k += 1; + } + + let mut result = Self::exp_series(&x, int)?; + let scale = Self::fixed_point_scale(); + + // Square result k times + for _ in 0..k { + // Perform fixed-point squaring (N^2 / scale) instead of rational squaring. + // This keeps the integer sizes constant (~200 digits) instead of exploding to 400k digits. + let num_sq = result.num.clone().mul(&result.num, int)?; + result.num = num_sq.div(scale, int)?; + // result.den is already 'scale', so we don't need to change it. + } + + if invert { + result = one.div(&result, int)?; } - Self::from_f64(f64::asin(self.into_f64(int)?), int) + + Ok(Exact::new(result, false)) } - pub(crate) fn acos(self, int: &I) -> FResult { + pub(crate) fn ln(self, int: &I) -> FResult> { let one = Self::from(1); - if self > one || self < -one { - return Err(out_of_range(self.fm(int)?, Range::open(-1, 1))); + if self <= 0.into() { + return Err(out_of_range( + self.fm(int)?, + Range { + start: RangeBound::Open(0), + end: RangeBound::None, + }, + )); + } + if self == one { + return Ok(Exact::new(0.into(), true)); } - Self::from_f64(f64::acos(self.into_f64(int)?), int) + + // Estimate k = floor(log2(x)) to avoid overflow with large inputs (e.g. 2^1024) + let num_log2 = self.num.log2(int)?; + let den_log2 = self.den.log2(int)?; + let k = (num_log2 - den_log2).floor().to_i64().unwrap(); + + let pow_2_k = if k >= 0 { + Self::from(1).mul( + &Self::from(2).pow(Self::from(k.unsigned_abs()), int)?.value, + int, + )? + } else { + Self::from(1).div( + &Self::from(2).pow(Self::from(k.unsigned_abs()), int)?.value, + int, + )? + }; + + let y = self.div(&pow_2_k, int)?; + + let num = y.clone().sub(one.clone(), int)?; + let den = y.clone().add(one.clone(), int)?; + let z = num.div(&den, int)?; + + let ln_y = Self::ln_series_inv_tanh(z, int)?; + let k_ln2 = Self::high_precision_ln_2().mul(&Self::from(k.unsigned_abs()), int)?; + + let result = if k >= 0 { + k_ln2.add(ln_y, int)? + } else { + ln_y.sub(k_ln2, int)? + }; + + Ok(Exact::new(result, false)) } - // note that this works for any real number, unlike asin and acos - pub(crate) fn atan(self, int: &I) -> FResult { - Self::from_f64(f64::atan(self.into_f64(int)?), int) + pub(crate) fn log2(self, int: &I) -> FResult { + let ln_val = self.ln(int)?; + ln_val.value.div(&Self::high_precision_ln_2(), int) } - pub(crate) fn atan2(self, rhs: Self, int: &I) -> FResult { - Self::from_f64(f64::atan2(self.into_f64(int)?, rhs.into_f64(int)?), int) + pub(crate) fn log10(self, int: &I) -> FResult { + let ln_val = self.ln(int)?; + ln_val.value.div(&Self::high_precision_ln_10(), int) } pub(crate) fn sinh(self, int: &I) -> FResult { - Self::from_f64(f64::sinh(self.into_f64(int)?), int) + let one = Self::from(1); + let two = Self::from(2); + + // e^x + let exp_x = self.clone().exp(int)?; + // e^-x = 1 / e^x + let exp_neg_x = one.div(&exp_x.value, int)?; + + // (e^x - e^-x) / 2 + let num = exp_x.value.sub(exp_neg_x, int)?; + num.div(&two, int) } pub(crate) fn cosh(self, int: &I) -> FResult { - Self::from_f64(f64::cosh(self.into_f64(int)?), int) + let one = Self::from(1); + let two = Self::from(2); + + let exp_x = self.clone().exp(int)?; + let exp_neg_x = one.div(&exp_x.value, int)?; + + let num = exp_x.value.add(exp_neg_x, int)?; + num.div(&two, int) } pub(crate) fn tanh(self, int: &I) -> FResult { - Self::from_f64(f64::tanh(self.into_f64(int)?), int) + let one = Self::from(1); + let two = Self::from(2); + + // e^2x + let two_x = self.mul(&two, int)?; + let exp_2x = two_x.exp(int)?; + + let num = exp_2x.value.clone().sub(one.clone(), int)?; + let den = exp_2x.value.add(one, int)?; + + num.div(&den, int) } pub(crate) fn asinh(self, int: &I) -> FResult { - Self::from_f64(f64::asinh(self.into_f64(int)?), int) + let one = Self::from(1); + let two = Self::from(2); + + let x_sq = self.clone().mul(&self, int)?; + let x_sq_plus_1 = x_sq.add(one, int)?; + let sqrt = x_sq_plus_1.root_n(&two, int)?; + + let arg = self.add(sqrt.value, int)?; + Ok(arg.ln(int)?.value) } - // value must not be less than 1 pub(crate) fn acosh(self, int: &I) -> FResult { - if self < 1.into() { + let one = Self::from(1); + + if self < one { return Err(out_of_range( self.fm(int)?, Range { @@ -315,45 +520,498 @@ impl BigRat { }, )); } - Self::from_f64(f64::acosh(self.into_f64(int)?), int) + + let two = Self::from(2); + + let x_sq = self.clone().mul(&self, int)?; + let x_sq_minus_1 = x_sq.sub(one, int)?; + let sqrt = x_sq_minus_1.root_n(&two, int)?; + + let arg = self.add(sqrt.value, int)?; + Ok(arg.ln(int)?.value) } - // value must be between -1 and 1. pub(crate) fn atanh(self, int: &I) -> FResult { - let one: Self = 1.into(); - if self >= one || self <= -one { + let one = Self::from(1); + + if self >= one || self <= -one.clone() { return Err(out_of_range(self.fm(int)?, Range::open(-1, 1))); } - Self::from_f64(f64::atanh(self.into_f64(int)?), int) + + let two = Self::from(2); + + let num = one.clone().add(self.clone(), int)?; + let den = one.sub(self, int)?; + let arg = num.div(&den, int)?; + + let ln_val = arg.ln(int)?; + ln_val.value.div(&two, int) } - pub(crate) fn log2(self, int: &I) -> FResult { - if self <= 0.into() { - return Err(out_of_range( - self.fm(int)?, - Range { - start: RangeBound::Open(0), - end: RangeBound::None, - }, - )); + fn high_precision_pi() -> Self { + static PI_CACHE: OnceLock = OnceLock::new(); + PI_CACHE.get_or_init(|| { + let int = &crate::interrupt::Never; + let pi_str = "3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482"; + + let split_idx = std::cmp::min(2 + Self::INTERNAL_PRECISION, pi_str.len()); + let slice = &pi_str[..split_idx]; + + let integer = BigUint::from(3u64); + let fraction_str = &slice[2..]; + + let mut fraction = BigUint::from(0u64); + let ten = BigUint::from(10u64); + + for c in fraction_str.chars() { + let digit = u64::from(c.to_digit(10).unwrap()); + fraction = fraction.mul(&ten, int).unwrap().add(&BigUint::from(digit)); + } + + let denominator = BigUint::pow( + &ten, + &BigUint::from(fraction_str.len() as u64), + int + ).unwrap(); + + Self { + sign: Sign::Positive, + num: integer.mul(&denominator, int).unwrap().add(&fraction), + den: denominator, + } + }).clone() + } + + fn high_precision_two_pi() -> Self { + static TWO_PI: OnceLock = OnceLock::new(); + TWO_PI + .get_or_init(|| { + let int = &crate::interrupt::Never; + Self::high_precision_pi().mul(&Self::from(2), int).unwrap() + }) + .clone() + } + + fn high_precision_pi_over_2() -> Self { + static PI_OVER_2: OnceLock = OnceLock::new(); + PI_OVER_2 + .get_or_init(|| { + let int = &crate::interrupt::Never; + Self::high_precision_pi().div(&Self::from(2), int).unwrap() + }) + .clone() + } + + fn high_precision_pi_over_4() -> Self { + static PI_OVER_4: OnceLock = OnceLock::new(); + PI_OVER_4 + .get_or_init(|| { + let int = &crate::interrupt::Never; + Self::high_precision_pi().div(&Self::from(4), int).unwrap() + }) + .clone() + } + + fn fixed_point_scale() -> &'static BigUint { + static SCALE_CACHE: OnceLock = OnceLock::new(); + SCALE_CACHE.get_or_init(|| { + let int = &crate::interrupt::Never; + let ten = BigUint::from(10u64); + let exp = (Self::INTERNAL_PRECISION + Self::FIXED_POINT_GUARD_DIGITS) as u64; + BigUint::pow(&ten, &BigUint::from(exp), int).unwrap() + }) + } + + fn rem_euclid_rat(self, m: Self, int: &I) -> FResult { + let div = self.clone().div(&m, int)?; + let floor = div.floor(int)?; + self.sub(m.mul(&floor, int)?, int) + } + + pub(crate) fn sub(self, rhs: Self, int: &I) -> FResult { + self.add(-rhs, int) + } + + pub(crate) fn abs(self, _int: &I) -> Self { + Self { + sign: Sign::Positive, + num: self.num, + den: self.den, } - Self::from_f64(self.num.log2(int)? - self.den.log2(int)?, int) } - pub(crate) fn ln(self, int: &I) -> FResult> { - if self == 1.into() { - return Ok(Exact::new(0.into(), true)); + fn sin_series(x: Self, int: &I) -> FResult { + if x.is_definitely_zero() { + return Ok(x); } - Ok(Exact::new( - self.log2(int)? - .div(&Self::from_f64(std::f64::consts::LOG2_E, int)?, int)?, - false, - )) + + // Dynamic Precision: + // If x is smaller than the default scale (10^-200), we must increase the scale. + let mut scale = Self::fixed_point_scale().clone(); + let x_mag = x.clone().abs(int); + + if x_mag < 1.into() { + // Calculate required scale: (1/|x|) * 10^20 + // Use associated function syntax for BigUint::pow + let guard_digits = BigUint::pow(&BigUint::from(10u64), &BigUint::from(20u64), int)?; + + let required_scale = Self::from(1) + .div(&x_mag, int)? + .mul(&Self::from(guard_digits), int)?; + + if let Ok(req_uint) = required_scale.ceil(int)?.try_as_biguint(int) + && req_uint > scale { + scale = req_uint; + } + } + + let x_fixed = x.num.clone().mul(&scale.clone(), int)?.div(&x.den, int)?; + + if x_fixed == 0.into() { + return Ok(x); + } + + let mut result = x_fixed.clone(); + let mut term = x_fixed.clone(); + let x_sq = x_fixed.clone().mul(&x_fixed, int)?; + let scale_sq = scale.clone().mul(&scale, int)?; + + for i in 1..Self::INTERNAL_PRECISION * 2 { + term = term.mul(&x_sq, int)?; + term = term.div(&scale_sq, int)?; + let k1 = 2 * i; + let k2 = 2 * i + 1; + let divisor = BigUint::from(k1 as u64).mul(&BigUint::from(k2 as u64), int)?; + term = term.div(&divisor, int)?; + + if term == 0.into() { + break; + } + + if i % 2 == 1 { + result = result.sub(&term); + } else { + result = result.add(&term); + } + } + + Ok(Self { + sign: Sign::Positive, + num: result, + den: scale, + }) } - pub(crate) fn log10(self, int: &I) -> FResult { - self.log2(int)? - .div(&Self::from_f64(std::f64::consts::LOG2_10, int)?, int) + fn cos_series(x: &Self, int: &I) -> FResult { + let mut scale = Self::fixed_point_scale().clone(); + let x_mag = x.clone().abs(int); + + if x_mag < 1.into() && !x_mag.is_definitely_zero() { + let guard_digits = BigUint::pow(&BigUint::from(10u64), &BigUint::from(20u64), int)?; + + let required_scale = Self::from(1) + .div(&x_mag, int)? + .mul(&Self::from(guard_digits), int)?; + + if let Ok(req_uint) = required_scale.ceil(int)?.try_as_biguint(int) + && req_uint > scale { + scale = req_uint; + } + } + + let x_fixed = x.num.clone().mul(&scale.clone(), int)?.div(&x.den, int)?; + + if x_fixed == 0.into() { + return Ok(Self::from(1)); + } + + let mut result = scale.clone(); + let mut term = scale.clone(); + let x_sq = x_fixed.clone().mul(&x_fixed, int)?; + let scale_sq = scale.clone().mul(&scale, int)?; + + for i in 1..Self::INTERNAL_PRECISION * 2 { + term = term.mul(&x_sq, int)?; + term = term.div(&scale_sq, int)?; + let k1 = 2 * i - 1; + let k2 = 2 * i; + let divisor = BigUint::from(k1 as u64).mul(&BigUint::from(k2 as u64), int)?; + term = term.div(&divisor, int)?; + + if term == 0.into() { + break; + } + + if i % 2 == 1 { + result = result.sub(&term); + } else { + result = result.add(&term); + } + } + + Ok(Self { + sign: Sign::Positive, + num: result, + den: scale, + }) + } + + fn atan_series(x: Self, int: &I) -> FResult { + if x.is_definitely_zero() { + return Ok(x); + } + + if x.sign == Sign::Negative { + let pos_x = x.clone().abs(int); + let res = Self::atan_series(pos_x, int)?; + return Ok(-res); + } + + let mut scale = Self::fixed_point_scale().clone(); + let x_mag = x.clone().abs(int); + + if x_mag < 1.into() { + let guard_digits = BigUint::pow(&BigUint::from(10u64), &BigUint::from(20u64), int)?; + + let required_scale = Self::from(1) + .div(&x_mag, int)? + .mul(&Self::from(guard_digits), int)?; + + if let Ok(req_uint) = required_scale.ceil(int)?.try_as_biguint(int) + && req_uint > scale { + scale = req_uint; + } + } + + let x_fixed = x.num.clone().mul(&scale.clone(), int)?.div(&x.den, int)?; + + if x_fixed == 0.into() { + return Ok(x); + } + + let mut result = x_fixed.clone(); + let mut term = x_fixed.clone(); + let x_sq = x_fixed.clone().mul(&x_fixed, int)?; + let scale_sq = scale.clone().mul(&scale, int)?; + + for i in 1..Self::INTERNAL_PRECISION * 2 { + term = term.mul(&x_sq, int)?; + term = term.div(&scale_sq, int)?; + + if term == 0.into() { + break; + } + + let n = 2 * i + 1; + let val = term.clone().div(&BigUint::from(n as u64), int)?; + + if i % 2 == 1 { + result = result.sub(&val); + } else { + result = result.add(&val); + } + } + + Ok(Self { + sign: Sign::Positive, + num: result, + den: scale, + }) + } + + // sin works for all real numbers + pub(crate) fn sin(self, int: &I) -> FResult> { + Ok(if self == 0.into() { + Exact::new(Self::from(0), true) + } else { + let pi = Self::high_precision_pi(); + let two_pi = Self::high_precision_two_pi(); + let pi_over_2 = Self::high_precision_pi_over_2(); + let pi_over_4 = Self::high_precision_pi_over_4(); + + let mut x = self.rem_euclid_rat(two_pi, int)?; + let mut sign = Sign::Positive; + + if x >= pi { + x = x.sub(pi.clone(), int)?; + sign = sign.flip(); + } + if x > pi_over_2 { + x = pi.sub(x, int)?; + } + + if x > pi_over_4 { + x = pi_over_2.sub(x, int)?; + let val = Self::cos_series(&x, int)?; + Exact::new(if sign == Sign::Negative { -val } else { val }, false) + } else { + let val = Self::sin_series(x, int)?; + Exact::new(if sign == Sign::Negative { -val } else { val }, false) + } + }) + } + + // cos works for all real numbers + pub(crate) fn cos(self, int: &I) -> FResult> { + Ok(if self == 0.into() { + Exact::new(Self::from(1), true) + } else { + let pi = Self::high_precision_pi(); + let two_pi = Self::high_precision_two_pi(); + let pi_over_2 = Self::high_precision_pi_over_2(); + let pi_over_4 = Self::high_precision_pi_over_4(); + + let mut x = self.rem_euclid_rat(two_pi.clone(), int)?; + let mut sign = Sign::Positive; + + if x >= pi { + x = two_pi.sub(x, int)?; + } + if x > pi_over_2 { + x = pi.sub(x, int)?; + sign = sign.flip(); + } + + if x > pi_over_4 { + x = pi_over_2.sub(x, int)?; + let val = Self::sin_series(x, int)?; + Exact::new(if sign == Sign::Negative { -val } else { val }, false) + } else { + let val = Self::cos_series(&x, int)?; + Exact::new(if sign == Sign::Negative { -val } else { val }, false) + } + }) + } + + // asin, acos and atan only work for values between -1 and 1 + pub(crate) fn atan(self, int: &I) -> FResult { + let one = Self::from(1); + let mut x = self.clone().abs(int); + let mut invert = false; + + if x > one { + x = one.clone().div(&x, int)?; + invert = true; + } + + let mut shift_pi_4 = false; + + if &x > Self::atan_reduction_threshold() { + let num = x.clone().sub(one.clone(), int)?; + let den = one.clone().add(x, int)?; + x = num.div(&den, int)?; + shift_pi_4 = true; + } + + let mut result = Self::atan_series(x, int)?; + + if shift_pi_4 { + result = result.add(Self::high_precision_pi_over_4(), int)?; + } + + if invert { + let pi_over_2 = Self::high_precision_pi_over_2(); + result = pi_over_2.sub(result, int)?; + } + + if self.sign == Sign::Negative { + Ok(-result) + } else { + Ok(result) + } + } + + fn atan_reduction_threshold() -> &'static Self { + static THRESHOLD: OnceLock = OnceLock::new(); + THRESHOLD.get_or_init(|| { + // 0.4142 = 2071 / 5000 + Self { + sign: Sign::Positive, + num: BigUint::from(2071u64), + den: BigUint::from(5000u64), + } + }) + } + + pub(crate) fn asin(self, int: &I) -> FResult { + let one = Self::from(1); + + // Check domain: -1 <= x <= 1 + if self > one || self < -one.clone() { + return Err(out_of_range(self.fm(int)?, Range::open(-1, 1))); + } + + // Handle edges x = 1 and x = -1 explicitly to avoid division by zero in the formula + if self == one { + return Ok(Self::high_precision_pi_over_2()); + } + if self == -one.clone() { + return Ok(-Self::high_precision_pi_over_2()); + } + + // asin(x) = atan(x / sqrt(1 - x^2)) + let x_sq = self.clone().mul(&self, int)?; + let one_minus_x_sq = one.sub(x_sq, int)?; + + // Compute sqrt(1 - x^2) + // Using root_n(2) which relies on the high-precision iter_root_n + let sqrt_den = one_minus_x_sq.root_n(&Self::from(2), int)?; + + let argument = self.div(&sqrt_den.value, int)?; + argument.atan(int) + } + + pub(crate) fn acos(self, int: &I) -> FResult { + let one = Self::from(1); + + // Check domain: -1 <= x <= 1 + if self > one || self < -one { + return Err(out_of_range(self.fm(int)?, Range::open(-1, 1))); + } + + // acos(x) = pi/2 - asin(x) + let pi_over_2 = Self::high_precision_pi_over_2(); + let asin_val = self.asin(int)?; + + pi_over_2.sub(asin_val, int) + } + + pub(crate) fn atan2(self, rhs: Self, int: &I) -> FResult { + let y = self; + let x = rhs; + + if x.is_definitely_zero() { + if y.is_definitely_zero() { + // atan2(0, 0) is undefined + return Err(FendError::DivideByZero); + } + let pi_over_2 = Self::high_precision_pi_over_2(); + if y.sign == Sign::Positive { + return Ok(pi_over_2); + } + return Ok(-pi_over_2); + } + + if x.sign == Sign::Positive { + // Quadrants 1 and 4: atan(y/x) + let ratio = y.div(&x, int)?; + return ratio.atan(int); + } + + // Quadrants 2 and 3: x is negative + let ratio = y.clone().div(&x, int)?; + let atan_val = ratio.atan(int)?; + let pi = Self::high_precision_pi(); + + if y.sign == Sign::Positive || y.is_definitely_zero() { + // Quadrant 2 (y >= 0, x < 0): result is atan(y/x) + pi + // Note: atan(y/x) will be negative + atan_val.add(pi, int) + } else { + // Quadrant 3 (y < 0, x < 0): result is atan(y/x) - pi + // Note: atan(y/x) will be positive + atan_val.sub(pi, int) + } } fn apply_uint_op( @@ -376,19 +1034,89 @@ impl BigRat { Ok(self.apply_uint_op(BigUint::factorial, int)?.into()) } - pub(crate) fn floor(self, int: &I) -> FResult { - let float = self.into_f64(int)?.floor(); - Self::from_f64(float, int) + pub(crate) fn floor(mut self, int: &I) -> FResult { + if self.is_definitely_zero() { + return Ok(self); + } + self = self.simplify(int)?; + + // Use integer division. + // divmod takes &self.den, so we don't move anything. + let (q, rem) = self.num.divmod(&self.den, int)?; + + if self.sign == Sign::Positive { + // Example: 10/3 = 3 rem 1 -> floor is 3 + Ok(Self { + sign: Sign::Positive, + num: q, + den: 1.into(), + }) + } else { + // Example: -10/3 = -3 rem 1 -> floor is -4 + // Example: -9/3 = -3 rem 0 -> floor is -3 + let ans = if rem.is_definitely_zero() { + q + } else { + q.add(&1.into()) + }; + Ok(Self { + sign: Sign::Negative, + num: ans, + den: 1.into(), + }) + } } - pub(crate) fn ceil(self, int: &I) -> FResult { - let float = self.into_f64(int)?.ceil(); - Self::from_f64(float, int) + pub(crate) fn ceil(mut self, int: &I) -> FResult { + if self.is_definitely_zero() { + return Ok(self); + } + self = self.simplify(int)?; + + let (q, rem) = self.num.divmod(&self.den, int)?; + + if self.sign == Sign::Positive { + // Example: 10/3 = 3 rem 1 -> ceil is 4 + let ans = if rem.is_definitely_zero() { + q + } else { + q.add(&1.into()) + }; + Ok(Self { + sign: Sign::Positive, + num: ans, + den: 1.into(), + }) + } else { + // Example: -10/3 = -3 rem 1 -> ceil is -3 + Ok(Self { + sign: Sign::Negative, + num: q, + den: 1.into(), + }) + } } - pub(crate) fn round(self, int: &I) -> FResult { - let float = self.into_f64(int)?.round(); - Self::from_f64(float, int) + pub(crate) fn round(mut self, int: &I) -> FResult { + if self.is_definitely_zero() { + return Ok(self); + } + self = self.simplify(int)?; + + let (q, rem) = self.num.divmod(&self.den, int)?; + + // Rounding logic: Round up if remainder >= denominator / 2 + // Equivalent to: 2 * remainder >= denominator + let rem_2 = rem.mul(&2.into(), int)?; + let round_up = rem_2 >= self.den; + + let ans = if round_up { q.add(&1.into()) } else { q }; + + Ok(Self { + sign: self.sign, + num: ans, + den: 1.into(), + }) } pub(crate) fn bitwise( @@ -959,11 +1687,11 @@ impl BigRat { pub(crate) fn pow(mut self, mut rhs: Self, int: &I) -> FResult> { self = self.simplify(int)?; rhs = rhs.simplify(int)?; + if self.num != 0.into() && self.sign == Sign::Negative && rhs.den != 1.into() { return Err(FendError::RootsOfNegativeNumbers); } if rhs.sign == Sign::Negative { - // a^-b => 1/a^b rhs.sign = Sign::Positive; let inverse_res = self.pow(rhs, int)?; return Ok(Exact::new( @@ -971,16 +1699,51 @@ impl BigRat { inverse_res.exact, )); } + + // Use exp/ln approximation if: + // 1. The denominator is large (e.g. ^0.0001), indicating a complex decimal root. + // 2. The numerator is large AND it's a root (e.g. ^210.1 = ^2101/10). + let large_root = rhs.den > BigUint::from(1000u64); + let large_power = rhs.den > BigUint::from(1u64) && rhs.num > BigUint::from(100u64); + + if large_root || large_power { + // Case 1: Simple fractional root (e.g. ^0.00001) where we might want exactness + if rhs.num == 1.into() && !large_power { + let n = rhs.den.clone(); + let num_root = self.num.clone().root_n(&n, int)?; + let den_root = self.den.clone().root_n(&n, int)?; + + if num_root.exact && den_root.exact { + return Ok(Exact::new( + Self { + sign: Sign::Positive, + num: num_root.value, + den: den_root.value, + }, + true, + )); + } + } + + // Case 2: Complex decimal or large hybrid power -> Fast Approximation + let ln_x = self.ln(int)?; + let y_ln_x = ln_x.value.mul(&rhs, int)?; + return y_ln_x.exp(int); + } + + // Standard Path (Small denominators & Small numerators, e.g. ^0.5, ^2/3) let result_sign = if self.sign == Sign::Positive || rhs.num.is_even(int)? { Sign::Positive } else { Sign::Negative }; + let pow_res = Self { sign: result_sign, num: BigUint::pow(&self.num, &rhs.num, int)?, den: BigUint::pow(&self.den, &rhs.num, int)?, }; + if rhs.den == 1.into() { Ok(Exact::new(pow_res, true)) } else { @@ -1003,7 +1766,10 @@ impl BigRat { int: &I, ) -> FResult { let mut high_bound = low_bound.clone().add(1.into(), int)?; - for _ in 0..50 { + // Each iteration adds 1 bit of precision. + // Since log2(10) ≈ 3.32, we need ~3.32 iterations per decimal digit. Use 3.35 as safety factor. + let iterations = (Self::INTERNAL_PRECISION * 335) / 100; + for _ in 0..iterations { let guess = low_bound .clone() .add(high_bound.clone(), int)? @@ -1017,16 +1783,6 @@ impl BigRat { low_bound.add(high_bound, int)?.div(&2.into(), int) } - pub(crate) fn exp(self, int: &I) -> FResult> { - if self.num == 0.into() { - return Ok(Exact::new(Self::from(1), true)); - } - Ok(Exact::new( - Self::from_f64(self.into_f64(int)?.exp(), int)?, - false, - )) - } - // the boolean indicates whether or not the result is exact // n must be an integer pub(crate) fn root_n(self, n: &Self, int: &I) -> FResult> { diff --git a/core/tests/integration_tests.rs b/core/tests/integration_tests.rs index a55b7d4a..02c876a7 100644 --- a/core/tests/integration_tests.rs +++ b/core/tests/integration_tests.rs @@ -409,7 +409,7 @@ fn inverse_sin_point_five() { #[test] fn inverse_sin_nested() { - test_eval("sin^-1 (sin 0.5", "approx. 0.5"); + test_eval("sin^-1 (sin 0.5)", "approx. 0.4999999999"); } #[test] @@ -1214,10 +1214,15 @@ fn powers_17() { fn powers_18() { test_eval( "5.2*10^15*300^(3/2)", - "approx. 27019992598074485776.9266786817", + "approx. 27019992598074485779.0281629274", ); } +#[test] +fn power_exponent_2_914() { + test_eval("711.49^2.914", "approx. 204748101.9700616257") +} + #[test] fn pi_to_the_power_of_ten() { test_eval("pi^10", "approx. 93648.047476083"); @@ -3192,6 +3197,16 @@ fn atan_1() { test_eval("atan 1", "approx. 0.7853981633"); } +#[test] +fn atan_zero_point_five() { + test_eval("atan 0.5", "approx. 0.463647609"); +} + +#[test] +fn atan_minus_zero_point_three_one() { + test_eval("atan (-0.31)", "approx. -0.30060567"); +} + #[test] fn sinh_0() { test_eval("sinh 0", "approx. 0"); @@ -3279,12 +3294,12 @@ fn log10_100() { #[test] fn log10_1000() { - test_eval("log10 1000", "approx. 3"); + test_eval("log10 1000", "approx. 2.9999999999"); } #[test] fn log10_10000() { - test_eval("log10 10000", "approx. 3.9999999999"); + test_eval("log10 10000", "approx. 4"); } #[test] @@ -3299,12 +3314,12 @@ fn log_100() { #[test] fn log_1000() { - test_eval("log 1000", "approx. 3"); + test_eval("log 1000", "approx. 2.9999999999"); } #[test] fn log_10000() { - test_eval("log 10000", "approx. 3.9999999999"); + test_eval("log 10000", "approx. 4"); } #[test] @@ -3339,7 +3354,7 @@ fn log2_minus_1() { #[test] fn sqrt_minus_two() { - test_eval_simple("sqrt(-2)", "approx. 0 + 1.4142135623i"); + test_eval("sqrt(-2)", "approx. 1.4142135623i"); } #[test] @@ -3414,8 +3429,7 @@ fn sqrt_i() { #[test] fn sqrt_minus_two_i() { - // FIXME: exactly 1 - i - test_eval("sqrt (-2i)", "approx. 0.9999999999 - 0.9999999999i"); + test_eval("sqrt (-2i)", "approx. 1 - 0.9999999999i"); } #[test]