11The Simplification System in Detail

This chapter explains in detail how simplification in Egison's computer algebra system cooperates with the type system. Chapter 9 introduced simplification as a built-in facility, but the simplification system is in fact assembled from several orthogonal, user-extensible mechanisms. This chapter walks through (1) declaring simplification rules, (2) generating rule sets mechanically by declaring ideals, (3) declaring mathematical functions, (4) selecting canonical forms with type annotations, (5) the promotion tower of types and its extension, and (6) quotient types, which live outside the tower.

11.1Declaring Simplification Rules

As we saw in Section 9.3 of Chapter 9, built-in simplification rules such as \(i^2 = -1\) are not embedded in the Haskell code of the interpreter; they are declared in the standard library with the declare rule statement. Users can declare new simplification rules with exactly the same declare rule statement. In other words, built-in rules and user rules rest on the same mechanism.

A declare rule statement is written in the following form.

declare rule (auto | name) (term | poly | frac) lhs = rhs

The first specifier determines how the rule is applied. A rule declared with auto (an automatic rule) is incorporated into the normalization procedure: it is applied automatically, repeatedly until it no longer applies, on every arithmetic operation (+, *, /, and so on) on mathematical expressions. Rules declared with a rule name (named rules) are described at the end of this section. The second specifier determines against which structure of an expression the left-hand side is matched. A term rule is applied to each term (monomial), a poly rule to each (sub-)polynomial, and a frac rule to each (sub-)fraction. In every case, the rule is also applied to subexpressions such as the arguments of function applications. The left-hand side can be a literal expression, or a pattern for mathematical expressions (Section 10.6) containing pattern variables $x and value patterns #x.

For example, the simplification of sqrt we saw in Section 9.3 of Chapter 9 (such as \(\sqrt{x} \cdot \sqrt{x} = x\)) is declared in the standard library as a rule that uses a pattern variable $x for the argument of sqrt. Conceptually, it is the following rule.

declare rule auto term (sqrt $x)^2 = x

Because the left-hand side's $x is a pattern variable, it matches whatever the argument of sqrt is, whether a concrete number or a symbol. So this single rule makes, for example, all of the following hold.

declare symbol a

(sqrt a)^2 -- a
(1 + sqrt 2)^2 -- 2 * 'sqrt 2 + 3

The actual library rule is a little more involved so that it handles not only \((\sqrt{x})^2\) but general powers \((\sqrt{x})^n\); its essence, however, is the single pattern-variable rule above.

The built-in simplification rules are declared with the same declare rule statement in the standard library file lib/math/normalize.egi. For example, the following is a part of the declarations realizing the simplifications introduced in Section 9.3 of Chapter 9.

declare rule auto term i^2 = -1
declare rule auto term log (exp $n) = n
declare rule auto term (exp $x)^$n = exp (n * x)

The first rule, i^2 = -1, is a rule about the imaginary unit \(i\). Mathematically, it introduces the algebraic extension \(\mathbb{Z}[i]\), i.e., the quotient \(\mathbb{Z}[X]/(X^2+1)\) of a polynomial ring, on the symbol \(i\). Simplification rules in Egison are global, but since symbols are globally unique, a rule can only act on expressions in which its symbol occurs. In other words, the uniqueness of symbols plays the role of scoping. The dual numbers (with \(\varepsilon^2 = 0\)) and roots of unity are algebraic extensions introduced on a symbol by the same pattern. The implementation also indexes simplification rules by the symbols occurring in their left-hand sides, and a rule is not even tried on expressions that do not contain those symbols. Therefore, declaring a rule hardly slows down computations unrelated to it. As in the second and third rules, rules whose left-hand sides contain pattern variables can also be declared.

A named rule, declared with a rule name, is not applied automatically; it is applied to an expression e only where one writes simplify e using name. For example, the addition formula for sine is a directional equation: applied from left to right it expands an expression, and applied from right to left it contracts one. Making such an equation an automatic rule would trigger the expansion in unintended places. Declaring it as a named rule instead lets us apply it explicitly, exactly where we want the expansion.

declare symbol a, b
declare rule sinAdd term sin ($x + $y) = sin x * cos y + cos x * sin y

