From d06075cc63ffc9b31a72d10b30f4a37a08aadd20 Mon Sep 17 00:00:00 2001 From: Ben Ruijl Date: Mon, 15 Apr 2024 11:45:37 +0200 Subject: [PATCH] Add together and apart for expressions - Improvements to apart function - Fix Hermite reduction - Handle constants better in rational integration algorithm - Fix printing of Mersennse prime field element --- src/api/cpp.rs | 118 ++++++++++++- src/api/python.rs | 55 ++++++ src/collect.rs | 78 ++++++++ src/domains/finite_field.rs | 4 +- src/domains/rational_polynomial.rs | 274 +++++++++++++++++++++-------- src/lib.rs | 21 +++ src/representations.rs | 3 +- symbolica.pyi | 23 +++ 8 files changed, 497 insertions(+), 79 deletions(-) diff --git a/src/api/cpp.rs b/src/api/cpp.rs index f7c577eb..99863630 100644 --- a/src/api/cpp.rs +++ b/src/api/cpp.rs @@ -422,12 +422,16 @@ unsafe extern "C" fn drop(symbolica: *mut Symbolica) { mod test { use std::ffi::{c_char, CStr}; - use super::{drop, init}; + use crate::domains::finite_field::Mersenne64; + + use super::{drop, init, set_options}; #[test] fn simplify() { let symbolica = unsafe { init() }; + unsafe { set_options(symbolica, true, false) }; + unsafe { super::set_vars(symbolica, b"d,y\0".as_ptr() as *const c_char) }; let input = "-(4096-4096*y^2)/(-3072+1024*d)*(1536-512*d)-(-8192+8192*y^2)/(2)*((-6+d)/2)-(-8192+8192*y^2)/(-2)*((-13+3*d)/2)-(-8192+8192*y^2)/(-4)*(-8+2*d)\0"; @@ -436,11 +440,15 @@ mod test { assert_eq!(result, "[32768-32768*y^2-8192*d+8192*d*y^2]"); + unsafe { set_options(symbolica, true, true) }; + let result = unsafe { super::simplify(symbolica, input.as_ptr() as *const i8, 0, false) }; let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); assert_eq!(result, "32768-32768*y^2-8192*d+8192*d*y^2"); + unsafe { set_options(symbolica, false, false) }; + let result = unsafe { super::simplify_factorized(symbolica, input.as_ptr() as *const i8, 0, true) }; let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); @@ -461,32 +469,132 @@ mod test { unsafe { super::set_vars(symbolica, b"d,y\0".as_ptr() as *const c_char) }; + let prime = 4293491017; + let input = "-(4096-4096*y^2)/(-3072+1024*d)*(1536-512*d)-(-8192+8192*y^2)/(2)*((-6+d)/2)-(-8192+8192*y^2)/(-2)*((-13+3*d)/2)-(-8192+8192*y^2)/(-4)*(-8+2*d)\0"; let result = - unsafe { super::simplify(symbolica, input.as_ptr() as *const i8, 4293491017, true) }; + unsafe { super::simplify(symbolica, input.as_ptr() as *const i8, prime, true) }; let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); assert_eq!(result, "[32768+4293458249*y^2+4293482825*d+8192*d*y^2]"); let result = - unsafe { super::simplify(symbolica, input.as_ptr() as *const i8, 4293491017, false) }; + unsafe { super::simplify(symbolica, input.as_ptr() as *const i8, prime, false) }; let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); assert_eq!(result, "32768+4293458249*y^2+4293482825*d+8192*d*y^2"); let result = unsafe { - super::simplify_factorized(symbolica, input.as_ptr() as *const i8, 4293491017, true) + super::simplify_factorized(symbolica, input.as_ptr() as *const i8, prime, true) }; let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); assert_eq!(result, "[32768+4293458249*y^2+4293482825*d+8192*d*y^2]"); let result = unsafe { - super::simplify_factorized(symbolica, input.as_ptr() as *const i8, 4293491017, false) + super::simplify_factorized(symbolica, input.as_ptr() as *const i8, prime, false) }; let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); unsafe { drop(symbolica) }; assert_eq!(result, "32768+4293458249*y^2+4293482825*d+8192*d*y^2"); } + + #[test] + fn simplify_mersenne() { + let symbolica = unsafe { init() }; + + unsafe { super::set_vars(symbolica, b"d,y\0".as_ptr() as *const c_char) }; + + let prime = Mersenne64::PRIME; + + let input = "-(4096-4096*y^2)/(-3072+1024*d)*(1536-512*d)-(-8192+8192*y^2)/(2)*((-6+d)/2)-(-8192+8192*y^2)/(-2)*((-13+3*d)/2)-(-8192+8192*y^2)/(-4)*(-8+2*d)\0"; + let result = + unsafe { super::simplify(symbolica, input.as_ptr() as *const i8, prime, true) }; + let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); + + assert_eq!( + result, + "[32768+2305843009213661183*y^2+2305843009213685759*d+8192*d*y^2]" + ); + + let result = + unsafe { super::simplify(symbolica, input.as_ptr() as *const i8, prime, false) }; + let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); + + assert_eq!( + result, + "32768+2305843009213661183*y^2+2305843009213685759*d+8192*d*y^2" + ); + + let result = unsafe { + super::simplify_factorized(symbolica, input.as_ptr() as *const i8, prime, true) + }; + let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); + + assert_eq!( + result, + "[32768+2305843009213661183*y^2+2305843009213685759*d+8192*d*y^2]" + ); + + let result = unsafe { + super::simplify_factorized(symbolica, input.as_ptr() as *const i8, prime, false) + }; + let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); + + unsafe { drop(symbolica) }; + assert_eq!( + result, + "32768+2305843009213661183*y^2+2305843009213685759*d+8192*d*y^2" + ); + } + + #[test] + fn simplify_u64_prime() { + let symbolica = unsafe { init() }; + + unsafe { super::set_vars(symbolica, b"d,y\0".as_ptr() as *const c_char) }; + + let prime = 18446744073709551163; + + let input = "-(4096-4096*y^2)/(-3072+1024*d)*(1536-512*d)-(-8192+8192*y^2)/(2)*((-6+d)/2)-(-8192+8192*y^2)/(-2)*((-13+3*d)/2)-(-8192+8192*y^2)/(-4)*(-8+2*d)\0"; + let result = + unsafe { super::simplify(symbolica, input.as_ptr() as *const i8, prime, true) }; + let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); + + assert_eq!( + result, + "[32768+18446744073709518395*y^2+18446744073709542971*d+8192*d*y^2]" + ); + + let result = + unsafe { super::simplify(symbolica, input.as_ptr() as *const i8, prime, false) }; + let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); + + assert_eq!( + result, + "32768+18446744073709518395*y^2+18446744073709542971*d+8192*d*y^2" + ); + + let result = unsafe { + super::simplify_factorized(symbolica, input.as_ptr() as *const i8, prime, true) + }; + let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); + + assert_eq!( + result, + "[32768+18446744073709518395*y^2+18446744073709542971*d+8192*d*y^2]" + ); + + let result = unsafe { + super::simplify_factorized(symbolica, input.as_ptr() as *const i8, prime, false) + }; + let result = unsafe { CStr::from_ptr(result).to_str().unwrap() }.to_owned(); + + unsafe { drop(symbolica) }; + assert_eq!( + result, + "32768+18446744073709518395*y^2+18446744073709542971*d+8192*d*y^2" + ); + } } diff --git a/src/api/python.rs b/src/api/python.rs index f8dd6527..1c9fa33f 100644 --- a/src/api/python.rs +++ b/src/api/python.rs @@ -2491,6 +2491,61 @@ impl PythonExpression { Ok(PythonExpression { expr: Arc::new(b) }) } + /// Compute the partial fraction decomposition in `x`. + /// + /// Examples + /// -------- + /// + /// >>> from symbolica import Expression + /// >>> x = Expression.var('x') + /// >>> p = Expression.parse('1/((x+y)*(x^2+x*y+1)(x+1))') + /// >>> print(p.apart(x)) + pub fn apart(&self, x: PythonExpression) -> PyResult { + let poly = self.expr.to_rational_polynomial::<_, _, u32>(&Q, &Z, None); + let x = poly + .get_variables() + .iter() + .position(|v| match (v, x.expr.as_view()) { + (Variable::Symbol(y), AtomView::Var(vv)) => *y == vv.get_symbol(), + (Variable::Function(_, f) | Variable::Other(f), a) => f.as_view() == a, + _ => false, + }) + .ok_or(exceptions::PyValueError::new_err(format!( + "Variable {} not found in polynomial", + x.__str__()? + )))?; + + let fs = poly.apart(x); + + let mut rn = Atom::new(); + Workspace::get_local().with(|ws| { + let mut res = ws.new_atom(); + let a = res.to_add(); + for f in fs { + a.extend(f.to_expression().as_view()); + } + + res.as_view().normalize(ws, &mut rn); + }); + + Ok(PythonExpression { expr: Arc::new(rn) }) + } + + /// Write the expression over a common denominator. + /// + /// Examples + /// -------- + /// + /// >>> from symbolica import Expression + /// >>> p = Expression.parse('v1^2/2+v1^3/v4*v2+v3/(1+v4)') + /// >>> print(p.together()) + pub fn together(&self) -> PyResult { + let poly = self.expr.to_rational_polynomial::<_, _, u32>(&Q, &Z, None); + Ok(PythonExpression { + expr: Arc::new(poly.to_expression()), + }) + } + /// Convert the expression to a polynomial, optionally, with the variables and the ordering specified in `vars`. /// All non-polynomial elements will be converted to new independent variables. pub fn to_polynomial(&self, vars: Option>) -> PyResult { diff --git a/src/collect.rs b/src/collect.rs index 09489627..bef4d5e6 100644 --- a/src/collect.rs +++ b/src/collect.rs @@ -1,6 +1,7 @@ use ahash::HashMap; use crate::{ + domains::{integer::Z, rational::Q}, representations::{Add, AsAtomView, Atom, AtomView, Symbol}, state::Workspace, }; @@ -51,6 +52,16 @@ impl Atom { pub fn coefficient<'a, T: AsAtomView<'a>>(&self, x: T) -> Atom { Workspace::get_local().with(|ws| self.as_view().coefficient_with_ws(x.as_atom_view(), ws)) } + + /// Write the expression over a common denominator. + pub fn together(&self) -> Atom { + self.as_view().together() + } + + /// Write the expression as a sum of terms with minimal denominators. + pub fn apart(&self, x: Symbol) -> Atom { + self.as_view().apart(x) + } } impl<'a> AtomView<'a> { @@ -341,6 +352,49 @@ impl<'a> AtomView<'a> { } } } + + /// Write the expression over a common denominator. + pub fn together(&self) -> Atom { + let mut out = Atom::new(); + self.together_into(&mut out); + out + } + + /// Write the expression over a common denominator. + pub fn together_into(&self, out: &mut Atom) { + self.to_rational_polynomial::<_, _, u32>(&Q, &Z, None) + .to_expression_into(out); + } + + /// Write the expression as a sum of terms with minimal denominators. + pub fn apart(&self, x: Symbol) -> Atom { + let mut out = Atom::new(); + + Workspace::get_local().with(|ws| { + self.apart_with_ws_into(x, ws, &mut out); + }); + + out + } + + /// Write the expression as a sum of terms with minimal denominators. + pub fn apart_with_ws_into(&self, x: Symbol, ws: &Workspace, out: &mut Atom) { + let poly = self.to_rational_polynomial::<_, _, u32>(&Q, &Z, None); + if let Some(v) = poly.get_variables().iter().position(|v| v == &x.into()) { + let mut a = ws.new_atom(); + let add = a.to_add(); + + let mut a = ws.new_atom(); + for x in poly.apart(v) { + x.to_expression_into(&mut a); + add.extend(a.as_view()); + } + + add.as_view().normalize(ws, out); + } else { + out.set_from_view(self); + } + } } #[cfg(test)] @@ -414,4 +468,28 @@ mod test { assert_eq!(out, ref_out); } + + #[test] + fn together() { + let input = Atom::parse("v1^2/2+v1^3/v4*v2+v3/(1+v4)").unwrap(); + let out = input.together(); + + let ref_out = + Atom::parse("(2*v4+2*v4^2)^-1*(2*v3*v4+v1^2*v4+v1^2*v4^2+2*v1^3*v2+2*v1^3*v2*v4)") + .unwrap(); + + assert_eq!(out, ref_out); + } + + #[test] + fn apart() { + let input = + Atom::parse("(2*v4+2*v4^2)^-1*(2*v3*v4+v1^2*v4+v1^2*v4^2+2*v1^3*v2+2*v1^3*v2*v4)") + .unwrap(); + let out = input.apart(State::get_symbol("v4")); + + let ref_out = Atom::parse("1/2*v1^2+v3*(v4+1)^-1+v1^3*v2*v4^-1").unwrap(); + + assert_eq!(out, ref_out); + } } diff --git a/src/domains/finite_field.rs b/src/domains/finite_field.rs index d202f418..aa5b6e81 100644 --- a/src/domains/finite_field.rs +++ b/src/domains/finite_field.rs @@ -685,13 +685,13 @@ impl Mersenne64 { impl std::fmt::Debug for Mersenne64 { fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { - f.write_fmt(format_args!("{}", Self::PRIME)) + std::fmt::Debug::fmt(&self.0, f) } } impl Display for Mersenne64 { fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { - f.write_fmt(format_args!("{}", Self::PRIME)) + self.0.fmt(f) } } diff --git a/src/domains/rational_polynomial.rs b/src/domains/rational_polynomial.rs index b3580f15..e154aa74 100644 --- a/src/domains/rational_polynomial.rs +++ b/src/domains/rational_polynomial.rs @@ -122,6 +122,10 @@ impl RationalPolynomial { self.denominator.unify_variables(&mut other.denominator); } + pub fn is_zero(&self) -> bool { + self.numerator.is_zero() + } + pub fn is_constant(&self) -> bool { self.numerator.is_constant() && self.denominator.is_constant() } @@ -743,72 +747,85 @@ where { /// Compute the partial fraction decomposition of the rational polynomial in `var`. pub fn apart(&self, var: usize) -> Vec { - let fs = self.denominator.factor(); - - if fs.len() == 1 { + if self.denominator.degree(var) == E::zero() { return vec![self.clone()]; } let rat_field = RationalPolynomialField::new_from_poly(&self.numerator); + let n = self + .numerator + .to_univariate(var) + .map_coeff(|c| c.clone().into(), rat_field.clone()); - let mut poly_univ = vec![]; - let mut exp = vec![E::zero(); self.numerator.nvars()]; - for (f, p) in &fs { - let f = f.clone().pow(*p); - - let l = f.to_univariate_polynomial_list(var); - let mut res: MultivariatePolynomial<_, E> = MultivariatePolynomial::new( - &rat_field, - Some(l.len()), - self.numerator.variables.clone().into(), - ); + let mut hs = vec![]; - for (p, e) in l { - exp[var] = e; - let one = p.one(); - res.append_monomial( - RationalPolynomial::from_num_den(p, one, &self.numerator.field, false), - &exp, - ); + let rem = if self.numerator.degree(var) >= self.denominator.degree(var) { + let d = self + .denominator + .to_univariate(var) + .map_coeff(|c| c.clone().into(), rat_field.clone()); + let (q, r) = n.quot_rem(&d); + if !q.is_zero() { + hs.push(Self::from_univariate(q)); } - poly_univ.push(res); - } + r + } else { + n + }; - let rhs = poly_univ[0].one(); - let deltas = MultivariatePolynomial::diophantine_univariate(&mut poly_univ, &rhs); + // partial fraction the denominator + let mut fs = self + .denominator + .factor() + .into_iter() + .map(|(x, p)| (x.to_univariate(var), p)) + .collect::>(); - let mut factors = Vec::with_capacity(deltas.len()); - for (d, (p, pe)) in deltas.into_iter().zip(fs.into_iter()) { - let mut unfold = rat_field.zero(); - for (c, e) in d - .coefficients - .into_iter() - .zip(d.exponents.chunks(self.numerator.nvars())) - { - unfold = &unfold - + &(&c - * &RationalPolynomial::from_num_den( - unfold - .numerator - .monomial(self.numerator.field.one(), e.to_vec()), - unfold.numerator.one(), - &self.numerator.field, - true, - )); + let mut constant = fs[0].0.one(); + fs.retain(|x| { + if x.0.is_constant() { + constant = &constant * &x.0.pow(x.1); + false + } else { + true } + }); - unfold = &unfold - * &RationalPolynomial::from_num_den( - self.numerator.clone(), - p.pow(pe), - &self.numerator.field, - false, - ); + let constant = + Self::from_univariate(constant.map_coeff(|c| c.clone().into(), rat_field.clone())) + .inv(); + + assert!(!fs.is_empty()); - factors.push(unfold.clone()); + let mut expanded = fs + .iter() + .map(|(f, p)| f.pow(*p).map_coeff(|c| c.clone().into(), rat_field.clone())) + .collect::>(); + + // perform partial fractioning + let deltas = if expanded.len() > 1 { + UnivariatePolynomial::diophantine(&mut expanded, &rem) + } else { + vec![rem] + }; + + let fs = fs + .into_iter() + .map(|(x, e)| (x.map_coeff(|c| c.clone().into(), rat_field.clone()), e)) + .collect::>(); + + for (d, (p, p_pow)) in deltas.into_iter().zip(fs) { + let exp = d.p_adic_expansion(&p); + let p_rat = Self::from_univariate(p); + for (pow, d_exp) in exp.into_iter().enumerate() { + hs.push( + &(&Self::from_univariate(d_exp) / &p_rat.pow(p_pow as u64 - pow as u64)) + * &constant, + ); + } } - factors + hs } } @@ -824,35 +841,54 @@ where /// else it is `r*log(a)`. pub fn integrate(&self, var: usize) -> (Vec, Vec<(Self, Self)>) { let rat_field = RationalPolynomialField::new_from_poly(&self.numerator); - let n = self.numerator.to_univariate(var); + let n = self + .numerator + .to_univariate(var) + .map_coeff(|c| c.clone().into(), rat_field.clone()); + + let d = self + .denominator + .to_univariate(var) + .map_coeff(|c| c.clone().into(), rat_field.clone()); - if self.denominator.is_one() { - // map coefficients to a rational polynomial - let n_conv = n.map_coeff(|c| c.clone().into(), rat_field); - let r = Self::from_univariate(n_conv.integrate()); + if d.is_constant() { + let r = &Self::from_univariate(n.integrate()) / &Self::from_univariate(d); return (vec![r], vec![]); } - let d = self.denominator.to_univariate(var); let (q, r) = n.quot_rem(&d); let mut v = if q.is_zero() { vec![] } else { - vec![Self::from_univariate( - q.map_coeff(|c| c.clone().into(), rat_field.clone()) - .integrate(), - )] + let n_conv = q.map_coeff(|c| c.clone().into(), rat_field.clone()); + vec![Self::from_univariate(n_conv.integrate())] }; // partial fraction the denominator - let fs = self + let mut fs = self .denominator .square_free_factorization() .into_iter() .map(|(x, p)| (x.to_univariate(var), p)) .collect::>(); + let mut constant = fs[0].0.one(); + fs.retain(|x| { + if x.0.is_constant() { + constant = &constant * &x.0.pow(x.1); + false + } else { + true + } + }); + + let mut constant = + Self::from_univariate(constant.map_coeff(|c| c.clone().into(), rat_field.clone())) + .inv(); + + assert!(!fs.is_empty()); + let mut expanded = fs .iter() .map(|(f, p)| f.pow(*p).map_coeff(|c| c.clone().into(), rat_field.clone())) @@ -869,17 +905,23 @@ where let fs = fs .into_iter() - .map(|(x, _)| x.map_coeff(|c| c.clone().into(), rat_field.clone())) + .map(|(x, e)| (x.map_coeff(|c| c.clone().into(), rat_field.clone()), e)) .collect::>(); let mut hs = vec![]; - for (d, p) in deltas.into_iter().zip(&fs) { + for (d, (p, p_pow)) in deltas.into_iter().zip(&fs) { let p_diff = p.derivative(); let (_, s, t) = p.eea(&p_diff); let p_full = Self::from_univariate(p.clone()); let mut d_exp = d.p_adic_expansion(&p); + + // grow to the same level as the pow + if d_exp.len() < *p_pow { + d_exp.resize(*p_pow, d.zero()); + } + // highest degree in 1/p last d_exp.reverse(); @@ -895,10 +937,15 @@ where let t_full = Self::from_univariate(t_cor); - v.push(-(&t_full / &(&rat_field.nth(i as u64) * &p_full.pow(i as u64)))); + let r = -(&t_full / &(&rat_field.nth(i as u64) * &p_full.pow(i as u64))); + if !r.is_zero() { + v.push(r); + } } - hs.push(d_exp.swap_remove(0)); + if !d_exp[0].is_zero() { + hs.push(d_exp.swap_remove(0)); + } } // create new temporary variable @@ -912,7 +959,7 @@ where let mut w = vec![]; - for (mut h, mut p) in hs.into_iter().zip(fs.into_iter()) { + for (mut h, (mut p, _)) in hs.into_iter().zip(fs.into_iter()) { for c in &mut p.coefficients { c.unify_variables(&mut t); } @@ -920,6 +967,8 @@ where c.unify_variables(&mut t); } + constant.unify_variables(&mut t); + // adding a variable changes the field p.field = RationalPolynomialField::new_from_poly(&p.coefficients[0].numerator); h.field = p.field.clone(); @@ -1040,9 +1089,9 @@ where res = &res + &(t.coefficient * &mm.into()); } - w.push((sol, res)); + w.push((&sol * &constant, res)); } else { - w.push((ff.clone().into(), res.clone())); + w.push((&ff.clone().into() * &constant, res.clone())); } } } @@ -1061,6 +1110,91 @@ mod test { state::State, }; + #[test] + fn hermite_reduction() { + use crate::representations::Atom; + let p: RationalPolynomial<_, _> = Atom::parse("1/(v1 + 1)^5") + .unwrap() + .to_rational_polynomial::<_, _, u8>(&Q, &Z, None); + + let (r, l) = p.integrate(0); + + assert_eq!( + r, + vec![Atom::parse("-1/(4+16*v1+24*v1^2+16*v1^3+4*v1^4)") + .unwrap() + .to_rational_polynomial::<_, _, u8>(&Q, &Z, r[0].get_variables().clone().into())] + ); + assert_eq!(l, vec![]); + } + + #[test] + fn constant() { + use crate::representations::Atom; + let p: RationalPolynomial<_, _> = Atom::parse("(v1^4+v2+v1*v2+2*v1)/(v2 + 1)") + .unwrap() + .to_rational_polynomial::<_, _, u8>( + &Q, + &Z, + Some(Arc::new(vec![ + State::get_symbol("v1").into(), + State::get_symbol("v2").into(), + ])), + ); + + let (r, l) = p.integrate(0); + assert_eq!( + r, + vec![ + Atom::parse("(10*v1*v2+10*v1^2+5*v1^2*v2+2*v1^5)/(10+10*v2)") + .unwrap() + .to_rational_polynomial::<_, _, u8>( + &Q, + &Z, + r[0].get_variables().clone().into() + ) + ] + ); + assert_eq!(l, vec![]); + } + + #[test] + fn mixed_denominator() { + use crate::representations::Atom; + let p: RationalPolynomial<_, _> = Atom::parse("(v1^4+v2+v1*v2+2*v1)/(v1)/(v2 + 1)") + .unwrap() + .to_rational_polynomial::<_, _, u8>( + &Q, + &Z, + Some(Arc::new(vec![ + State::get_symbol("v1").into(), + State::get_symbol("v2").into(), + ])), + ); + + let (r, l) = p.integrate(0); + + let v = l[0].0.get_variables().clone(); + + assert_eq!( + r, + vec![Atom::parse("(8*v1+4*v1*v2+v1^4)/(4+4*v2)") + .unwrap() + .to_rational_polynomial::<_, _, u8>(&Q, &Z, r[0].get_variables().clone().into())] + ); + assert_eq!( + l, + vec![( + Atom::parse("v2/(1+v2)") + .unwrap() + .to_rational_polynomial::<_, _, u8>(&Q, &Z, v.clone().into()), + Atom::parse("v1") + .unwrap() + .to_rational_polynomial::<_, _, u8>(&Q, &Z, v.clone().into()), + ),] + ); + } + #[test] fn three_factors() { use crate::representations::Atom; diff --git a/src/lib.rs b/src/lib.rs index d3742490..2969d74e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,3 +1,24 @@ +//! Symbolica is a blazing fast computer algebra system. +//! +//! It can be used to perform mathematical operations, +//! such as symbolic differentiation, integration, simplification, +//! pattern matching and solving equations. +//! +//! For example: +//! +//! ``` +//! use symbolica::{representations::Atom, state::State}; +//! +//! fn main() { +//! let input = Atom::parse("x^2*log(2*x + y) + exp(3*x)").unwrap(); +//! let a = input.derivative(State::get_symbol("x")); +//! println!("d({})/dx = {}:", input, a); +//! } +//! ``` +//! +//! Check out the [guide](https://symbolica.io/docs/get_started.html) for more information, examples, +//! and additional documentation. + use std::{ collections::HashMap, env, diff --git a/src/representations.rs b/src/representations.rs index 492445dc..450e4621 100644 --- a/src/representations.rs +++ b/src/representations.rs @@ -568,8 +568,7 @@ impl Atom { /// # state::{FunctionAttribute, State}, /// # }; /// # fn main() { -/// /// -/// let f_id = State::get_or_insert_fn("f", Some(vec![FunctionAttribute::Symmetric])).unwrap(); +/// let f_id = State::get_symbol_with_attributes("f", vec![FunctionAttribute::Symmetric]).unwrap(); /// let fb = FunctionBuilder::new(f_id); /// let a = fb /// .add_arg(&Atom::new_num(3)) diff --git a/symbolica.pyi b/symbolica.pyi index be7f0905..4782fedc 100644 --- a/symbolica.pyi +++ b/symbolica.pyi @@ -697,6 +697,29 @@ class Expression: ) -> Expression: """Taylor expand in `x` around `expansion_point` to depth `depth`.""" + def apart(self, x: Expression) -> Expression: + """Compute the partial fraction decomposition in `x`. + + Examples + -------- + + >>> from symbolica import Expression + >>> x = Expression.var('x') + >>> p = Expression.parse('1/((x+y)*(x^2+x*y+1)(x+1))') + >>> print(p.apart(x)) + """ + + def together(self) -> Expression: + """Write the expression over a common denominator. + + Examples + -------- + + >>> from symbolica import Expression + >>> p = Expression.parse('v1^2/2+v1^3/v4*v2+v3/(1+v4)') + >>> print(p.together()) + """ + def to_polynomial(self, vars: Optional[Sequence[Expression]] = None) -> Polynomial: """Convert the expression to a polynomial, optionally, with the variable ordering specified in `vars`. All non-polynomial parts will be converted to new, independent variables.