arithmetic 1.2 → 1.3

raw patch · 7 files changed

+526/−22 lines, 7 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• arithmetic.cabal +3/−1
• src/Arithmetic/Modular.hs +29/−15
• src/Arithmetic/Polynomial.hs +237/−0
• src/Arithmetic/Prime/Factor.hs +12/−0
• src/Arithmetic/Quadratic.hs +2/−2
• src/Arithmetic/Ring.hs +53/−0
• src/Test.hs +190/−4

@@ -1,5 +1,5 @@ name: arithmetic-version: 1.2+version: 1.3 category: Number Theory synopsis: Natural number arithmetic license: MIT@@ -32,11 +32,13 @@     Arithmetic.Lucas,     Arithmetic.Modular,     Arithmetic.Montgomery,+    Arithmetic.Polynomial,     Arithmetic.Prime,     Arithmetic.Prime.Factor,     Arithmetic.Prime.Sieve,     Arithmetic.Quadratic,     Arithmetic.Random,+    Arithmetic.Ring,     Arithmetic.Utility,     Arithmetic.Utility.Heap,     Arithmetic.Williams

@@ -11,9 +11,10 @@ where  import OpenTheory.Primitive.Natural-import OpenTheory.Natural.Divides+import qualified OpenTheory.Natural.Divides as Divides  import Arithmetic.Utility+import qualified Arithmetic.Ring as Ring  normalize :: Natural -> Natural -> Natural normalize n x = x `mod` n@@ -21,28 +22,44 @@ add :: Natural -> Natural -> Natural -> Natural add n x y = normalize n (x + y) -double :: Natural -> Natural -> Natural-double n x = add n x x- negate :: Natural -> Natural -> Natural negate n x =     if y == 0 then y else n - y   where     y = normalize n x +multiply :: Natural -> Natural -> Natural -> Natural+multiply n x y = normalize n (x * y)++divide :: Natural -> Natural -> Natural -> Maybe Natural+divide n x y =+    if g == n then if Divides.divides n x then Just 0 else Nothing+    else if Divides.divides g x then Just (multiply n (x `div` g) s)+    else Nothing+  where+    (g,(s,_)) = Divides.egcd y n  -- s * y == g (mod n)++ring :: Natural -> Ring.Ring Natural+ring n =+    Ring.Ring {Ring.fromNatural = normalize n,+               Ring.add = add n,+               Ring.negate = Arithmetic.Modular.negate n,+               Ring.multiply = multiply n,+               Ring.divide = divide n}++double :: Natural -> Natural -> Natural+double = Ring.double . ring+ subtract :: Natural -> Natural -> Natural -> Natural subtract n x y =     if y <= x then normalize n (x - y)     else Arithmetic.Modular.negate n (y - x) -multiply :: Natural -> Natural -> Natural -> Natural-multiply n x y = normalize n (x * y)- square :: Natural -> Natural -> Natural-square n x = multiply n x x+square = Ring.square . ring  exp :: Natural -> Natural -> Natural -> Natural-exp n = multiplyExponential (multiply n) 1+exp = Ring.exp . ring  exp2 :: Natural -> Natural -> Natural -> Natural exp2 n x k = if k == 0 then normalize n x else functionPower (square n) k x@@ -53,10 +70,7 @@     else if g == 1 then Just s     else Nothing   where-    (g,(s,_)) = egcd x n+    (g,(s,_)) = Divides.egcd x n -divide :: Natural -> Natural -> Natural -> Maybe Natural-divide n x y =-    case invert n y of-      Nothing -> Nothing-      Just z -> Just (multiply n x z)+divides :: Natural -> Natural -> Natural -> Bool+divides n a b = Divides.divides (gcd n b) (gcd n a)