sin (a + b) -- 'sin (b + a)
simplify (sin (a + b)) using sinAdd
-- ('cos b) * ('sin a) + ('cos a) * ('sin b)

11.2Declaring Ideals

As we saw in the previous section, declaring a simplification rule such as \(i^2 = -1\) amounts to introducing an algebraic extension — a quotient of a polynomial ring by an ideal — on top of a symbol. Pushing this view further, the rule set itself can be generated mechanically from the generators of the ideal. This section introduces the declare ideal statement and the Gröbner-basis library functions underneath it.

Hand-written rules have a completeness problem. Introduce symbols s2, s3, s6 standing for \(\sqrt{2}\), \(\sqrt{3}\), \(\sqrt{6}\); the relations evident from the definitions are the three \(s_2^2 = 2\), \(s_3^2 = 3\), and \(s_6 = s_2 s_3\). Declaring just these three as rules, however, leaves \(s_2 s_6\) (which is \(\sqrt{2}\cdot\sqrt{6} = 2\sqrt{3}\)) unreduced. Writing every needed rule by hand is hard, and even when it succeeds, termination and confluence — reaching the same normal form regardless of application order — can only be checked by inspection.

Gröbner basis theory gives the mathematical answer to this problem. From a set of generators one can compute a complete rule set (the reduced Gröbner basis) for which normal forms modulo the ideal are unique, and the termination and confluence of the resulting rules hold as theorems. In Egison, Buchberger's algorithm is implemented as the library function groebnerBasis, written in Egison itself with pattern matching (lib/math/algebra/groebner.egi).

declare symbol s2, s3, s6

groebnerBasis [s2^2 - 2, s3^2 - 3, s6 - s2 * s3]
-- [s2^2 - 2, s2 s3 - s6, s2 s6 - 2 * s3,
--  s3^2 - 3, s3 s6 - 3 * s2, s6^2 - 6]

The user supplied the three obvious relations; completion yields the six rules that form the multiplication table of \(\{1, \sqrt{2}, \sqrt{3}, \sqrt{6}\}\). Each element reads as a rewrite rule from its leading monomial to the (negated) remaining terms, such as \(s_2 s_6 \to 2 s_3\).

A computed basis can be applied to an expression explicitly with the function polyNF. polyNF gb e returns the normal form of e modulo the ideal. In particular, a result of 0 means that e equals zero as a consequence of the relations — membership in the ideal.

def gb := groebnerBasis [s2^2 - 2, s3^2 - 3, s6 - s2 * s3]

polyNF gb ((s2 + s3)^2) -- 2 * s6 + 5
polyNF gb (s2 * s3 - s6) -- 0

Note that polyNF uses the polynomials it is given as the basis as-is (it does not complete them); to get unique normal forms, pass a basis that went through groebnerBasis, as above.

To apply the rule set all the time, use the declare ideal statement. The Gröbner basis is computed from the generators once, at declaration time, and each element is registered as a term-level auto rule. From then on the rules apply automatically on every arithmetic operation, like ordinary simplification rules.

declare symbol s2, s3, s6

declare ideal [s2^2 - 2, s3^2 - 3, s6 - s2 * s3]

s2 * s6 -- 2 * s3
(s2 + s3)^2 -- 2 * s6 + 5
(s2 + s3) * (s2 - s3) -- -1

In fact, the standard library's rules for the primitive cube root of unity \(w\) have been replaced: instead of the two hand-written rules (\(w^3 = 1\) and \(w^2 = -1-w\)), the library now contains the single declaration declare ideal [w^2 + w + 1]. The generated rule \(w^2 \to -1-w\) also derives the reduction of \(w^3\) in two steps, so the completeness of the rule set becomes a theorem instead of a convention.

Compound atoms work as-is. Since sin θ and cos θ are factors (atoms) of the mathematical expression data (Section 10.1), the Pythagorean relation can be declared without any change of variables.

declare symbol `$\theta$`

declare ideal [(sin `$\theta$`)^2 + (cos `$\theta$`)^2 - 1]

