diff --git a/HaskellForMaths.cabal b/HaskellForMaths.cabal --- a/HaskellForMaths.cabal +++ b/HaskellForMaths.cabal @@ -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, diff --git a/Math/Algebras/Octonions.hs b/Math/Algebras/Octonions.hs new file mode 100644 --- /dev/null +++ b/Math/Algebras/Octonions.hs @@ -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] diff --git a/Math/Algebras/Quaternions.hs b/Math/Algebras/Quaternions.hs --- a/Math/Algebras/Quaternions.hs +++ b/Math/Algebras/Quaternions.hs @@ -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) diff --git a/Math/Algebras/Structures.hs b/Math/Algebras/Structures.hs --- a/Math/Algebras/Structures.hs +++ b/Math/Algebras/Structures.hs @@ -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) diff --git a/Math/CommutativeAlgebra/GroebnerBasis.hs b/Math/CommutativeAlgebra/GroebnerBasis.hs --- a/Math/CommutativeAlgebra/GroebnerBasis.hs +++ b/Math/CommutativeAlgebra/GroebnerBasis.hs @@ -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]/\. +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) diff --git a/Math/CommutativeAlgebra/Polynomial.hs b/Math/CommutativeAlgebra/Polynomial.hs --- a/Math/CommutativeAlgebra/Polynomial.hs +++ b/Math/CommutativeAlgebra/Polynomial.hs @@ -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. diff --git a/Math/NumberTheory/Factor.hs b/Math/NumberTheory/Factor.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Factor.hs @@ -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. + diff --git a/Math/NumberTheory/Prime.hs b/Math/NumberTheory/Prime.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Prime.hs @@ -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) + diff --git a/Math/Projects/ChevalleyGroup/Exceptional.hs b/Math/Projects/ChevalleyGroup/Exceptional.hs --- a/Math/Projects/ChevalleyGroup/Exceptional.hs +++ b/Math/Projects/ChevalleyGroup/Exceptional.hs @@ -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] diff --git a/Math/Test/TCommutativeAlgebra/TPolynomial.hs b/Math/Test/TCommutativeAlgebra/TPolynomial.hs new file mode 100644 --- /dev/null +++ b/Math/Test/TCommutativeAlgebra/TPolynomial.hs @@ -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 diff --git a/Math/Test/TNumberTheory/TPrimeFactor.hs b/Math/Test/TNumberTheory/TPrimeFactor.hs new file mode 100644 --- /dev/null +++ b/Math/Test/TNumberTheory/TPrimeFactor.hs @@ -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))) + ] diff --git a/Math/Test/TestAll.hs b/Math/Test/TestAll.hs --- a/Math/Test/TestAll.hs +++ b/Math/Test/TestAll.hs @@ -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 + ]