HaskellForMaths 0.4.0 → 0.4.1

raw patch · 12 files changed

+751/−49 lines, 12 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• HaskellForMaths.cabal +6/−2
• Math/Algebras/Octonions.hs +69/−0
• Math/Algebras/Quaternions.hs +115/−3
• Math/Algebras/Structures.hs +1/−1
• Math/CommutativeAlgebra/GroebnerBasis.hs +14/−1
• Math/CommutativeAlgebra/Polynomial.hs +27/−35
• Math/NumberTheory/Factor.hs +190/−0
• Math/NumberTheory/Prime.hs +118/−0
• Math/Projects/ChevalleyGroup/Exceptional.hs +6/−5
• Math/Test/TCommutativeAlgebra/TPolynomial.hs +80/−0
• Math/Test/TNumberTheory/TPrimeFactor.hs +119/−0
• Math/Test/TestAll.hs +6/−2

@@ -1,5 +1,5 @@    Name:                HaskellForMaths
-   Version:             0.4.0
+   Version:             0.4.1
    Category:            Math
    Description:         A library of maths code in the areas of combinatorics, group theory, commutative algebra, and non-commutative algebra. The library is mainly intended for educational purposes, but does have efficient implementations of several fundamental algorithms.
    Synopsis:            Combinatorics, group theory, commutative algebra, non-commutative algebra
@@ -31,8 +31,10 @@         Math/Test/TCombinatorics/TDigraph.hs
         Math/Test/TCombinatorics/TIncidenceAlgebra.hs
         Math/Test/TCombinatorics/TMatroid.hs
+        Math/Test/TCommutativeAlgebra/TPolynomial.hs
         Math/Test/TCommutativeAlgebra/TGroebnerBasis.hs
         Math/Test/TCore/TField.hs
+        Math/Test/TNumberTheory/TPrimeFactor.hs
         Math/Test/TProjects/TMiniquaternionGeometry.hs
 
    Library
@@ -46,7 +48,8 @@         Math.Algebra.NonCommutative.NCPoly, Math.Algebra.NonCommutative.GSBasis, Math.Algebra.NonCommutative.TensorAlgebra,
         Math.Algebras.AffinePlane, Math.Algebras.Commutative, Math.Algebras.GroupAlgebra,
         Math.Algebras.LaurentPoly, Math.Algebras.Matrix, Math.Algebras.NonCommutative,
-        Math.Algebras.Quaternions, Math.Algebras.Structures, Math.Algebras.TensorAlgebra,
+        Math.Algebras.Octonions, Math.Algebras.Quaternions,
+        Math.Algebras.Structures, Math.Algebras.TensorAlgebra,
         Math.Algebras.TensorProduct, Math.Algebras.VectorSpace,
         Math.Combinatorics.Graph, Math.Combinatorics.GraphAuts, Math.Combinatorics.StronglyRegularGraph,
         Math.Combinatorics.Design, Math.Combinatorics.FiniteGeometry, Math.Combinatorics.Hypergraph,
@@ -55,6 +58,7 @@         Math.Common.IntegerAsType, Math.Common.ListSet,
         Math.CommutativeAlgebra.Polynomial, Math.CommutativeAlgebra.GroebnerBasis,
         Math.Core.Utils, Math.Core.Field,
+        Math.NumberTheory.Prime, Math.NumberTheory.Factor,
         Math.Projects.RootSystem,
         Math.Projects.Rubik, Math.Projects.MiniquaternionGeometry,
         Math.Projects.ChevalleyGroup.Classical, Math.Projects.ChevalleyGroup.Exceptional,

