Skip to content

Commit

Permalink
Add support for floating point coefficients
Browse files Browse the repository at this point in the history
  • Loading branch information
@benruijl
benruijl committed Jun 6, 2024
1 parent 2f28123 commit c39dca7
Show file tree
Hide file tree
Showing 6 changed files with 174 additions and 20 deletions.
33 changes: 32 additions & 1 deletion src/atom/coefficient.rs
View file
Original file line number Diff line number Diff line change
@@ -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},
Expand All @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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")
}
Expand All @@ -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")
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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)
}
Expand Down
115 changes: 104 additions & 11 deletions src/coefficient.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down Expand Up @@ -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<u64>, FiniteFieldIndex),
RationalPolynomial(RationalPolynomial<IntegerRing, u16>),
}

impl Eq for Coefficient {}

impl From<i64> for Coefficient {
fn from(value: i64) -> Self {
Coefficient::Rational(value.into())
Expand Down Expand Up @@ -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(),
}
Expand All @@ -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))
Expand All @@ -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");
}
}
}
}
Expand All @@ -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))
Expand Down Expand Up @@ -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");
}
}
}
}
Expand Down Expand Up @@ -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<u64>, FiniteFieldIndex),
RationalPolynomial(SerializedRationalPolynomial<'a>),
Expand All @@ -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")
Expand All @@ -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")
}
Expand All @@ -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")
Expand All @@ -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")
}
Expand Down Expand Up @@ -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")
Expand All @@ -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")
}
Expand All @@ -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(),
}
Expand All @@ -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())
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -611,6 +653,26 @@ impl Add<CoefficientView<'_>> 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");
}
}
}
}
Expand Down Expand Up @@ -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");
}
}
}
}
Expand All @@ -698,6 +780,7 @@ impl Add<i64> 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);
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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());
}
}
3 changes: 3 additions & 0 deletions src/evaluate.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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")
}
Expand Down
Loading

0 comments on commit c39dca7

Please sign in to comment.