@@ -0,0 +1,237 @@+{- |+module: Arithmetic.Polynomial+description: Polynomial arithmetic+license: MIT++maintainer: Joe Leslie-Hurd <joe@gilith.com>+stability: provisional+portability: portable+-}+module Arithmetic.Polynomial+where++import OpenTheory.Primitive.Natural+import OpenTheory.List+import Data.List as List+import Data.Maybe as Maybe++import qualified Arithmetic.Ring as Ring++data Polynomial a =+     Polynomial+       {carrier :: Ring.Ring a,+        coefficients :: [a]}++instance (Eq a, Show a) => Show (Polynomial a) where+  show p =+      if null ps then "0"+      else List.intercalate " + " ps+    where+      r = carrier p+      z = Ring.zero r+      o = Ring.one r+      ps = showC (0 :: Natural) (coefficients p)++      showC _ [] = []+      showC k (x : xs) = showM x k ++ showC (k + 1) xs++      showM x k =+          if x == z then []+          else [(if k /= 0 && x == o then "" else show x) +++                (if k == 0 then ""+                 else ("x" ++ (if k == 1 then "" else "^" ++ show k)))]++fromCoefficients :: Eq a => Ring.Ring a -> [a] -> Polynomial a+fromCoefficients r cs =+    Polynomial+      {carrier = r,+       coefficients = norm cs}+  where+    z = Ring.zero r++    zcons x xs = if null xs && x == z then [] else x : xs++    norm [] = []+    norm (x : xs) = zcons x (norm xs)++zero :: Ring.Ring a -> Polynomial a+zero r =+    Polynomial+      {carrier = r,+       coefficients = []}++isZero :: Polynomial a -> Bool+isZero = null . coefficients++constant :: Eq a => Ring.Ring a -> a -> Polynomial a+constant r x = fromCoefficients r [x]++destConstant :: Polynomial a -> Maybe a+destConstant p =+    case coefficients p of+      [] -> Just (Ring.zero r)+      [c] -> Just c+      _ -> Nothing+  where+    r = carrier p++isConstant :: Polynomial a -> Bool+isConstant = Maybe.isJust . destConstant++fromNatural :: Eq a => Ring.Ring a -> Natural -> Polynomial a+fromNatural r = constant r . Ring.fromNatural r++one :: Eq a => Ring.Ring a -> Polynomial a+one r = constant r (Ring.one r)++multiplyByPower :: Polynomial a -> Natural -> Polynomial a+multiplyByPower p k =+    if k == 0 || null cs then p+    else p {coefficients = replicate (fromIntegral k) z ++ cs}+  where+    r = carrier p+    z = Ring.zero r+    cs = coefficients p++monomial :: Eq a => Ring.Ring a -> a -> Natural -> Polynomial a+monomial r x = multiplyByPower (constant r x)++variablePower :: Eq a => Ring.Ring a -> Natural -> Polynomial a+variablePower r = monomial r (Ring.one r)++variable :: Eq a => Ring.Ring a -> Polynomial a+variable r = variablePower r 1++degree :: Polynomial a -> Natural+degree = naturalLength . coefficients++leadingCoefficient :: Polynomial a -> Maybe a+leadingCoefficient p =+    case coefficients p of+      [] -> Nothing+      cs -> Just (last cs)++nthCoefficient :: Polynomial a -> Natural -> a+nthCoefficient p k =+    if k < degree p then coefficients p !! (fromIntegral k)+    else Ring.zero (carrier p)++isMonic :: Eq a => Polynomial a -> Bool+isMonic p =+    case leadingCoefficient p of+      Nothing -> False+      Just c -> c == Ring.one (carrier p)++-- Horner's method+evaluate :: Polynomial a -> a -> a+evaluate p x =+    foldr eval (Ring.zero r) (coefficients p)+  where+    r = carrier p+    eval c z = Ring.add r c (Ring.multiply r x z)++addCoefficients :: Ring.Ring a -> [a] -> [a] -> [a]+addCoefficients r =+    addc+  where+    addc [] ys = ys+    addc xs [] = xs+    addc (x : xs) (y : ys) = Ring.add r x y : addc xs ys++add :: Eq a => Polynomial a -> Polynomial a -> Polynomial a+add p q =+    fromCoefficients r (addCoefficients r ps qs)+  where+    r = carrier p+    ps = coefficients p+    qs = coefficients q++negate :: Polynomial a -> Polynomial a+negate p =+    Polynomial+      {carrier = r,+       coefficients = map (Ring.negate r) pl}+  where+    r = carrier p+    pl = coefficients p++multiply :: Eq a => Polynomial a -> Polynomial a -> Polynomial a+multiply p q =+    case coefficients q of+      [] -> zero r+      qh : qt ->+          fromCoefficients r (foldr multc [] (coefficients p))+        where+          z = Ring.zero r++          madd pc cs = addCoefficients r (map (Ring.multiply r pc) qt) cs++          multc pc cs =+              if pc == z then z : cs+              else Ring.multiply r pc qh : madd pc cs+  where+    r = carrier p++multiplyByScalar :: Eq a => Polynomial a -> a -> Polynomial a+multiplyByScalar p x = multiply p (constant (carrier p) x)++invert :: Polynomial a -> Maybe (Polynomial a)+invert p =+    case coefficients p of+      [x] -> case Ring.invert r x of+               Nothing -> Nothing+               Just y -> Just (Polynomial {carrier = r, coefficients = [y]})+      _ -> Nothing+  where+    r = carrier p++subtract :: Eq a => Polynomial a -> Polynomial a -> Polynomial a+subtract p = Ring.subtract (ring (carrier p)) p++quotientRemainder :: Eq a => Polynomial a -> Polynomial a ->+                     Maybe (Polynomial a, Polynomial a)+quotientRemainder p q =+    if d_p < d_q then Just (zero r, p)+    else case leadingCoefficient q of+           Nothing -> Nothing+           Just q_m ->+               go [] p (d_p - d_q)+             where+               sub f k =+                   if f_m == z then Just (z,f)+                   else case Ring.divide r f_m q_m of+                          Nothing -> Nothing+                          Just c ->+                              Just (c, Arithmetic.Polynomial.subtract f g)+                            where+                              g = multiplyByPower (multiplyByScalar q c) k+                 where+                   f_m = nthCoefficient f (d_q + k)++               go cs f k =+                   case sub f k of+                     Nothing -> Nothing+                     Just (c,g) -> go' (c : cs) g k++               go' cs f k =+                   if k == 0 then Just (fromCoefficients r cs, f)+                   else go cs f (k - 1)+  where+    r = carrier p+    z = Ring.zero r+    d_p = degree p+    d_q = degree q++divide :: Eq a => Polynomial a -> Polynomial a -> Maybe (Polynomial a)+divide p q =+    case quotientRemainder p q of+      Nothing -> Nothing+      Just (x,y) -> if isZero y then Just x else Nothing++ring :: Eq a => Ring.Ring a -> Ring.Ring (Polynomial a)+ring r =+    Ring.Ring {Ring.fromNatural = fromNatural r,+               Ring.add = add,+               Ring.negate = Arithmetic.Polynomial.negate,+               Ring.multiply = multiply,+               Ring.divide = divide}