(sin `$\theta$`)^4 - (cos `$\theta$`)^4 -- 2 * 'sin `$\theta$`^2 - 1
(sin `$\theta$`)^2 + 2 * (cos `$\theta$`)^2 -- - 'sin `$\theta$`^2 + 2

Which atoms survive in normal forms is determined by an ordering: symbols follow their declare symbol declaration order, and compound atoms follow their order of appearance in the generators, with earlier atoms surviving. In the example above, sin θ was written first in the generator, so the rules are generated in the direction that keeps \(\sin\) and eliminates even powers of \(\cos\).

declare ideal and polyNF divide the labor between applying always and applying once. An oriented normal form is not always what one wants, and applying heavy relations on every operation can slow the whole computation down. For a relation used only as a final check — such as the defining polynomial of a coordinate — applying it once at the point of comparison with polyNF is the better fit.

11.3Declaring Mathematical Functions

Mathematical functions such as sqrt and exp, introduced in Section 9.2, are also defined with dedicated declaration statements. The declare mathfunc statement declares a mathematical function whose applications remain as symbolic values. An application of a function that has only been declared is not simplified; it simply becomes a symbolic value. Combining it with the declare apply statement defines a simplification that is executed once, at application time. For example, the logarithm function log is declared in the standard library as follows.

declare mathfunc log
declare apply log x :=
  match x as mathValue with
    | #1 -> 0
    | #e -> 1
    | _ -> 'log x

The body of a declare apply can be an arbitrary Egison expression; for inputs that cannot be simplified, it returns a symbolic value such as 'log x using the single quote (Section 10.3). The simplification \(\sqrt{8} \rightarrow 2\sqrt{2}\) we saw in Section 9.2 is also defined by the declare apply of sqrt.

While the simplification of a declare apply runs only once at function application time, the simplification rules of declare rule are applied repeatedly, until convergence, on every normalization. The division of labor is as follows: algorithmic simplifications that inspect the argument's value and compute, such as \(\sqrt{8} \rightarrow 2\sqrt{2}\), are defined with declare apply, while rewrites on the shape of expressions, such as \(i^2 = -1\), are defined with declare rule.

Furthermore, the declare derivative statement extends the differentiation operator d/d (Section 9.5) to new mathematical functions. The built-in differentiation formulas of the standard library are also declared with this statement.

declare derivative sin = cos
declare derivative cos = \z -> - sin z
declare derivative exp = exp
declare derivative log = \z -> 1 / z
declare derivative sqrt = \z -> 1 / (2 * sqrt z)

Combining these declarations, a new mathematical function can be added to the computer algebra system in a few lines. For example, the hyperbolic secant function \(\mathrm{sech}\) can be introduced as follows. The derivative of composite functions (the chain rule) is handled by d/d itself, so only the shape of the derivative needs to be declared.

declare symbol x

declare mathfunc sech
declare apply sech x :=
  match x as mathValue with
    | #0 -> 1
    | _ -> 'sech x
declare derivative sech = \z -> - sech z * tanh z

sech 0 -- 1
sech x -- 'sech x
d/d (sech (x^2)) x -- -2 * ('sech (x^2)) * ('tanh (x^2)) * x

11.4Selecting Canonical Forms with Types

The same mathematical expression can be represented by several different data structures. For example,

can be held as a flat polynomial over the atoms \(\{i, x\}\), or as a polynomial in \(x\) whose coefficients are Gaussian integers in \(\mathbb{Z}[i]\). Which canonical form is appropriate depends on what one intends to do next (differentiate with respect to \(x\), or exploit properties of \(\mathbb{Z}[i]\)), so it cannot be determined from the value alone. In Egison, type annotations select the canonical form.

declare symbol i, x

def v := (2 + 3 * i) + (1 - i) * x + 4 * x^2

def flat : Poly Integer [i, x] := v
inspect flat
-- "4 * x^2 - i x + x + 3 * i + 2 : Poly Integer [i, x]"

def byX : Poly (Poly Integer [i]) [x] := v
inspect byX
-- "4 * x^2 + (- i + 1) * x + 3 * i + 2 : Poly (Poly Integer [i]) [x]"

