Polynomials

Polynomials are an important sub-type of general mathematical expression and lie at the heart of many mathematical problems. Since polynomials have a constrained form, Symbolica has its own data structures for univariate and multivariate polynomials, which allow 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:

Python

from symbolica import Polynomial
p = Polynomial.parse('x^2*y + 1 + 2*x + y^4', ['x', 'y'])
print(p)

Rust

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:

Python

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)

Rust

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:

Python

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))

Rust

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(&IntegerRing::new(), None);

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

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

Prime fields and Galois fields

The coefficient ring of the polynomial can be a prime field or Galois field. For example, below we define a polynomial over \(\mathbb{Z}_3\):

Python

from symbolica import *
p = Expression.parse('x^2+7x+2').to_polynomial(modulus=3)
print('poly =', p)

Rust

use symbolica::domains::finite_field::Zp;
use symbolica::poly::polynomial::MultivariatePolynomial;
use symbolica::atom::Atom;

fn main() {
    let p = Atom::parse("x^2+7x+2").to_polynomial::<_, u8>(&Zp::new(3));
    println!("poly = {}", p);
}

Note that in Rust, \(\mathbb{Z_2}\) is represented with the Z2 type.

which yields x^2+x+2.

Below we construct the Galois field \(GF(3^2)\):

Python

from symbolica import *
p = Expression.parse('x^2+7x+2').to_polynomial(modulus=3, extension=(2, Expression.symbol('t')))
print('poly =', p)

Rust

use symbolica::domains::finite_field::Zp;
use symbolica::domains::algebraic_numbers::AlgebraicExtension;
use symbolica::poly::polynomial::MultivariatePolynomial;
use symbolica::atom::Atom;

fn main() {
    let g = AlgebraicExtension::galois_field(3, 2, State::get_symbol("t"));
    let p = Atom::parse("x^2+7x+2").to_polynomial::<_, u8>(&g);
    println!("poly = {}", p);
}

where the minimal polynomial is determined automatically.

Algebraic numbers

Algebraic numbers can be defined by providing their minimal polynomial. For example, the extension \(\mathbb{Q}[\sqrt{2}]\) has the minimal polynomial \(t^2-2\). In the example below we define a polynomial in \(\mathbb{Q}[\sqrt{2}]\) and factor it:

Python

p = Expression.parse('x^4 - 10x^2 + 1').to_polynomial(minimal_poly=Expression.parse('t^2-2'))
for f, e in p.factor():
    print(f)

Rust

use symbolica::domains::finite_field::Zp;
use symbolica::domains::algebraic_numbers::AlgebraicExtension;
use symbolica::poly::polynomial::MultivariatePolynomial;
use symbolica::atom::Atom;

fn main() {
    let g = AlgebraicExtension::new(Atom::parse("t^2-2"));
    let p = Atom::parse("x^4 - 10x^2 + 1").to_polynomial::<_, u8>(&g).  to_number_field(&g);

    for (f, _) in p.factor() {
        println!("{}", f);
    }
}

we obtain:

-1+(-2*t)*x+x^2 % -2+t^2
-1+(2*t)*x+x^2 % -2+t^2

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.

Python

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)

Rust

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:

Python

from symbolica import RationalPolynomial
p = RationalPolynomial.parse('1/x+(1+x)^2/(1+x+y)', ['x', 'y'])
print(p)

Rust

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.

Python

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)

Rust

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:

Python

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))

Rust

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]
    );
}