Introduction

What is Algebraeon?

Algebraeon is a computer algebra system written purely in Rust. It implements algorithms for working with matrices, polynomials, algebraic numbers, factorizations, etc. The focus is on exact algebraic computations over approximate numerical solutions.

Source Code and Contributing

The project is open source and is hosted on GitHub here. Contributions are welcome and should be made via GitHub.

Stability

The API is subject to large breaking changes at this time. I hope to stabilize things more in the not too distant future.

Crates

Algebraeon is published to crates.io under an umbrella crate algebraeon and re-exports the sub-crates:

Algorithms

To give a flavour of what Algebraeon can do some of the algorithms it implements are listed:

  • Euclids algorithm for GCD and the extended version for obtaining Bezout coefficients.
  • Polynomial GCD computations using subresultant pseudo-remainder sequences.
  • AKS algorithm for natural number primality testing.
  • Matrix algorithms including:
    • Putting a matrix into Hermite normal form. In particular putting it into echelon form.
    • Putting a matrix into Smith normal form.
    • Gram–Schmidt algorithm for orthogonalization and orthonormalization.
    • Putting a matrix into Jordan normal.
    • Finding the general solution to a linear or affine system of equations.
  • Polynomial factoring algorithms including:
    • Kronecker's method for factoring polynomials over the integers (slow).
    • Berlekamp-Zassenhaus algorithm for factoring polynomials over the integers.
    • Berlekamp's algorithm for factoring polynomials over finite fields.
    • Cantor–Zassenhaus algorithm for factoring polynomials over finite fields.
    • Trager's algorithm for factoring polynomials over algebraic number fields.
  • Expressing symmetric polynomials in terms of elementary symmetric polynomials.
  • Computations with algebraic numbers:
    • Real root isolation and arithmetic.
    • Complex root isolation and arithmetic.
  • Computations with multiplication tables for small finite groups.
  • Todd-Coxeter algorithm for the enumeration of finite index cosets of a finitely generated groups.

Getting Started

As a Library

Run cargo add algebraeon in the root of your rust project to add the latest version of Algebraeon to your dependencies in your Cargo.toml

[dependencies]
algebraeon = "0.0.11"

Copy one of the example in the next section to get started.

Examples

Factoring Integers

To factor large integers using Algebraeon

#![allow(unused)]
fn main() {
use std::str::FromStr;
use algebraeon::{nzq::Natural, rings::rings::natural::factorization::factor};

let n = Natural::from_str("706000565581575429997696139445280900").unwrap();
let f = factor(n.clone()).unwrap();
println!("{} = {}", n, f);
/*
Output:
    706000565581575429997696139445280900 = 2^2 × 5^2 × 6988699669998001 × 1010203040506070809
*/
}

Algebraeon implements Lenstra elliptic-curve factorization for quickly finding prime factors with around 20 digits.

Factoring Polynomials

Factor the polynomials \(x^2 - 5x + 6\) and \(x^{15} - 1\).

#![allow(unused)]
fn main() {
use algebraeon::rings::{polynomial::*, structure::*};
use algebraeon::nzq::Integer;

let x = &Polynomial::<Integer>::var().into_ergonomic();
let f = (x.pow(2) - 5*x + 6).into_verbose();
println!("f(λ) = {}", f.factor().unwrap());
/*
Output:
    f(λ) = 1 * ((-2)+λ) * ((-3)+λ)
*/

let f = (x.pow(15) - 1).into_verbose();
println!("f(λ) = {}", f.factor().unwrap());
/*
Output:
    f(λ) = 1 * ((-1)+λ) * (1+λ+λ^2) * (1+λ+λ^2+λ^3+λ^4) * (1+(-1)λ+λ^3+(-1)λ^4+λ^5+(-1)λ^7+λ^8)
*/
}

so

\[x^2 - 5x + 6 = (x-2)(x-3)\]

\[x^{15}-1 = (x-1)(x^2+x+1)(x^4+x^3+x^2+x+1)(x^8-x^7+x^5-x^4+x^3-x+1)\]

