From 9de3aa22b2d1d7864a19277603b14b5b7da6aece Mon Sep 17 00:00:00 2001 From: Ben Ruijl Date: Mon, 14 Oct 2024 10:45:30 +0200 Subject: [PATCH] Evaluate built-in mathematical constants - Fix digit count for negative Python decimals - Make sure state is initialized when querying built-in constants --- src/api/python.rs | 2 +- src/domains/dual.rs | 28 +++++++++ src/domains/float.rs | 137 +++++++++++++++++++++++++++++++++++++++++++ src/evaluate.rs | 12 ++-- src/state.rs | 19 ++++++ symbolica.pyi | 2 +- 6 files changed, 194 insertions(+), 6 deletions(-) diff --git a/src/api/python.rs b/src/api/python.rs index 0343d95..6e138a1 100644 --- a/src/api/python.rs +++ b/src/api/python.rs @@ -1753,7 +1753,7 @@ impl<'a> FromPyObject<'a> for PythonMultiPrecisionFloat { // get the number of accurate digits let digits = a .chars() - .skip_while(|x| *x == '.' || *x == '0') + .skip_while(|x| *x == '.' || *x == '0' || *x == '-') .filter(|x| *x != '.') .take_while(|x| x.is_ascii_digit()) .count(); diff --git a/src/domains/dual.rs b/src/domains/dual.rs index 937ac81..15491df 100644 --- a/src/domains/dual.rs +++ b/src/domains/dual.rs @@ -731,6 +731,34 @@ macro_rules! create_hyperdual_from_components { } impl $crate::domains::float::Real for $t { + #[inline(always)] + fn pi(&self) -> Self { + let mut res = self.zero(); + res.values[0] = self.values[0].pi(); + res + } + + #[inline(always)] + fn e(&self) -> Self { + let mut res = self.zero(); + res.values[0] = self.values[0].e(); + res + } + + #[inline(always)] + fn euler(&self) -> Self { + let mut res = self.zero(); + res.values[0] = self.values[0].euler(); + res + } + + #[inline(always)] + fn phi(&self) -> Self { + let mut res = self.zero(); + res.values[0] = self.values[0].phi(); + res + } + #[inline(always)] fn norm(&self) -> Self { let n = self.values[0].norm(); diff --git a/src/domains/float.rs b/src/domains/float.rs index 54dcacc..1f58473 100644 --- a/src/domains/float.rs +++ b/src/domains/float.rs @@ -275,6 +275,15 @@ pub trait ConstructibleFloat: NumericalFloatLike { } pub trait Real: NumericalFloatLike { + /// The constant π, 3.1415926535... + fn pi(&self) -> Self; + /// Euler's number, 2.7182818... + fn e(&self) -> Self; + /// The Euler-Mascheroni constant, 0.5772156649... + fn euler(&self) -> Self; + /// The golden ratio, 1.6180339887... + fn phi(&self) -> Self; + fn norm(&self) -> Self; fn sqrt(&self) -> Self; fn log(&self) -> Self; @@ -426,6 +435,26 @@ impl ConstructibleFloat for f64 { } impl Real for f64 { + #[inline(always)] + fn pi(&self) -> Self { + std::f64::consts::PI + } + + #[inline(always)] + fn e(&self) -> Self { + std::f64::consts::E + } + + #[inline(always)] + fn euler(&self) -> Self { + 0.57721566490153286 + } + + #[inline(always)] + fn phi(&self) -> Self { + 1.6180339887498948 + } + #[inline(always)] fn norm(&self) -> Self { f64::abs(*self) @@ -777,6 +806,26 @@ impl RealNumberLike for F64 { } impl Real for F64 { + #[inline(always)] + fn pi(&self) -> Self { + std::f64::consts::PI.into() + } + + #[inline(always)] + fn e(&self) -> Self { + std::f64::consts::E.into() + } + + #[inline(always)] + fn euler(&self) -> Self { + 0.57721566490153286.into() + } + + #[inline(always)] + fn phi(&self) -> Self { + 1.6180339887498948.into() + } + #[inline(always)] fn norm(&self) -> Self { self.0.norm().into() @@ -1638,6 +1687,26 @@ impl RealNumberLike for Float { } impl Real for Float { + #[inline(always)] + fn pi(&self) -> Self { + MultiPrecisionFloat::with_val(self.prec(), rug::float::Constant::Pi).into() + } + + #[inline(always)] + fn e(&self) -> Self { + self.one().exp() + } + + #[inline(always)] + fn euler(&self) -> Self { + MultiPrecisionFloat::with_val(self.prec(), rug::float::Constant::Euler).into() + } + + #[inline(always)] + fn phi(&self) -> Self { + (self.one() + self.from_i64(5).sqrt()) / 2 + } + #[inline(always)] fn norm(&self) -> Self { self.0.clone().abs().into() @@ -2188,6 +2257,34 @@ impl RealNumberLike for ErrorPropagatingFloat { } impl Real for ErrorPropagatingFloat { + fn pi(&self) -> Self { + ErrorPropagatingFloat { + value: self.value.pi(), + prec: self.prec, + } + } + + fn e(&self) -> Self { + ErrorPropagatingFloat { + value: self.value.e(), + prec: self.prec, + } + } + + fn euler(&self) -> Self { + ErrorPropagatingFloat { + value: self.value.euler(), + prec: self.prec, + } + } + + fn phi(&self) -> Self { + ErrorPropagatingFloat { + value: self.value.phi(), + prec: self.prec, + } + } + fn norm(&self) -> Self { ErrorPropagatingFloat { value: self.value.norm(), @@ -2422,6 +2519,26 @@ macro_rules! simd_impl { } impl Real for $t { + #[inline(always)] + fn pi(&self) -> Self { + std::f64::consts::PI.into() + } + + #[inline(always)] + fn e(&self) -> Self { + std::f64::consts::E.into() + } + + #[inline(always)] + fn euler(&self) -> Self { + 0.57721566490153286.into() + } + + #[inline(always)] + fn phi(&self) -> Self { + 1.6180339887498948.into() + } + #[inline(always)] fn norm(&self) -> Self { (*self).abs() @@ -3344,6 +3461,26 @@ impl NumericalFloatLike for Complex { /// Following the same conventions and formulas as num::Complex. impl Real for Complex { + #[inline] + fn pi(&self) -> Self { + Complex::new(self.re.pi(), self.im.zero()) + } + + #[inline] + fn e(&self) -> Self { + Complex::new(self.re.e(), self.im.zero()) + } + + #[inline] + fn euler(&self) -> Self { + Complex::new(self.re.euler(), self.im.zero()) + } + + #[inline] + fn phi(&self) -> Self { + Complex::new(self.re.phi(), self.im.zero()) + } + #[inline] fn norm(&self) -> Self { Complex::new(self.norm_squared().sqrt(), self.im.zero()) diff --git a/src/evaluate.rs b/src/evaluate.rs index a03658c..8f11af2 100644 --- a/src/evaluate.rs +++ b/src/evaluate.rs @@ -4086,10 +4086,14 @@ impl<'a> AtomView<'a> { "Rational polynomial coefficient not yet supported for evaluation".to_string(), ), }, - AtomView::Var(v) => Err(format!( - "Variable {} not in constant map", - State::get_name(v.get_symbol()) - )), + AtomView::Var(v) => match v.get_symbol() { + State::E => Ok(coeff_map(&1.into()).e()), + State::PI => Ok(coeff_map(&1.into()).pi()), + _ => Err(format!( + "Variable {} not in constant map", + State::get_name(v.get_symbol()) + )), + }, AtomView::Fun(f) => { let name = f.get_symbol(); if [State::EXP, State::LOG, State::SIN, State::COS, State::SQRT].contains(&name) { diff --git a/src/state.rs b/src/state.rs index 7c26165..4c15a16 100644 --- a/src/state.rs +++ b/src/state.rs @@ -199,12 +199,21 @@ impl State { } } + #[inline(always)] pub(crate) unsafe fn symbol_from_id(id: u32) -> Symbol { + if ID_TO_STR.len() == 0 { + let _ = *STATE; // initialize the state + } + ID_TO_STR[id as usize].0 } /// Iterate over all defined symbols. pub fn symbol_iter() -> impl Iterator { + if ID_TO_STR.len() == 0 { + let _ = *STATE; // initialize the state + } + ID_TO_STR .iter() .skip(SYMBOL_OFFSET.load(Ordering::Relaxed)) @@ -410,14 +419,24 @@ impl State { } /// Get the name for a given symbol. + #[inline] pub fn get_name(id: Symbol) -> &'static str { + if ID_TO_STR.len() == 0 { + let _ = *STATE; // initialize the state + } + &ID_TO_STR[id.get_id() as usize + SYMBOL_OFFSET.load(Ordering::Relaxed)] .1 .name } /// Get the user-specified normalization function for the symbol. + #[inline] pub fn get_normalization_function(id: Symbol) -> Option<&'static NormalizationFunction> { + if ID_TO_STR.len() == 0 { + let _ = *STATE; // initialize the state + } + ID_TO_STR[id.get_id() as usize + SYMBOL_OFFSET.load(Ordering::Relaxed)] .1 .function diff --git a/symbolica.pyi b/symbolica.pyi index 30b374c..185f908 100644 --- a/symbolica.pyi +++ b/symbolica.pyi @@ -3590,7 +3590,7 @@ class Integer: """Check if the 64-bit number `n` is a prime number.""" @classmethod - def solve_integer_relation(_cls, x: Sequence[int | float | Decimal], tolerance: float, max_coeff: Optional[int] = None, gamma: Optional[float | Decimal] = None) -> Sequence[int]: + def solve_integer_relation(_cls, x: Sequence[int | float | Decimal], tolerance: float | Decimal, max_coeff: Optional[int] = None, gamma: Optional[float | Decimal] = None) -> Sequence[int]: """Use the PSLQ algorithm to find a vector of integers `a` that satisfies `a.x = 0`, where every element of `a` is less than `max_coeff`, using a specified tolerance and number of iterations. The parameter `gamma` must be more than or equal to `2/sqrt(3)`.