@@ -0,0 +1,69 @@+-- Copyright (c) 2011, David Amos. All rights reserved.++{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, TypeSynonymInstances, NoMonomorphismRestriction #-}++-- |A module defining the (non-associative) algebra of octonions over an arbitrary field.+--+-- The octonions are the algebra defined by the basis {1,i0,i1,i2,i3,i4,i5,i6},+-- where each i_n^2 = -1, and i_n+1*i_n+2 = i_n+4 (where the indices are modulo 7).+module Math.Algebras.Octonions where++import Math.Core.Field+import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct hiding (i1,i2)+import Math.Algebras.Structures+import Math.Algebras.Quaternions++import Math.Combinatorics.FiniteGeometry (ptsAG)+++-- Conway & Smith, On Quaternions and Octonions++-- OCTONIONS++data OBasis = O Int deriving (Eq,Ord)+-- map (return . O) [-1..6] -> [1,i0,i1,i2,i3,i4,i5,i6]++type Octonion k = Vect k OBasis++instance Show OBasis where+    show (O n) | n == -1 = "1"+               | 0 <= n && n <= 6 = "i" ++ show n+               | otherwise = error "Octonion: invalid basis element"++i0, i1, i2, i3, i4, i5, i6 :: Octonion Q+i0 = return (O 0)+i1 = return (O 1)+i2 = return (O 2)+i3 = return (O 3)+i4 = return (O 4)+i5 = return (O 5)+i6 = return (O 6)++i_ :: Num k => Int -> Octonion k+i_ n = return (O n)++instance (Num k) => Algebra k OBasis where+    unit x = x *> return (O (-1))+    mult = linear m where+        m (O (-1), O n) = return (O n)+        m (O n, O (-1)) = return (O n)+        m (O a, O b) = case (b-a) `mod` 7 of+                       0 -> -1+                       1 -> i_ ((a+3) `mod` 7)       -- i_n+1 * i_n+2 == i_n+4+                       2 -> i_ ((a+6) `mod` 7)       -- i_n+2 * i_n+4 == i_n+1+                       3 -> -1 *> i_ ((a+1) `mod` 7) -- i_n+1 * i_n+4 == -i_n+2+                       4 -> i_ ((a+5) `mod` 7)       -- i_n+4 * i_n+1 == i_n+2+                       5 -> -1 *> i_ ((a+4) `mod` 7) -- i_n+4 * i_n+2 == -i_n+1+                       6 -> -1 *> i_ ((a+2) `mod` 7) -- i_n+2 * i_n+1 == -i_n+4++instance Num k => HasConjugation k OBasis where+    conj = (>>= conj') where+        conj' (O n) = (if n == -1 then 1 else -1) *> return (O n)+    -- ie conj = linear conj', but avoiding unnecessary nf call+    sqnorm x = sum $ map ((^2) . snd) $ terms x+    -- sqnorm x = scalarPart (x * conj x)+            +-- Hence, the octonions inherit a Fractional instance++-- octonions fq = [sum $ zipWith (\x n -> x *> i_ n) xs [-1..6] | xs <- ptsAG 8 fq]

@@ -1,16 +1,20 @@ -- Copyright (c) 2010, David Amos. All rights reserved. -{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, NoMonomorphismRestriction #-}-+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, TypeSynonymInstances, NoMonomorphismRestriction #-} +-- |A module defining the algebra of quaternions over an arbitrary field.+--+-- The quaternions are the algebra defined by the basis {1,i,j,k}, where i^2 = j^2 = k^2 = ijk = -1 module Math.Algebras.Quaternions where -import Math.Algebra.Field.Base+import Math.Core.Field import Math.Algebras.VectorSpace import Math.Algebras.TensorProduct import Math.Algebras.Structures  +-- Conway & Smith, On Quaternions and Octonions+ -- QUATERNIONS  data HBasis = One | I | J | K deriving (Eq,Ord)@@ -38,11 +42,119 @@                mult' (K,I) = return J                mult' (I,K) = -1 *> return J +-- |The quaternions have {1,i,j,k} as basis, where i^2 = j^2 = k^2 = ijk = -1. i,j,k :: Num k => Quaternion k i = return I j = return J k = return K ++class Algebra k a => HasConjugation k a where+    -- |A conjugation operation is required to satisfy the following laws:+    --+    -- * conj (x+y) = conj x + conj y+    --+    -- * conj (x*y) = conj y * conj x  (note the order-reversal)+    --+    -- * conj (conj x) = x+    --+    -- * conj x = x if and only if x in k+    conj :: Vect k a -> Vect k a+    -- |The squared norm is defined as sqnorm x = x * conj x. It satisfies:+    --+    -- * sqnorm (x*y) = sqnorm x * sqnorm y+    --+    -- * sqnorm (unit k) = k^2, for k a scalar+    sqnorm :: Vect k a -> k++-- |If an algebra has a conjugation operation, then it has multiplicative inverses,+-- via 1/x = conj x / sqnorm x+instance (Fractional k, Ord a, Show a, HasConjugation k a) => Fractional (Vect k a) where+    recip 0 = error "recip 0"+    recip x = (1 / sqnorm x) *> conj x+    fromRational q = fromRational q *> 1++-- |The scalar part of the quaternion w+xi+yj+zk is w. Also called the real part.+scalarPart :: (Num k) => Quaternion k -> k+scalarPart = coeff One++-- |The vector part of the quaternion w+xi+yj+zk is xi+yj+zk. Also called the pure part.+vectorPart :: (Num k) => Quaternion k -> Quaternion k+vectorPart q = q - scalarPart q *> 1++instance Num k => HasConjugation k HBasis where+    conj = (>>= conj') where+        conj' One = return One+        conj' imag = -1 *> return imag+    -- ie conj = linear conj', but avoiding unnecessary nf call+    sqnorm x = sum $ map ((^2) . snd) $ terms x+    -- sqnorm x = scalarPart (conj x * x) -- the vector part will be zero anyway+    -- sqnorm x = x <.> x+{-+instance Fractional k => Fractional (Quaternion k) where+    recip 0 = error "Quaternion.recip 0"+    recip x = (1 / sqnorm x) *> conj x+    fromRational q = fromRational q *> 1+-}++x <.> y = scalarPart (conj x * y)+-- x <..> y = 1/2 * (sqnorm x + sqnorm y - sqnorm (x-y))++++x^-1 = recip x++-- Conway p40+refl q = \x -> -q * conj x * q++-- Given a linear function f on the quaternions, return the matrix representing it,+-- relative to a given basis. The matrix is considered as acting on the right.+asMatrix f bs = [ let fi = f ei in [ej <.> fi | ej <- bs] | ei <- bs  ]+-- It is possible to write this function using coeff, instead of <.>,+-- but then you have to pass in I,J,K, instead of i,j,k, which is uglier.++-- Conway p24+-- A homomorphism from H\0 to SO3+-- if q is restricted to unit quaternions, this is a double cover of SO3 (since q, -q induce same rotation)+-- The unit quaternions form the group Spin3+reprSO3' q = \x -> q^-1 * x * q++-- |Given a non-zero quaternion q in H, the map x -> q^-1 * x * q defines an action on the 3-dimensional vector space+-- of pure quaternions X (ie linear combinations of i,j,k). It turns out that this action is a rotation of X,+-- and this is a surjective group homomorphism from H* onto SO3. If we restrict q to the group of unit quaternions+-- (those of norm 1), then this homomorphism is 2-to-1 (since q and -q give the same rotation).+-- This shows that the multiplicative group of unit quaternions is isomorphic to Spin3, the double cover of SO3.+--+-- @reprSO3 q@ returns the 3*3 matrix representing this map.+reprSO3 :: (Fractional k) => Quaternion k -> [[k]]+reprSO3 q = reprSO3' q `asMatrix` [i,j,k]+-- It's clear from the definition that repr3' q leaves scalars invariant++-- for achiral elts, ie GO3\SO3, we compose the above with conj++-- For unit quaternions, this is a double cover of SO4 (since (l,r), (-l,-r) induce same rotation)+-- Ordered pairs of unit quaternions form the group Spin4+reprSO4' (l,r) = \x -> l^-1 * x * r+-- then (l1,r1) * (l2,r2) -> (l1*l2,r1*r2)+-- having l^-1 is required for this to work++-- |Given a pair of unit quaternions (l,r), the map x -> l^-1 * x * r defines an action on the 4-dimensional space+-- of quaternions. It turns out that this action is a rotation, and this is a surjective group homomorphism+-- onto SO4. The homomorphism is 2-to-1 (since (l,r) and (-l,-r) give the same map).+-- This shows that the multiplicative group of pairs of unit quaternions (with pointwise multiplication)+-- is isomorphic to Spin4, the double cover of SO4.+--+-- @reprSO4 (l,r)@ returns the 4*4 matrix representing this map.+reprSO4 :: (Fractional k) => (Quaternion k, Quaternion k) -> [[k]]+reprSO4 (l,r) = reprSO4' (l,r) `asMatrix` [1,i,j,k]+-- could consider checking that l,r are unit length - except that this is hard to achieve working over Q++reprSO4d lr = reprSO4 (p1 lr, p2 lr)++-- for achiral elts, GO4\SO4, we compose the above with conj+++-- DUAL SPACE OF QUATERNIONS AS COALGEBRA  one',i',j',k' :: Num k => Vect k (Dual HBasis) one' = return (Dual One)

@@ -29,7 +29,7 @@     unit :: k -> Vect k b     mult :: Vect k (Tensor b b) -> Vect k b --- |An instance declaration for Coalgebra k b is saying that the vector space Vect k b is a k-algebra.+-- |An instance declaration for Coalgebra k b is saying that the vector space Vect k b is a k-coalgebra. class Coalgebra k b where     counit :: Vect k b -> k     comult :: Vect k b -> Vect k (Tensor b b)

@@ -227,7 +227,7 @@ -- ie V(I+J) = V(I) intersect V(J) sumI :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>      [Vect k m] -> [Vect k m] -> [Vect k m]-sumI fs gs = gb $ fs ++ gs+sumI fs gs = gb (fs ++ gs)  -- Cox et al, p183 -- |Given ideals I and J, their product I.J is the ideal generated by all products {f.g | f \<- I, g \<- J}.@@ -363,6 +363,13 @@           as <-> [] = as           [] <-> bs = map negate bs +-- |In the case where every variable v occurs in some generator g of the homogeneous ideal (the usual case),+-- then the vs can be inferred from the gs.+-- @hilbertSeriesQA' gs@ returns the Hilbert series for the quotient algebra k[vs]/\<gs\>.+hilbertSeriesQA' :: (Fractional k, Ord k, MonomialConstructor m, Ord (m v), Monomial (m v), Algebra k (m v)) =>+    [Vect k (m v)] -> [Integer]+hilbertSeriesQA' gs = hilbertSeriesQA vs gs where vs = toSet (concatMap vars gs)+ -- |For i \>\> 0, the Hilbert function becomes a polynomial in i, called the Hilbert polynomial. hilbertPolyQA :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) =>      [Vect k m] -> [Vect k m] -> GlexPoly Q String@@ -372,5 +379,11 @@           hilbertPolyQA' [] x = product [ x + fromInteger j | j <- [1..n-1] ] / (fromInteger $ product [1..n-1])           hilbertPolyQA' (m:ms) x = hilbertPolyQA' ms x - hilbertPolyQA' (ms `quotientP` m) (x - fromIntegral (deg m)) +hilbertPolyQA' :: (Fractional k, Ord k, MonomialConstructor m, Ord (m v), Monomial (m v), Algebra k (m v)) =>+    [Vect k (m v)] -> GlexPoly Q String+hilbertPolyQA' gs = hilbertPolyQA vs gs where vs = toSet (concatMap vars gs)+ -- The dimension of a variety dim vs gs = 1 + deg (hilbertPolyQA vs gs)++dim' gs = 1 + deg (hilbertPolyQA' gs)

@@ -1,6 +1,6 @@ -- Copyright (c) 2011, David Amos. All rights reserved. -{-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses, FlexibleInstances #-}+{-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses, FlexibleInstances, DeriveFunctor #-}  -- |A module defining the algebra of commutative polynomials over a field k. -- Polynomials are represented as the free k-vector space with the monomials as basis.@@ -14,6 +14,7 @@ module Math.CommutativeAlgebra.Polynomial where  import Math.Core.Field+import Math.Core.Utils (toSet) import Math.Algebras.VectorSpace import Math.Algebras.TensorProduct import Math.Algebras.Structures@@ -22,7 +23,6 @@ -- Multiplication is defined by the Mon (monoid) class. Division is defined in this class. -- The functions here are primarily intended for internal use only. class (Eq m, Show m, Mon m) => Monomial m where-    -- mdivmaybe :: m -> m -> Maybe m     mdivides :: m -> m -> Bool     mdiv :: m -> m -> m     mgcd :: m -> m -> m@@ -32,8 +32,6 @@  -- mlcm m1 m2 = let m = mgcd m1 m2 in mmult m1 (mdiv m2 m) --- mdividesne m1 m2 = m1 /= m2 && mdivides m1 m2- mproperlydivides m1 m2 = m1 /= m2 && mdivides m1 m2  @@ -65,7 +63,7 @@ -- -- No Ord instance is defined for MonImpl v, so it cannot be used as the basis for a free vector space of polynomials. -- Instead, several different newtype wrappers are provided, corresponding to different monomial orderings.-data MonImpl v = M Int [(v,Int)] deriving (Eq)+data MonImpl v = M Int [(v,Int)] deriving (Eq, Functor) -- The initial Int is the degree of the monomial. Storing it speeds up equality tests and comparisons  instance Show v => Show (MonImpl v) where@@ -78,13 +76,6 @@     mmult (M si xis) (M sj yjs) = M (si+sj) $ addmerge xis yjs  instance (Ord v, Show v) => Monomial (MonImpl v) where-{--    mdivmaybe m1 m2 =-        let m@(M s xis) = mdiv m1 m2-        in if s < 0 || any (< 0) (map snd xis)-           then Nothing-           else Just m--}     mdivides (M si xis) (M sj yjs) = si <= sj && mdivides' xis yjs where         mdivides' ((x,i):xis) ((y,j):yjs) =             case compare x y of@@ -131,7 +122,7 @@ -- -- Lex stands for lexicographic ordering. -- For example, in Lex ordering, monomials up to degree two would be ordered as follows: x^2+xy+xz+x+y^2+yz+y+z^2+z+1.-newtype Lex v = Lex (MonImpl v) deriving (Eq, Mon, Monomial, MonomialConstructor) -- GeneralizedNewtypeDeriving+newtype Lex v = Lex (MonImpl v) deriving (Eq, Functor, Mon, Monomial, MonomialConstructor) -- GeneralizedNewtypeDeriving  instance Show v => Show (Lex v) where     show (Lex m) = show m@@ -177,7 +168,7 @@ -- -- Glex stands for graded lexicographic. Thus monomials are ordered first by degree, then by lexicographic order. -- For example, in Glex ordering, monomials up to degree two would be ordered as follows: x^2+xy+xz+y^2+yz+z^2+x+y+z+1.-newtype Glex v = Glex (MonImpl v) deriving (Eq, Mon, Monomial, MonomialConstructor) -- GeneralizedNewtypeDeriving+newtype Glex v = Glex (MonImpl v) deriving (Eq, Functor, Mon, Monomial, MonomialConstructor) -- GeneralizedNewtypeDeriving  instance Show v => Show (Glex v) where     show (Glex m) = show m@@ -213,7 +204,7 @@ -- For example, in Grevlex ordering, monomials up to degree two would be ordered as follows: x^2+xy+y^2+xz+yz+z^2+x+y+z+1. -- -- In general, Grevlex leads to the smallest Groebner bases.-newtype Grevlex v = Grevlex (MonImpl v) deriving (Eq, Mon, Monomial, MonomialConstructor) -- GeneralizedNewtypeDeriving+newtype Grevlex v = Grevlex (MonImpl v) deriving (Eq, Functor, Mon, Monomial, MonomialConstructor) -- GeneralizedNewtypeDeriving  instance Show v => Show (Grevlex v) where     show (Grevlex m) = show m@@ -243,7 +234,7 @@  -- ELIMINATION ORDER -data Elim2 a b = Elim2 !a !b deriving (Eq)+data Elim2 a b = Elim2 !a !b deriving (Eq, Functor)  instance (Ord a, Ord b) => Ord (Elim2 a b) where     compare (Elim2 a1 b1) (Elim2 a2 b2) = compare (a1,b1) (a2,b2)@@ -260,7 +251,6 @@     mmult (Elim2 a1 b1) (Elim2 a2 b2) = Elim2 (mmult a1 a2) (mmult b1 b2)  instance (Monomial a, Monomial b) => Monomial (Elim2 a b) where-    -- mdivmaybe :: Ord v => m v -> m v -> Maybe (m v)     mdivides (Elim2 a1 b1) (Elim2 a2 b2) = mdivides a1 a2 && mdivides b1 b2     mdiv (Elim2 a1 b1) (Elim2 a2 b2) = Elim2 (mdiv a1 a2) (mdiv b1 b2)     mgcd (Elim2 a1 b1) (Elim2 a2 b2) = Elim2 (mgcd a1 a2) (mgcd b1 b2)@@ -325,7 +315,12 @@                   Nothing -> error ("eval: no binding given for " ++ show x) -- The type could be more general than this, but haven't so far found a use case +-- |List the variables used in a polynomial+vars :: (Num k, Ord k, MonomialConstructor m, Ord (m v)) =>+     Vect k (m v) -> [Vect k (m v)]+vars f = toSet [ var v | (m,_) <- terms f, v <- map fst (mindices m) ] + -- DIVISION ALGORITHM FOR POLYNOMIALS  lt (V (t:ts)) = t -- leading term@@ -377,24 +372,20 @@         in quotRemMP' h' (reverse us', r+lth)     splitlt (V (t:ts)) = (V [t], V ts) -{---- version using mdivmaybe - surprisingly, this seems to be slower-quotRemMP2 f gs = quotRemMP' f (replicate n 0, 0) where-    n = length gs-    quotRemMP' 0 (us,r) = (us,r)-    quotRemMP' h (us,r) = divisionStep h (gs,[],us,r)-    divisionStep h (g:gs,us',u:us,r) =-        case tdivmaybe (lt h) (lt g) of-        Nothing -> divisionStep h (gs,u:us',us,r)-        Just t -> let h' = h - V [t] * g-                      u' = u + V [t]-                  in quotRemMP' h' (reverse us' ++ u':us, r)-    divisionStep h ([],us',[],r) =-        let (lth,h') = splitlt h-        in quotRemMP' h' (reverse us', r+lth)-    splitlt (V (t:ts)) = (V [t], V ts)--} +rewrite f gs = rewrite' (f,0) gs where+    rewrite' (0,r) _ = r+    rewrite' (l,r) (h:hs) =+        if lt h `tdivides` lt l -- if lhs of "rewrite rule" h matches+        then let l' = l - V [lt l `tdiv` lt h] * h -- apply rewrite rule to eliminate leading term+             in rewrite' (l',r) gs -- then start again and try to eliminate the new lt.+        else rewrite' (l,r) hs -- else try the next potential divisor+    rewrite' (l,r) [] = -- none of the rewrite rules matches lt l+        let (h,t) = split l+        in rewrite' (t, r + h) gs -- so move it into the remainder r, and try to rewrite the other terms+    split (V (t:ts)) = (V [t], V ts)++ infixl 7 %%  -- |@f %% gs@ is the reduction of a polynomial f with respect to a list of polynomials gs.@@ -403,7 +394,8 @@ -- and is zero if and only if f is in I. (%%) :: (Fractional k, Monomial m, Ord m, Algebra k m) =>      Vect k m -> [Vect k m] -> Vect k m-f %% gs = r where (_,r) = quotRemMP f gs+f %% gs = rewrite f gs+-- f %% gs = r where (_,r) = quotRemMP f gs   -- |As a convenience, a partial instance of Fractional is defined for polynomials.

@@ -0,0 +1,190 @@+-- Copyright (c) 2006-2011, David Amos. All rights reserved.++module Math.NumberTheory.Factor where++import Math.NumberTheory.Prime+import Data.Either (lefts)+import Data.List (zip4)++-- Cohen, A Course in Computational Algebraic Number Theory, p488+++-- return (u,v,d) s.t ua+vb = d, with d = gcd a b+extendedEuclid a b+    | b == 0 = (1,0,a)+    | otherwise = let (q,r) = a `quotRem` b        -- a == qb + r+                      (s,t,d) = extendedEuclid b r -- s*b+t*r == d+                  in (t,s-q*t,d)                   -- s*b+t*(a-q*b) == d+++-- ELLIPTIC CURVE ARITHMETIC++data EllipticCurve = EC Integer Integer Integer deriving (Eq, Show)+-- EC p a b represents the curve y^2 == x^3+ax+b over Fp++data EllipticCurvePt = Inf | P Integer Integer deriving (Eq, Show)+-- P x y++isEltEC _ Inf = True+isEltEC (EC n a b) (P x y) = (y*y - x*x*x - a*x - b) `mod` n == 0+++-- Koblitz p34++-- Addition in an elliptic curve, with bailout if the arithmetic fails (giving a potential factor of n)+ecAdd _ Inf pt = Right pt+ecAdd _ pt Inf = Right pt+ecAdd (EC n a b) (P x1 y1) (P x2 y2)+    | x1 /= x2 = let (_,v,d) = extendedEuclid n ((x1-x2) `mod` n)  -- we're expecting d == 1, v == 1/(x1-x2) (mod n)+                     m = (y1-y2) * v `mod` n+                     x3 = (m*m - x1 - x2) `mod` n+                     y3 = (-y1 + m*(x1 - x3)) `mod` n+                  in if d == 1 then Right (P x3 y3) else Left d+    | x1 == x2 = if (y1 + y2) `mod` n == 0  -- includes the case y1 == y2 == 0+                 then Right Inf+                 else let (_,v,d) = extendedEuclid n ((2*y1) `mod` n)  -- we're expecting d == 1, v == 1/(2*y1) (mod n)+                          m = (3*x1*x1 + a) * v `mod` n+                          x3 = (m*m - 2*x1) `mod` n+                          y3 = (-y1 + m*(x1 - x3)) `mod` n+                      in if d == 1 then Right (P x3 y3) else Left d+-- Note that b isn't actually used anywhere++-- Note, only the final `mod` n calls when calculating x3, y3 are necessary+-- and the code is faster if the others are removed++-- Scalar multiplication in an elliptic curve+ecSmult _ 0 _ = Right Inf+ecSmult ec k pt | k > 0 = ecSmult' k pt Inf+    where -- ecSmult' k q p = k * q + p+          ecSmult' 0 _ p = Right p+          ecSmult' i q p = let p' = if odd i then ecAdd ec p q else Right p+                               q' = ecAdd ec q q+                           in case (p',q') of+                              (Right p'', Right q'') -> ecSmult' (i `div` 2) q'' p''+                              (Left _, _) -> p'+                              (_, Left _) -> q'+++-- ELLIPTIC CURVE FACTORISATION++-- We choose an elliptic curve E over Zn, and a point P on the curve+-- We then try to calculate kP, for suitably chosen k+-- What we are hoping is that at some stage we will fail because we can't invert an element in Zn+-- This will lead to finding a non-trivial factor of n+++discriminantEC a b = 4 * a * a * a + 27 * b * b++-- perform a sequence of scalar multiplications in the elliptic curve, hoping for a bailout+ecTrial ec@(EC n a b) ms pt+    | d == 1 = ecTrial' ms pt+    | otherwise = Left d+    where d = gcd n (discriminantEC a b)+          ecTrial' [] pt = Right pt+          ecTrial' (m:ms) pt = case ecSmult ec m pt of+                               Right pt' -> ecTrial' ms pt'+                               Left d -> Left d+-- In effect, we're calculating ecSmult ec (product ms) pt, but an m at a time++l n = exp (sqrt (log n * log (log n)))+-- L(n) is some sort of measure of the average smoothness of numbers up to n+-- # [x <= n | x is L(n)^a-smooth] = n L(n)^(-1/2a+o(1))  -- Cohen p 482++-- q is the largest prime we're looking for - normally sqrt n+-- the b figure here is from Cohen p488+multipliers q = [p' | p <- takeWhile (<= b) primes, let p' = last (takeWhile (<= b) (powers p))]+    where b = round ((l q) ** (1/sqrt 2))+          powers x = iterate (*x) x++findFactorECM n | gcd n 6 == 1 =+    let ms = multipliers (sqrt $ fromInteger n)+    in head $ filter ( (/= 0) . (`mod` n) ) $+       lefts [ecTrial (EC n a 1) ms (P 0 1) | a <- [1..] ]+    -- the filter is because d might be a multiple of n,+    -- for example if the problem was that the discriminant was divisible by n+++-- |List the prime factors of n (with multiplicity).+-- The algorithm uses trial division, followed by the elliptic curve method if necessary.+-- The running time increases with the size of the second largest prime factor of n.+-- It can find 10-digit prime factors in seconds, but can struggle with 20-digit prime factors.+pfactors :: Integer -> [Integer]+pfactors n | n > 0 = pfactors' n $ takeWhile (< 10000) primes+           | n < 0 = -1 : pfactors' (-n) (takeWhile (< 10000) primes)+    where pfactors' n (d:ds) | n == 1 = []+                             | n < d*d = [n]+                             | r == 0 = d : pfactors' q (d:ds)+                             | otherwise = pfactors' n ds+                             where (q,r) = quotRem n d+          pfactors' n [] = pfactors'' n+          pfactors'' n = if isMillerRabinPrime n then [n]+                         else let d = findFactorParallelECM n -- findFactorECM n+                              in merge (pfactors'' d) (pfactors'' (n `div` d))++merge (x:xs) (y:ys) =+    case compare x y of+    LT -> x : merge xs (y:ys)+    EQ -> x : y : merge xs ys+    GT -> y : merge (x:xs) ys+merge xs ys = xs ++ ys+++-- Cohen p489+-- find inverse of as mod n in parallel, or a non-trivial factor of n+parallelInverse n as = if d == 1 then Right bs else Left $ head [d' | a <- as, let d' = gcd a n, d' /= 1]+    where c:cs = reverse $ scanl (\x y -> x*y `mod` n) 1 as+          ds = scanl (\x y -> x*y `mod` n) 1 (reverse as)+          (u,_,d) = extendedEuclid c n+          bs = reverse [ u*nota `mod` n | nota <- zipWith (*) cs ds]+-- let m = length as+-- then the above code requires O(m) mod calls - in fact 3m-3 calls (?)++parallelEcAdd n ecs ps1 ps2 =+    case parallelInverse n (zipWith f ps1 ps2) of+    Right invs -> Right [g ec p1 p2 inv | (ec,p1,p2,inv) <- zip4 ecs ps1 ps2 invs]+    Left d -> Left d+    where f Inf pt = 1+          f pt Inf = 1+          f (P x1 y1) (P x2 y2) | x1 /= x2 = x1-x2 -- slightly faster not to `mod` n here+                                | x1 == x2 = 2*y1 -- slightly faster not to `mod` n here+          -- inverses = parallelInverse n $ zipWith f ps1 ps2+          g _ Inf pt _ = pt+          g _ pt Inf _ = pt+          g (EC n a b) (P x1 y1) (P x2 y2) inv+              | x1 /= x2 = let m = (y1-y2) * inv -- slightly faster not to `mod` n here+                               x3 = (m*m - x1 - x2) `mod` n+                               y3 = (-y1 + m*(x1 - x3)) `mod` n+                           in P x3 y3+              | x1 == x2 = if (y1 + y2) `elem` [0,n] -- `mod` n == 0  -- includes the case y1 == y2 == 0+                           then Inf+                           else let m = (3*x1*x1 + a) * inv -- slightly faster not to `mod` n here+                                    x3 = (m*m - 2*x1) `mod` n+                                    y3 = (-y1 + m*(x1 - x3)) `mod` n+                                 in P x3 y3++parallelEcSmult _ _ 0 pts = Right $ map (const Inf) pts+parallelEcSmult n ecs k pts | k > 0 = ecSmult' k pts (map (const Inf) pts)+    where -- ecSmult' k qs ps = k * qs + ps+          ecSmult' 0 _ ps = Right ps+          ecSmult' k qs ps = let ps' = if odd k then parallelEcAdd n ecs ps qs else Right ps+                                 qs' = parallelEcAdd n ecs qs qs+                             in case (ps',qs') of+                                (Right ps'', Right qs'') -> ecSmult' (k `div` 2) qs'' ps''+                                (Left _, _) -> ps'+                                (_, Left _) -> qs'++parallelEcTrial n ecs ms pts+    | all (==1) ds = ecTrial' ms pts+    | otherwise = Left $ head $ filter (/=1) ds+    where ds = [gcd n (discriminantEC a b) | EC n a b <- ecs]+          ecTrial' [] pts = Right pts+          ecTrial' (m:ms) pts = case parallelEcSmult n ecs m pts of+                                Right pts' -> ecTrial' ms pts'+                                Left d -> Left d++findFactorParallelECM n | gcd n 6 == 1 =+    let ms = multipliers (sqrt $ fromInteger n)+    in head $ filter ( (/= 0) . (`mod` n) ) $+       lefts [parallelEcTrial n [EC n (a+i) 1 | i <- [1..100]] ms (replicate 100 (P 0 1)) | a <- [0,100..] ]+-- 100 at a time is chosen heuristically.+

@@ -0,0 +1,118 @@+-- Copyright (c) 2006-2011, David Amos. All rights reserved.++{-# LANGUAGE NoMonomorphismRestriction #-}+++module Math.NumberTheory.Prime where++import System.Random+import System.IO.Unsafe+++isTrialDivisionPrime n+    | n > 1 = isNotDivisibleBy primes+    | otherwise = False+    where isNotDivisibleBy (d:ds) | d*d > n         = True+                                  | n `mod` d == 0  = False+                                  | otherwise       = isNotDivisibleBy ds++-- |A (lazy) list of the primes+primes :: [Integer]+primes = 2 : 3 : filter isTrialDivisionPrime (concat [ [m6-1,m6+1] | m6 <- [6,12..] ])++-- initial version. This isn't going to be very good if n has any "large" prime factors (eg > 10000)+pfactors1 n | n > 0 = pfactors' n primes+            | n < 0 = -1 : pfactors' (-n) primes+    where pfactors' n (d:ds) | n == 1 = []+                             | n < d*d = [n]+                             | r == 0 = d : pfactors' q (d:ds)+                             | otherwise = pfactors' n ds+                             where (q,r) = quotRem n d+++-- MILLER-RABIN TEST+-- Cohen, A Course in Computational Algebraic Number Theory, p422+-- Koblitz, A Course in Number Theory and Cryptography+++-- Let n-1 = 2^s * q, q odd+-- Then n is a strong pseudoprime to base b if+-- either b^q == 1 (mod n)+-- or b^(2^r * q) == -1 (mod n) for some 0 <= r < s+-- (For we know that if n is prime, then b^(n-1) == 1 (mod n)++isStrongPseudoPrime n b =+    let (s,q) = split2s 0 (n-1)  -- n-1 == 2^s * q, with q odd+    in isStrongPseudoPrime' n (s,q) b++isStrongPseudoPrime' n (s,q) b+    | b' == 1 = True+    | otherwise = n-1 `elem` squarings+    where b' = power_mod b q n     -- b' = b^q `mod` n+          squarings = take s $ iterate (\x -> x*x `mod` n) b' -- b^(2^r *q) for 0 <= r < s++-- split2s 0 m returns (s,t) such that 2^s * t == m, t odd+split2s s t = let (q,r) = t `quotRem` 2+              in if r == 0 then split2s (s+1) q else (s,t)++-- power_mod b t n == b^t mod n+power_mod b t n = powerMod' b 1 t+    where powerMod' x y 0 = y+          powerMod' x y t = powerMod' (x*x `mod` n) (if even t then y else x*y `mod` n) (t `div` 2)++isMillerRabinPrime' n+    | n >= 4 =+        let (s,q) = split2s 0 (n-1) -- n-1 == 2^s * q, with q odd+        in do g <- getStdGen+              let rs = randomRs (2,n-1) g+              return $ all (isStrongPseudoPrime' n (s,q)) (take 25 rs)+    | n >= 2 = return True+    | otherwise = return False+-- Cohen states that if we restrict our rs to single word numbers, we can use a more efficient powering algorithm++-- isMillerRabinPrime :: Integer -> Bool+isMillerRabinPrime n = unsafePerformIO (isMillerRabinPrime' n)+++-- |Is this number prime? The algorithm consists of using trial division to test for very small factors,+-- followed if necessary by the Miller-Rabin probabilistic test.+isPrime :: Integer -> Bool+isPrime n | n > 1 = isPrime' $ takeWhile (< 100) primes+          | otherwise = False+    where isPrime' (d:ds) | n < d*d = True+                          | otherwise = let (q,r) = quotRem n d+                                        in if r == 0 then False else isPrime' ds+          isPrime' [] = isMillerRabinPrime n+-- the < 100 is found heuristically to be about the point at which trial division stops being worthwhile++notPrime :: Integer -> Bool+notPrime = not . isPrime++-- |Given n, @prevPrime n@ returns the greatest p, p < n, such that p is prime+prevPrime :: Integer -> Integer+prevPrime n | n > 5 = head $ filter isPrime $ candidates+            | n < 3 = error "prevPrime: no previous primes"+            | n == 3 = 2+            | otherwise = 3+            where n6 = (n `div` 6) * 6+                  candidates = dropWhile (>= n) $ concat [ [m6+5,m6+1] | m6 <- [n6, n6-6..] ]++-- |Given n, @nextPrime n@ returns the least p, p > n, such that p is prime+nextPrime :: Integer -> Integer+nextPrime n | n < 2 = 2+            | n < 3 = 3+            | otherwise = head $ filter isPrime $ candidates+            where n6 = (n `div` 6) * 6+                  candidates = dropWhile (<= n) $ concat [ [m6+1,m6+5] | m6 <- [n6, n6+6..] ]++-- slightly better version. This is okay so long as n has at most one "large" prime factor (> 10000)+-- if it has more, it does at least tell you, via an error message, that it has run into difficulties+pfactors2 n | n > 0 = pfactors' n $ takeWhile (< 10000) primes+            | n < 0 = -1 : pfactors' (-n) (takeWhile (< 10000) primes)+    where pfactors' n (d:ds) | n == 1 = []+                             | n < d*d = [n]+                             | r == 0 = d : pfactors' q (d:ds)+                             | otherwise = pfactors' n ds+                             where (q,r) = quotRem n d+          pfactors' n [] = if isMillerRabinPrime n then [n] else error ("pfactors2: can't factor " ++ show n)+

@@ -4,8 +4,9 @@ 
 import Data.List as L
 
-import Math.Algebra.Field.Base
-import Math.Algebra.Field.Extension hiding ( (<+>), (<*>) )
+-- import Math.Algebra.Field.Base
+-- import Math.Algebra.Field.Extension hiding ( (<+>), (<*>) )
+import Math.Core.Field
 import Math.Algebra.LinearAlgebra
 
 import Math.Algebra.Group.PermutationGroup hiding (fromList)
@@ -13,14 +14,13 @@ import Math.Algebra.Group.RandomSchreierSims as RSS
 
 import Math.Combinatorics.FiniteGeometry (ptsAG)
--- import ClassicalChevalleyGroup (ptsAG)
 
 
 -- Follows Conway's notation
 
 -- The octonion xinf + x0 i0 + x1 i1 + ... + x6 i6
 -- is represented as O [(-1,xinf),(0,x0),(1,x1),...,(6,x6)]
--- where a list element may be omitted if one if the coefficient is zero
+-- where a list element may be omitted if the coefficient is zero
 
 
 newtype Octonion k = O [(Int,k)] deriving (Eq, Ord)
@@ -176,7 +176,8 @@ 
 gamma4s = [x | x <- unitImagOctonions f4, isOrthogonal (O [(0,1)]) x, isOrthogonal (O [(1,1)]) x, isOrthogonal (O [(3,1)]) x]
 
-gamma4 = autFrom (O [(0,1::F4)]) (O [(1,1)]) (O [(5,embed x),(6,embed $ 1+x)])
+-- gamma4 = autFrom (O [(0,1::F4)]) (O [(1,1)]) (O [(5,embed x),(6,embed $ 1+x)])
+gamma4 = autFrom (O [(0,1::F4)]) (O [(1,1)]) (O [(5,a4),(6,1+a4)])
 
 alpha4' = fromPairs [(x, x %^ alpha4) | x <- unitImagOctonions f4]
 beta4' = fromPairs [(x, x %^ beta4) | x <- unitImagOctonions f4]

@@ -0,0 +1,80 @@+-- Copyright (c) 2011, David Amos. All rights reserved.++module Math.Test.TCommutativeAlgebra.TPolynomial where++import Test.HUnit++import Data.List as L+import Math.Core.Field+-- import Math.Algebras.VectorSpace+import Math.CommutativeAlgebra.Polynomial+++testlistPolynomial = TestList [+    testlistMonomialOrders1,+    testlistMonomialOrders2,+    testlistOverFiniteField,+    testlistEval,+    testlistSubst+    ]+++testcaseMonomialOrder1 desc monomials expected =+    TestCase $ assertEqual desc expected (L.sort monomials)++testlistMonomialOrders1 = TestList [+   let [x,y,z] = map lexvar ["x","y","z"]+   in testcaseMonomialOrder1 "Lex" [1,x,y,z,x^2,x*y,x*z,y^2,y*z,z^2] [x^2,x*y,x*z,x,y^2,y*z,y,z^2,z,1],+   let [x,y,z] = map glexvar ["x","y","z"]+   in testcaseMonomialOrder1 "Glex" [1,x,y,z,x^2,x*y,x*z,y^2,y*z,z^2] [x^2,x*y,x*z,y^2,y*z,z^2,x,y,z,1],+   let [x,y,z] = map grevlexvar ["x","y","z"]+   in testcaseMonomialOrder1 "Grevlex" [1,x,y,z,x^2,x*y,x*z,y^2,y*z,z^2] [x^2,x*y,y^2,x*z,y*z,z^2,x,y,z,1]+   ]+++testcaseMonomialOrder2 desc poly string =+    TestCase $ assertEqual desc string (show poly)++testlistMonomialOrders2 = TestList [+    let [x,y,z] = map lexvar ["x","y","z"]+    in testcaseMonomialOrder2 "Lex" ((x+y+z+1)^3)+       "x^3+3x^2y+3x^2z+3x^2+3xy^2+6xyz+6xy+3xz^2+6xz+3x+y^3+3y^2z+3y^2+3yz^2+6yz+3y+z^3+3z^2+3z+1",+    let [x,y,z] = map glexvar ["x","y","z"]+    in testcaseMonomialOrder2 "Glex" ((x+y+z+1)^3)+       "x^3+3x^2y+3x^2z+3xy^2+6xyz+3xz^2+y^3+3y^2z+3yz^2+z^3+3x^2+6xy+6xz+3y^2+6yz+3z^2+3x+3y+3z+1",+    let [x,y,z] = map grevlexvar ["x","y","z"]+    in testcaseMonomialOrder2 "Grevlex" ((x+y+z+1)^3)+       "x^3+3x^2y+3xy^2+y^3+3x^2z+6xyz+3y^2z+3xz^2+3yz^2+z^3+3x^2+6xy+3y^2+6xz+6yz+3z^2+3x+3y+3z+1"+    ]+++testcaseOverFiniteField desc input expected = TestCase $ assertEqual desc expected input++testlistOverFiniteField = TestList [+    let [x,y,z] = map var ["x","y","z"] :: [GlexPoly F3 String] in+    testcaseOverFiniteField "F3" ((x+y+z)^3) (x^3+y^3+z^3),+    let [x,y,z] = map var ["x","y","z"] :: [GlexPoly F5 String] in+    testcaseOverFiniteField "F5" ((x+y+z)^5) (x^5+y^5+z^5)+    ]+++testcaseEval poly point value =+    TestCase $ assertEqual "Eval" value (eval poly point)++testlistEval =+    let [x,y,z] = map glexvar ["x","y","z"] in TestList [+        testcaseEval (x^2+y^2-z^2) [(x,3),(y,4),(z,5)] 0,+        testcaseEval (z-1) [(x,100),(y,-100),(z,1)] 0+        ]+++testcaseSubst poly subs poly' =+    TestCase $ assertEqual "Subst" poly' (subst poly subs)++testlistSubst =+    let [x,y,z,x',y',z'] = map glexvar ["x","y","z","x'","y'","z'"] in TestList [+        testcaseSubst (x^2+y^2-z^2) [(x,(x'-z')/2),(y,y'),(z,(x'+z')/2)] (y'^2-x'*z')+        ]+++-- The division algorithm is adequately tested by the GroebnerBasis tests

@@ -0,0 +1,119 @@+-- Copyright (c) 2011, David Amos. All rights reserved.++module Math.Test.TNumberTheory.TPrimeFactor where++import Data.List ( (\\) )++import Math.NumberTheory.Prime+import Math.NumberTheory.Factor++import Test.HUnit+++testlistPrimeFactor = TestList [+    testlistSmallPrimes,+    testlistMillerRabin,+    testlistMersennePrimes,+    testlistMersenneNonPrimes,+    testlistFermatPrimes,+    testlistFermatNonPrimes,+    testlistCullenPrimes,+    testlistCullenNonPrimes,+    testlistWoodallPrimes,+    testlistWoodallNonPrimes,+    testlistWagstaffPrimes,+    testlistWagstaffNonPrimes,+    testlistFermatFactors,+    testlistNextPrime,+    testlistPrevPrime,+    testlistConsistentFactors,+    testlistFactorOrder+    ]+++testlistSmallPrimes = TestList [+    TestCase $ assertEqual "small primes"+               [False,True,True,False,True,False,True,False,False,False]+               (map isPrime [1..10]),+    TestCase $ assertBool "negative primes" (all notPrime [-10..0])+    ]++testcaseMillerRabin n = TestCase $ assertEqual ("MillerRabin " ++ show n)+                                   (isTrialDivisionPrime n) (isMillerRabinPrime n)++testlistMillerRabin = TestList $ map testcaseMillerRabin $ [1 :: Integer ..1000] ++ [10^6..10^6+10^3]++-- Source: http://en.wikipedia.org/wiki/Mersenne_prime+testlistMersennePrimes = TestList+    [TestCase $ assertBool ("Mersenne " ++ show p) (isPrime (2^p-1)) |+     p <- [2,3,5,7,13,17,19,31,61,89,107,127,521] ]++testlistMersenneNonPrimes = TestList+    [TestCase $ assertBool ("Mersenne " ++ show p) (notPrime (2^p-1)) |+     p <- [11,23,29,37,41,43,47,53,59,67,71,73,79,83,97,101,103,109,113] ]++-- http://en.wikipedia.org/wiki/Fermat_prime+testlistFermatPrimes = TestList+    [TestCase $ assertBool ("Fermat " ++ show n) (isPrime (2^2^n + 1)) | n <- [0..4] ]++testlistFermatNonPrimes = TestList+    [TestCase $ assertBool ("Fermat " ++ show n) (notPrime (2^2^n + 1)) | n <- [5..10] ]++-- http://en.wikipedia.org/wiki/Cullen_number+testlistCullenPrimes = TestList+    [TestCase $ assertBool ("Cullen " ++ show n) (isPrime (n * 2^n + 1)) | n <- [141] ]++testlistCullenNonPrimes = TestList+    [TestCase $ assertBool ("Cullen " ++ show n) (notPrime (n * 2^n + 1)) | n <- [2..100] ]++-- http://en.wikipedia.org/wiki/Woodall_number+testlistWoodallPrimes = TestList+    [TestCase $ assertBool ("Woodall " ++ show n) (isPrime (n * 2^n - 1)) |+     n <- [2,3,6,30,75,81,115,123,249,362,384] ]++testlistWoodallNonPrimes = TestList+    [TestCase $ assertBool ("Woodall " ++ show n) (notPrime (n * 2^n - 1)) |+     n <- [2..100] \\ [2,3,6,30,75,81,115,123,249,362,384] ]++-- http://en.wikipedia.org/wiki/Wagstaff_prime+testlistWagstaffPrimes = TestList+    [TestCase $ assertBool ("Wagstaff " ++ show n) (isPrime ((2^n + 1) `div` 3)) |+     n <- [3,5,7,11,13,17,19,23,31,43,61,79,101,127,167,191,199] ]++testlistWagstaffNonPrimes = TestList+    [TestCase $ assertBool ("Wagstaff " ++ show n) (notPrime ((2^n + 1) `div` 3)) |+     n <- takeWhile (<200) primes \\ [3,5,7,11,13,17,19,23,31,43,61,79,101,127,167,191,199] ]+++testcaseKnownFactors n ps = TestCase $ assertEqual (show n) ps (pfactors n)++testlistFermatFactors = TestList [+    testcaseKnownFactors (2^2^5+1) [641, 6700417],+    testcaseKnownFactors (2^2^6+1) [274177, 67280421310721]+    ]+++testlistNextPrime = TestList+    [TestCase $ assertEqual (show n) p (nextPrime n) |+     (n,p) <- [(0,2),(1,2),(2,3),(3,5),(4,5),(5,7),(6,7),(7,11),(8,11),(9,11),+               (10,11),(11,13),(12,13),(13,17),(14,17),(15,17),(16,17),(17,19),(18,19),(19,23),+               (20,23),(21,23),(22,23),(23,29),(24,29),(25,29),(26,29),(27,29),(28,29),(29,31),(30,31)] ]++testlistPrevPrime = TestList+    [TestCase $ assertEqual (show n) p (prevPrime n) |+     (n,p) <- [(3,2),(4,3),(5,3),(6,5),(7,5),(8,7),(9,7),+               (10,7),(11,7),(12,11),(13,11),(14,13),(15,13),(16,13),(17,13),(18,17),(19,17),+               (20,19),(21,19),(22,19),(23,19),(24,23),(25,23),(26,23),(27,23),(28,23),(29,23),(30,29)] ]+++testcaseConsistentFactors n = TestCase $ assertBool (show n) (product (pfactors n) == n)++testlistConsistentFactors = TestList $ map testcaseConsistentFactors $ [10^6..10^6+10^2] ++ [10^16..10^16+10^2]+++testlistFactorOrder =+    let f1 = nextPrime 50000; f2 = nextPrime 70000 in TestList [+    TestCase (assertEqual "" [2,2,2,3,3,5] (pfactors (2^3*3^2*5))),+    TestCase (assertEqual "" [f1,f1,f1,f2,f2] (pfactors (f1^3*f2^2))),+    TestCase (assertEqual "" [2,2,2,3,3,5,f1,f1,f1,f2,f2] (pfactors (2^3*3^2*5*f1^3*f2^2)))+    ]

@@ -15,7 +15,9 @@ import Math.Test.TCombinatorics.TIncidenceAlgebra
 import Math.Test.TCombinatorics.TMatroid
 import Math.Test.TCombinatorics.TPoset
+import Math.Test.TCommutativeAlgebra.TPolynomial
 import Math.Test.TCommutativeAlgebra.TGroebnerBasis
+import Math.Test.TNumberTheory.TPrimeFactor
 import Math.Test.TProjects.TMiniquaternionGeometry
 
 
@@ -46,5 +48,7 @@     testlistIncidenceAlgebra,
     testlistMatroid,
     testlistPoset,
-    testlistGroebnerBasis
-    ]+    testlistPolynomial,
+    testlistGroebnerBasis,
+    testlistPrimeFactor
+    ]