Linear Systems of Equations

Find the general solution to the linear system

\[a \begin{pmatrix}3 \\ 4 \\ 1\end{pmatrix} + b \begin{pmatrix}2 \\ 1 \\ 2\end{pmatrix} + c \begin{pmatrix}1 \\ 3 \\ -1\end{pmatrix} = \begin{pmatrix}5 \\ 5 \\ 3\end{pmatrix}\]

for integers \(a\), \(b\) and \(c\).

#![allow(unused)]
fn main() {
use algebraeon::rings::linear::matrix::Matrix;
use algebraeon::nzq::Integer;
let x = Matrix::<Integer>::from_rows(vec![vec![3, 4, 1], vec![2, 1, 2], vec![1, 3, -1]]);
let y = Matrix::<Integer>::from_rows(vec![vec![5, 5, 3]]);
let s = x.row_solution_lattice(y);
s.pprint();
/*
Output:
    Start Affine Lattice
    Offset
    ( 2    0    -1 )
    Start Linear Lattice
    ( 1    -1    -1 )
    End Linear Lattice
    End Affine Lattice
*/
}

so the general solution is all \(a\), \(b\), \(c\) such that

\[\begin{pmatrix}a \\ b \\ c\end{pmatrix} = \begin{pmatrix}2 \\ 0 \\ -1\end{pmatrix} + t\begin{pmatrix}1 \\ -1 \\ -1\end{pmatrix}\]

for some integer \(t\).

Complex Root Isolation

Find all complex roots of the polynomial \[f(x) = x^5 + x^2 - x + 1\]

#![allow(unused)]
fn main() {
use algebraeon::rings::{polynomial::*, structure::*};
use algebraeon::nzq::Integer;

let x = &Polynomial::<Integer>::var().into_ergonomic();
let f = (x.pow(5) + x.pow(2) - x + 1).into_verbose();
// Find the complex roots of f(x)
for root in f.all_complex_roots() {
    println!("root {} of degree {}", root, root.degree());
}
/*
Output:
    root ≈-1.328 of degree 3
    root ≈0.662-0.559i of degree 3
    root ≈0.662+0.559i of degree 3
    root -i of degree 2
    root i of degree 2
*/
}

Despite the output, the roots found are not numerical approximations. Rather, they are stored internally as exact algebraic numbers by using isolating boxes in the complex plane.

Factoring Multivariable Polynomials

Factor the following multivariable polynomial with integer coefficients

\[f(x, y) = 6x^4 - 6x^3y^2 + 6xy - 6x - 6y^3 + 6y^2\]

#![allow(unused)]
fn main() {
use algebraeon::{nzq::Integer, rings::{polynomial::*, structure::*}};

let x = &MultiPolynomial::<Integer>::var(Variable::new("x")).into_ergonomic();
let y = &MultiPolynomial::<Integer>::var(Variable::new("y")).into_ergonomic();

let f = (6 * (x.pow(4) - x.pow(3) * y.pow(2) + x * y - x - y.pow(3) + y.pow(2))).into_verbose();
println!("f(x, y) = {}", f.factor().unwrap());

/*
Output:
    f(x, y) = 1 * ((3)1) * ((2)1) * (x+(-1)y^2) * (x^3+y+(-1)1)
*/
}

so the factorization of \(f(x, y)\) is

\[f(x, y) = 2 \times 3 \times (x^3 + y - 1) \times (y^2 - x)\]

P-adic Root Finding

Find the \(2\)-adic square roots of \(17\).

#![allow(unused)]
fn main() {
use algebraeon::nzq::{Natural, Integer};
use algebraeon::rings::{polynomial::*, structure::*};
let x = Polynomial::<Integer>::var().into_ergonomic();
let f = (x.pow(2) - 17).into_verbose();
for mut root in f.all_padic_roots(&Natural::from(2u32)) {
    println!("{}", root.truncate(&20.into()).string_repr()); // Show 20 2-adic digits
}
/*
Output:
    ...00110010011011101001
    ...11001101100100010111
*/
}

