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 and if it occurs in expanded form:

from symbolica import Polynomial
p = Polynomial.parse('x^2*y + 1 + 2*x + y^4', ['x', 'y'])
print(p)
use symbolica::{
    parser::Token,
    poly::polynomial::MultivariatePolynomial,
    printer::{PolynomialPrinter, PrintOptions},
    rings::integer::IntegerRing,
    state::State,
};

fn main() {
    let var_name_map = ["x".into(), "y".into()];
    let vars: Vec<_> = var_name_map
        .iter()
        .map(|x| State::get_symbol(x))
        .collect();

    let p: MultivariatePolynomial<IntegerRing, u8> = Token::parse("x^2*y + 1 + 2*x + y^4")
        .unwrap()
        .to_polynomial(IntegerRing::new(), &vars, &var_name_map)
        .unwrap();

    println!(
        "{}",
        PolynomialPrinter::new(&p, &state, PrintOptions::default())
    );
}

or can be converted from any Symbolica expression:

from symbolica import Expression
x, y = Expression.symbols('x', 'y')
e = x**2*y + 1 + 2*x + y**4
p = e.to_polynomial()
print(p)
use symbolica::{
    parser::Token,
    poly::polynomial::MultivariatePolynomial,
    printer::{PolynomialPrinter, PrintOptions},
    rings::integer::IntegerRing,
};

fn main() {
    let a = Atom::parse("x^2*y + 1 + 2*x + y^4").unwrap();

    let p: MultivariatePolynomial<IntegerRing, u8> =
        a.to_polynomial(IntegerRing::new(), None).unwrap();

    println!("{}", p);
}

Converting from a Symbolica expressions will automatically define non-polynomial parts as new independent variables:

from symbolica import Expression
e = Expression.parse("x^2 + x + f(z + 1) / y + f(z + 1)^2")
p = e.to_polynomial()

print('poly =', p)
for i, var in enumerate(p.get_var_list()):
    print('var {} = {}'.format(i, var))
use symbolica::domains::integer::IntegerRing;
use symbolica::poly::polynomial::MultivariatePolynomial;
use symbolica::atom::Atom;

fn main() {
    let a = Atom::parse("x^2 + x + f(z + 1) / y + f(z + 1)^2").unwrap();
    let p: MultivariatePolynomial<_, u8> = a
        .to_polynomial_with_conversion(&IntegerRing::new());

    println!("poly = {}", p);
    for (i, a) in p.var_map.unwrap().iter().enumerate() {
        println!("var {} = {}", i, a.to_string());
    }
}

where f(z + 1) and 1/y are newly created variables.

Groebner basis

A Groebner basis of a polynomial system can be computed using lexicographical ordering or reverse graded lexicographical order (grevlex). The latter is often much faster.

from symbolica import Polynomial
basis = Polynomial.groebner_basis(
    [Expression.parse("a b c d - 1").to_polynomial(),
     Expression.parse("a b c + a b d + a c d + b c d").to_polynomial(),
     Expression.parse("a b + b c + a d + c d").to_polynomial(),
     Expression.parse("a + b + c + d").to_polynomial()],
    grevlex=False,
    print_stats=True
)
for p in basis:
    print(p)
use symbolica::{
    domains::finite_field::{FiniteField, FiniteFieldCore},
    poly::{groebner::GroebnerBasis, polynomial::MultivariatePolynomial, GrevLexOrder},
    atom::Atom,
    state::State,
};

fn main() {
    for x in 'a'..='z' {
        State::get_symbol(x.to_string());
    }

    // cyclic-4
    let polys = [
        "a b c d - 1",
        "a b c + a b d + a c d + b c d",
        "a b + b c + a d + c d",
        "a + b + c + d",
    ];

    let ideal: Vec<MultivariatePolynomial<_, u16>> = polys
        .iter()
        .map(|x| {
            Atom::parse(x)
                .unwrap()
                .to_polynomial(&FiniteField::<u32>::new(13), None)
                .unwrap()
        })
        .collect();

    // compute the Groebner basis with lex ordering
    let gb = GroebnerBasis::new(&ideal, true);

    println!("Lex order basis:");
    for g in &gb.system {
        println!("\t{}", g);
    }

    // compute the Groebner basis with grevlex ordering by converting the polynomials
    let grevlex_ideal: Vec<_> = ideal.iter().map(|p| p.reorder::<GrevLexOrder>()).collect();
    let gb = GroebnerBasis::new(&grevlex_ideal, true);
    println!("Grevlex order basis:");
    for g in &gb.system {
        println!("\t{}", g);
    }
}

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)
use symbolica::{
    parser::Token,
    printer::{PrintOptions, RationalPolynomialPrinter},
    rings::{
        integer::IntegerRing, rational::RationalField, rational_polynomial::RationalPolynomial,
    },
    state::{State, Workspace},
};

fn main() {
    let var_name_map = ["x".into(), "y".into()];
    let vars: Vec<_> = var_name_map
        .iter()
        .map(|x| State::get_symbol(x))
        .collect();

    let p: RationalPolynomial<IntegerRing, u8> = Token::parse("x^2*y + 1 + 2*x + y^4")
        .unwrap()
        .to_rational_polynomial(
            RationalField::new(),
            IntegerRing::new(),
            &vars,
            &var_name_map,
        )
        .unwrap();

    println!("{}", p);
}

or can be converted from a Symbolica expression.

from symbolica import Expression
x, y = Expression.symbols('x', 'y')
e = 1/x+(1+x)^2/(1+x+y)
p = e.to_rational_polynomial()
print(p)
fn main() {
    let a = Atom::parse("x^2*y + 1 + 2*x + y^4").unwrap();

    let p: RationalPolynomial<IntegerRing, u8> = a
        .to_rational_polynomial(
            RationalField::new(),
            IntegerRing::new(),
            None,
        )
        .unwrap();

    println!("{}", p);
}

Contrary to polynomials conversion, all rational polynomial conversion routines attempt to do expansions internally.

If the input contains non-rational polynomial parts, these will be considered as new independent variables. For example:

from symbolica import Expression
e = Expression.parse("x^2 + x + f(z + 1) / y + f(z + 1)^2")
p = e.to_rational_polynomial()

print('poly =', p)
for i, var in enumerate(p.get_var_list()):
    print('var {} = {}'.format(i, var))
use symbolica::domains::integer::IntegerRing;
use symbolica::domains::rational_polynomial::RationalPolynomial;
use symbolica::atom::Atom;

fn main() {
    let a = Atom::parse("x^2 + x + f(z + 1) / y + f(z + 1)^2").unwrap();

    let p: RationalPolynomial<_, u8> = a
        .to_rational_polynomial(&IntegerRing::new(), &IntegerRing::new(), None);

    println!("poly = {}", p);
    for (i, a) in p.numerator.variables.unwrap().iter().enumerate() {
        println!("var {} = {}", i, a.to_string());
    }
}

where f(z + 1), is a newly created variable.

Small-exponent versions

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},
    atom::Atom,
    rings::rational::RationalField,
};

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 poly: MultivariatePolynomial<_, u8> = Atom::parse(SIGMA)
        .unwrap()
        .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]
    );
}