A Purely Functional Computer Algebra System Embedded in Haskell
aa r X i v : . [ c s . S C ] S e p A Purely Functional Computer Algebra SystemEmbedded in Haskell
Hiromi Ishii
Doctoral Program in Mathematics,University of Tsukuba, Tsukuba, Ibaraki 305-8571, Japan [email protected]
Abstract.
We demonstrate how methods in
Functional Programming can be used to implement a computer algebra system. As a proof-of-concept, we present the computational-algebra package. It is a com-puter algebra system implemented as an embedded domain-specific lan-guage in
Haskell , a purely functional programming language. Utilisingmethods in functional programming and prominent features of Haskell,this library achieves safety, composability, and correctness at the sametime. To demonstrate the advantages of our approach, we have imple-mented advanced Gr¨obner basis algorithms, such as Faug`ere’s F and F , in a composable way. Keywords:
Gr¨obner basis; signature-based algorithms; computationalalgebra; functional programming; Haskell; type system; formal methods;property-based testing; implementation report.
In the last few decades, the area of computational algebra has grown larger. Manyalgorithms have been proposed, and there have emerged plenty of computeralgebra systems. Such systems must achieve correctness , composability and safety so that one can implement and examine new algorithms within them. Morespecifically, we want to achieve the following goals: Composability means that users can easily implement algorithms or mathe-matical objects so that they work seamlessly with existing features.
Safety prevents users and implementors from writing “wrong” code. For ex-ample, elements in different rings, e.g. Q [ x, y, z ] and Q [ w, x, y ], should betreated differently and must not directly be added. Also, it is convenient tohave handy ways to convert, inject, or coerce such values. Correctness of algorithms, with respect to prescribed formal specifications,should be guaranteed with a high assurance.We apply methods in the area of functional programming to achieve thesegoals. As a proof-of-concept, we present the computational-algebra pack-age [12]. It is implemented as an embedded domain-specific language in the able 1.
Symbols in Code FragmentsSymbol Code Symbol Code Symbol Code Symbol Code N Nat Z Integer Q Rational F p F p :: :: = == = /= λ ~ x → e \ ~ x -> e × * ⋉ !* ⌢ ++ ⊖ %- ≃ :~: ∼ ~ → -> ← <- = ⇒ ==> ⇒ => := .= : ⇐ .%= ⊆ ‘Subset‘ ≤ <= ◦ . ∧ .&&. • .* < $ > <$> < ∗ > <*> ∀ forall Haskell
Language [10]. More precisely, we adopt the
Glasgow Haskell Compiler (GHC) [7] as our hosting language. We use GHC because: its type-system allowsus to build a safe and composable interface for computer algebra; lazy evalua-tion enables us to treat infinite objects intuitively; declarative style sometimesreduces a burden of writing mathematical programs; purity permits a wide rangeof equational optimisation; and there is a plenty of libraries for functional meth-ods, especially property-based testing . These methods are not widely adopted inthis area; an exception is
DoCon [23], a pioneering work combining Haskell andcomputer algebra. Our system is designed with more emphasis on safety andcorrectness than DoCon, adding more ingredients. Although we use a functionallanguage, some methods in this paper are applicable in imperative languages.This paper is organised as follows. In Section 2, we discuss how the progressivetype-system of GHC enables us to build a safe and expressive type-system fora computer algebra. Then, in Section 3, we see how the method of property-based testing can be applied to verify the correctness of algebraic programs in alightweight and top-down manner. To demonstrate the practical advantages ofHaskell, Section 4 gives a brief description of the current implementations of theHilbert-driven, F and F algorithms. We also take a simple benchmark there.We summarise the paper and discuss related and future works in Section 5.In what follows, we use symbols in Table 1 in code fragments for readability. In this section, we will see how the progressive type-level functionalities of GHCcan be exploited to construct a safe, composable and flexible type-system fora computer algebra system. There are several existing works on type-systemsfor computer algebra, such as in Java and Scala [18, 15], and DoCon. However,none of them achieves the same level of safety and composability as our approach,which utilises the power of dependent types and type-level functions . We use type-classes , an ad-hoc polymorphism mechanism in Haskell, to en-code an algebraic hierarchy. This idea is not particularly new (for example, seeechveliani [23] or Jolly [15]), and we build our system on top of the existing algebra package [17], which provides a fine-grained abstract algebraic hierarchy.
Code 1
Group structure, coded in the algebra package class Additive a where ( + ) :: a → a → a class Additive a ⇒ Monoidal a where zero :: a class Monoidal a ⇒ Group a where negate :: a → a Code 1 illustrates a simplified version of the algebraic hierarchy up to
Group provided by the algebra package. Each statement between class or ⇒ and where , such as Additive a or Monoidal a , expresses the constraint for types.For example, Lines 1 and 2 express “a type a is Additive if it is endowed witha binary operation +”, and Lines 3 and 4 that “a type a is Monoidal if it is
Additive and has a distinguished element called zero ”.Note that, none of these requires the “proof” of algebraic axioms. Hence, onecan accidentally write a non-associative
Additive -instance, or non-distributive
Ring -instance . This sounds rather “unsafe”, and we will see how this could beaddressed reasonably in Section 3. Expressing algebraic hierarchy using type-class hierarchy, or class inheritance,is not so new and they are already implemented in DoCon or JAS. However,these systems lack a functionality to distinguish the arity of polynomials or thedenominator of a quotient ring. In particular, DoCon uses sample argumentsto indicate such parameters, and they cannot be checked at compile-time. Toovercome these restrictions, we use
Dependent Types .For example, Code 2 presents the simplified definition of the class
IsOrdPoly for polynomials. We provide an abstract class for polynomials, not just an imple-mentation, to enable users to choose appropriate internal representations fittingtheir use-cases.The class definition includes not only functions, but also associated types , or type-level functions : Arity , MOrder and
Coeff . Respectively, they correspond tothe number of variables, the monomial ordering and the coefficient ring.Note that liftMap corresponds to the universality of the polynomial ring R [ X , . . . , X n ]; i.e. the free associative commutative R -algebra over { , . . . , n } . Indeed, one can use dependent types , described in the next subsection, to require suchproofs. However, this is too heavy for the small outcome, and does not currently workfor primitive types. ode 2
A type-class for polynomials class ( Module ( Coeff poly ) poly , Commutative poly , Ring poly , CoeffRing ( Coeff poly ), IsMonomialOrder ( MOrder poly )) ⇒ IsOrdPoly poly where type Arity poly :: N type MOrder poly :: Type type Coeff poly :: Type liftMap :: ( Module ( Scalar ( Coeff poly )) alg , Ring alg ) ⇒ ( N < Arity poly → alg ) → poly → alg leadTerm :: poly → ( Coeff poly , OrdMonom ( MOrder poly ) n) ... Code 3
Examples for polynomial instances instance ( IsMonomialOrder ord , CoeffRing r) ⇒ IsOrdPoly ( OrdPoly r ord n) where type Arity ( OrdPoly r ord n) = n type MOrder ( OrdPoly r ord n) = ord type Coeff ( OrdPoly r ord n) = r ... f :: OrdPoly Q Grevlex 3 f = let [x ,y ,z] = vars in x ^ 2 × y + × x + z + instance ( CoeffRing r) ⇒ IsOrdPoly ( Unipol r) where type Arity ( OrdPoly r ord n) = 1 type MOrder ( OrdPoly r ord n) = Lex type Coeff ( OrdPoly r ord n) = r ... In theory, this function suffices to characterise the polynomial ring. However, forthe sake of efficiency, we also include some other operations in the definition.Code 3 shows example instance definitions for the standard multivariate andunivariate polynomial ring types. Note that, in Lines 8 and 12, number literal expressions type contexts. Types depending on expressions arecalled
Dependent Types in type theory. GHC supports them via the
PromotedData-types language extension [27] since version 7.4. Our library heavily uses thisfunctionality, and achieves the type-safety preventing users from unintendedlyconfusing elements from different rings.
In theory, we can use liftMap to cast between any elements of “compatible”polynomial rings. To reduce the burden to write boilerplate casting functions, ourlibrary comes with smart functions, as shown in Code 4. The convPoly function ode 4
Various casting function, with simplified type-signatures convPoly :: ( Coeff r ∼ Coeff r ’, MOrder r ∼ MOrder r ’, Arity r ∼ Arity r ’) ⇒ r → r ’ injVars :: ( Arity r ≤ Arity r ’, Coeff r ∼ Coeff r ’) ⇒ r → r ’ injVarsOffset :: (n + Arity r ≤ Arity r ’, Coeff r ∼ Coeff r ’) ⇒ Sing n → r → r ’ maps a polynomial into one with the same setting but different representation;e.g. OrdPoly Q Lex 1 into
Unipol Q . The next injVars function maps anelement of R [ X , . . . , X n ] into another polynomial ring with the same coefficientring, but with more number of variables, e.g. R [ X , . . . , X n + m ], regardless ofordering. For example, it maps Unipol Q into OrdPoly Q Grevelx 3 . Then, injVarsOffset is a variant of injVars which maps variables with offset; forexample,1 i n j V a r s O f f s e t [ sn | | ] :: Unipol Q → P o l y n o m i a l Q maps Q [ X ] into Q [ X , . . . , X ] with X X . Here, [sn | | ] is called a singleton for the type-level natural number , first introduced by Eisenberg et al. [4]. Moreprecisely, for any type-level natural n , there is the unique expression sing ::Sing n and we can use it as a tag for type-level arguments.To work with type-level naturals, we sometimes have to prove some con-straints. For example, suppose we want to write a variant of injVars mappingvariables to the end of those of the target polynomial ring, instead of the begin-ning . We might first write it as follows:1 i n j V a r s A t E n d :: ( Arity r ≤ Arity r ’ , Coeff r ∼ Coeff r ’) ⇒ r → r ’ i n j V a r s A t E n d = let sn = sing :: Sing ( Arity r ) sm = sing :: Sing ( Arity r ’) in i n j V a r s O f f s e t ( sm ⊖ sn ) -- Errors ! However, GHC cannot see
Arity r’ - Arity r + Arity r ≤ Arity r’ . Al-though this constraint is rather clear to us, we have to give the compiler itsproof. We have developed the type-natural package [14] which includes typical“lemmas”. For example, we can use the minusPlus lemma to fix this:1 -- From type - n a t u r a l: m i n u s P l u s :: Sing n → Sing m → IsTrue ( m ≤ n ) → (( n - m ) + m ) ≃ n i n j V a r s A t E n d :: ( Arity r ≤ Arity r ’ , Coeff r ∼ Coeff r ’) ⇒ r → r ’ i n j V a r s A t E n d = let sn = sing :: Sing ( Arity r ) sm = sing :: Sing ( Arity r ’) in w i t h R e f l ( m i n u s P l u s sm sn W i t n e s s) $ i n j V a r s O f f s e t ( sm ⊖ sn ) Since giving such a proof each time is rather tedious, we can use type-checkerplugins to let the compiler try to prove constraints automatically. In particu-lar, the author developed the ghc-typelits-presburger plugin [13] to resolvepropositions in Presburger arithmetic at compile time.Our library also provides the
LabPoly type, which converts existing polyno-mial types into “ labelled ” ones. For example, one can write as follows:1 f :: L a b P o l y ( P o l y n o m i a l Q
3) ’[ " x " , " y " , " z " ] f = 5 × × × + This relies on the
DataKinds and
OverloadedLabels language extensions ofGHC. GHC’s type system is strong enough to reject illegal terms and types,such as Q ) ’[ "a" ] ( w is not listed as a variable) or LabPoly (Polynomial Q
3) ’[ "x" , "y" , "x" ] (the variable x occurs twice).Using the type-level information, one can invoke the canonical inclusion mapsnaturally as follows:1 f :: LabPoly ’ Q G r e v l e x ’[ " x " , " y " , " z " ] f = × × + × × × + g :: LabPoly ’ Q Lex ’[ " w " , " z " , " y " , " u " , " x " ] g = c a n o n i c a l M a p f -- Where : c a n o n i c a l M a p :: ( xs ⊆ ys , Wraps xs poly , Wraps ys poly ’ , I s P o l y n o m i a l poly , I s P o l y n o m i a l poly ’ , Coeff poly ∼ Coeff poly ’) ⇒ L a b P o l y poly xs → L a b P o l y poly ’ ys
Since the casting functions are implemented generically, they sometimes intro-duce unnecessary overhead. For example, if one uses injVars with the same source and target types, it should just be the identity function. Fortunately, wecan use the type-safe
Rewriting Rule functionality of GHC to achieve this:1 { -
Each rewriting rule fires at compile-time, if there is a term matching the left-handside of the rule and having the same type as the right-hand side.In Haskell, it suffices just to consider algebraic laws to write down customrewriting rules. This is due to the purity of Haskell. That is, every expression inaskell is pure, in a sense that they evaluate to the same result when given thesame arguments. Note that this does not mean that Haskell cannot treat valueswith side-effects; indeed, the type-system of Haskell distinguishes pure and im-pure values at type-level, and one can treat impure operations without violatingpurity as a whole. The trick behind this situation is to describe side-effects assome kind of abstract instructions, instead of treating impure values directly.Hence, for example, duplicating the same term does not make any difference inits meaning, provided that it is algebraically correct. Such a rewriting rule isused extensively in Haskell. For example, Stream Fusion [3] uses them to elim-inate unnecessary intermediate expressions and fuse complicated functions intoefficient one-path constructions. Yet, DoCon did not do any optimisation usingrewriting rules.In our library, we also use rewriting rules to remove idempotent applicationssuch as “grading” a monomial ordering twice, e.g:1 { - ∀ ord . graded ( graded ord ) = graded ord The safety we achieved in this section cannot be achieved at compile-time with-out dependent types and type-level functions. Existing works using type-classesor class inheritance to encode algebraic hierarchy, such as JAS or DoCon, lackthis level of safety. In theory, one can achieve the same level of safety even in astatically-typed imperative language, if it supports a kind of dependent types.For example, in C ++ , templates with non-type arguments can be used to simulatedependent types. On the other hand, in Java, Generics do not allow non-typearguments and we need to mimic Peano numerals with classes. In either case,it requires much effort to prove the properties of naturals within them, becausethey lack dedicated support for type-level naturals or type-checker plugins.On the other hand, to make use of rewriting rules, we need purity as discussedabove. In this section, we will address the correctness issue, in a top-down, or lightweight manner. Especially, we apply the method of property-based testing [1] to verifythe correctness of our implementation. The idea is that one specifies the formalproperties that the implemented algorithms and types must satisfy, and checksif they hold by testing them against randomly or exhaustively generated inputs.Although it is not as rigorous as a theorem proving, it still gives a guarantee ofthe correctness at high assurance, after repeating tests time after time. ode 5
Formal Specification of Algebraic Programs prop_division :: Q → Property prop_division q = q = = ⇒ ( recip q × q = ∧ q × recip q = ∧ q × = q ∧ × q = q prop_passesSTest n = forAll ( idealOfArity n) $ λ ideal → let gs = calcGroebnerBasis ( toIdeal ideal ) in all ( isZero ◦ (‘ modPoly ‘ gs )) [ sPoly f g | f ← gs , g ← gs , f = g] Code 5 presents the example specifications for algebraic programs. In Lines1 through 4, prop_division states that the implementation of Q must satisfythe axioms of division ring. The prop_passesSTest function demand the resultof calcGroebnerBasis to pass the S -test. The tester accepts the specificationsabove, generates a specified number of inputs (default: 100) and tests againstthem. If all the inputs satisfy the specifications, it successfully halts; otherwise,it reports counterexamples, which is useful while debugging. There are several libraries for property-based testing adopting different strategiesto generate inputs. For example, QuickCheck [1] generates inputs randomly,while SmallCheck [26] exhaustively enumerates inputs in the depth-increasingorder. Even though there are other implementations of property-based testers inlanguages other than Haskell [11], it does not seem that it is applied in existingsystems, such as Singular [9], JAS or DoCon.By its generative nature, property-based testing has several drawbacks andpitfalls. First, evidently, it cannot assure the validity as rigorously as the for-mal theorem proving , unless the input space is finite. There are several piecesof research that combine formal theorem proving and computational algebrato rigorously certify correctness of implementations (for example, [24, 2]). Thesefirst formalise the theory of Gr¨obner basis in the constructive type-theory. Then,execute them within the host theorem proving language, or extract the programinto other languages. However, by its nature, this approach requires everythingto be proven formally. It is not so easy a task to prove the correctness of everypart of a program, even with help from automatic provers. Even if one managesto finish the proof of the validity of some algorithm, when one wants to optimiseit afterwards, then one must prove the “equivalence” or validity of that optimisa-tion. Moreover, it is sometimes the case that the validity, or even termination, ofthe algorithm remains unknown when it is implemented; e.g. the correctness andtermination of Fauger`e’s F [6] are proven very recently [25]. Furthermore, thereis an obvious restriction that we can extract programs only into the languagesupported by the theorem prover. We consider these conditions too restrictive,and decided to adopt theorem proving only in trivial arity arithmetic.Secondly, if the algorithm has a bad time complexity, property-based testscan easily explode. Specifically, since Gr¨obner bases have double-exponentialworst time complexity, randomly generated input can take much time to be pro-cessed. One might reduce the burden by combining randomised and enumerativegeneration strategies carefully, but there is still a possibility that there are smallinputs which take much time. To avoid such a circumstance, one can reduce thenumber of inputs, however it also reduces the assurance of validity.Finally, they are not so good at treating existential properties . AlthoughSmallCheck provides the existential quantifier in its vocabulary, it just triesto find solutions up to a prescribed depth. If solutions are relatively “larger”than its inputs, this results in false-negative failures. For example, one canwrite the following specification that demands each element of the result of calcGroebnerBasis to be a member of the original ideal, however it does notwork as expected:1 p r o p _ g b I n c ideal = let j = c a l c G r o e b n e r B a s i s ideal in exists $ λ cs → and ( z i p W i t h ( λ f gs → f = dot ideal gs ) j cs ) In the above, dot i g denotes the “dot-product”. As a workaround, we cur-rently combine inter-process communication with property-based testing. Morespecifically, we invoke a reliable existing implementation, such as SINGULAR,inside the spec as follows:1 p r o p _ g b I n c = forAll a r b i t r a r y $ λ i → m o n a d i c I O $ do let gs = c a l c G r o e b n e r B a s i s i is ← e v a l S i n g u l a r I d e a l W i t h [] [] $ funE " reduce " [ idealE gs , funE " g r o e b n e r" [ idealE i ]] return $ all isZero is Thus, if the existential property in question is decidable and has an existingreliable implementation, then it might be better to call it inside specifications. F and F algorithms for calculating Gr¨obner bases In this section, we will focus on three algorithms as case-studies: the Hilbert-driven, F and F algorithms. Firstly, we demonstrate the power of lazinessand parallelism by the Hilbert-driven algorithm. Then by the F interface, weillustrate the practical example of composability. Finally, we skim through thesimplified version of the main routine of F and see how imperative programmingwith mutable states can be written purely in Haskell. For our purpose, we williscuss only a fragment of implementations that elucidates the advantages ofHaskell, rather than the entire implementation and theoretical details. Basic API for homogenisation data Homogenised poly instance IsOrdPoly poly ⇒ IsOrdPoly ( Homogenised poly ) where type Arity ( Homogenised poly ) = 1 + Arity poly type MOrder ( Homogenised poly ) = HomogOrder ( MOrder poly ) type Coeff ( Homogenised poly ) = Coeff poly ... homogenise :: IsOrdPoly poly ⇒ poly → Homogenised poly unhomogenise :: IsOrdPoly poly ⇒ Homogenised poly → poly calcGBViaHomog :: ( Field ( Coeff poly ), IsOrdPoly poly ) ⇒ ( ∀ r. ( Field ( Coeff r), IsOrdPoly r) ⇒ Ideal r → [r ]) → Ideal poly → [ poly ] calcGBViaHomog calc i | all isHomogeneous i = calc i | otherwise = map unhomogenise ( calc ( fmap homogenise i )) Homogenisation is a powerful tool in Gr¨obner basis computation. If I ⊆ k [ X ]is a non-homogeneous ideal and ¯ I ⊆ k [ x, X ] its homogenisation, then one canget a Gr¨obner basis for I by unhomogenising the Gr¨obner basis ¯ G for ¯ I w.r.t. asuitably induced monomial ordering. In this way, any Gr¨obner basis algorithmfor homogeneous ideals can be converted into one for non-homogeneous ones.Code 6 is an API for these operations. The type Homogenised poly repre-sents polynomials obtained by homogenising polynomials of type poly . Then calcGBViaHomog calc i first checks if the input i is homogeneous. If it is so,then it applies the argument calc to its input directly (Line 15); otherwise,it first homogenises the input, applies calc , and then unhomogenises it to getthe final result (Line 16). Note that, though it uses the same term calc inboth cases, they have different types. In the first case, since it just feeds an in-put directly, calc has type Ideal poly → [poly] . On the other hand, in thenon-homogeneous case, it is applied after homogenisation, hence it is of type Ideal (Homogenised poly) → [Homogenised poly] . Thus, calcGBViaHomog takes a polymorphic function as its first argument and this is why we have ∀ inside the type of the first argument. Such a nested polymorphic type is calleda rank n polymorphic type , and it is supported by GHC’s RankNTypes languageextension . This can be achieved in object-oriented language with subtyping and Generics. ode 7
Data-type of and operations on Hilbert–Poincar´e series data HPS n = HPS { taylor :: [ Z ], hpsNumerator :: Unipol Z } instance Eq ( HPS a) where ( = ) = ( = ) ‘on ‘ hpsNumerator instance Additive ( HPS n) where HPS cs f + HPS ds g = HPS ( zipWith ( + ) cs ds ) ( f + g) instance LeftModule ( Unipol Z ) ( HPS n) where f • HPS cs g = HPS ( conv ( taylor f ⌢ repeat 0) cs ) (f × g) conv :: [ Z ] → [ Z ] → [ Z ] conv (x : xs ) (y : ys ) = let parSum a b c = a ‘par‘ b ‘par‘ c ‘seq‘ (a + b + c) in x × y : zipWith3 parSum ( map (x × ) ys ) ( map (y × ) xs ) (0 : conv xs ys ) For example, one can use the so-called
Hilbert-driven algorithm as the firstargument to calcGBViaHomog . It first computes a Gr¨obner basis w.r.t. a lightermonomial ordering, compute the Hilbert–Poincar´e series (HPS) with it and use itto compute Gr¨obner basis w.r.t. the heavier ordering. In this procedure, we needthe following operations on HPS: Equality test on HPS’s, n th Taylor coefficientof the given HPS, and the Z [ X ]-module operation on HPS. Code 7 illustratessuch an interface for HPS. For equality test, we use the numerator hpsNumerator of the closed form, and an infinite list taylor maintains Taylor coefficients. Bythe lazy nature of Haskell, we can intuitively treat infinite lists and write aconvolution on them. In Line 12, par and seq specify the evaluation strategy .In brief, expressions x and y in “ x ‘par‘ y ” (resp. seq ) are evaluated parallelly (resp. sequentially ). Since every expression is pure in Haskell, we can safely takeadvantage of parallelism, without a possibility of changing results. F F is one of the most efficient algorithms for Gr¨obner basis computation andintroduced by Faug`ere [5]. Briefly, F reduces more than two polynomials atonce, replacing S -polynomial remaindering in the Buchberger Algorithm withthe Gaussian elimination of the matrices. This means that the efficiency of F reduces to that of Gaussian elimination and the internal representation of matri-ces. Thus, it is useful if we can easily switch internal representations and elim-ination algorithms. For this purpose, we provide type-classes for mutable andimmutable matrices which admit row operations and a dedicated Gaussian elim-ination. Code 8 demonstrates the interface for immutable and mutable matrices( Matrix and
MMatrix ) and the type signature of our F implementation ( f4 ). InLines 1 and 6, the last type argument a of Matrix and
MMatrix corresponds tothe type of coefficients. Note that, one can give different instance definitions for ode 8
Matrix classes and the F function class MMatrix mat a where fromRows :: [ Vector a] → ST s ( mat s a) scaleRow :: Multiplicative a ⇒ Int → a → mat s a → ST s () ... class MMatrix ( Mutable mat ) a ⇒ Matrix mat a where type Mutable mat :: ⋆ → ⋆ freeze :: 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 , IsOrdPoly poly , Field ( Coeff poly ), Matrix mat ( Coeff poly )) ⇒ proxy mat → Strategy poly w → Ideal poly → [ poly ] the same mat but different coefficient types a . For example, one can implementefficient Gaussian elimination on F p for Matrix Mat F p , and then use it in thedefinition of Matrix Mat Q , with the Hensel lifting or Chinese remaindering.In Line 15, the first argument of f4 of type proxy mat specifies the internalrepresentation mat of matrices. In addition, f4 takes a selection strategy as thesecond argument. Here, the selection strategy is abstracted as a weighting func-tion to some ordered types, and we store intermediate polynomials in a heapand select all the polynomials with the minimum weight at each iteration. F algorithm Finally, we present the simplified version of the main routine of Faug`ere’s F [6](Code 9). Readers may be surprised that the code looks much imperative. This ismade possible by the ST monad [19], which encapsulates side-effects introducedby mutable states and prevents them from leaking outside. We use a functionalheap to choose the polynomial vectors with the least signature, demonstratingthe fusion of functional and imperative styles.
We also take a simple benchmark and the result is shown in Table 2 (examples aretaken from Giovini et al. [8]). This compares the algorithms implemented in our computational-algebra package and Singular. The first four rows correspondto the alrorithms implemented in our library; i.e. the Buchberger algorithm op-timised with syzygy and sugar strategy (B), the degree-by-degree algorithm forhomogeneous ideals (DbyD), the Hilbert-driven algorithm (Hilb), and F . S(gr)and S(sba) stand for the groebner and sba functions in the Singular computer ode 9 Main Routine of the F Algorithm f5 :: ( Field ( Coeff pol ), IsOrdPoly pol ) ⇒ Vector pol → [( Vector pol , pol )] f5 ( map monoize → i0 ) = runST $ do let n = length i0 gs ← newSTRef [] ps ← newSTRef $ H. fromList [ basis n i | i ← [0.. n -1]] syzs ← newSTRef [ sVec ( i0 ! m) ( i0 ! n) | m ← [0.. n -1] , n ← [0.. j -1] ] whileJust_ (H. viewMin < $ > readSTRef ps ) $ λ ( Entry sig g , ps ’) → do ps := ps ’ ( gs0 , ss0 ) ← ( ,) < $ > readSTRef gs < ∗ > readSTRef syzs unless ( standardCriterion sig ss0 ) $ do let (h , ph ) = reduceSignature i0 g gs0 h ’ = map ( × injectCoeff ( recip $ leadingCoeff ph )) h if isZero ph then syzs : ⇐ ( mkEntry h : ) else do let adds = fromList $ mapMaybe ( regSVec ( ph , h ’)) gs0 ps : ⇐ H. union adds gs : ⇐ (( monoize ph , mkEntry h ’) :) map ( λ (p , Entry _ a) → (a , p )) < $ > readSTRef gs algebra system 4.0.3. The complete source-code is available on GitHub [12] .The benchmark program is compiled with GHC 8.2.2 with flags -O2 -threaded-rtsopts -with-rtsopts=-N , and ran on an Intel Xeon E5-2690 at 2.90 GHz,RAM 128GB, Linux 3.16.0-4 (SMP), using 10 cores in parallel. We used theGauge framework to report the run-time of our library, and the rtimer primi-tive for Singular. For actual benchmark codes, see http://bit.ly/hbench1 and hbench2 . Unfortunately, in our system, F takes much more computing time,hence we did not include the result. The results show that, among the algo-rithms implemented in our system, F works fine in general, though it takesmuch time in some specific cases. Nevertheless, there remains much room forimprovement to compete with the state-of-the-art implementations such as Sin-gular, although there is one case where our implementation is slightly faster thanSingular’s groebner function. In this paper, we have demonstrated how we can adopt the methods developedin the area of functional programming to build a computer algebra system. Someof these methods are also applicable in imperative languages. More specifically, we used the implementation in commit . able 2. Benchmark results (ms) I (Lex) I (Grevlex) I (Lex) I (Grevlex) I (Grevlex)B 1 . × . × . × . × . × DbyD 6 . × . × . × . × . × Hilb 1 . × . × . × . × . × F . × . × . × . × . × S(gr) 2 . × . × − . × . × − . × − S(sba) 2 . × − . × − . × − . × − . × − I := h y − xy − y z + 3 x + 30 xz − z + 140 yt − u, xy − y z − x y + 45 xyz − yz + 210 y t − xt + 70 zt + 126 yu i I := h w + x + y + z, wx + xy + yz + zw, wxy + xyz + yzw + zwx, wxyz − i I := h x − x − x − y, x − z, x − t i In Section 2, we presented a type-system strong enough to detect algebraicerrors at compile-time. For example, our system can distinguish number of vari-ables of polynomial rings at type-level thanks to dependent types. It also enablesus to automatically generate casting functions and we saw how their overheadcan be reduced using rewriting rules. As for type-systems for a computer algebrasystem, there are several existing works [18, 23]. However, these systems are notsafe enough for discriminating variable arity at type-level and don’t make use ofrewriting rules.In Section 3, we successfully applied the method of property-based testing forverification of the implementation, which is lightweight compared to the existingtheorem-prover based approach [2, 24]. Although property-based testing is notas rigorous as theorem proving, it is lightweight and can be applied to algo-rithms not yet proven to be valid or terminate and available also for imperativelanguages.We have seen that, in Section 4, other features of Haskell, such as higher-orderpolymorphism, parallelism and laziness, can also be easily applied to computeralgebra by actual examples. Even though they are shown as fragments of code,we expect them to be convincing.Since some of the methods in this paper, such as dependent types or property-based testing, are not limited to the functional paradigm, it might be interestingto investigate their applicability in the imperative settings.From the viewpoint of efficiency, there are much to be done. For example,efficiency of our current F implementation is far inferior to that of the na¨ıveBuchberger algorithm, and other algorithms are far much slower than state-of-the-art implementations such as Singular. To optimise implementations, we canmake more use of Rewriting Rules and efficient data structures. Also, the par-allelism must undoubtedly play an important role. Fortunately, there are plentyof the parallel computation functionalities in Haskell, such as Regular ParallelArrays [16] and parallel package [22], and another book by Marlow [21] on gen-eral topics in parallelism in Haskell. Also, there is an existing work by Lobachevet al. [20] on parallel symbolic computation in Eden, a dialect of Haskell witharallelism support. Although Eden is retired, the methods introduced theremight be helpful. Acknowledgements
The author would like to thank my supervisor, Prof. Akira Terui, for discussions,and to anonymous reviewers for helpful comments. This research is supportedby Grant-in-Aid for JSPS Research Fellow Number 17J00479, and partially byGrants-in-Aid for Scientific Research 16K05035. This is a pre-print of an articlepublished in “Computer Algebra in Scientific Computing” (2018). The final au-thenticated version is available online at: https://doi.org/10.1007/978-3-319-99639-4
References
1. Claessen, K., Hughes, J.: QuickCheck: a lightweight tool for random testing ofHaskell programs. In: Proceedings of the Fifth ACM SIGPLAN InternationalConference on Functional Programming. ICFP ’00, pp. 268–279. ACM, New York,NY, USA (2000). http://doi.acm.org/10.1145/351240.351266
2. Coquand, T., Persson, H.: Gr¨obner bases in type theory. In: Types for Proofs andPrograms, pp. 33–46. Springer, Berlin, Heidelberg (1999)3. Coutts, D., Leshchinskiy, R., Stewart, D.: Stream Fusion. From lists to streamsto nothing at all. In: Proceedings of the 12th ACM SIGPLAN International Con-ference on Functional Programming. ICFP ’07. (2007)4. Eisenberg, R. A., Weirich, S.: Dependently typed programming with singletons.ACM SIGPLAN Notices - Haskell ’12 (12), 117–130 (2012)5. Faug`ere, J.-C.: A new efficient algorithm for computing Gr¨obner bases ( F ). Jour-nal of Pure and Applied Algebra (1), 61–88 (1999)6. Faug`ere, J.-C.: A new efficient algorithm for computing Gr¨obner bases withoutreduction to zero ( F ). In: Proceedings of the 2002 International Symposium onSymbolic and Algebraic Computation, pp. 75–83. ACM, Lille, France (2002)7. GHC Team: The Glasgow Haskell Compiler. (2018). Accessed 20188. Giovini, A., Mora, T., Niesi, G., Robbiano, L., Traverso, C.: “One sugar cube,please” or selection strategies in the Buchberger algorithm. In: Proceedings ofthe 1991 International Symposium on Symbolic and Algebraic Computation. IS-SAC’91, pp. 5–4. ACM, (1991)9. Greuel, G.-M., Pfister, G.: A Singular Introduction to Commutative Algebra. 2nd.Springer, (2007)10. Haskell Committee: The Haskell Programming Language. http://haskell.org/
11. Hypothesis: Most testing is ineffective - Hypothesis. https://hypothesis.works (2018). Accessed 06/05/201812. Ishii, H.: The computational-algebra package. 2018. https://konn.github.io/computational-algebra
13. Ishii, H.: The ghc-typelits-presburger package. http://hackage.haskell.org/package/ghc-typelits-presburger (2017)14. Ishii, H.: The type-natural package. 2013. http://hackage.haskell.org/package/type-natural
5. Jolly, R.: Categories as type classes in the Scala Algebra System. In: ComputerAlgebra in Scientific Computing, pp. 209–218. Springer, Cham (2013)16. Keller, G., Chakravarty, M. M., Leshchinskiy, R., Peyton Jones, S., Lippmeier,B.: Regular, shape-polymorphic, parallel arrays in Haskell. In: Proceedings ofthe 15th ACM SIGPLAN International Conference on Functional Programming.ICFP ’10, pp. 261–272. ACM, Baltimore, Maryland, USA (2010)17. Kmett, E. A.: The algebra package. http://hackage.haskell.org/package/algebra (2011). Accessed 201818. Kredel, H., Jolly, R.: Generic, type-safe and object oriented computer algebrasoftware. In: Computer Algebra in Scientific Computing, pp. 162–177. Springer,Berlin, Heidelberg (2010)19. Launchbury, J., Peyton Jones, S. L.: Lazy functional state threads. In: Proceed-ings of the ACM SIGPLAN 1994 Conference on Programming Language Designand Implementation. PLDI ’94, pp. 24–35. ACM, Orlando, Florida, USA (1994)20. Lobachev, O., Loogen, R.: Implementing data parallel rational multiple-residuearithmetic in Eden. In: Computer Algebra in Scientific Computing, pp. 178–193.Springer, Berlin, Heidelberg (2010)21. Marlow, S.: Parallel and Concurrent Programming in Haskell: Techniques forMulticore and Multithreaded Programming. O’Reilly Media, (2013)22. Marlow, S., Maier, P., Loidl, H.-W., Aswad, M. K., Trinder, P.: Seq no more:better strategies for Parallel Haskell. In: Proceedings of the Third ACM HaskellSymposium on Haskell. Haskell ’10, pp. 91–102. ACM, Baltimore, Maryland, USA(2010). http://doi.acm.org/10.1145/1863523.1863535
23. Mechveliani, S. D.: Computer algebra with Haskell: applying functional–categorial–“lazy” programming. In: Proceedings of International WorkshopCAAP, pp. 203–211. (2001)24. Mechveliani, S. D.: DoCon-A a Provable Algebraic Domain Construc-tor. (2018). Accessed 06/05/201825. Pan, S., Hu, Y., Wang, B.: The termination of the F5 algorithm revisited. In:Proceedings of the 38th International Symposium on Symbolic and AlgebraicComputation. ISSAC ’13, pp. 291–298. ACM, Boston, Maine, USA (2013)26. Runciman, C., Naylor, M., Lindblad, F.: SmallCheck and Lazy SmallCheck: au-tomatic exhaustive testing for small values. In: Proceedings of the First ACMSIGPLAN Symposium on Haskell. Haskell ’08, pp. 37–48. ACM, Victoria, BC,Canada (2008). http://doi.acm.org/10.1145/1411286.1411292http://doi.acm.org/10.1145/1411286.1411292