Truncating to the last 16 bits it can be verified that, modulo \(2^{16}\), the square of these values is \(17\).

#![allow(unused)]
fn main() {
let a = 0b0010011011101001u16;
assert_eq!(a.wrapping_mul(a), 17u16);
let b = 0b1101100100010111u16;
assert_eq!(b.wrapping_mul(b), 17u16);
}

Enumerating a Finitely Generated Group

Let \(G\) be the finitely generated group generated by \(3\) generators \(a\), \(b\), \(c\) subject to the relations \(a^2 = b^2 = c^2 = (ab)^3 = (bc)^5 = (ac)^2 = e\).

\[G = \langle a, b, c : a^2 = b^2 = c^2 = (ab)^3 = (bc)^5 = (ac)^2 = e \rangle\]

Using Algebraeon, \(G\) is found to be a finite group of order \(120\):

#![allow(unused)]
fn main() {
use algebraeon::groups::free_group::todd_coxeter::*;
let mut g = FinitelyGeneratedGroupPresentation::new();
// Add the 3 generators
let a = g.add_generator();
let b = g.add_generator();
let c = g.add_generator();
// Add the relations
g.add_relation(a.pow(2));
g.add_relation(b.pow(2));
g.add_relation(c.pow(2));
g.add_relation((&a * &b).pow(3));
g.add_relation((&b * &c).pow(5));
g.add_relation((&a * &c).pow(2));
// Count elements
let (n, _) = g.enumerate_elements();
assert_eq!(n, 120);
}

Jordan Normal Form of a Matrix

#![allow(unused)]
fn main() {
use algebraeon::nzq::{Rational};
use algebraeon::rings::{linear::matrix::*, rings::isolated_algebraic::complex::*};
use algebraeon::sets::structure::*;
// Construct a matrix
let a = Matrix::<Rational>::from_rows(vec![
    vec![5, 4, 2, 1],
    vec![0, 1, -1, -1],
    vec![-1, -1, 3, 0],
    vec![1, 1, -1, 2],
]);
// Put it into Jordan Normal Form
let j = MatrixStructure::new(ComplexAlgebraic::structure()).jordan_normal_form(&a);
j.pprint();
/*
Output:
    / 2    0    0    0 \
    | 0    1    0    0 |
    | 0    0    4    1 |
    \ 0    0    0    4 /
*/
}

Computing Discriminants

Algebraeon can find an expression for the discriminant of a polynomial in terms of the polynomials coefficients.

#![allow(unused)]
fn main() {
use algebraeon::rings::polynomial::*;
use algebraeon::nzq::Integer;

let a_var = Variable::new("a");
let b_var = Variable::new("b");
let c_var = Variable::new("c");
let d_var = Variable::new("d");
let e_var = Variable::new("e");

let a = MultiPolynomial::<Integer>::var(a_var);
let b = MultiPolynomial::<Integer>::var(b_var);
let c = MultiPolynomial::<Integer>::var(c_var);
let d = MultiPolynomial::<Integer>::var(d_var);
let e = MultiPolynomial::<Integer>::var(e_var);

let p =
    Polynomial::<MultiPolynomial<Integer>>::from_coeffs(vec![c.clone(), b.clone(), a.clone()]);
println!("p(λ) = {}", p);
println!("disc(p) = {}", p.discriminant().unwrap());

println!();

let p = Polynomial::<MultiPolynomial<Integer>>::from_coeffs(vec![
    d.clone(),
    c.clone(),
    b.clone(),
    a.clone(),
]);
println!("p(λ) = {}", p);
println!("disc(p) = {}", p.discriminant().unwrap());

println!();

let p = Polynomial::<MultiPolynomial<Integer>>::from_coeffs(vec![
    e.clone(),
    d.clone(),
    c.clone(),
    b.clone(),
    a.clone(),
]);
println!("p(λ) = {}", p);
println!("disc(p) = {}", p.discriminant().unwrap());

/*
Output:
    p(λ) = (c)+(b)λ+(a)λ^2
    disc(p) = (-4)ac+b^2

    p(λ) = (d)+(c)λ+(b)λ^2+(a)λ^3
    disc(p) = (-27)a^2d^2+(18)abcd+(-4)ac^3+(-4)b^3d+b^2c^2
*/
}