The type Poly Integer [i, x] denotes integer-coefficient polynomials over the atom set \(\{i, x\}\). Nesting a polynomial type in the coefficient position, as in Poly (Poly Integer [i]) [x], denotes the canonical form “polynomials in \(x\) with coefficients in \(\mathbb{Z}[i]\)”.

An atom set can be written by naming its symbols explicitly, as in [i, x] or [x] (a closed atom set), or left open by writing [..]. An open atom set [..] means “read the atoms off the value”: the type does not fix which symbols become the polynomial's variables. The outer atom set of a nested polynomial type can likewise be [..] (fixing the inner coefficients to \(\mathbb{Z}[i]\) while leaving the outer variables to the value).

def byX3 : Poly (Poly Integer [i]) [..] := v
inspect byX3
-- "4 * x^2 + (- i + 1) * x + 3 * i + 2 : Poly (Poly Integer [i]) [x]"

Here the outer atom set is filled with the x occurring in the value v, so the observed type is Poly (Poly Integer [i]) [x] (it would become [x, y] for a value containing \(x\) and \(y\)). The promotion tower of the next section is likewise written with open atom sets such as Poly Integer [..], to denote general polynomials not tied to specific symbols.

Conversely, the inner atom set can also be left open. Fixing only the outer set means “keep exactly the named symbols as the polynomial's variables, and fold every other atom into the coefficients”.

def byI : Poly (Poly Integer [..]) [i] := v
inspect byI
-- "(- x + 3) * i + 4 * x^2 + x + 2 : Poly (Poly Integer [x]) [i]"

Here only \(i\) remains as the outer variable, and the same value is now regrouped as the polynomial \((3 - x)i + (4x^2 + x + 2)\) in \(i\) with coefficients in \(\mathbb{Z}[x]\). However, [..] may appear at most once within one nested polynomial type (a tower). With a single open atom set, every atom's destination is determined uniquely (atoms listed in a closed set go to that level; all remaining atoms go to the open level), whereas opening two positions (as in Poly (Poly Integer [..]) [..]) would make the destination of the remaining atoms ambiguous, so such a type is rejected as a type error.

A type annotation does more than classify the value: it reorganizes the value's internal structure into the shape of the type (reshape). Reshaping never changes the value. Indeed, the difference of the two canonical forms computes to \(0\).

flat - byX -- 0

Arithmetic between differently shaped representations simplifies correctly because of the convention that operations always return results in the default flat form. Nested canonical forms appear exactly where annotations request them and never propagate silently across operations.

The inspect used above displays a value together with the type read off the value (its observed type). Frequently used types can be given names with the declare cas-type statement. Aliases are transparent: displayed types use the expanded original form.

declare cas-type GaussianPolyX := Poly (Poly Integer [i]) [x]

def byX2 : GaussianPolyX := v
inspect byX2
-- "4 * x^2 + (- i + 1) * x + 3 * i + 2 : Poly (Poly Integer [i]) [x]"

11.5The Promotion Tower and Its Extension

The CAS types are ordered by value-preserving embeddings. The backbone of this order is the promotion tower, which turns the algebraic inclusions \(\mathbb{Z} \subset \mathbb{Q} \subset \mathbb{Z}[x] \subset \mathbb{Q}[x] \subset \mathbb{Q}(x)\) directly into a staircase of types.

Integer                        -- integers
  <: Frac Integer              -- rationals
  <: Poly Integer [..]         -- integer-coefficient polynomials
  <: Poly (Frac Integer) [..]  -- rational-coefficient polynomials
  <: Frac (Poly Integer [..])  -- rational functions
  <: MathValue                 -- all mathematical expressions

Each step of the tower is an embedding as a genuine subset; no value ever changes. In addition, atom-set inclusion (Poly Integer [x] <: Poly Integer [x, y]) and embeddings through the coefficient position belong to the order as well. When values of different types are mixed in one operation, they climb the tower and are automatically promoted to the least type containing both (their join).

def p : Poly Integer [x] := x + 1
def q : Frac Integer := 1 / 2

def r := p + q
inspect r -- "x + 3 / 2 : Poly (Frac Integer) [x]"