@@ -22,6 +22,9 @@  newtype Factor = Factor {unFactor :: Map.Map Natural Natural} +primePowers :: Factor -> [(Natural,Natural)]+primePowers = Map.toList . unFactor+ one :: Factor one = Factor {unFactor = Map.empty} @@ -51,6 +54,15 @@  isPrime :: Factor -> Bool isPrime = Maybe.isJust . destPrime++destRSA :: Factor -> Maybe (Natural,Natural)+destRSA f =+    case primePowers f of+      [(p,1),(q,1)] -> Just (p,q)+      _ -> Nothing++isRSA :: Factor -> Bool+isRSA = Maybe.isJust . destRSA  multiply :: Factor -> Factor -> Factor multiply f1 f2 =

@@ -23,8 +23,8 @@   where     bisect l u =         if m == l then l-	else if m * m <= n then bisect m u-	else bisect l m+        else if m * m <= n then bisect m u+        else bisect l m       where         m = (l + u) `div` 2 

@@ -0,0 +1,53 @@+{- |+module: Arithmetic.Ring+description: An abstract ring type+license: MIT++maintainer: Joe Leslie-Hurd <joe@gilith.com>+stability: provisional+portability: portable+-}+module Arithmetic.Ring+where++import OpenTheory.Primitive.Natural+import qualified Data.Maybe as Maybe++import Arithmetic.Utility++data Ring a = Ring+    {fromNatural :: Natural -> a,+     add :: a -> a -> a,+     negate :: a -> a,+     multiply :: a -> a -> a,+     divide :: a -> a -> Maybe a}++zero :: Ring a -> a+zero r = fromNatural r 0++one :: Ring a -> a+one r = fromNatural r 1++two :: Ring a -> a+two r = fromNatural r 2++double :: Ring a -> a -> a+double r x = add r x x++subtract :: Ring a -> a -> a -> a+subtract r x y = add r x (Arithmetic.Ring.negate r y)++square :: Ring a -> a -> a+square r x = multiply r x x++exp :: Ring a -> a -> Natural -> a+exp r = multiplyExponential (multiply r) (one r)++exp2 :: Ring a -> a -> Natural -> a+exp2 r x k = functionPower (square r) k x++divides :: Ring a -> a -> a -> Bool+divides r x y = Maybe.isJust (divide r x y)++invert :: Ring a -> a -> Maybe a+invert r x = divide r (one r) x

