diff --git a/src/api/python.rs b/src/api/python.rs index 43117c7..22a9a83 100644 --- a/src/api/python.rs +++ b/src/api/python.rs @@ -1157,24 +1157,13 @@ impl<'a> FromPyObject<'a> for ConvertibleToExpression { let a = format!("{}", num); let i = Integer::from_large(rug::Integer::parse(&a).unwrap().complete()); Ok(ConvertibleToExpression(Atom::new_num(i).into())) - } else if let Ok(f) = ob.extract::() { - if !f.is_finite() { - return Err(exceptions::PyValueError::new_err("Number must be finite")); - } - - Ok(ConvertibleToExpression( - Atom::new_num(rug::Float::with_val(53, f)).into(), - )) - } else if ob.is_instance(get_decimal(ob.py()).as_ref(ob.py()))? { - let a = ob.call_method0("__str__").unwrap().extract::<&str>()?; - Ok(ConvertibleToExpression( - Atom::new_num( - Float::parse(a) - .unwrap() - .complete((a.len() as f64 * LOG2_10).ceil() as u32), - ) - .into(), + } else if let Ok(_) = ob.extract::<&str>() { + // disallow direct string conversion + Err(exceptions::PyValueError::new_err( + "Cannot convert to expression", )) + } else if let Ok(f) = ob.extract::() { + Ok(ConvertibleToExpression(Atom::new_num(f.0).into())) } else { Err(exceptions::PyValueError::new_err( "Cannot convert to expression", @@ -1263,11 +1252,32 @@ impl<'a> FromPyObject<'a> for PythonMultiPrecisionFloat { fn extract(ob: &'a pyo3::PyAny) -> PyResult { if ob.is_instance(get_decimal(ob.py()).as_ref(ob.py()))? { let a = ob.call_method0("__str__").unwrap().extract::<&str>()?; - Ok(Float::parse(a, Some(a.len() as u32)).unwrap().into()) + + // get the number of accurate digits + let digits = a + .chars() + .skip_while(|x| *x == '.' || *x == '0') + .take_while(|x| x.is_ascii_digit()) + .count(); + + Ok(Float::parse( + a, + Some((digits as f64 * std::f64::consts::LOG2_10).ceil() as u32), + ) + .map_err(|_| exceptions::PyValueError::new_err("Not a floating point number"))? + .into()) } else if let Ok(a) = ob.extract::<&str>() { - Ok(Float::parse(a, None).unwrap().into()) + Ok(Float::parse(a, None) + .map_err(|_| exceptions::PyValueError::new_err("Not a floating point number"))? + .into()) } else if let Ok(a) = ob.extract::() { - Ok(Float::with_val(53, a).into()) + if a.is_finite() { + Ok(Float::with_val(53, a).into()) + } else { + Err(exceptions::PyValueError::new_err( + "Floating point number is not finite", + )) + } } else { Err(exceptions::PyValueError::new_err( "Not a valid multi-precision float", @@ -1522,9 +1532,9 @@ impl PythonExpression { Ok(result) } - /// Create a new Symbolica number from an int, a float, or a string. + /// Create a new Symbolica number from an int, a float, a Decimal, or a string. /// A floating point number is kept as a float with the same precision as the input, - /// but it can also be truncated by specifying the maximal denominator value. + /// but it can also be converted to the smallest rational number given a `relative_error`. /// /// Examples /// -------- @@ -1533,37 +1543,33 @@ impl PythonExpression { /// 1/2 /// /// >>> print(Expression.num(1/3)) - /// >>> print(Expression.num(0.33, 5)) + /// >>> print(Expression.num(0.33, 0.1)) + /// >>> print(Expression.num('0.333`3')) + /// >>> print(Expression.num(Decimal('0.1234'))) /// 3.3333333333333331e-1 /// 1/3 + /// 3.33e-1 + /// 1.2340e-1 #[classmethod] pub fn num( _cls: &PyType, py: Python, num: PyObject, - max_denom: Option, + relative_error: Option, ) -> PyResult { if let Ok(num) = num.extract::(py) { Ok(Atom::new_num(num).into()) } else if let Ok(num) = num.extract::<&PyLong>(py) { let a = format!("{}", num); PythonExpression::parse(_cls, &a) - } else if let Ok(f) = num.extract::(py) { - if !f.is_finite() { - return Err(exceptions::PyValueError::new_err("Number must be finite")); - } - - if let Some(max_denom) = max_denom { - let mut r: Rational = f.into(); - r = r.truncate_denominator(&(max_denom as u64).into()); + } else if let Ok(f) = num.extract::(py) { + if let Some(relative_error) = relative_error { + let mut r: Rational = f.0.into(); + r = r.round(&relative_error.into()).into(); Ok(Atom::new_num(r).into()) } else { - Ok(Atom::new_num(rug::Float::with_val(53, f)).into()) + Ok(Atom::new_num(f.0).into()) } - } else if let Ok(f) = num.extract::(py) { - Ok(Atom::new_num(f.0.clone()).into()) - } else if let Ok(a) = num.extract::<&str>(py) { - PythonExpression::parse(_cls, a) } else { Err(exceptions::PyValueError::new_err("Not a valid number")) } @@ -2008,14 +2014,25 @@ impl PythonExpression { } } - /// Convert all coefficients to floats, with a given decimal precision. - pub fn to_float(&self, decimal_prec: u32) -> PythonExpression { - self.expr.to_float(decimal_prec).into() + /// Convert all coefficients to floats with a given precision `decimal_prec``. + /// The precision of floating point coefficients in the input will be truncated to `decimal_prec`. + pub fn coefficients_to_float(&self, decimal_prec: u32) -> PythonExpression { + self.expr.coefficients_to_float(decimal_prec).into() } - /// Convert all floating point coefficients to rationals, with a given maximal denominator. - pub fn float_to_rat(&self, max_denominator: Integer) -> PythonExpression { - self.expr.float_to_rat(&max_denominator).into() + /// Map all floating point and rational coefficients to the best rational approximation + /// in the interval `[self*(1-relative_error),self*(1+relative_error)]`. + pub fn rationalize_coefficients(&self, relative_error: f64) -> PyResult { + if relative_error <= 0. || relative_error > 1. { + return Err(exceptions::PyValueError::new_err( + "Relative error must be between 0 and 1", + )); + } + + Ok(self + .expr + .rationalize_coefficients(&relative_error.into()) + .into()) } /// Create a pattern restriction based on the wildcard length before downcasting. diff --git a/src/atom/representation.rs b/src/atom/representation.rs index f85c8d8..0ffd6e2 100644 --- a/src/atom/representation.rs +++ b/src/atom/representation.rs @@ -169,6 +169,7 @@ impl Atom { Ok(a.as_view().rename(state_map)) } + #[allow(dead_code)] pub(crate) unsafe fn from_raw(raw: RawAtom) -> Self { match raw[0] & TYPE_MASK { NUM_ID => Atom::Num(Num::from_raw(raw)), diff --git a/src/coefficient.rs b/src/coefficient.rs index 9ce5f35..6765d8d 100644 --- a/src/coefficient.rs +++ b/src/coefficient.rs @@ -859,21 +859,39 @@ impl Atom { self.as_view().set_coefficient_ring(vars) } - /// Convert all coefficients to floats with a given precision. - pub fn to_float(&self, decimal_prec: u32) -> Atom { + /// Convert all coefficients to floats with a given precision `decimal_prec``. + /// The precision of floating point coefficients in the input will be truncated to `decimal_prec`. + pub fn coefficients_to_float(&self, decimal_prec: u32) -> Atom { let mut a = Atom::new(); - self.as_view().to_float_into(decimal_prec, &mut a); + self.as_view() + .coefficients_to_float_into(decimal_prec, &mut a); a } - /// Convert all coefficients to floats with a given precision. - pub fn to_float_into(&self, decimal_prec: u32, out: &mut Atom) { - self.as_view().to_float_into(decimal_prec, out); + /// Convert all coefficients to floats with a given precision `decimal_prec``. + /// The precision of floating point coefficients in the input will be truncated to `decimal_prec`. + pub fn coefficients_to_float_into(&self, decimal_prec: u32, out: &mut Atom) { + self.as_view().coefficients_to_float_into(decimal_prec, out); } - /// Map all floating point coefficients to rational numbers, using a given maximum denominator. - pub fn float_to_rat(&self, max_denominator: &Integer) -> Atom { - self.as_view().float_to_rat(max_denominator) + /// Map all coefficients using a given function. + pub fn map_coefficient Coefficient + Copy>(&self, f: F) -> Atom { + self.as_view().map_coefficient(f) + } + + /// Map all coefficients using a given function. + pub fn map_coefficient_into Coefficient + Copy>( + &self, + f: F, + out: &mut Atom, + ) { + self.as_view().map_coefficient_into(f, out); + } + + /// Map all floating point and rational coefficients to the best rational approximation + /// in the interval `[self*(1-relative_error),self*(1+relative_error)]`. + pub fn rationalize_coefficients(&self, relative_error: &Rational) -> Atom { + self.as_view().rationalize_coefficients(relative_error) } } @@ -1052,15 +1070,16 @@ impl<'a> AtomView<'a> { } } - /// Convert all coefficients to floats with a given precision. - pub fn to_float(&self, decimal_prec: u32) -> Atom { + /// Convert all coefficients to floats with a given precision `decimal_prec``. + /// The precision of floating point coefficients in the input will be truncated to `decimal_prec`. + pub fn coefficients_to_float(&self, decimal_prec: u32) -> Atom { let mut a = Atom::new(); - self.to_float_into(decimal_prec, &mut a); + self.coefficients_to_float_into(decimal_prec, &mut a); a } - - /// Convert all coefficients to floats with a given precision. - pub fn to_float_into(&self, decimal_prec: u32, out: &mut Atom) { + /// Convert all coefficients to floats with a given precision `decimal_prec``. + /// The precision of floating point coefficients in the input will be truncated to `decimal_prec`. + pub fn coefficients_to_float_into(&self, decimal_prec: u32, out: &mut Atom) { let binary_prec = (decimal_prec as f64 * LOG2_10).ceil() as u32; Workspace::get_local().with(|ws| self.to_float_impl(binary_prec, true, false, ws, out)) @@ -1083,7 +1102,7 @@ impl<'a> AtomView<'a> { } CoefficientView::Float(f) => { let mut f = f.to_float(); - if f.prec() != binary_prec { + if f.prec() > binary_prec { f.set_prec(binary_prec); out.to_num(Coefficient::Float(f)); } else { @@ -1181,20 +1200,31 @@ impl<'a> AtomView<'a> { } } - /// Map all floating point coefficients to rational numbers, using a given maximum denominator. - pub fn float_to_rat(&self, max_denominator: &Integer) -> Atom { + /// Map all floating point and rational coefficients to the best rational approximation + /// in the interval `[self*(1-relative_error),self*(1+relative_error)]`. + pub fn rationalize_coefficients(&self, relative_error: &Rational) -> Atom { let mut a = Atom::new(); - self.map_coefficient_into( - |c| match c { - CoefficientView::Float(f) => f - .to_float() - .to_rational() - .truncate_denominator(max_denominator) - .into(), - _ => c.to_owned(), - }, - &mut a, - ); + Workspace::get_local().with(|ws| { + self.map_coefficient_impl( + |c| match c { + CoefficientView::Float(f) => { + f.to_float().to_rational().round(relative_error).into() + } + CoefficientView::Natural(n, d) => { + let r = Rational::new(n, d); + r.round(relative_error).into() + } + CoefficientView::Large(r) => Rational::from_large(r.to_rat()) + .round(relative_error) + .into(), + _ => c.to_owned(), + }, + true, + false, + ws, + &mut a, + ) + }); a } @@ -1375,26 +1405,27 @@ mod test { ); assert_eq!( r, - "1.5707963267948966192313216916397514420985846996875529104874722*x+21.472450421034925" + "1.57079632679489661923132169163975144209858469968755291048747*x+21.4724504210349" ); } #[test] fn float_convert() { let expr = Atom::parse("1/2 x + 238947/128903718927 + sin(3/4)").unwrap(); - let expr = expr.to_float(60); + let expr = expr.coefficients_to_float(60); let r = format!( "{}", AtomPrinter::new_with_options(expr.as_view(), PrintOptions::file()) ); - assert_eq!(r, "5.0000000000000000000000000000000000000000000000000000000000000e-1*x+6.8164061370918581635917066956651198726148569775622233288512875e-1"); + assert_eq!(r, "5.00000000000000000000000000000000000000000000000000000000000e-1*x+6.81640613709185816359170669566511987261485697756222332885129e-1"); } #[test] fn float_to_rat() { let expr = Atom::parse("1/2 x + 238947/128903718927 + sin(3/4)").unwrap(); - let expr = expr.to_float(60); - let expr = expr.float_to_rat(&1000.into()); - assert_eq!(expr, Atom::parse("1/2*x+349/512").unwrap()); + let expr = expr.coefficients_to_float(60); + let expr = expr.rationalize_coefficients(&(1, 10000).into()); + println!("expr {}", expr); + assert_eq!(expr, Atom::parse("1/2*x+137/201").unwrap()); } } diff --git a/src/domains/float.rs b/src/domains/float.rs index e608b85..c43462a 100644 --- a/src/domains/float.rs +++ b/src/domains/float.rs @@ -1,5 +1,5 @@ use std::{ - f64::consts::LOG2_10, + f64::consts::{LOG10_2, LOG2_10}, fmt::{self, Debug, Display, Formatter, LowerExp, Write}, ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign}, }; @@ -338,13 +338,32 @@ impl Debug for Float { impl Display for Float { fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result { - Display::fmt(&self.0, f) + // print only the significant digits + // the original float value may not be reconstructible + // from this output + if f.precision().is_none() { + f.write_fmt(format_args!( + "{0:.1$}", + self.0, + (self.0.prec() as f64 * LOG10_2).floor() as usize + )) + } else { + Display::fmt(&self.0, f) + } } } impl LowerExp for Float { fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result { - LowerExp::fmt(&self.0, f) + if f.precision().is_none() { + f.write_fmt(format_args!( + "{0:.1$e}", + self.0, + (self.0.prec() as f64 * LOG10_2).floor() as usize + )) + } else { + LowerExp::fmt(&self.0, f) + } } } @@ -799,6 +818,10 @@ impl Float { self.0.set_prec(prec); } + pub fn is_finite(&self) -> bool { + self.0.is_finite() + } + /// Parse a float from a string. /// Precision can be specified by a trailing backtick followed by the precision. /// For example: ```1.234`20``` for a precision of 20 decimal digits. @@ -824,7 +847,14 @@ impl Float { .complete(prec), )) } else { - let prec = ((s.len() as f64 * LOG2_10).ceil() as u32).max(53); + // get the number of accurate digits + let digits = s + .chars() + .skip_while(|x| *x == '.' || *x == '0') + .take_while(|x| x.is_ascii_digit()) + .count(); + + let prec = ((digits as f64 * LOG2_10).ceil() as u32).max(53); Ok(Float( MultiPrecisionFloat::parse(s) .map_err(|e| e.to_string())? @@ -1856,6 +1886,12 @@ macro_rules! simd_impl { simd_impl!(f64x2, pow_f64x2); simd_impl!(f64x4, pow_f64x4); +impl From for Rational { + fn from(value: Float) -> Self { + value.to_rational() + } +} + impl LowerExp for Rational { fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { // lower-exp is ignored for integers diff --git a/src/domains/rational.rs b/src/domains/rational.rs index aabc698..75f7b92 100644 --- a/src/domains/rational.rs +++ b/src/domains/rational.rs @@ -399,6 +399,61 @@ impl Rational { } } + /// Round the rational to the one with the smallest numerator or denominator in the interval + /// `[self * (1-relative_error), self * (1+relative_error)]`, where + /// `0 < relative_error < 1`. + pub fn round(&self, relative_error: &Rational) -> Rational { + if self.is_negative() { + self.round_in_interval( + self.clone() * (Rational::one() + relative_error), + self.clone() * (Rational::one() - relative_error), + ) + } else { + self.round_in_interval( + self.clone() * (Rational::one() - relative_error), + self.clone() * (Rational::one() + relative_error), + ) + } + } + + /// Round the rational to the one with the smallest numerator or denominator in the interval + /// `[l, u]`, where `l < u`. + pub fn round_in_interval(&self, mut l: Rational, mut u: Rational) -> Rational { + assert!(l < u); + + let mut flip = false; + if l.is_negative() && u.is_negative() { + flip = true; + (l, u) = (-u, -l); + } else if l.is_negative() { + return Rational::zero(); + } + + // use continued fractions to find the best approximation in an interval + let (mut ln, mut ld) = (l.numerator(), l.denominator()); + let (mut un, mut ud) = (u.numerator(), u.denominator()); + + // h1/k1 accumulates the shared continued fraction terms of l and u + let (mut h1, mut h0, mut k1, mut k0): (Integer, Integer, Integer, Integer) = + (1.into(), 0.into(), 0.into(), 1.into()); + + loop { + let a = &(&ln - &1.into()) / &ld; // get next term in continued fraction + (ld, ud, ln, un) = (&un - &a * &ud, &ln - &a * &ld, ud, ld); // subtract and invert + (h1, h0) = (&a * &h1 + &h0, h1); + (k1, k0) = (&a * &k1 + &k0, k1); + if ln <= ld { + let res: Rational = (h1 + &h0, k1 + &k0).into(); + + if flip { + return -res; + } else { + return res; + } + } + } + } + /// Reconstruct a rational number `q` from a value `v` in a prime field `p`, /// such that `q ≡ v mod p`. /// @@ -968,3 +1023,35 @@ impl<'a> std::iter::Sum<&'a Self> for Rational { iter.fold(Rational::zero(), |a, b| a + b) } } + +#[cfg(test)] +mod test { + use crate::domains::rational::Rational; + + #[test] + fn rounding() { + let r = Rational::new(11, 10); + let res = r.round_in_interval(Rational::new(1, 1), Rational::new(12, 10)); + assert_eq!(res, Rational::new(1, 1)); + + let r = Rational::new(11, 10); + let res = r.round_in_interval(Rational::new(2, 1), Rational::new(3, 1)); + assert_eq!(res, Rational::new(2, 1)); + + let r = Rational::new(503, 1500); + let res = r.round(&Rational::new(1, 10)); + assert_eq!(res, Rational::new(1, 3)); + + let r = Rational::new(-503, 1500); + let res = r.round(&Rational::new(1, 10)); + assert_eq!(res, Rational::new(-1, 3)); + + let r = crate::domains::float::Float::from(rug::Float::with_val( + 1000, + rug::float::Constant::Pi, + )) + .to_rational(); + let res = r.round(&Rational::new(1, 100000000)); + assert_eq!(res, (93343, 29712).into()); + } +} diff --git a/src/evaluate.rs b/src/evaluate.rs index 4617f5c..65947cb 100644 --- a/src/evaluate.rs +++ b/src/evaluate.rs @@ -224,7 +224,7 @@ mod test { assert_eq!( format!("{}", r), - "6.0000000099840062521194578624390895167558285149387196915810785" + "6.00000000998400625211945786243908951675582851493871969158108" ); } } diff --git a/src/parser.rs b/src/parser.rs index 11c5be7..7efe9bb 100644 --- a/src/parser.rs +++ b/src/parser.rs @@ -1330,10 +1330,7 @@ mod test { "{}", AtomPrinter::new_with_options(input.as_view(), PrintOptions::file()) ); - assert_eq!( - r, - "1.200000000000000000003*x+2*exp(5)+1.12340000000000002e28" - ); + assert_eq!(r, "1.2000000000000000000*x+2*exp(5)+1.123400000000000e28"); } #[test] diff --git a/symbolica.pyi b/symbolica.pyi index 25e30cc..29882d3 100644 --- a/symbolica.pyi +++ b/symbolica.pyi @@ -193,10 +193,10 @@ class Expression: """ @classmethod - def num(_cls, num: int | float | str | Decimal, max_denom: Optional[int] = None) -> Expression: + def num(_cls, num: int | float | str | Decimal, relative_error: Optional[float] = None) -> Expression: """Create a new Symbolica number from an int, a float, or a string. A floating point number is kept as a float with the same precision as the input, - but it can also be truncated by specifying the maximal denominator value. + but it can also be converted to the smallest rational number given a `relative_error`. Examples -------- @@ -205,9 +205,13 @@ class Expression: 1/2 >>> print(Expression.num(1/3)) - >>> print(Expression.num(0.33, 5)) + >>> print(Expression.num(0.33, 0.1)) + >>> print(Expression.num('0.333`3')) + >>> print(Expression.num(Decimal('0.1234'))) 3.3333333333333331e-1 1/3 + 3.33e-1 + 1.2340e-1 """ @classmethod @@ -382,11 +386,13 @@ class Expression: transformations can be applied. """ - def to_float(self, decimal_prec: int) -> Expression: - """Convert all coefficients to floats, with a given decimal precision.""" + def coefficients_to_float(self, decimal_prec: int) -> Expression: + """Convert all coefficients to floats with a given precision `decimal_prec``. + The precision of floating point coefficients in the input will be truncated to `decimal_prec`.""" - def float_to_rat(self, max_denominator: int) -> Expression: - """Convert all floating point coefficients to rationals, with a given maximal denominator.""" + def rationalize_coefficients(self, relative_error: float) -> Expression: + """Map all floating point and rational coefficients to the best rational approximation + in the interval `[self*(1-relative_error),self*(1+relative_error)]`.""" def req_len(self, min_length: int, max_length: int | None) -> PatternRestriction: """ @@ -1305,7 +1311,7 @@ class Transformer: def series( self, x: Expression, - expansion_point: Expression, + expansion_point: Expression | int | float | Decimal, depth: int, depth_denom: int = 1, depth_is_absolute: bool = True