When the atom sets differ, values are promoted to the type over the union of the atom sets.

declare cas-type GaussianInt := Poly Integer [i]
declare cas-type Zsqrt2 := Poly Integer [sqrt2]

def a : GaussianInt := 1 + i
def b : Zsqrt2 := 1 + sqrt2

def c := a + b
inspect c -- "sqrt2 + i + 2 : Poly Integer [i, sqrt2]"

The sum of a \(\mathbb{Z}[i]\) value and a \(\mathbb{Z}[\sqrt{2}]\) value lands in \(\mathbb{Z}[i, \sqrt{2}]\) without any annotation.

The tower shown so far is built in, but it is not fixed: users can extend the tower. The declare cas-type statement of the previous section names canonical forms, and declare cas-subtype adds new edges (embeddings) to the tower. For example, the built-in order deliberately does not relate the flat form (Poly Integer [i, x]) and the nested form (Poly (Poly Integer [i]) [x]) of the previous section. They are different canonical forms of the same set of values, and which one to promote to is the user's choice—and precisely this choice is the information a user writes into the tower.

declare cas-subtype Poly Integer [i, x] <: Poly (Poly Integer [i]) [x]

After this declaration, operations mixing the flat and the nested form promote to the nested form.

So that arbitrary edges cannot break the tower, every declared edge is checked to keep the join of every pair of types unique over the whole order (the join-semilattice property). A declaration that would break uniqueness is rejected together with a suggested completing edge to declare first.

declare cas-subtype Poly Integer [x] <: Poly Integer [i, x]
-- Warning: this edge is already derivable (redundant edge);
--   its endpoints are still registered as nodes for the check

declare cas-subtype Poly Integer [i] <: Poly (Poly Integer [i]) [x]
-- Error: join would become ambiguous (D1 semilattice check).
--   pair (Poly Integer [i], Poly Integer [x]) would get
--   minimal upper bounds
--   {Poly Integer [i, x], Poly (Poly Integer [i]) [x]}
--   hint: declare the completing edge first:
--     declare cas-subtype Poly Integer [i, x]
--       <: Poly (Poly Integer [i]) [x]

Placing \(\mathbb{Z}[i]\) inside the nested form would split the common upper bounds of \(\mathbb{Z}[i]\) and \(\mathbb{Z}[x]\) into two—the flat and the nested form—so the declaration is rejected. The suggested completing edge is always a mathematically valid embedding (the flat form regroups into the nested one), and declaring it first makes the original edge acceptable. In this way joins remain unique over the extended tower, and implicit promotion stays deterministic. Note that the check runs over the types that have appeared in declarations so far (the redundant edge on the first line above is written to register Poly Integer [x] for the check).

11.6Quotient Types

Quotients such as \(\mathbb{Z}/7\mathbb{Z}\) (the world of remainders modulo 7) behave differently from the types so far. \(12\) and \(5\) are different values in \(\mathbb{Z}\) but equal in \(\mathbb{Z}/7\mathbb{Z}\). Since equality itself depends on the type, quotients cannot join the order of value-preserving embeddings of the previous section. In Egison, quotient types are declared outside the order with the declare cas-quotient statement.

declare cas-quotient Mod7 := Integer by (\n -> modulo n 7)

def a : Mod7 := projMod7 5
def b : Mod7 := projMod7 4

a + b -- 2
projMod7 12 == projMod7 5 -- True

The declaration automatically derives the projection projMod7 from integers into Mod7, the representative extraction reprMod7, the ring operations, and the equality test. Every operation reduces its representative, so values always stay in the range \(0\) to \(6\). At declaration time, the system also tests, on sample points, that the reduction function is idempotent and compatible with addition and multiplication (a congruence).

Mixing a Mod7 value with a raw integer (as in a + 10) is a type error. This prevents silently confusing equality in \(\mathbb{Z}\) with equality in \(\mathbb{Z}/7\mathbb{Z}\); crossing between the two worlds is explicit, via projMod7 and reprMod7. Functions can also be defined directly on the quotient. For example, the inverse in \(\mathbb{F}_7\) is one line by Fermat's little theorem (\(z^{-1} = z^5\)).

