From c39dca736ef6379a0ddca272445b5a14a19fc9a7 Mon Sep 17 00:00:00 2001 From: Ben Ruijl Date: Thu, 6 Jun 2024 10:48:30 +0200 Subject: [PATCH] Add support for floating point coefficients --- src/atom/coefficient.rs | 33 +++++++++++- src/coefficient.rs | 115 ++++++++++++++++++++++++++++++++++++---- src/evaluate.rs | 3 ++ src/parser.rs | 35 +++++++++--- src/poly.rs | 3 ++ src/printer.rs | 5 ++ 6 files changed, 174 insertions(+), 20 deletions(-) diff --git a/src/atom/coefficient.rs b/src/atom/coefficient.rs index 1c54b9a2..a3e89388 100644 --- a/src/atom/coefficient.rs +++ b/src/atom/coefficient.rs @@ -1,8 +1,13 @@ +use std::io::Write; + use bytes::{Buf, BufMut}; use rug::{integer::Order, Integer as MultiPrecisionInteger}; use crate::{ - coefficient::{Coefficient, CoefficientView, SerializedRational, SerializedRationalPolynomial}, + coefficient::{ + Coefficient, CoefficientView, SerializedFloat, SerializedRational, + SerializedRationalPolynomial, + }, domains::{ finite_field::FiniteFieldElement, integer::{Integer, IntegerRing, Z}, @@ -19,6 +24,7 @@ const U64_NUM: u8 = 0b00000100; const FIN_NUM: u8 = 0b00000101; const ARB_NUM: u8 = 0b00000111; const RAT_POLY: u8 = 0b00001000; +const FLOAT: u8 = 0b00001001; const U8_DEN: u8 = 0b00010000; const U16_DEN: u8 = 0b00100000; const U32_DEN: u8 = 0b00110000; @@ -152,6 +158,15 @@ impl PackedRationalNumberWriter for Coefficient { .write_digits(&mut dest[old_len + num_digits..], Order::Lsf); } }, + Coefficient::Float(f) => { + dest.put_u8(FLOAT); + + // TODO: improve serialization + let s = f.to_string_radix(16, None); + dest.put_u64_le(s.len() as u64 + 4); + dest.put_u32_le(f.prec()); + dest.write_all(s.as_bytes()).unwrap(); + } Coefficient::FiniteField(num, f) => { dest.put_u8(FIN_NUM); (num.0, f.0 as u64).write_packed(dest); // this adds an extra tag @@ -216,6 +231,7 @@ impl PackedRationalNumberWriter for Coefficient { Rational::Natural(num, den) => (*num, *den).write_packed_fixed(dest), Rational::Large(_) => todo!("Writing large packed rational not implemented"), }, + Coefficient::Float(_) => todo!("Writing large packed rational not implemented"), Coefficient::RationalPolynomial(_) => { todo!("Writing packed rational polynomial not implemented") } @@ -236,6 +252,10 @@ impl PackedRationalNumberWriter for Coefficient { 1 + (n, d).get_packed_size() + n as u64 + d as u64 } }, + Coefficient::Float(f) => { + let s = f.to_string_radix(16, None); + 1 + 8 + 4 + s.len() as u64 + } Coefficient::FiniteField(m, i) => 2 + (m.0, i.0 as u64).get_packed_size(), Coefficient::RationalPolynomial(_) => { unimplemented!("Cannot get the packed size of a rational polynomial") @@ -267,6 +287,14 @@ impl PackedRationalNumberReader for [u8] { CoefficientView::RationalPolynomial(SerializedRationalPolynomial(&start[..len])), source, ) + } else if disc == FLOAT { + let len = source.get_u64_le() as usize; + let start = source; + source.advance(len); + ( + CoefficientView::Float(SerializedFloat(&start[..len])), + source, + ) } else if (disc & NUM_MASK) == ARB_NUM { let (num, den); (num, den, source) = source.get_frac_i64(); @@ -519,6 +547,9 @@ impl PackedRationalNumberReader for [u8] { let size = get_size_of_natural(var_size & NUM_MASK) + get_size_of_natural((var_size & DEN_MASK) >> 4); dest.advance(size as usize); + } else if v_num == FLOAT { + let size = dest.get_u64_le() as usize; + dest.advance(size); } else { unreachable!("Unsupported numerator/denominator type {}", disc) } diff --git a/src/coefficient.rs b/src/coefficient.rs index fefd578b..9fba675a 100644 --- a/src/coefficient.rs +++ b/src/coefficient.rs @@ -5,9 +5,10 @@ use std::{ }; use ahash::HashMap; +use bytes::Buf; use rug::{ integer::Order, - ops::{NegAssign, Pow as RPow}, + ops::{CompleteRound, NegAssign, Pow as RPow}, Integer as MultiPrecisionInteger, Rational as MultiPrecisionRational, }; use smallvec::{smallvec, SmallVec}; @@ -41,13 +42,16 @@ pub trait ConvertToRing: Ring { /// A coefficient that can appear in a Symbolica expression. /// In most cases, this is a rational number but it can also be a finite field element or /// a rational polynomial. -#[derive(Debug, Clone, PartialEq, Eq)] +#[derive(Debug, Clone, PartialEq)] pub enum Coefficient { Rational(Rational), + Float(rug::Float), FiniteField(FiniteFieldElement, FiniteFieldIndex), RationalPolynomial(RationalPolynomial), } +impl Eq for Coefficient {} + impl From for Coefficient { fn from(value: i64) -> Self { Coefficient::Rational(value.into()) @@ -119,6 +123,7 @@ impl Coefficient { pub fn is_zero(&self) -> bool { match self { Coefficient::Rational(r) => r.is_zero(), + Coefficient::Float(f) => f.is_zero(), Coefficient::FiniteField(num, _field) => num.0 == 0, Coefficient::RationalPolynomial(r) => r.numerator.is_zero(), } @@ -144,10 +149,7 @@ impl Add for Coefficient { let f = State::get_finite_field(i1); Coefficient::FiniteField(f.add(&n1, &n2), i1) } - (Coefficient::FiniteField(_, _), _) => { - panic!("Cannot add finite field to non-finite number. Convert other number first?"); - } - (_, Coefficient::FiniteField(_, _)) => { + (Coefficient::FiniteField(_, _), _) | (_, Coefficient::FiniteField(_, _)) => { panic!("Cannot add finite field to non-finite number. Convert other number first?"); } (Coefficient::Rational(r), Coefficient::RationalPolynomial(rp)) @@ -171,6 +173,14 @@ impl Add for Coefficient { Coefficient::RationalPolynomial(r) } } + (Coefficient::Rational(r), Coefficient::Float(f)) + | (Coefficient::Float(f), Coefficient::Rational(r)) => { + Coefficient::Float(r.to_multi_prec_float(f.prec()) + f) + } + (Coefficient::Float(f1), Coefficient::Float(f2)) => Coefficient::Float(f1 + f2), + (Coefficient::Float(_), _) | (_, Coefficient::Float(_)) => { + panic!("Cannot add float to finite-field number or rational polynomial"); + } } } } @@ -194,10 +204,7 @@ impl Mul for Coefficient { let f = State::get_finite_field(i1); Coefficient::FiniteField(f.mul(&n1, &n2), i1) } - (Coefficient::FiniteField(_, _), _) => { - panic!("Cannot multiply finite field to non-finite number. Convert other number first?"); - } - (_, Coefficient::FiniteField(_, _)) => { + (Coefficient::FiniteField(_, _), _) | (_, Coefficient::FiniteField(_, _)) => { panic!("Cannot multiply finite field to non-finite number. Convert other number first?"); } (Coefficient::Rational(r), Coefficient::RationalPolynomial(mut rp)) @@ -227,6 +234,14 @@ impl Mul for Coefficient { Coefficient::RationalPolynomial(r) } } + (Coefficient::Rational(r), Coefficient::Float(f)) + | (Coefficient::Float(f), Coefficient::Rational(r)) => { + Coefficient::Float(r.to_multi_prec_float(f.prec()) * f) + } + (Coefficient::Float(f1), Coefficient::Float(f2)) => Coefficient::Float(f1 * f2), + (Coefficient::Float(_), _) | (_, Coefficient::Float(_)) => { + panic!("Cannot add float to finite-field number or rational polynomial"); + } } } } @@ -257,10 +272,22 @@ impl<'a> SerializedRational<'a> { #[derive(Debug, Copy, Clone, PartialEq, Eq)] pub struct SerializedRationalPolynomial<'a>(pub &'a [u8]); +#[derive(Debug, Copy, Clone, PartialEq, Eq)] +pub struct SerializedFloat<'a>(pub &'a [u8]); + +impl<'a> SerializedFloat<'a> { + pub fn to_float(&self) -> rug::Float { + let mut d = self.0; + let prec = d.get_u32_le(); + rug::Float::parse_radix(&d, 16).unwrap().complete(prec) + } +} + /// A view of a coefficient that keeps GMP rationals serialized. #[derive(Debug, Copy, Clone, PartialEq, Eq)] pub enum CoefficientView<'a> { Natural(i64, i64), + Float(SerializedFloat<'a>), Large(SerializedRational<'a>), FiniteField(FiniteFieldElement, FiniteFieldIndex), RationalPolynomial(SerializedRationalPolynomial<'a>), @@ -276,6 +303,7 @@ impl ConvertToRing for RationalField { fn element_from_coefficient(&self, number: Coefficient) -> Self::Element { match number { Coefficient::Rational(r) => r, + Coefficient::Float(_) => panic!("Cannot convert float to rational"), Coefficient::FiniteField(_, _) => panic!("Cannot convert finite field to rational"), Coefficient::RationalPolynomial(_) => { panic!("Cannot convert rational polynomial to rational") @@ -288,6 +316,9 @@ impl ConvertToRing for RationalField { match number { CoefficientView::Natural(r, d) => Rational::Natural(r, d), CoefficientView::Large(r) => Rational::Large(r.to_rat()), + CoefficientView::Float(_) => { + panic!("Cannot convert float to rational") + } CoefficientView::FiniteField(_, _) => { panic!("Cannot convert finite field to rational") } @@ -311,6 +342,7 @@ impl ConvertToRing for IntegerRing { debug_assert!(r.is_integer()); r.numerator() } + Coefficient::Float(_) => panic!("Cannot convert float to integer"), Coefficient::FiniteField(_, _) => panic!("Cannot convert finite field to integer"), Coefficient::RationalPolynomial(_) => { panic!("Cannot convert rational polynomial to rational") @@ -330,6 +362,9 @@ impl ConvertToRing for IntegerRing { debug_assert!(r.denom() == &1); Integer::from_large(r.numer().clone()) } + CoefficientView::Float(_) => { + panic!("Cannot convert float to integer") + } CoefficientView::FiniteField(_, _) => { panic!("Cannot convert finite field to integer") } @@ -360,6 +395,7 @@ where &r.numerator().to_finite_field(self), &r.denominator().to_finite_field(self), ), + Coefficient::Float(_) => panic!("Cannot convert float to finite field"), Coefficient::FiniteField(_, _) => panic!("Cannot convert finite field to other one"), Coefficient::RationalPolynomial(_) => { panic!("Cannot convert rational polynomial to finite field") @@ -384,6 +420,9 @@ where &Integer::Large(d).to_finite_field(self), ) } + CoefficientView::Float(_) => { + panic!("Cannot convert float to finite field") + } CoefficientView::FiniteField(_, _) => { panic!("Cannot convert finite field to other one") } @@ -401,7 +440,8 @@ impl CoefficientView<'_> { Rational::Natural(n, d) => Coefficient::Rational((n, d).into()), Rational::Large(l) => Coefficient::Rational(l.into()), }, - CoefficientView::Large(_) + CoefficientView::Float(_) + | CoefficientView::Large(_) | CoefficientView::FiniteField(_, _) | CoefficientView::RationalPolynomial(_) => self.to_owned(), } @@ -411,6 +451,7 @@ impl CoefficientView<'_> { match self { CoefficientView::Natural(num, den) => Coefficient::Rational((*num, *den).into()), CoefficientView::Large(r) => Coefficient::Rational(r.to_rat().into()), + CoefficientView::Float(f) => Coefficient::Float(f.to_float().into()), CoefficientView::FiniteField(num, field) => Coefficient::FiniteField(*num, *field), CoefficientView::RationalPolynomial(p) => { Coefficient::RationalPolynomial(p.deserialize()) @@ -490,6 +531,7 @@ impl CoefficientView<'_> { pub fn is_integer(&self) -> bool { match self { CoefficientView::Natural(_, d) => *d == 1, + CoefficientView::Float(_) => false, CoefficientView::Large(r) => r.to_rat().is_integer(), CoefficientView::FiniteField(_, _) => true, CoefficientView::RationalPolynomial(_) => false, @@ -611,6 +653,26 @@ impl Add> for CoefficientView<'_> { Coefficient::RationalPolynomial(r) } } + (CoefficientView::Natural(n, d), CoefficientView::Float(f)) + | (CoefficientView::Float(f), CoefficientView::Natural(n, d)) => { + let f = f.to_float(); + Coefficient::Float(Rational::from((n, d)).to_multi_prec_float(f.prec()) + f) + } + (CoefficientView::Large(r), CoefficientView::Float(f)) + | (CoefficientView::Float(f), CoefficientView::Large(r)) => { + let r = r.to_rat(); + let f = f.to_float(); + Coefficient::Float(Rational::from_large(r).to_multi_prec_float(f.prec()) + f) + } + (CoefficientView::Float(f1), CoefficientView::Float(f2)) => { + Coefficient::Float(f1.to_float() + f2.to_float()) + } + (CoefficientView::Float(_), CoefficientView::RationalPolynomial(_)) => { + panic!("Cannot add float to rational polynomial"); + } + (CoefficientView::RationalPolynomial(_), CoefficientView::Float(_)) => { + panic!("Cannot add float to rational polynomial"); + } } } } @@ -686,6 +748,26 @@ impl Mul for CoefficientView<'_> { Coefficient::RationalPolynomial(r) } } + (CoefficientView::Natural(n, d), CoefficientView::Float(f)) + | (CoefficientView::Float(f), CoefficientView::Natural(n, d)) => { + let f = f.to_float(); + Coefficient::Float(Rational::from((n, d)).to_multi_prec_float(f.prec()) * f) + } + (CoefficientView::Large(r), CoefficientView::Float(f)) + | (CoefficientView::Float(f), CoefficientView::Large(r)) => { + let r = r.to_rat(); + let f = f.to_float(); + Coefficient::Float(Rational::from_large(r).to_multi_prec_float(f.prec()) * f) + } + (CoefficientView::Float(f1), CoefficientView::Float(f2)) => { + Coefficient::Float(f1.to_float() * f2.to_float()) + } + (CoefficientView::Float(_), CoefficientView::RationalPolynomial(_)) => { + panic!("Cannot multiply float to rational polynomial"); + } + (CoefficientView::RationalPolynomial(_), CoefficientView::Float(_)) => { + panic!("Cannot multiply float to rational polynomial"); + } } } } @@ -698,6 +780,7 @@ impl Add for CoefficientView<'_> { CoefficientView::Natural(n1, d1) => { Coefficient::Rational(Rational::Natural(n1, d1) + &other.into()) } + CoefficientView::Float(f) => Coefficient::Float(f.to_float() + other), CoefficientView::Large(r1) => (r1.to_rat() + other).into(), CoefficientView::FiniteField(n1, i1) => { let f = State::get_finite_field(i1); @@ -902,6 +985,8 @@ mod test { use crate::{atom::Atom, domains::rational::Rational, state::State}; + use super::Coefficient; + #[test] fn coeff_conversion() { let expr = Atom::parse("v1*coeff(v2+v3/v4)+v1*coeff(v2)").unwrap(); @@ -951,4 +1036,12 @@ mod test { assert_eq!(a, expr); } + + #[test] + fn float() { + let expr = Atom::parse("1/2 x + 5.8912734891723").unwrap(); + let c = Coefficient::Float(rug::Float::with_val(200, rug::float::Constant::Pi)); + let expr = expr * &Atom::new_num(c); + println!("{}", expr.expand()); + } } diff --git a/src/evaluate.rs b/src/evaluate.rs index 62eaa17c..27bd0760 100644 --- a/src/evaluate.rs +++ b/src/evaluate.rs @@ -68,6 +68,9 @@ impl<'a> AtomView<'a> { AtomView::Num(n) => match n.get_coeff_view() { CoefficientView::Natural(n, d) => coeff_map(&Rational::Natural(n, d)), CoefficientView::Large(l) => coeff_map(&Rational::Large(l.to_rat())), + CoefficientView::Float(_) => { + unimplemented!("Float not yet supported for evaluation") + } CoefficientView::FiniteField(_, _) => { unimplemented!("Finite field not yet supported for evaluation") } diff --git a/src/parser.rs b/src/parser.rs index 1bbb5650..cd5b1ee8 100644 --- a/src/parser.rs +++ b/src/parser.rs @@ -1,14 +1,14 @@ -use std::{fmt::Write, string::String, sync::Arc}; +use std::{f64::consts::LOG2_10, fmt::Write, string::String, sync::Arc}; use bytes::Buf; -use rug::Integer as MultiPrecisionInteger; +use rug::{ops::CompleteRound, Integer as MultiPrecisionInteger}; use smallvec::SmallVec; use smartstring::{LazyCompact, SmartString}; use crate::{ atom::Atom, - coefficient::ConvertToRing, + coefficient::{Coefficient, ConvertToRing}, domains::{integer::Integer, Ring}, poly::{polynomial::MultivariatePolynomial, Exponent, Variable}, state::{State, Workspace}, @@ -384,7 +384,15 @@ impl Token { Ok(x) => { out.to_num(x.into()); } - Err(e) => return Err(format!("Could not parse number: {}", e)), + Err(_) => match rug::Float::parse(n) { + Ok(f) => { + // derive precision from string length, should be overestimate + out.to_num(Coefficient::Float( + f.complete((n.len() as f64 * LOG2_10).ceil() as u32), + )); + } + Err(e) => Err(format!("Error parsing number: {}", e))?, + }, }, Token::ID(x) => { out.to_var(state.get_symbol_impl(x)); @@ -478,9 +486,20 @@ impl Token { out: &mut Atom, ) -> Result<(), String> { match self { - Token::Number(n) => { - out.to_num(n.parse::()?.into()); - } + Token::Number(n) => match n.parse::() { + Ok(x) => { + out.to_num(x.into()); + } + Err(_) => match rug::Float::parse(n) { + Ok(f) => { + // derive precision from string length, should be overestimate + out.to_num(Coefficient::Float( + f.complete((n.len() as f64 * LOG2_10).ceil() as u32), + )); + } + Err(e) => Err(format!("Error parsing number: {}", e))?, + }, + }, Token::ID(name) => { let index = var_name_map .iter() @@ -661,7 +680,7 @@ impl Token { } } ParseState::Number => { - if c != '_' && c != ' ' && !c.is_ascii_digit() { + if c != '_' && c != ' ' && c != '.' && !c.is_ascii_digit() { state = ParseState::Any; stack.push(Token::Number(id_buffer.as_str().into())); id_buffer.clear(); diff --git a/src/poly.rs b/src/poly.rs index 9b7f7c38..20540f1f 100644 --- a/src/poly.rs +++ b/src/poly.rs @@ -559,6 +559,9 @@ impl<'a> AtomView<'a> { Err("Exponent too large or negative or a fraction") } } + CoefficientView::Float(_) => { + Err("Float is not supported in conversion routine") + } CoefficientView::FiniteField(_, _) => { Err("Finite field not supported in conversion routine") } diff --git a/src/printer.rs b/src/printer.rs index 29633e18..d8815fa5 100644 --- a/src/printer.rs +++ b/src/printer.rs @@ -339,6 +339,10 @@ impl<'a> FormattedPrintNum for NumView<'a> { f.write_fmt(format_args!("{}", num.unsigned_abs())) } } + CoefficientView::Float(fl) => { + let float = fl.to_float(); + f.write_fmt(format_args!("{}", float)) + } CoefficientView::Large(r) => { let rat = r.to_rat().abs(); if !opts.latex @@ -565,6 +569,7 @@ impl<'a> FormattedPrintPow for PowView<'a> { || if let AtomView::Num(n) = b { match n.get_coeff_view() { CoefficientView::Natural(n, d) => n < 0 || d != 1, + CoefficientView::Float(_) => true, // TODO CoefficientView::Large(r) => r.is_negative() || !r.to_rat().is_integer(), CoefficientView::FiniteField(n, i) => { opts.symmetric_representation_for_finite_field