A Purely Functional
Computer Algebra System
Embedded in Haskell
Hiromi ISHII
University of Tsukuba, Japan
19th September 2018
CASC 2018 @ LilleMotivation
We apply methods in Functional
Programming such as:
Dependent Types, Property Based Testing,Purity and Rewriting Rules, ...
... to implement a computer algebra system
with:
Safety, Correctness, and Composability.

Our System
Embedded Domain Specific Language (EDSL) in
Haskell.
Haskell: A Statically-typed lazy functional
programming language
We take advantages of powerful extensions of
Glasgow Haskell Compiler (GHC) to design acomputer algebra system
Spoiler: Some methods are applicable also in otherlanguages or paradigms!

Table of Contents
Type System for Safety andComposability
Property-based Testing: Lightweightcorrectness
Current Status & Examples
Future Works & Conclusions

Table of Contents
Type System for Safety andComposability
Property-based Testing: Lightweightcorrectness
Current Status & Examples
Future Works & Conclusions

Type System for Safety
and Composability

Type-systems foralgebraType System: system to decide how values must be typed.
Types: Tags, or Invariants on values to enforce safety;"Typed terms never get stuck"
There are some existing works; e.g:Less-typed: DoCon[1]: Haskell,Java Algebra System[2].
Dependently-typed: DoCon-A[3], Coquand et al.[4]Our system sits between of the above two
We utilise a weak form of "Dependent Types"

Dependent Types?
Types depending on expressions
FULL Dep Types can simulate the higherorder logic, used in proof assistants
We use WEAK dependent types dependingon naturals and list of strings to...
distinguish the # of vars andlabel variables with unique name

"Safety" wewant here!

Example:
Polynomial arity
Suppose we have two polyn rings: R[X1, ... Xn]
and R[Y1, ..., Ym], possibly n ≠ m.
Less-typed approach: both are
represented by the same type Poly r.
Compiler should refuse such a
confusion of different rings, since it's
unclear how vars must be mapped!
Let's make Poly dependent on n or m!

Arity parametrised
polynomials
Old Poly has "kind" (type of type):
Poly ! Type! TypeOur type: Poly r n
Poly ! Type ! ℕ! TypeNow depends on Nat, not only Types!
Such types are NOT directly available in Javaor Plain Haskell.
We can still simulate type nats by phantom
types, but it adds burden w/o native support.

Throwing errors at
Compile-time
f1 ! Poly ℕ 1
f1 = 3 * #x ^ 2 + 2 * #x + 1
f2 ! Poly ℚ 1
f2 = 3 * #x ^ 2 + 2 * #x + 1
λ> f1 + f2
Different
Coeff.
Couldn't match type ‘ℚ’ with ‘ℕ’(so what?)
Type-level
g1 ! Poly ℚ 2naturals!
g1 = let [x,y] = vars in x * 2 + y
g2 ! Poly ℚ 3
g2 = let [x,y,z] = vars in z * y + x
λ> g1 + g2
Couldn't match type ‘2’ with ‘3’

Different Arity!

Generic I/F for
PolynsWe also provide a genericinterface with type-classes.
➡ MoreMaking library morecomposable

polyns optimised for
univariate case or
homogenisation, ...
We use type-level
functions to repr. theirarities, monomial
orderings and coeffs.

class(Module (Coeff poly) poly,
Ring poly, Ring (Coeff poly),
IsMonomialOrder (MOrder poly))! IsOrdPoly poly wheretype Arity poly ! ℕ
type MOrder poly ! Type
type Coeff poly ! Type
liftMap! (Module (CoeffType-levelpoly) alg,
Ring alg)
functions!
! (ℕ<Arity poly ! alg) ! poly! alg!!➡ Examples

Casting functions• We cannot add directly polyns with exactly thesame setting but with different types, by design.
• (f ! Unipol ℚ)+(g ! OrdPol ℚ Lex 1) → Error!
• We provide various casters for explicit casting!