def inv7 (z : Mod7) : Mod7 := z * z * z * z * z

reprMod7 (inv7 (projMod7 3)) -- 5

The base of a quotient need not be the integers; it can be a polynomial ring. This lets us declare, as a single type, the quotient \(\mathbb{Z}[x]/(x^3)\) that treats \(x^3 = 0\) and drops every term of degree \(3\) or higher — power series truncated below degree \(3\). The reduction function scans the terms and keeps only those whose \(x\)-degree is below \(3\).

declare symbol x

-- keep a term only when its x-degree is below 3
def dropHigh (t : MathValue) : MathValue :=
  match t as mathValue with
    | term _ ((#x, $k) :: []) -> if k < 3 then t else 0
    | _ -> t

declare cas-quotient Series3 := Poly Integer [x] by
  (\p -> match p as mathValue with
           | poly $ts -> sum (map dropHigh ts)
           | _ -> p)

def s : Series3 := projSeries3 (1 + x)

reprSeries3 (s * s) -- x^2 + 2 * x + 1
reprSeries3 (projSeries3 ((1 + x)^4)) -- 6 * x^2 + 4 * x + 1

The \(x^3\) and \(x^4\) terms of \((1 + x)^4 = 1 + 4x + 6x^2 + 4x^3 + x^4\) vanish, leaving \(6x^2 + 4x + 1\). Apart from the base being \(\mathbb{Z}[x]\), the mechanism is exactly the same as for Mod7. This \(x^3 = 0\) could also be expressed by the global simplification rule declare rule auto term x^3 = 0, but then the truncation would happen everywhere x occurs. Declaring it as a quotient confines the truncation to values of type Series3, while ordinary Poly Integer [x] values keep their terms of degree \(3\) and above. This is the benefit of treating a coefficient-domain quotient as an independent type outside the tower.

Finally, as an example where the machinery of this chapter comes together, let us build the finite field of order four, \(\mathrm{GF}(4) = \mathbb{F}_2[\alpha]/(\alpha^2+\alpha+1)\). Two disciplines are needed: counting coefficients modulo \(2\) (a coefficient quotient), and reducing powers of \(\alpha\) by the relation \(\alpha^2 + \alpha + 1 = 0\) (an ideal, Section 11.2). The library function finiteFieldReduce builds the reduce that composes the two, from a prime \(p\) and the generators of the minimal polynomial. The Gröbner basis is computed once, at declaration time, and the field divisions use Fermat inverses (\(p\) must be prime).

declare symbol `$\alpha$`

declare cas-quotient GF4 := MathValue by finiteFieldReduce 2 [`$\alpha$`^2 + `$\alpha$` + 1]

def a : GF4 := projGF4 `$\alpha$`
def u : GF4 := projGF4 1

reprGF4 ((a + u) * (a + u)) -- `$\alpha$`
reprGF4 (a * (a + u)) -- 1
reprGF4 (a + a) -- 0

The field structure is maintained on every operation: \((\alpha+1)^2 = \alpha\), \(\alpha(\alpha+1) = 1\) (they are inverses of each other), and \(a + a = 0\) (characteristic \(2\)).

That the base is MathValue matters. The type contains not just the field scalars but polynomials with \(\mathrm{GF}(4)\) coefficients from the start, because the reduce applies the coefficient discipline inside every term.

declare symbol x, y

reprGF4 (projGF4 (x^2 + 3*x + `$\alpha$`^2)) -- x^2 + x + `$\alpha$` + 1

(projGF4 (`$\alpha$`*x + y)) * (projGF4 (`$\alpha$`*x + y))
  == projGF4 ((`$\alpha$`*x)^2 + y^2) -- True

In the first, the coefficient \(3\) reduces to \(1\) and \(\alpha^2\) to \(\alpha + 1\), inside the terms of the polynomial. The second is the identity \((u+v)^2 = u^2 + v^2\), which holds exactly in fields of characteristic \(2\). Two orthogonal mechanisms — coefficient quotients (this section) and rule generation from ideals (Section 11.2) — assemble into a finite field by nothing more than composing reduce functions.