so

\[\mathop{\text{disc}}(ax^2 + bx + c) = b^2 - 4ac\]

\[\mathop{\text{disc}}(ax^3 + bx^2 + cx + d) = b^2c^2 - 4ac^3 - 4b^3d - 27a^2d^2 + 18abcd\]

The Structure System

This section describes how Algebraeon represents the mathematical idea of sets with some additional structure such as groups, rings, fields, and so on.

Structure Types

Motivation

In mathematics there are many instances of sets with some additional structure, for example:

  • The set of integers with its ordered ring structure.
  • The set of rational numbers with its ordered field structure.
  • For each natural number \(n \in \mathbb{N}\), the finite set of integers modulo \(n\) with its ring structure.
  • The set of all algebraic numbers in \(\mathbb{C}\) with its field structure.
  • The set of all ideals in a ring with the operations of ideal addition, ideal intersection, and ideal multiplication. Depending on the ring, ideals may be uniquely factorable as a product of prime ideals.

The approach taken by Algebraeon to represent such sets with additional structure is well illustrated by the example of the ring of integers modulo \(n \in \mathbb{N}\). This would be done as follows:

  • Define a structure type IntegersModuloN whose objects shall represent the ring of integers modulo \(n\) for different values of \(n\).
  • Implement the desired structure by implementing signature traits on the structure type. In the case of IntegersModuloN the required signature traits are:
    • Signature the base signature trait.
    • SetSignature<Set = Integer> so that instances of Integer shall be used to represent elements of the set of integers modulo \(n\).
    • EqSignature so that pairs of Integers can be tested for equality modulo \(n\).
    • FiniteSetSignature so that a list of all integers modulo \(n\) can be produced.
    • SemiRingSignature so that Integers can be added and multiplied modulo \(n\).
    • RingSignature so that Integers can be subtracted modulo \(n\). In practice, this could look like
#![allow(unused)]
fn main() {
use algebraeon::nzq::traits::AbsDiff;
  use algebraeon::nzq::{Integer, Natural};
  use algebraeon::sets::structure::*;

  #[derive(Debug, Clone, PartialEq, Eq)]
  struct IntegersModuloN {
      n: Natural,
  }

  impl Signature for IntegersModuloN {}

  impl SetSignature for IntegersModuloN {
      type Set = Integer;

      fn is_element(&self, x: &Integer) -> bool {
          x < &self.n
      }
  }

  impl EqSignature for IntegersModuloN {
      fn equal(&self, a: &Integer, b: &Integer) -> bool {
          a == b
      }
  }

  use algebraeon::rings::structure::{RingSignature, SemiRingSignature};

  impl SemiRingSignature for IntegersModuloN {
      fn zero(&self) -> Self::Set {
          Integer::ZERO
      }

      fn one(&self) -> Self::Set {
          (Integer::ONE % &self.n).into()
      }

      fn add(&self, a: &Self::Set, b: &Self::Set) -> Self::Set {
          ((a + b) % &self.n).into()
      }

      fn mul(&self, a: &Self::Set, b: &Self::Set) -> Self::Set {
          ((a * b) % &self.n).into()
      }
  }

  impl RingSignature for IntegersModuloN {
      fn neg(&self, a: &Self::Set) -> Self::Set {
          (-a % &self.n).into()
      }
  }

  let mod_6 = IntegersModuloN { n: 6u32.into() };
  // Since we've given `mod_6` the structure of a ring, Algebraeon implements
  // the repeated squaring algorithm for taking very large powers modulo `n`.
  assert!(mod_6.equal(
      &mod_6.nat_pow(&2.into(), &1000000000000u64.into()),
      &4.into()
  ));
}