convPoly ! (Coeff r ~ Coeff r', MOrder r ~ MOrder r',
Arity r ~ Arity r')
Exactly the same! r ! r'
settings

injVars! (Arity r ≤ Arity r', Coeff r ~ Coeff r')! r ! r'
Cast into"more" variables

Labeled Polyns
We want more flexible control of reordering of vars!
LabPoly converts any polyn type into "labelled" one, eachvariables with the unique name (LabPoly' is a synonym).
canonicalMap does exactly what we expect!
data LabPoly poly (vs ! [Symbol])
type LabPoly' r ord vs = LabPoly (OrdPoly r ord (Len vs)) vs
Specify variable intuitively

f ! LabPoly' Rational Grevlex '["x", "y"]f = 5 * #x ^ 2 * #y ^ 3 - #y + 1
f' ! LabPoly' Rationalf' = canonicalMap f

Type-level
strings!

Lex '["a", "y", "b", "x", "z"]λ> canonicalMap f ! LabPoly' Rational Grevlex '["y"]error: Couldn't match type ‘'False’ with ‘'True’Why not full
Dependent Types?
Encoding everything in Dependent Types means
proving everything (including termination)
We sometimes want to implement algorithms
whose termination is remain unknown!
We require proofs only for arity arithmetic
We developed lemma collection andcompiler plugin for Presburger arith tominimize burden.
➡ Example
Floating Point Numbers doesn't form a ring; it
can't be treated directly in such settings!

Type-system:
Summary
We use weak dependent-types for type-safety:
distinguish polynomials with different # of vars
Type-naturals are simulatable in other langs.
We save "manual proofs" by compiler plugin.
Automatically induces maps between polyns.
Rewriting Rules to reduce overheads, thanks to
the Purity (difficult in impure langs).

➡ MoreOmitted: safer quotient rings, without full dep
types but with higher polymorphism!

➡ MoreTable of Contents
Type System for Safety andComposability
Property-based Testing: Lightweightcorrectness
Current Status & Examples
Future Works & Conclusions

Table of Contents
Type System for Safety andComposability
Property-based Testing: Lightweightcorrectness
Current Status & Examples
Future Works & Conclusions

Property-based Testing:Lightweight correctness

Property-based[6]
Testing (PBT)
Tests program against formal spec, feeding some# of randomly generated or enumerated inputs.
More robust than unit tests w/ fixed inputs.
Multiple strategies for generating inputs
Not limited to Haskell; e.g. Hypothesis[7] in Python
No proof required; we can even verify algorithms
whose validity is unknown.

Example
Every non-zero

prop_division ! ℚ ! Property
rational has inverse
prop_division q =q ! 0 !! (recip q * q ! 1 && q * recip q ! 1)&& q * 1 ! q && 1 * q = q
For all n-

1 is mult unit

variate ideal
prop_passesSTestn =S-poly of any two
forAll (idealOfArity n) $ \ i !distinct elements
let gs = calcGroebnerBasis i
reduces to zeroin all (isZero . (`modPoly` gs))
[sPoly f g | f ! gs, g ! gs, f ! g]Drawbacks
Not as rigorous as formal theorem proving as in DoCon-A;
trade-off for flexibility.
Testing may take much time
Since G.b. comp has doubly-exponential worst timecomplexity, tests may explode.
We can reduce # of inputs, at the sacrifice of theconfidence.
By its nature, not so good at treating existential props.
Invoke external decision proc in such cases if available➡ MorePBT: Summary
Checks if formal specs are satisfied by
testing against generated inputs
More rigorous than fixed-input unit
tests, but less than theorem proving
Applicable to experimental algorithms
Available in many other languages
Worst complexity and existential properties
are bottlenecks for PBT.

Table of Contents
Type System for Safety andComposability
Property-based Testing: Lightweightcorrectness
Current Status & Examples
Future Works & Conclusions

Table of Contents
Type System for Safety andComposability
Property-based Testing: Lightweightcorrectness
Current Status & Examples
Future Works & Conclusions

Current Status &ExamplesImplemented
Algorithms
Groebner Basis Computation
Buchberger (naive, syzygy, sugar)
Basis conversion (FGLM, deg-by-deg, Hilbert)
Faugère's F5 and F4
Quotients by Zero-dimensional ideals
Cantor-Zassenhaus factorisation
Fields: ℚ, 𝔽p, Galois Fields, ℚ̅ (naïve)We will present a signature-based algorithm following an outline very much likethat of Buchberger’s algorithm to make connections with other approaches we havestudied more apparent. Faugère’s presentation of the original F5 algorithm looked
quite different.

F5: Pseudocode in[8]
CLO
Input: F = ( f1 , . . . , fs ), fi ∈ R
Output: φ(G), a Gröbner basis for I = ⟨ f1 , . . . , fs ⟩G := ∅

P := {e1 , . . . , es }S := {− fj ei + fi ej | 1 ≤ i < j ≤ s}
WHILE P ̸= ∅ DO

g := the element of smallest signature in P
P := P \ {g}IF Criterion(g, G ∪ S) = false THEN
h := a regular s-reduction of g by G
IF φ(h) = 0 THEN
S := S ∪ {h}

ELSE

1
h
LC(φ(h))P := P ∪ {S(k, h) | k ∈ G and S(k, h) is regular}h :=G := G ∪ {h}RETURN φ(G)Example: F5 impl
f5 ! (Field (Coeff p), IsOrdPoly pol)! Vector p ! [(Vector p, p)]
f5 i = runST $ do
let n = length i
gs ! newSTRef []ps ! newSTRef $ fromList [ basis n k | k ! [0!n-1]]
syzs ! newSTRef [ sVec im in | m ! [0!n-1], n ! [0!j-1] ]whileJust_ (H.viewMin !! readSTRef ps) $ \ (Entry sig g, ps') ! do
ps .= ps'(gs', ss') ! (,) !! readSTRef gs !! readSTRef syzs
unless (standardCriterion sig ss') $ dolet (h, ph) = reduceSignature i g gs'
h' = map (* injCoeff (recip $ leadCoeff ph)) h
if isZero ph then syzs .%= (mkEntry h : )else dolet adds = fromList $ mapMaybe (regSVec (ph, h')) gs'
ps .%= H.union adds
gs .%= ((monoize ph, mkEntry h') :)map (\ (p, Entry _ a) ! (a, p)) !! readSTRef gs

Other Impls
Generic Matrix I/F for F4➡ MorePluggable Gaussian Elimination;
users can use custom matrices.
Using laziness and parallelism
combinators in Hilbert-driven alg.

➡ MorePower series as infinite list,
computing convolutions parallellyBenchmarks (ms)Faster thanSingular's gb!
Our lib

SingularI1 (lex)

I1 (grevlex)

I2 (lex)

Buch

1.861

13.59

14.28

4.204

800.3

Hilb

104.4

160.2

25.64

16.76

7785

F5

0.5623

3.869

2.992

1.389

7.173

gb

2.55501.05542.50370.89040.9090sba

0.27170.37680.24030.25920.4221

Fastest
2nd

I2 (grevlex) I3 (grevlex)

Intel Xeon E5-2690 at 2.90 GHz, RAM 128GB, Linux 3.16.0-4 (SMP), using 10 cores

I1 := ⟨35y 4 − 30xy 2 − 210y 2z + 3x 2 + 30xz − 105z 2 + 140yt − 21u,

5xy 3 − 140y 3z − 3x 2y + 45xyz − 420yz 2 + 210y 2t − 25xt + 70zt + 126yu⟩
I2 := ⟨w + x + y + z, wx + xy + yz + zw, wxy + xyz + yzw + zwx, wxyz − 1⟩
I3 := ⟨x 31 − x 6 − x − y, x 8 − z, x 10 − t⟩

Our F4 impl took much execution time and not included
Slower than state-of-the-art impl in most cases

More on F5

Fastest
2nd
3rd

SignatureOrder

I1 (lex)

I1 (grevlex)

I2 (lex)POT

0.4138

4.262

2.837

1.286

17.62

TOP

0.5977

4.288

5.728

3.461

6.781

Our

t-POT

0.5623

3.869

2.992

1.389

7.173

F5

t-TOP

0.4860

3.879

3.100

1.360

7.319

d-POT

2.986

3.764

3.297

1.342

7.040

d-TOP

3.631

4.138

5.178

3.521

6.709

Singu

gb

2.55501.05542.50370.89040.9090larsba

0.27170.37680.24030.25920.4221I2 (grevlex) I3 (grevlex)➡ Different Env

Some heuristics can help? (TOP for high degree.... etc.)Table of Contents
Type System for Safety andComposability
Property-based Testing: Lightweightcorrectness
Current Status & Examples
Future Works & Conclusions

Table of Contents
Type System for Safety andComposability
Property-based Testing: Lightweightcorrectness
Current Status & Examples
Future Works & Conclusions

Future Works &Conclusions

Future Works
More performance tuning needed for F4
Hensel lifting or Chinese Remaindering for
Matrices...
Parallelism and GPU
Haskell's Purity empowers parallelism
There are several parallel matrices [9, 10]More Aggressive Heuristics and Rewriting?
Mixture of theorem-proving, automated proving andproperty-based testingConclusions
With weak dependent-types and higher
polymorphism, we can achieve more type-safety,
retaining flexibility as an experimental env.
Type-class and type naturals enables us to makegeneric and composable interface
Rewriting Rules can reduce the overhead
Property-based testing enables us to verify the
impl. in a lightweight manner.
Some methods are applicable in other paradigms!

References
1. Mechveliani, S. D.: Computer algebra with Haskell: applying functional–
categorial–“lazy” programming. In: Proceedings of International
Workshop CAAP, pp. 203–211. (2001)
2. Kredel, H., Jolly, R.: Generic, type-safe and object oriented computer algebra
software. In: Computer Algebra in Scientific Computing, pp. 162–
177. Springer, Berlin, Heidelberg (2010)
3. Mechveliani, S. D.: DoCon-A a Provable Algebraic Domain Constructor.
http: //www.botik.ru/pub/local/Mechveliani/docon-A/2.02/
manual.pdf (2018). Accessed 06/05/2018
4. Coquand, T., Persson, H.: Gröbner bases in type theory. In: Types for
Proofs and Programs, pp. 33–46. Springer, Berlin, Heidelberg (1999)

References
5. Kiselyov, O., Shan, C.: Functional Pearl: implicit configurations -- or, typeclass reflect the values of types. In: Proceedings of the 2004 ACM
SIGPLAN workshop on Haskell, pp.33-44. ACM (2004).
6. Claessen, K., Hughes, J.: QuickCheck: a lightweight tool for random testingof Haskell programs. In: Proceedings of the Fifth ACM SIGPLAN
International Conference on Functional Programming. ICFP ’00,
pp. 268–279. ACM (2000).
7. Hypothesis: Most testing is ineffective - Hypothesis. https://hypothesis.works (2018).
8. Cox, D.A., Little, J., O'Shea, D.: Ideals, Varieties and Algorithms.
Springer (2015).

References
9. Keller, G., Chakravarty, M. M., Leshchinskiy, R., Peyton Jones, S.,
Lippmeier, B.: Regular, shape-polymorphic, parallel arrays in Haskell. In:
Proceedings of the 15th ACM SIGPLAN International Conference
on Functional Programming. ICFP ’10, pp. 261–272. ACM,
Baltimore, Maryland, USA (2010)
10. Chakravarty, M.M.T., Keller, G., Lee, S., McDonell, T.L., Grover, V.:
Accelerating Haskell array codes with multicore GPUs. In Proceedings of thesixth workshop on Declarative aspects of multicore programming(DAMP '11). ACM (2011).

Thank You!

Any Questions?
With weak dependent-types and higher
polymorphism, we can achieve more type-safety,
retaining flexibility as an experimental env.
Type-class and type naturals enables us to makegeneric and composable interface
Rewriting Rules can reduce the overhead
Property-based testing enables us to verify the
impl. in a lightweight manner.
Some methods are applicable in other paradigms!

Appendix

Algebraic Hierarchy
as Type-classes
Expressing algebraic
hierarchy intuitively
Not so new idea
NO proofs ofalgebraic laws
required
Property-based
testing can fix thisclass Additive a where(+) ! a ! a ! aclass Additive a! Monoidal a wherezero ! a!!
class (Multiplicative a,
Monoidal a)! Ring a where
fromInteger ! ℕ ! a

Examples of
Polynomial Rings
Univariate poly,
impl. as a coeff list

data Unipol r = Unipol [r]instance CoeffRing r ! IsOrdPoly (Unipol r) wheretype Arity (Unipol r) = 1
type Coeff (Unipol r) = r
Multivar poly,
type MOrder (Unipol r) = Lex
as a fin map from!!
monomials
data OrdPoly r ord n = Unipol (Map (OMonom ord n) r)
instance (IsMonomialOrder ord, CoeffRing r)! IsOrdPoly (OrdPoly r ord n) wheretype Arity (Unipol r ord n) = n
type Coeff (Unipol r ord n) = r
type MOrder (Unipol r ord n) = ord!!Example: Easy arity
proofsExample: mapping vars to the end!
We want: injVarsEnd ! (Arity r ≤ Arity r') ! r ! r'with X1, …, Xn ↦ Xm - n + 1, …, Xm.

Singleton: type-level

We use:argument

injVarsOffset!(k + Arity r ≤ Arity r') ! Sing k ! r ! r'with X1, …, Xn ↦ Xk + 1, …, Xk + n.
Solution? injVarsOffset (sing ! Sing (m - n))
GHC cannot see m - n + n ≤ m !

Convincing GHC
with Proofs
We developed the type-natural package withmany proofs on natural numbers
Answer::: (m - n) + n = m

withRefl (minusPlus m n Witness) $injVarsOffset (sing ! Sing (m - n))

We also devised a type-checker plugin toautomatically proof props within Presburger
arithmetic.
With helps from these, much less effort isneeded to convince GHC.

Purity and Rewriting
Rules: Reducing Overhead
Casting functions are implemented genericaly;
imposes extra overhead if the mapping is trivial
We can use Rewriting Rules to reduce overhead!
LHS will be replace by RHS if the type matches.
Since every expr in Haskell is pure (w/o sideeffect), we can concentrate on algebraic validity!
convPoly b/w same types
must be identity{-# RULES
More efficient var shifting"convPoly/id" convPoly = id
within specific impl.
"injVars/monotone"
injVars (Poly dic) = Poly (mapKeysMonotone padZero dic) #-}Quotient Ring Example(omitted in the paper)
We can still achieve type-safe quotient rings!
Enabled by "Implicit Configurations" [5]tech, which utilises Rank-N Polymorphism.
Quotient of ring r by an ideal i
We cannot
data Quot r i
add Quot's
modIdeal ! Reifies i (Ideal r) ! r ! Quot r i
withQuot ! Ideal polywith! (∀i. Reifies i (Ideal poly) ! Quot poly i) different! polyi's.
This ∀ prevents ideal from leaking
instance (Reifies i (Ideal r), IsOrdPoly r) ! Ring (Quot r i)

Reading: "i carries info of an ideal / R"Existential
Properties
Naive spec for gb(I) ⊆ I:
∀ f ⃗ ∀g ∈ gb(⟨f1, ... fn⟩) ∃(c1, ..., cn) s.t. f = c1 f1 + ... + cn fn
There is no guarantee that the tester canfind c's by its generation strategy,
resulting in false-negatives!
Solution: resort to existing, externalreliable decision procedure
Our package calls SINGULAR in the spec.

Matrix I/FM for F4class MMatrix mat a where
scaleRow ! Mult a ! Int ! a ! mat s a ! ST s ()!!class MMatrix (Mutable mat) a ! Matrix mat a wheretype Mutable mat ! Type ! Typefreeze ! Mutable mat s a ! ST s (mat a)!!gaussReduction ! Field a ! mat a ! mat a
type Strategy f w = f ! f ! w
f4 ! (Ord w, !!, Matrix mat (Coeff p))! proxy mat ! Strategy p w ! Ideal p ! [p]F4 algorithm reduces g.b. to Gauss elimination
We provide im/mutable matrices classes to make algorithm
composable
We also abstract selection strategies as weighting functionHilbert Driven
data HPS n = HPS { taylor ! [ℕ], numerator ! Unipol ℕ }instance Eq (HPS a) where(!) = (!) `on` numerator
instance Additive (HPS n) where
HPS cs f + HPS ds g = HPS (zipWith (+) cs ds) (f + g)
instance LeftModule (Unipol Integer) (HPS n) wheref .* HPS cs g = HPS (conv (taylor f ! repeat 0) cs) (f * g)conv ! [ℕ] ! [ℕ] ! [ℕ]conv (x : xs) (y : ys) =let parSum a b c = a `par` b `par` c `seq` (a + b + c) inx * y :
zipWith3 parSum (map (x*) ys) (map (y*) xs) (0 : conv xs ys)

Power series as inf list, and numerator for equality
We use parallelism combinators in convolution!

Benchmark in Other
Env
SignatureOrder

I1 (lex)

I1 (grevlex)

I2 (lex)POT

0.3543

2.391

1.889

0.932211.14

TOP0.3315

2.641

4.342

2.774

5.328

Our

t-POT

0.7545

3.175

2.091

1.071

5.611

F5

t-TOP

0.6171

3.458

2.148

1.026

8.779

d-POT

3.523

2.304

1.694

0.9248

4.740

d-TOP

3.164

2.822

4.204

2.697

5.519

FastestSingu

gb

3.19551.32312.9396

1.0662

1.2709

2nd

lar

sba

0.46700.45020.37420.37580.48013rd

I2 (grevlex) I3 (grevlex)

NB: Benchmarked on my Laptop (2.8 GHz Intel
Core i7, 16GB RAM; different than other results)