Overview

The computational-algebra is the computational algebra system, implemented as a Embedded Domain Specific Language (EDSL) in Haskell. This library provides many functionality for computational algebra, especially ideal computation such as Groebner basis calculation.

Thanks to Haskell’s powerful language features, this library achieves the following goals:

Type-Safety

Haskell’s static type system enforces static correctness and prevents you from violating invariants.

Flexibility

With the powerful type-system of Haskell, we can write highly abstract program resulted in easy-to-extend system.

Efficiency

Haskell comes with many aggressive optimization mechanism and parallel computation features, which enables us to write efficient program.

This package currently provides the following functionalities:

• Groebner basis calculation w.r.t. arbitrary monomial ordering
  • Currently using Buchberger’s algorithm with some optimization
  • Faugere’s $F_4$ algorithms is experimentally implemented, but currently not as fast as Buchberger’s algorithm
• Computation in the (multivariate) polynomial ring over arbitarary field and its quotient ring
  • Ideal membership problem
  • Ideal operations such as intersection, saturation and so on.
  • Zero-dimensional ideal operation and conversion via FGLM algorithm
  • Variable elimination
• Find numeric solutions for polynomial system with real coefficient

Requirements and Installation

Old version of this package is uploaded on Hackage, but it’s rather outdated. Most recent version of computational-algebra is developed on GitHub.

It uses the most agressive language features recently implemented in Glasgow Haskell Compiler, so it requires at least GHC 8.0.1 and also it depends on many packages currently not available on Hackage, but you can install it fairly easily with help of The Haskell Tool Stack.

$ curl -sSL https://get.haskellstack.org/ | sh
  # if you haven't install Stack yet
$ git clone https://github.com/konn/computational-algebra
$ cd computational-algebra
$ stack build

In addition, you may need to install GSL and LAPACK (for matrix computation) beforehand. You can install them via Homebrew (OS X), apt-get, or other major package management systems.

Example