@@ -11,12 +11,13 @@   ( main ) where -import qualified Test.QuickCheck as QuickCheck import OpenTheory.Primitive.Natural import OpenTheory.Natural import OpenTheory.Natural.Divides import qualified OpenTheory.Natural.Bits as Bits+import qualified Data.Maybe as Maybe import qualified OpenTheory.Natural.Prime as Prime+import qualified Test.QuickCheck as QuickCheck import qualified OpenTheory.Primitive.Random as Random import qualified OpenTheory.Natural.Uniform as Uniform @@ -26,7 +27,9 @@ import qualified Arithmetic.Prime.Factor as Factor import qualified Arithmetic.Modular as Modular import qualified Arithmetic.Montgomery as Montgomery+import qualified Arithmetic.Polynomial as Polynomial import qualified Arithmetic.Quadratic as Quadratic+import qualified Arithmetic.Ring as Ring import qualified Arithmetic.Williams as Williams  propPrimes :: Natural -> Bool@@ -70,13 +73,43 @@     a = Uniform.random n rnd     b = Modular.negate n a +propModularSubtract :: Natural -> Natural -> Natural -> Bool+propModularSubtract np a b =+    Modular.subtract n a b ==+    Ring.subtract r (Ring.fromNatural r a) (Ring.fromNatural r b)+  where+    n = np + 1+    r = Modular.ring n++propModularExp2 :: Natural -> Natural -> Natural -> Bool+propModularExp2 np x k =+    Modular.exp2 n x k == Ring.exp2 r (Ring.fromNatural r x) k+  where+    n = np + 1+    r = Modular.ring n++propModularDivide :: Natural -> Natural -> Natural -> Bool+propModularDivide np a b =+    case Modular.divide n a b of+      Nothing -> not (Modular.divides n a b)+      Just c -> Modular.multiply n b c == Modular.normalize n a && c < n+  where+    n = np + 1++propModularDivides :: Natural -> Natural -> Natural -> Bool+propModularDivides np a b =+    Modular.divides n a b ==+    Ring.divides r (Ring.fromNatural r a) (Ring.fromNatural r b)+  where+    n = np + 1+    r = Modular.ring n+ propModularInvert :: Natural -> Natural -> Bool propModularInvert np a =-    case Modular.invert n a of-      Nothing -> gcd n a /= 1-      Just b -> Modular.multiply n a b == Modular.normalize n 1 && b < n+    Modular.invert n a == Ring.invert r (Ring.fromNatural r a)   where     n = np + 1+    r = Modular.ring n  propFermat :: Natural -> Random.Random -> Bool propFermat pp rnd =@@ -352,6 +385,142 @@   where     n = 2 * np + 5 +propPolynomialConstantDegree :: Natural -> Natural -> Bool+propPolynomialConstantDegree np cp =+    Polynomial.degree (Polynomial.constant r c) == if c == 0 then 0 else 1+  where+    n = np + 1+    r = Modular.ring n+    c = Ring.fromNatural r cp++propPolynomialFromNaturalDegree :: Natural -> Natural -> Bool+propPolynomialFromNaturalDegree np c =+    Polynomial.degree (Polynomial.fromNatural r c) ==+    if divides n c then 0 else 1+  where+    n = np + 1+    r = Modular.ring n++propPolynomialAddDegree :: Natural -> [Natural] -> [Natural] -> Bool+propPolynomialAddDegree np ps qs =+    if d_p == d_q then d_pq <= d_p+    else d_pq == max d_p d_q+  where+    n = np + 1+    r = Modular.ring n+    p = Polynomial.fromCoefficients r (map (Ring.fromNatural r) ps)+    q = Polynomial.fromCoefficients r (map (Ring.fromNatural r) qs)+    d_p = Polynomial.degree p+    d_q = Polynomial.degree q+    pq = Polynomial.add p q+    d_pq = Polynomial.degree pq++propPolynomialNegateDegree :: Natural -> [Natural] -> Bool+propPolynomialNegateDegree np ps =+    Polynomial.degree (Polynomial.negate p) == Polynomial.degree p+  where+    n = np + 1+    r = Modular.ring n+    p = Polynomial.fromCoefficients r (map (Ring.fromNatural r) ps)++propPolynomialMultiplyDegree :: Natural -> [Natural] -> [Natural] -> Bool+propPolynomialMultiplyDegree np ps qs =+    if d_p == 0 || d_q == 0 then d_pq == 0+    else d_pq + 1 <= d_p + d_q+  where+    n = np + 1+    r = Modular.ring n+    p = Polynomial.fromCoefficients r (map (Ring.fromNatural r) ps)+    q = Polynomial.fromCoefficients r (map (Ring.fromNatural r) qs)+    d_p = Polynomial.degree p+    d_q = Polynomial.degree q+    pq = Polynomial.multiply p q+    d_pq = Polynomial.degree pq++propPolynomialConstantEvaluate :: Natural -> Natural -> Natural -> Bool+propPolynomialConstantEvaluate np cp xp =+    Polynomial.evaluate (Polynomial.constant r c) x == c+  where+    n = np + 1+    r = Modular.ring n+    c = Ring.fromNatural r cp+    x = Ring.fromNatural r xp++propPolynomialFromNaturalEvaluate :: Natural -> Natural -> Natural -> Bool+propPolynomialFromNaturalEvaluate np c xp =+    Polynomial.evaluate (Polynomial.fromNatural r c) x ==+    Ring.fromNatural r c+  where+    n = np + 1+    r = Modular.ring n+    x = Ring.fromNatural r xp++propPolynomialAddEvaluate ::+    Natural -> [Natural] -> [Natural] -> Natural -> Bool+propPolynomialAddEvaluate np ps qs xp =+    Polynomial.evaluate (Polynomial.add p q) x ==+    Ring.add r (Polynomial.evaluate p x) (Polynomial.evaluate q x)+  where+    n = np + 1+    r = Modular.ring n+    p = Polynomial.fromCoefficients r (map (Ring.fromNatural r) ps)+    q = Polynomial.fromCoefficients r (map (Ring.fromNatural r) qs)+    x = Ring.fromNatural r xp++propPolynomialNegateEvaluate :: Natural -> [Natural] -> Natural -> Bool+propPolynomialNegateEvaluate np ps xp =+    Polynomial.evaluate (Polynomial.negate p) x ==+    Ring.negate r (Polynomial.evaluate p x)+  where+    n = np + 1+    r = Modular.ring n+    p = Polynomial.fromCoefficients r (map (Ring.fromNatural r) ps)+    x = Ring.fromNatural r xp++propPolynomialMultiplyEvaluate ::+    Natural -> [Natural] -> [Natural] -> Natural -> Bool+propPolynomialMultiplyEvaluate np ps qs xp =+    Polynomial.evaluate (Polynomial.multiply p q) x ==+    Ring.multiply r (Polynomial.evaluate p x) (Polynomial.evaluate q x)+  where+    n = np + 1+    r = Modular.ring n+    p = Polynomial.fromCoefficients r (map (Ring.fromNatural r) ps)+    q = Polynomial.fromCoefficients r (map (Ring.fromNatural r) qs)+    x = Ring.fromNatural r xp++propPolynomialQuotientRemainder :: Natural -> [Natural] -> [Natural] -> Bool+propPolynomialQuotientRemainder np ps qs =+    case Polynomial.quotientRemainder p q of+      Nothing -> True+      Just (a,b) -> Polynomial.coefficients p ==+                    Polynomial.coefficients (Polynomial.add (Polynomial.multiply a q) b)+  where+    n = np + 1+    r = Modular.ring n+    p = Polynomial.fromCoefficients r (map (Ring.fromNatural r) ps)+    q = Polynomial.fromCoefficients r (map (Ring.fromNatural r) qs)++propPolynomialQuotientRemainderMonic ::+    Natural -> [Natural] -> [Natural] -> Bool+propPolynomialQuotientRemainderMonic np ps qs =+    Maybe.isJust (Polynomial.quotientRemainder p q)+  where+    n = np + 2+    r = Modular.ring n+    p = Polynomial.fromCoefficients r (map (Ring.fromNatural r) ps)+    q = Polynomial.fromCoefficients r (map (Ring.fromNatural r) (qs ++ [1]))++{-+np = (0 :: Natural)+ps = ([] :: [Natural])+qs = ([] :: [Natural])+n = np + 2+r = Modular.ring n+p = Polynomial.fromCoefficients r (map (Ring.fromNatural r) ps)+q = Polynomial.fromCoefficients r (map (Ring.fromNatural r) (qs ++ [1]))+-}+ check :: QuickCheck.Testable prop => String -> prop -> IO () check desc prop =     do putStr (desc ++ "\n  ")@@ -369,6 +538,10 @@        check "Generating random RSA moduli" propRandomRSA        check "Trial division" propTrialDivision        check "Modular negate" propModularNegate+       check "Modular subtract" propModularSubtract+       check "Modular exp2" propModularExp2+       check "Modular divide" propModularDivide+       check "Modular divides" propModularDivides        check "Modular invert" propModularInvert        check "Fermat's little theorem" propFermat        check "Montgomery invariant" propMontgomeryInvariant@@ -399,4 +572,17 @@        check "Williams sequence exponential" propWilliamsNthExp        check "Williams sequence equals two" propWilliamsNthEqTwo        check "Williams factorization works" propWilliamsFactor+       check "Polynomial constant degree" propPolynomialConstantDegree+       check "Polynomial fromNatural degree" propPolynomialFromNaturalDegree+       check "Polynomial add degree" propPolynomialAddDegree+       check "Polynomial negate degree" propPolynomialNegateDegree+       check "Polynomial multiply degree" propPolynomialMultiplyDegree+       check "Polynomial constant evaluate" propPolynomialConstantEvaluate+       check "Polynomial fromNatural evaluate" propPolynomialFromNaturalEvaluate+       check "Polynomial add evaluate" propPolynomialAddEvaluate+       check "Polynomial negate evaluate" propPolynomialNegateEvaluate+       check "Polynomial multiply evaluate" propPolynomialMultiplyEvaluate+       check "Polynomial quotient remainder" propPolynomialQuotientRemainder+       check "Polynomial quotient remainder monic"+         propPolynomialQuotientRemainderMonic        return ()