Canonical Structures

Sometimes the situation is simple and we only want to define one set with structure rather than a family of sets, for example, the set of all rational numbers. Since sets with structure are represented in Algebraeon objects of structure tyes we will need a structure type with exactly once instance. This can be done explicitly like so

#![allow(unused)]
fn main() {
use algebraeon::{nzq::Rational, rings::structure::*, sets::structure::*};

#[derive(Debug, Clone, PartialEq, Eq)]
pub struct MyRational {
    value: Rational,
}

#[derive(Debug, Clone, PartialEq, Eq)]
pub struct MyRationalCanonicalStructure {}

impl Signature for MyRationalCanonicalStructure {}

impl SetSignature for MyRationalCanonicalStructure {
    type Set = MyRational;

    fn is_element(&self, _x: &Self::Set) -> bool {
        true
    }
}

impl EqSignature for MyRationalCanonicalStructure {
    fn equal(&self, x: &Self::Set, y: &Self::Set) -> bool {
        x == y
    }
}
}

However, Algebraeon provides a derive macro CanonicalStructure which reduces the boilerplate above to

#![allow(unused)]
fn main() {
use algebraeon::{nzq::Rational, rings::structure::*, sets::structure::*};

#[derive(Debug, Clone, PartialEq, Eq, CanonicalStructure)]
pub struct MyRational {
    value: Rational,
}
}

In any case, once we have the structure type MyRationalCanonicalStructure implementing Signature + SetSignature<Set = MyRational> we can go on to implement more structure traits like RingSignature and FieldSignature for MyRationalCanonicalStructure to give the set of instances of MyRational the structure of a ring or a field.

Multivariate Polynomials

Example

#![allow(unused)]
fn main() {
use algebraeon::nzq::Integer;
use algebraeon::rings::polynomial::*;
use algebraeon::rings::structure::*;
use std::collections::HashMap;

// Define symbols for the variables we wish to use
let x_var = Variable::new("x");
let y_var = Variable::new("y");

// Construct the polynomials x and y from the variables
let x = &MultiPolynomial::<Integer>::var(x_var.clone());
let y = &MultiPolynomial::<Integer>::var(y_var.clone());

// x and y can now be used just like other ring elements in Algebraeon
let f = MultiPolynomial::add(x, y).nat_pow(&3u32.into());

// f = x^3+(3)x^2y+(3)xy^2+y^3
println!("f = {}", f);

// For evaluating f the inputs are specified using the variable symbols
let v = f.evaluate(HashMap::from([(x_var, &1.into()), (y_var, &2.into())]));

// v = f(1, 2) = 27
println!("v = {}", v);
}

Ideals in Algebraic Rings of Integers

Factoring example

#![allow(unused)]
fn main() {
 use algebraeon::{
    nzq::*,
    rings::{
        rings::{
            algebraic_number_fields::ring_of_integer_extensions::RingOfIntegersExtension,
            finite_fields::modulo::Modulo,
        },
        polynomial::Polynomial,
        structure::*,
    },
    sets::structure::*,
};

// Construct the ring of integers Z[i]
let x = &Polynomial::<Rational>::var().into_ergonomic();
let anf = (x.pow(2) + 1).into_verbose().algebraic_number_field();
let roi = anf.ring_of_integers();

// The ideal (27i - 9) in Z[i]
let ideal = roi.principal_ideal(&roi.try_anf_to_roi(&(27 * x - 9).into_verbose()).unwrap());

// Factor the ideal
for (prime_ideal, power) in roi
    .factor_ideal(&ideal)
    .unwrap()
    .into_factor_powers()
{
    println!("power = {power} prime_ideal_factor = {:?}", prime_ideal);
}

// There's not yet a nice way to print ideals so the output is messy
// But it prints the following factorization into primes
// ideal = (1+i) * (1+2i) * (3)^2
}