Evaluate built-in mathematical constants

- Fix digit count for negative Python decimals - Make sure state is initialized when querying built-in constants
benruijl · Oct 14, 2024 · 9de3aa2 · 9de3aa2
1 parent 7e52bd3
commit 9de3aa2
Show file tree

Hide file tree

Showing 6 changed files with 194 additions and 6 deletions.
diff --git 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
@@ -731,6 +731,34 @@ macro_rules! create_hyperdual_from_components {
         }
 
         impl<T: $crate::domains::float::Real> $crate::domains::float::Real for $t<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
@@ -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<T: RealNumberLike> RealNumberLike for ErrorPropagatingFloat<T> {
 }
 
 impl<T: Real + RealNumberLike> Real for ErrorPropagatingFloat<T> {
+    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<T: NumericalFloatLike> NumericalFloatLike for Complex<T> {
 
 /// Following the same conventions and formulas as num::Complex.
 impl<T: Real> Real for Complex<T> {
+    #[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
@@ -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
@@ -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<Item = (Symbol, &'static str)> {
+        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
@@ -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)`.