Polynomials
Polynomials are an important sub-type of general mathematical expression and lie at the heart of many mathematical problems. Since polynomials have the constrained form, Symbolica has its own datatype for multivariate polynomials. This data structure allows for state-of-the-art polynomial arithmetic.
Construction
Polynomials can be efficiently constructed directly from a string if the variables are provided
Python
Rustfrom symbolica import Polynomial
p = Polynomial.parse('x**2*y + 1 + 2*x + y**4', ['x', 'y'])
print(p)TODOor can be converted from a Symbolica expression.
from symbolica import Expression
x, y = Expression.vars('x', 'y')
e = x**2*y + 1 + 2*x + y**4
p = e.to_polynomial()
print(p)TODOIn both cases, the input expression already needs to be in an expanded polynomial form, so an expansion may have to be performed first by the user.
Rational polynomials
Rational polynomials can also be efficiently constructed directly from a string:
from symbolica import RationalPolynomial
p = RationalPolynomial.parse('1/x+(1+x)^2/(1+x+y)', ['x', 'y'])
print(p)TODOor can be converted from a Symbolica expression.
from symbolica import Expression
x, y = Expression.vars('x', 'y')
e = 1/x+(1+x)^2/(1+x+y)
p = e.to_rational_polynomial()
print(p)TODOContrary to polynomials conversion, all rational polynomial conversion routines attempt to do expansions internally.
Both polynomials and rational polynomials have variants that limit the power of the variables, which increases performance.
Symbolica notation
The fastest way to parse a rational polynomial from a string is to provide it in Symbolica notation, which is the following format:
[numerator,denominator]with square brackets. The greatest common divisor of the numerator and the denominator must be 1. Each term in the two polynomials must be in one the following format:
- The product operator is
* - The coefficient must be placed first. If the coefficient is one, it can be omitted.
- The variable with exponent is given as
x^n. If the power is one, it can be omitted
For example:
[2*x+x^2+x^3*y^2,1+x+y]The terms do not have to be sorted in Symbolica ordering, but the such a conversion will impact performance.
Evaluation optimization
Evaluating a polynomial is a commonplace operation in computer algebra. Naively evaluating a polynomial is prohibitively slow for large expressions, especially when it has to be evaluated many times, such as for Monte Carlo simulations/integrations. Symbolica provides state-of-the-art ways to write the polynomial in a form that makes it very fast to evaluate, often 10-100 times faster than a naive evaluation.
It tries to find the best multivariate Horner scheme, performs common subexpression elimination and recycles intermediate variables. It comes close to what a programming language compiler would do.
In the example below, the (small) expression SIGMA is optimized for evaluation. Then, a C++ program is written that can evaluate the polynomial efficiently (with templated floating point type). Additionally, the expression can be evaluated by Symbolica itself, using a floating point type of choice. IN the example below, the polynomial is evaluated with the f64 type and the SIMD type f64x4. Evaluation with f64 is about 4.5 times slower than the C++ code compiled with O3.
use symbolica::{
poly::evaluate::{BorrowedHornerScheme, InstructionSetPrinter},
representations::Atom,
rings::rational::RationalField,
state::{State, Workspace},
};
use symbolica::poly::polynomial::MultivariatePolynomial;
use wide::f64x4;
// a small for this example
const SIGMA: &str = "-a5^3*b0^5+a4*a5^2*b0^4*b1-a4^2*a5*b0^4*b2+a4^3*b0^4*b3-a3*a5^2*
b0^3*b1^2+2*a3*a5^2*b0^4*b2+a3*a4*a5*b0^3*b1*b2-3*a3*a4*a5*b0^4*
b3-a3*a4^2*b0^3*b1*b3-a3^2*a5*b0^3*b2^2+2*a3^2*a5*b0^3*b1*b3+a3^2
*a4*b0^3*b2*b3-a3^3*b0^3*b3^2+a2*a5^2*b0^2*b1^3-3*a2*a5^2*b0^3*b1
*b2+3*a2*a5^2*b0^4*b3-a2*a4*a5*b0^2*b1^2*b2+2*a2*a4*a5*b0^3*b2^2+
a2*a4*a5*b0^3*b1*b3+a2*a4^2*b0^2*b1^2*b3-2*a2*a4^2*b0^3*b2*b3+a2*
a3*a5*b0^2*b1*b2^2-2*a2*a3*a5*b0^2*b1^2*b3-a2*a3*a5*b0^3*b2*b3-a2
*a3*a4*b0^2*b1*b2*b3+3*a2*a3*a4*b0^3*b3^2+a2*a3^2*b0^2*b1*b3^2-
a2^2*a5*b0^2*b2^3+3*a2^2*a5*b0^2*b1*b2*b3-3*a2^2*a5*b0^3*b3^2+
a2^2*a4*b0^2*b2^2*b3-2*a2^2*a4*b0^2*b1*b3^2-a2^2*a3*b0^2*b2*b3^2+
a2^3*b0^2*b3^3-a1*a5^2*b0*b1^4+4*a1*a5^2*b0^2*b1^2*b2-2*a1*a5^2*
b0^3*b2^2-4*a1*a5^2*b0^3*b1*b3+a1*a4*a5*b0*b1^3*b2-3*a1*a4*a5*
b0^2*b1*b2^2-a1*a4*a5*b0^2*b1^2*b3+5*a1*a4*a5*b0^3*b2*b3-a1*a4^2*
b0*b1^3*b3+3*a1*a4^2*b0^2*b1*b2*b3-3*a1*a4^2*b0^3*b3^2-a1*a3*a5*
b0*b1^2*b2^2+2*a1*a3*a5*b0*b1^3*b3+2*a1*a3*a5*b0^2*b2^3-4*a1*a3*
a5*b0^2*b1*b2*b3+3*a1*a3*a5*b0^3*b3^2+a1*a3*a4*b0*b1^2*b2*b3-2*a1
*a3*a4*b0^2*b2^2*b3-a1*a3*a4*b0^2*b1*b3^2-a1*a3^2*b0*b1^2*b3^2+2*
a1*a3^2*b0^2*b2*b3^2+a1*a2*a5*b0*b1*b2^3-3*a1*a2*a5*b0*b1^2*b2*b3
-a1*a2*a5*b0^2*b2^2*b3+5*a1*a2*a5*b0^2*b1*b3^2-a1*a2*a4*b0*b1*
b2^2*b3+2*a1*a2*a4*b0*b1^2*b3^2+a1*a2*a4*b0^2*b2*b3^2+a1*a2*a3*b0
*b1*b2*b3^2-3*a1*a2*a3*b0^2*b3^3-a1*a2^2*b0*b1*b3^3-a1^2*a5*b0*
b2^4+4*a1^2*a5*b0*b1*b2^2*b3-2*a1^2*a5*b0*b1^2*b3^2-4*a1^2*a5*
b0^2*b2*b3^2+a1^2*a4*b0*b2^3*b3-3*a1^2*a4*b0*b1*b2*b3^2+3*a1^2*a4
*b0^2*b3^3-a1^2*a3*b0*b2^2*b3^2+2*a1^2*a3*b0*b1*b3^3+a1^2*a2*b0*
b2*b3^3-a1^3*b0*b3^4+a0*a5^2*b1^5-5*a0*a5^2*b0*b1^3*b2+5*a0*a5^2*
b0^2*b1*b2^2+5*a0*a5^2*b0^2*b1^2*b3-5*a0*a5^2*b0^3*b2*b3-a0*a4*a5
*b1^4*b2+4*a0*a4*a5*b0*b1^2*b2^2+a0*a4*a5*b0*b1^3*b3-2*a0*a4*a5*
b0^2*b2^3-7*a0*a4*a5*b0^2*b1*b2*b3+3*a0*a4*a5*b0^3*b3^2+a0*a4^2*
b1^4*b3-4*a0*a4^2*b0*b1^2*b2*b3+2*a0*a4^2*b0^2*b2^2*b3+4*a0*a4^2*
b0^2*b1*b3^2+a0*a3*a5*b1^3*b2^2-2*a0*a3*a5*b1^4*b3-3*a0*a3*a5*b0*
b1*b2^3+6*a0*a3*a5*b0*b1^2*b2*b3+3*a0*a3*a5*b0^2*b2^2*b3-7*a0*a3*
a5*b0^2*b1*b3^2-a0*a3*a4*b1^3*b2*b3+3*a0*a3*a4*b0*b1*b2^2*b3+a0*
a3*a4*b0*b1^2*b3^2-5*a0*a3*a4*b0^2*b2*b3^2+a0*a3^2*b1^3*b3^2-3*a0
*a3^2*b0*b1*b2*b3^2+3*a0*a3^2*b0^2*b3^3-a0*a2*a5*b1^2*b2^3+3*a0*
a2*a5*b1^3*b2*b3+2*a0*a2*a5*b0*b2^4-6*a0*a2*a5*b0*b1*b2^2*b3-3*a0
*a2*a5*b0*b1^2*b3^2+7*a0*a2*a5*b0^2*b2*b3^2+a0*a2*a4*b1^2*b2^2*b3
-2*a0*a2*a4*b1^3*b3^2-2*a0*a2*a4*b0*b2^3*b3+4*a0*a2*a4*b0*b1*b2*
b3^2-3*a0*a2*a4*b0^2*b3^3-a0*a2*a3*b1^2*b2*b3^2+2*a0*a2*a3*b0*
b2^2*b3^2+a0*a2*a3*b0*b1*b3^3+a0*a2^2*b1^2*b3^3-2*a0*a2^2*b0*b2*
b3^3+a0*a1*a5*b1*b2^4-4*a0*a1*a5*b1^2*b2^2*b3+2*a0*a1*a5*b1^3*
b3^2-a0*a1*a5*b0*b2^3*b3+7*a0*a1*a5*b0*b1*b2*b3^2-3*a0*a1*a5*b0^2
*b3^3-a0*a1*a4*b1*b2^3*b3+3*a0*a1*a4*b1^2*b2*b3^2+a0*a1*a4*b0*
b2^2*b3^2-5*a0*a1*a4*b0*b1*b3^3+a0*a1*a3*b1*b2^2*b3^2-2*a0*a1*a3*
b1^2*b3^3-a0*a1*a3*b0*b2*b3^3-a0*a1*a2*b1*b2*b3^3+3*a0*a1*a2*b0*
b3^4+a0*a1^2*b1*b3^4-a0^2*a5*b2^5+5*a0^2*a5*b1*b2^3*b3-5*a0^2*a5*
b1^2*b2*b3^2-5*a0^2*a5*b0*b2^2*b3^2+5*a0^2*a5*b0*b1*b3^3+a0^2*a4*
b2^4*b3-4*a0^2*a4*b1*b2^2*b3^2+2*a0^2*a4*b1^2*b3^3+4*a0^2*a4*b0*
b2*b3^3-a0^2*a3*b2^3*b3^2+3*a0^2*a3*b1*b2*b3^3-3*a0^2*a3*b0*b3^4+
a0^2*a2*b2^2*b3^3-2*a0^2*a2*b1*b3^4-a0^2*a1*b2*b3^4+a0^3*b3^5";
fn main() {
let mut state = State::new();
let workspace = Workspace::default();
let poly: MultivariatePolynomial<_, u8> = Atom::parse(SIGMA, &mut state, &workspace)
.unwrap()
.as_view()
.to_polynomial(RationalField::new(), None)
.unwrap();
let (h, _ops, scheme) = poly.optimize_horner_scheme(4000);
let mut i = h.to_instr(poly.nvars);
println!(
"Number of operations={}, with scheme={:?}",
BorrowedHornerScheme::from(&h).op_count_cse(),
scheme,
);
i.fuse_operations();
for _ in 0..100_000 {
if !i.common_pair_elimination() {
break;
}
i.fuse_operations();
}
let op_count = i.op_count();
let o = i.to_output(true);
let o_f64 = o.convert::<f64>();
println!("Writing output to evaluate.cpp");
std::fs::write(
"evaluate.cpp",
format!(
"{}",
InstructionSetPrinter {
instr: &o,
var_map: poly.var_map.as_ref().unwrap(),
state: &state,
mode: symbolica::poly::evaluate::InstructionSetMode::CPP(
symbolica::poly::evaluate::InstructionSetModeCPPSettings {
write_header_and_test: true,
always_pass_output_array: false,
}
)
}
),
)
.unwrap();
let mut evaluator = o_f64.evaluator();
println!("Final number of operations={}", op_count);
println!(
"Evaluation = {}",
evaluator.evaluate(&(0..poly.nvars).map(|x| x as f64 + 1.).collect::<Vec<_>>())[0]
);
// evaluate with simd
let o_f64x4 = o.convert::<f64x4>();
let mut evaluator = o_f64x4.evaluator();
println!(
"Evaluation with simd = {:?}",
evaluator.evaluate(
&(0..poly.nvars)
.map(|x| f64x4::new([x as f64 + 1., x as f64 + 2., x as f64 + 3., x as f64 + 4.]))
.collect::<Vec<_>>()
)[0]
);
}