{-# LANGUAGE ConstraintKinds, DataKinds, GADTs, KindSignatures     #-}
{-# LANGUAGE MultiParamTypeClasses, NoImplicitPrelude              #-}
{-# LANGUAGE NoMonomorphismRestriction, QuasiQuotes, TypeOperators #-}
module Main where
import Algebra.Algorithms.Groebner
import Algebra.Field.Finite
import Algebra.Prelude
import Data.Type.Ordinal.Builtin

-- | 0-th variable of polynomial ring with at least one variable.
-- Variables are 0-origin.
x :: (KnownNat n, CoeffRing r, IsMonomialOrder n order, (0 :< n) ~ 'True)
  => OrderedPolynomial r order n
x = var [od|0|]

-- | 1-st variable of polynomial ring with at least two variable.
y :: (KnownNat n, CoeffRing r, IsMonomialOrder n order, (1 :< n) ~ 'True)
  => OrderedPolynomial r order n
y = var [od|1|]

-- | The last variable of
z :: Polynomial Rational 3
z = var [od|2|]

-- | f in QQ[x,y,z]
f :: OrderedPolynomial Rational Grevlex 3
f = 1%2*x*y^2 + y^2

-- | map f to the F_5[x,y,z], where F_5 = ZZ/5ZZ
f' :: Polynomial (F 5) 3
f' = mapCoeff (\r -> fromInteger (numerator r) / fromInteger (denominator r) ) f

-- | ideal of QQ[x,y,a,b,c,s]
heron :: Ideal (OrderedPolynomial Rational Lex 6)
heron = toIdeal [ 2 * s - a * y
                , b^2 - (x^2 + y^2)
                , c^2 - ((a - x)^2 + y^2)
                ]
  where
    -- | Get the last four variables of QQ[x,y,a,b,c,s]
    [_, _, a, b, c, s] = vars

main :: IO ()
main = do
  print f
  print f'
  print $ x * f'^2
  print $ calcGroebnerBasis heron
  -- print $ f' * 5 + f -- ^ Type error!

Type Interface for Polynomials and Algebraic Structures

computational-algebra provides well-typed interface. In this section, we will see how this package represents mathematical objects by type.

Type-level natural numbers and singletons

As we will see below, we use type-level natural number to indicate the number of variables. So, let’s see how we express natural numbers as type.

The type-natural package provides the functionality to treat type-level natural numbers seamlesly. That library also provides Peano numerals, but it is enough to use *.Builtin module for our purposes.

Sometimes we have to specify the type-level natural as function argument explicitly. We use so-called singletons for the type natural in such case. To generate singletons for type-level naturals, we can use snat quasiquoter from Data.Type.Natural.Builtin in type-natural package. Furthermore, the singletons package provides unified way to do with singletons. For more detail, please read the original paper of singletons.

For technical reason, the compiler must know the information of specific type-level natural number. This constraint is expressed as the type-classKnownNat n from GHC.TypeLits module, where $n$ is type-level natural. It seems rather noisy to have these constraints all around, but if we have singleton value sn :: SNat n for some n, then we can give such information to the compiler by withKnownNat from Data.Singletons.TypeLits of singletons package:

import Data.Singletons.TypeLits

func :: KnownNat n => a -> b -> ...

caller :: SNat n -> ...
caller sn = withKnownNat n $ func ...

Algebraic structures and operations

To express algebraic structures, we use the type classes from algebra package. For example, if we say “k is a field”, it means that k is an instance of Field class from algebra package. As mentioned above, we can compute the Groebner basis for ideals in polynomial rings over arbitary field. This means we can compute bases for those with coefficient field an instance of Field.

The ring and field operations for objects implemented in this package is provided as the instance function of Ring and Field classes. Of course, this package also provides instances for the standard type classes such as Num and Fractional, but we recommend to use the operation from algebra with NoImplicitPrlude option. We provide the convenient module Algebra.Prelude to use with NoImplicitPrlude option.

Polynomial

The type for the polynomials and operations are defined in Algebra.Ring.Polynomial module.

OrderedPolynomial r ord n represents

• the $n$-variate polynomial ring,
• over the coefficient ring r,
• with terms sorted w.r.t. the monomial ordering ord.

In the above, n should have kind Nat and r should be at least an instance of CoeffRing, which is essentially equivalent to “Commutative Ring with decidable equality”, but usually the Field for practical usage. The monomial ordering ord should be the instance of IsMonomialOrder. More precisely, ord must have instance for IsMonomial Order n ord if and only if ord stands for some monomial ordering on n-ary polynomial.

Let’s see example. $\mathbb{Q}[x,y,z]$ (the trivariate polynomial ring over the rational number) with Lex ordering is represented by the type OrderedPolynomial Rational Lex 3. Polynomial Rational 3 is short for OrderedPolynomial Rational Grevlex 3.

Abstract type-classes for Polynomials

Sometimes, one might want to use different implementation for polynomials optimized to specific task. For such a purpose, we provide two abstract type-classes IsPolynomial and IsOrderedPolynomial, defined in Algebra.Ring.Polynomial.Class module. Indeed, many algebraic operations for OrderedPolynomial are provided via these classes.

The first IsPolynomial class abstracts polynomial rings by the universality offree commutative algebra over commutative ring. The instance IsPolynomial poly means “poly is a polynomial ring”. The class comes with two associated types: Arity and Coefficient. For example, we have the following instance for OrderedPolynomial:

instance (CoeffRing r, IsMonomialOrder n ord)
      => IsPolynomial (OrderedPolynomial r ord n) where
  type Arity (OrderedPolynomial r ord n) = n
  type Coefficient (OrderedPolynomial r ord n) = r
  ...

As their name indicates, Arity poly stands for the arity of poly, that is, the number of variables of poly and Coefficient poly stands for the coefficient ring of poly. The essential class functions of it is var and liftMap:

class IsPolynomial poly where
  var     :: Ordinal (Arity poly) -> poly
  liftMap :: (Module (Scalar (Coefficient poly)) alg,
              Ring alg, Commutative alg)
          => (Ordinal (Arity poly) -> alg) -> poly -> alg

var n stands for the $n$-th variable of poly. The type Ordinal n is provided in Data.Type.Ordinal.Builtin of type-natural package, and it stands for the natural numbers less than n. So, in the context of polynomial, you can think Ordinal n as “variable of $n$-variate polynomial”. One can construct ordinals safely by the quasiquoter od provided in Data.Type.Ordinal.Builtin, when we use QuasiQutoes language extension. For example, [od|3|] stands for the third ordinal. [od|3|] :: Ordinal 4 typechecks, but [od|3|] :: Ordinal 2 is rejected in compile-time¹.

The latter function liftMap seems to have odd type, but it is just an algebraic substitution mapping. First, f :: Ordinal n -> A can be seen as “$A$-valued assignment for each variable”. Then liftMap f p extends f onto entire polynomial ring $R[X_1,\ldots,X_n]$, just substituting each variables in p using f and taking products in $A$. These are what we have calld “the universality of free algebra over commutative rings”, as pictured the following diagram:

\[\begin{xy} \xymatrix @C=10ex @R=15ex { R[X_1, \ldots, X_n] \ar @{.>}[r]^-{\mathop{\mathtt{liftMap}} f} & A\\ \{X_1, \ldots, X_n\} \ar[u]^{\mathtt{var}} \ar[ur]_{f} } \end{xy}\]

Although, we can derive other algebraic operations from these two functions in theory, but for the practical purpose, IsPolynomial class have other algebraic operations as its member functions, which can be overridden by instance-specific optimized implementation.

Polynomials and Monomial Orderings

IsPolynomial class doesn’t incorporate any information on monomial orderings. Polynomial rings with operations related monomial orderings is abstracted in IsOrderedPolynomial. This class comes with associated type MOrder, which stands for the monomial ordering of given polynomial type:

instance (...) => IsOrderedPolynomial (OrderedPolynomial r ord n) where
  type MOrder (OrderedPolynomial r ord n) = ord
  ...

This class provide the interfaces to retrieve information related to monomial orderings, such as leadingTerm, leadingMonomial and leadingCoeff.

By default, computational-algebra provides the following monomial orderings:

• Lex, the lexicographical order,
• Grlex, the graded lex order, and
• Grevlex, the graded reversed lex order.

In addition to the basic monomial orders listed above, we can construct new monomial orderings from existing ones with following:

• Graded ord, the graded order which first compares the grade (i.e. total degree) and break the tie with ord,
• ProductOrder n m ord ord', the product order which compares first n variables with ord, then the rest $m$ variables with ord', and
• WeightOrder ws ord, weighted order which compares the dot-product with ws first and then break the tie with ord.

We provide the Revlex, the reversed lex order. Revlex is not the monomial order, but we can construct monomial orderings from it with above constructors. For example, Graded Revlex is equivalent to Grevlex.

Other utility functions and related type-classes are defined in the module Algebra.Ring.Polynomial.Monomial.

How to write IsMonomialOrder instances?

If you should use the monomial ordering which cannot constructed from the above, and you have proven that ordering is really a monomial ordering, you can just implement an instance for the IsMonomialOrder.

In computational-algebra, monomials are essentially represented as Sized n Int, the Int vector of size $n$ and each $k$-th element stands for the power of $k$-th variable.

More precisely, there are two types representing monomials: Monomial and OrderedMonomial. The type Monomial n is just a synonym of Sized n Int, which is mathematically equal to $\mathbb{N}^n$. You can manipulate the value of Monomial n with functions provided by Data.Sized.Builtin from sized package.

OrderedMonomial is just a newtype wrapping Monomial tagged with additional monomial ordering information:

newtype OrderedMonomial ord n =
        OrderedMonomial { getMonomial :: Monomial n }

Note that type parameter ord doesn’t appear in the right hand side of its definition. Such type-parameters are called phantom types. The type OrderedMonomial itself doesn’t incorporate any implementation of monomial ordering, but its phantom type paramter ord carries such information.

We use IsOrder classs to retrieve ordering infromation from such pahntom types:

class IsOrder n ord where
  cmpMonomial :: Proxy ord -> Monomial n -> Monomial n -> Ordering

That is, IsOrder n ord stands for the “ord is ordering on $\mathbb{N}^n$” and cmpMonomial is the function to compare two monomials. The first argument Proxy ord is just to indicate “which order to use”, otherwise cmpMonomial can be ambiguous. For example, we have following instance for Lex ²:

instance KnownNat n => IsOrder n Lex where
  cmpMonomial _ NilL      NilL      = EQ
  cmpMonomial _ (n :< ns) (m :< ms)
    | n < m     = LT
    | n > m     = GT
    | otherwise = cmpMonomial ns ms

The type Ordering is one of the Haskell’s standard data-type which stands for the “comparison result” of two values; that is, compare n m returns LT if $n < m$, GT if $n > m$ and EQ if they are equal. Haskell’s Monoid type-class and Data.Ord module provides more convenient way to write such a comparison function. For example, we can rewrite above definition as follows:

cmpMonomial _ ns ms = mconcat (zipWithSame compare ns ms)

where zipWithSame is imported from Data.Sized.Builtin from sized package. Monoid opertions for Ordering can be considered as left-biased “breaking tie” operator.

The Ord instance for Monomial ord n is defined if IsOrder n ord is defined. But the above definition only requires ord to be “total order”; it should be monomial ordering to treat do polynomials. So, if one have proven that some ord is actually a monomial order, one should declare the instance for IsMonomialOrder as below:

instance IsMonomialOrder n Lex

IsMonomialOrder doesn’t provide any additional member function, but it is defined to distinguish mere ordering with monomial ordering. It is instance-implementor’s responsibility to assure that it is really a monomial ordering³.

So in this way, we can define the custom monomial ordering.

There is yet another type-class for monomial orderings: IsStrongMonomialOrder. IsOrder and IsMonomialOrder takes fixed arity as its parameter, but sometimes we require orderings to work with arbitarary many variables. If some specific oreder ord has IsMonomialOrder n ord for each $n$, then GHC automatically generates the instance IsStrongMonomialOrder ord. One can use cmpAnyMonomial function to compare monomials with different arity for such an ordering.

Variants of polynomial types

There are several polynomial types shipped with this library other than OrderedPolynomial:

• Unipol r defined in Algebra.Ring.Polynomial.Univariate, which stands for univariate polynomial $R[x]$ over some commutative ring $R$. It comes with operations optimized to univariate polynomials, such as efficient substitution using Horner’s rule and fast multplication using Karatsuba algorithm.

• LabPolynomial poly vars defined in Algebra.Ring.Polynomial.Labeled. It wraps existing polynomial type poly and have the same operation with it, but it labels each variables in poly by vars. Parameter vars is the type-level list of unique symbols which have the length equal to the arity of poly and. For example: haskell LabPolynomial (OrderedPolynomial Rational Grevlex 3) '["x", "y", "z"] is essentially the same as OrderedPolynomial Rational Grevlex 3 , but each variable is “labeled” with names $x$, $y$ and $z$ when we prrety-print values. By default, the following type-synonym is provided for convenience:

  type LabPolynomial' r ord '[x] = LabPolynomial (Unipol r) '[x]
  type LabPolynomial' r ord vs   = LabPolynomial (OrderedPolynomial r ord (Length vs)) vs
  type LabUnipol r x             = LabPolynomial (Unipol r) '[x]

  It also provides strongly-typed inclusion mapping. For exmaple, compiler can statically generate inclusion mapping from LabPolynomial poly '["x", "y", "z"] to LabPolynomial poly '["z", "a", "x", "b", "c"]. Furthermore, with GHC’s OverloadedLabels extension, one can use #<var> syntax to represent variables safely. For example the following type-checks and we can get what we wanted:

  #x * #y - 5 * #a^2 :: LabPolynomial' Rational Grevlex '["a", "x", "y"]

  And #z :: LabUnipol Rational "x" is statically rejected by compiler at compile-time. One limitation is that we can only use #<var> syntax only for variables starting with small alphabet and whithout any white-spaces.

Of course, users can define their custom polynomial types and made them instance of IsOrdredPolynomial. The module Algebra.Ring.Polynomial.Class provides the function injectVars, which converts between different polynomial type with the same coefficient, just mapping each variable to corresponding one with the same index in the target. Sometimes (e.g. variable elimination) one might want to permute variables. In such a case, you can just use liftMap, subst or their variants.

Other polynomial operations

• The module Algebra.Ring.Polynomial.Factorise implements the factorisation algorithm for integer-coefficient univariate polynomials.
• The module Algebra.Algorithms.ZeroDim provides various algorithms to work with zero-dimensional ideals.

Provided Algebraic Structures

This package comes with the following types for algebraic structures:

• Integer for the ring of integers,
• Rational for the rational field,
  • N.B. computational-algebra’s Rational is different from the default Rational type of Haskell. This is so because Haskell’s Ratio a type requires superfluous constraints for some algebraic instances.
• Indeed, Rational is a type synonym for Fraction Integer, where Fraction r stands for the field of fractions of an integral domain r.

Aside from the basic structurs above, we have the following structures: finite fields, quotient rings of polynomial rings, and the field of algebraic reals, which we will describe below.

Finite Fields

Algebra.Field.Finite provides the type-class for finite fields FiniteField and concrete types for prime field F p which corresponds to $\mathbb{F}_p = \mathbb{Z}/p\mathbb{Z}$. Note that, this type doesn’t check primarity of type parameter $p$ (too expensive!).

For other general finite fields other than prime fields (Galois Field), you can use Algebra.Field.Galois module provides types GF p n, which corresponds to $\mathbb{F}_{p^n}$. We use Conway polynomial for internal representation of Galois Fields. As a default, computational-algebra comes with the information of Conway polynomials for 10th power of 2,3,5,7,11. Users can easily add the information by just defining ConwayPolynomial p n instace for specific $p$ an $n$ as follows:

instance ConwayPolynomial 19 1 where
  conwayPolynomial _ _ = x ^2 + 18 * x + 2
    where x = var 0 :: Unipol (F 19)

Although we are planning to implement the functionality to automatically calculate Conway Polynomial, it is recomended to provide concrete value for each specific $p$ and $n$ to gain the efficiency. The primitive constant(s) stands for a primitive element of $\mathbb{F}_{p^n}$, i.e. a generator of the multiplicative group $\mathbb{F}_{p^n}^\times$ of units.

Galois Field computation with arbitrary irreducible polynomials

Although Conway polynomials provides systematic way to treat field extensions, it takes some computational overhead to compute Conway polynomial. So if one doesn’t need to treat field extension, it is enough to chose arbitrary irreducible polynomial of degree $n$ with coeffcients in $\mathbb{F}_p$ to do computation.

Internally, the type GF p n is synonym for GF' p n (Conway p n); here, Conway p n is a placeholder to retrieve the information of conway polynomial for $\mathbb{F}_{p^n}$. Actual computation algorithm for Galios fields is defined for GF' p n f for f carrying information of such an irreducible polynomial. So if we have some irreducible $p \in \mathbb{F}_p[x]$ with $\deg(p) = n$, one can compute in $\mathbb{F}_{p^n}$ by reflecting the information of $p$ to parameter f. The reflection package provides general way to do such a type-level reflection. Based on that, Algebra.Field.Galois provides utility function to reflect given irreducible polynomial to type-level: withIrreducible. Suppose $p \in \mathbb{F}_5$ is irreducible and $\deg(p) = 7$. Then we can do computation in $\mathbb{F}_{5^7}$ as follows:

withIrreducible p $ \pxy ->
  show (sqrt (3 `asProxyTypeOf` pxy))

In above, pxy is Proxy type to carry the information of reflected field and asProxyTypeOf forces literals to be interpreted as an element of the reflected field. One thing to note is that the type variable f dynamically reflecting polynomial cannot leak outside of given functions. For example, the value GF' p n f itself cannot be taken out from withIrreducible :

withIrreducible p $ \pxy ->
  primitive * (2 * primivite - 1) `asProxyTypeOf` pxy -- type error!

In such a situation, one cannot “take out” the reulst directly, but one can still extract the linear representation of it:

withIrreducible p $ \pxy ->
  linearRepGF (primitive * (2 * primivite - 1) `asProxyTypeOf` pxy) -- OK!

On the other hand, if we adopt Conway polynomials as a representation, one can do any computation without any scope restriction as this. This is because Conway p n carries information of an irreducible polynomial statically. So you can define Reifies instance for your custom placeholder type and store the information of some specific irreducible polynomial, then you can do such a calculation without any scoping problem:

data MyPoly = MyPoly -- ^ Just for placeholder

instance Reifies MyPoly (Unipol (F 5)) where
  reflect _ = x^2 + 2 x + 4

type MyGF5'2 = GF' 5 2 MyPoly
...

Also, Algebra.Field.Galois comes with monadic function generateIrreducible to find irreducible polynomials and reifyGF' combining these two functions. There is another function withGF' to retrieve linear representation of elements of Galois Field. See documents for more information.

Quotient rings

The type Quotient k ord n ideal stands for the quotient ring of n-variate polynomial ring over the field k. In order to distinguish the quotient ring over different ideals, we parametrize ideals in type. We use the functionalities provided by reflection package here, again.

Algebraic Reals

Algebra.Field.AlgebraicReal module provides the type Algebraic for the field of algebraic reals, i.e. real roots of real coefficient polynomials.

Internally, every algebraic real is represented by the real-coefficient polynomial $f \in \mathbb{R}[X]$ and the interval $I \subseteq \mathbb{R}$ which contains exactly one real root of $f$.

Aside the basic field operations, we currently provide the following operations on algebraic reals:

• nthRoot, where nthRoot n r calculates an $n$^th real root of the given algebraic real r. If there is no real root, it returns Nothing.
• approximate: approximate e r returns an approximating rational number $q$ with $\lvert q - r \rvert < e$. There is also a type-generic variant approxFractional, which returns any Fractional number (such as Double or Float ) instead of Rational.
• realRoots: for the univariate polynomial $f$, realRoots f computes all the real roots of $f$.
  • There is also complexRoots which computes all the complex roots of $f$, but it comes with really naive implementation and not ready for the practical usage.

Links

• API Documents

Publication

• Hiromi ISHII, “A Purely Functional Computer Algebra System Embedded in Haskell”, preprint (to appear in the Proceedings of The 20th International Workshop on Computer Algebra in Scientific Computing (CASC 2018)).