arithmetic 1.0 → 1.1

raw patch · 8 files changed

+375/−142 lines, 8 filesdep +opentheory-primePVP ok

version bump matches the API change (PVP)

Dependencies added: opentheory-prime

API changes (from Hackage documentation)

Files

• arithmetic.cabal +12/−9
• src/Arithmetic/ContinuedFraction.hs +102/−0
• src/Arithmetic/Modular.hs +13/−0
• src/Arithmetic/Smooth.hs +86/−0
• src/Arithmetic/SquareRoot.hs +72/−0
• src/IntegerDivides.hs +0/−31
• src/NaturalDivides.hs +0/−37
• src/Test.hs +90/−65

@@ -1,5 +1,5 @@ name: arithmetic-version: 1.0+version: 1.1 category: Number Theory synopsis: Natural number arithmetic license: MIT@@ -10,7 +10,7 @@ maintainer: Joe Leslie-Hurd <joe@gilith.com> description:   This package implements a library of natural number arithmetic functions,-  including Montgomery multiplication.+  including Montgomery multiplication and continued fractions.  Library   build-depends:@@ -20,14 +20,18 @@     opentheory-primitive >= 1.0 && < 2.0,     opentheory >= 1.0 && < 2.0,     opentheory-bits >= 1.0 && < 2.0,-    opentheory-divides >= 1.0 && < 2.0+    opentheory-divides >= 1.0 && < 2.0,+    opentheory-prime >= 1.0 && < 2.0   hs-source-dirs: src   ghc-options: -Wall   exposed-modules:+    Arithmetic.ContinuedFraction,     Arithmetic.Modular,     Arithmetic.Montgomery,     Arithmetic.Prime,-    Arithmetic.Random+    Arithmetic.Random,+    Arithmetic.Smooth,+    Arithmetic.SquareRoot  executable arithmetic   build-depends:@@ -37,7 +41,8 @@     opentheory-primitive >= 1.0 && < 2.0,     opentheory >= 1.0 && < 2.0,     opentheory-bits >= 1.0 && < 2.0,-    opentheory-divides >= 1.0 && < 2.0+    opentheory-divides >= 1.0 && < 2.0,+    opentheory-prime >= 1.0 && < 2.0   hs-source-dirs: src   ghc-options: -Wall   main-is: Main.hs@@ -51,10 +56,8 @@     opentheory-primitive >= 1.0 && < 2.0,     opentheory >= 1.0 && < 2.0,     opentheory-bits >= 1.0 && < 2.0,-    opentheory-divides >= 1.0 && < 2.0+    opentheory-divides >= 1.0 && < 2.0,+    opentheory-prime >= 1.0 && < 2.0   hs-source-dirs: src   ghc-options: -Wall   main-is: Test.hs-  other-modules:-    IntegerDivides,-    NaturalDivides

@@ -0,0 +1,102 @@+{- |+module: Arithmetic.ContinuedFraction+description: Continued fractions+license: MIT++maintainer: Joe Leslie-Hurd <joe@gilith.com>+stability: provisional+portability: portable+-}+module Arithmetic.ContinuedFraction+where++import OpenTheory.Primitive.Natural++newtype ContinuedFraction =+    ContinuedFraction {unContinuedFraction :: (Natural,[Natural])}+  deriving Eq++fromNatural :: Natural -> ContinuedFraction+fromNatural n = ContinuedFraction (n,[])++goldenRatio :: ContinuedFraction+goldenRatio = ContinuedFraction (1, repeat 1)++naturalLogarithmBase :: ContinuedFraction+naturalLogarithmBase =+    ContinuedFraction (2, go 2)+  where+    go n = 1 : n : 1 : go (n + 2)++convergentsFn :: (Natural -> a) -> (a -> a -> a) -> (a -> a -> a) ->+                 [Natural] -> a -> a -> [a]+convergentsFn lift add mult =+    go+  where+    go [] _ _ = []+    go (q : qs) y x =+        z : go qs z y+      where+        z = add (mult (lift q) y) x++numerators :: (Natural -> a) -> (a -> a -> a) -> (a -> a -> a) ->+              ContinuedFraction -> [a]+numerators lift add mult (ContinuedFraction (q0,qs)) =+    x : convergentsFn lift add mult qs x one+  where+    x = lift q0+    one = lift 1++denominators :: (Natural -> a) -> (a -> a -> a) -> (a -> a -> a) ->+                ContinuedFraction -> [a]+denominators lift add mult (ContinuedFraction (_,qs)) =+    one : convergentsFn lift add mult qs one zero+  where+    one = lift 1+    zero = lift 0++convergents :: (Natural -> a) -> (a -> a -> a) -> (a -> a -> a) ->+               (a -> a -> a) -> ContinuedFraction -> [a]+convergents lift add mult divf cf =+    zipWith divf nums dens+  where+    nums = numerators lift add mult cf+    dens = denominators lift add mult cf++unstableConvergents :: Eq a => [a] -> [a]+unstableConvergents [] = error "empty convergents"+unstableConvergents (q0 : qs) =+    q0 : go q0 qs+  where+    go _ [] = []+    go x (h : t) = if x == h then [] else h : go h t++fractionalConvergents :: Fractional a => ContinuedFraction -> [a]+fractionalConvergents = convergents fromIntegral (+) (*) (/)++rationalConvergents :: ContinuedFraction -> [Rational]+rationalConvergents = convergents fromIntegral (+) (*) (/)++toDouble :: ContinuedFraction -> Double+toDouble = last . unstableConvergents . fractionalConvergents++instance Show ContinuedFraction where+  show = show . toDouble++fromRealFrac :: RealFrac a => a -> ContinuedFraction+fromRealFrac x =+    ContinuedFraction (q0, go y)+  where+    go s =+      if s == 0.0 then []+      else let (q,t) = properFraction (1.0 / s) in q : go t++    (q0,y) = properFraction x++invert :: ContinuedFraction -> Maybe ContinuedFraction+invert (ContinuedFraction (q0,qs)) =+    if q0 /= 0 then Just (ContinuedFraction (0, q0 : qs))+    else+      case qs of+        [] -> Nothing+        h : t -> Just (ContinuedFraction (h,t))

@@ -11,6 +11,7 @@ where  import OpenTheory.Primitive.Natural+import OpenTheory.Natural.Divides import qualified OpenTheory.Natural.Bits as Bits  multiplyExponential :: (a -> a -> a) -> a -> a -> Natural -> a@@ -60,3 +61,15 @@  exp2 :: Natural -> Natural -> Natural -> Natural exp2 n x k = functionPower (square n) k x++invert :: Natural -> Natural -> Maybe Natural+invert n x =+    if g == 1 then Just s else Nothing+  where+    (g,(s,_)) = 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)

@@ -0,0 +1,86 @@+{- |+module: Arithmetic.Smooth+description: Smooth numbers+license: MIT++maintainer: Joe Leslie-Hurd <joe@gilith.com>+stability: provisional+portability: portable+-}+module Arithmetic.Smooth+where++import qualified Data.List as List+import OpenTheory.Primitive.Natural+import qualified OpenTheory.Natural.Bits as Bits+import OpenTheory.Natural.Divides+import qualified OpenTheory.Natural.Prime as Prime++factorOut :: Natural -> Natural -> Maybe (Natural,Natural)+factorOut p =+    go 0+  where+    go k n =+      if divides p n then go (k + 1) (n `div` p)+      else if k == 0 then Nothing+      else Just (k,n)++factorList :: [Natural] -> Natural -> ([(Natural,Natural)],Natural)+factorList ps n =+    case ps of+      [] -> ([],n)+      p : pt ->+        case factorOut p n of+	  Nothing -> factorList pt n+	  Just (k,m) ->+            let (pks,q) = factorList pt m in+            ((p,k) : pks, q)++factorBase :: Natural -> Natural -> ([(Natural,Natural)],Natural)+factorBase k = factorList (take (fromIntegral k) Prime.primes)++multiplyBase :: ([(Natural,Natural)],Natural) -> Natural+multiplyBase =+    \(pks,m) -> foldr mult m pks+  where+    mult (p,k) m = (p ^ k) * m++newtype Smooth =+    Smooth {unSmooth :: ([(Natural,Natural)],Natural)}+  deriving (Eq,Ord)++instance Show Smooth where+  show s =+      if null factors then "1" else List.intercalate "*" factors+    where+      factors = map showPk pks ++ showRest+      showRest = if n == 1 then [] else [showWidth]+      showWidth = if w < 20 then show n+                  else "[" ++ show w ++ "]"+      showPk (p,k) = show p ++ showExp k+      showExp k = if k == 1 then "" else "^" ++ show k+      (pks,n) = unSmooth s+      w = Bits.width n++fromNatural :: Natural -> Natural -> Smooth+fromNatural k = Smooth . factorBase k++toNatural :: Smooth -> Natural+toNatural = multiplyBase . unSmooth++factoring :: Smooth -> Maybe [(Natural,Natural)]+factoring s =+    if n == 1 then Just pks else Nothing+  where+    (pks,n) = unSmooth s++next :: Natural -> Natural -> Smooth+next k =+    go+  where+    go n =+        case factoring s of+          Nothing -> go (n + 1)+          Just _ -> s+      where+        s = fromNatural k n

@@ -0,0 +1,72 @@+{- |+module: Arithmetic.SquareRoot+description: Natural number square root+license: MIT++maintainer: Joe Leslie-Hurd <joe@gilith.com>+stability: provisional+portability: portable+-}+module Arithmetic.SquareRoot+where++import OpenTheory.Primitive.Natural+import qualified Data.List as List++import qualified Arithmetic.ContinuedFraction as ContinuedFraction++floor :: Natural -> Natural+floor n =+    if n < 2 then n else bisect 0 n+  where+    bisect l u =+        if m == l then l+	else if m * m <= n then bisect m u+	else bisect l m+      where+        m = (l + u) `div` 2++ceiling :: Natural -> Natural+ceiling n =+    if sqrtn * sqrtn == n then sqrtn else sqrtn + 1+  where+    sqrtn = Arithmetic.SquareRoot.floor n++continuedFraction :: Natural -> ContinuedFraction.ContinuedFraction+continuedFraction n =+    ContinuedFraction.ContinuedFraction (sqrtn,qs)+  where+    sqrtn = Arithmetic.SquareRoot.floor n++    ps = continuedFractionPeriodicTail n sqrtn++    qs = if null ps then [] else cycle ps++continuedFractionPeriodic :: Natural -> [Natural]+continuedFractionPeriodic n =+    continuedFractionPeriodicTail n sqrtn+  where+    sqrtn = Arithmetic.SquareRoot.floor n++continuedFractionPeriodicTail :: Natural -> Natural -> [Natural]+continuedFractionPeriodicTail n sqrtn =+    List.unfoldr go (sqrtn,sqrtd)+  where+    sqrtd = n - sqrtn * sqrtn++-- (sqrt(n) + a) / b = c + 1 / x ==>+-- x = b / (sqrt(n) + a - c * b)+--   = b / (sqrt(n) - (c * b - a))+--   = (b * (sqrt(n) + (c * b - a))) / (n - (c * b - a)^2)+    advance (a,b) =+        (c,(d,e))+      where+        c = (sqrtn + a) `div` b+        d = c * b - a+        e = (n - d * d) `div` b++    go (a,b) =+        case b of+          0 -> Nothing+          1 -> Just (2 * a, (0,0))+          _ -> Just (advance (a,b))

@@ -1,31 +0,0 @@-{- |-module: IntegerDivides-description: Integer division algorithms-license: MIT--maintainer: Joe Leslie-Hurd <joe@gilith.com>-stability: provisional-portability: portable--}-module IntegerDivides-where--divides :: Integer -> Integer -> Bool-divides 0 b = b == 0-divides a b = abs b `mod` abs a == 0--egcd :: Integer -> Integer -> (Integer,(Integer,Integer))-egcd a 0 = (a,(1,0))-egcd a b =-    (g, (t, s - (a `div` b) * t))-  where-    (g,(s,t)) = egcd b (a `mod` b)--chineseRemainder :: Integer -> Integer -> Integer -> Integer -> Integer-chineseRemainder a b =-    \x y -> (x * tb + y * sa) `mod` ab-  where-    (_,(s,t)) = egcd a b-    ab = a * b-    sa = s * a-    tb = t * b

@@ -1,37 +0,0 @@-{- |-module: NaturalDivides-description: Natural number division algorithms-license: MIT--maintainer: Joe Leslie-Hurd <joe@gilith.com>-stability: provisional-portability: portable--}-module NaturalDivides-where--import OpenTheory.Primitive.Natural--divides :: Natural -> Natural -> Bool-divides 0 b = b == 0-divides a b = b `mod` a == 0--egcd :: Natural -> Natural -> (Natural,(Natural,Natural))-egcd a 0 = (a,(1,0))-egcd a b =-    if c == 0-    then (b, (1, a `div` b - 1))-    else (g, (u, t + (a `div` b) * u))-  where-    c = a `mod` b-    (g,(s,t)) = egcd c (b `mod` c)-    u = s + (b `div` c) * t--chineseRemainder :: Natural -> Natural -> Natural -> Natural -> Natural-chineseRemainder a b =-    \x y -> (x * tb + y * sa) `mod` ab-  where-    (_,(s,t)) = egcd a b-    ab = a * b-    sa = s * a-    tb = (a - t) * b

@@ -14,72 +14,97 @@ import qualified Test.QuickCheck as QuickCheck import OpenTheory.Primitive.Natural import OpenTheory.Natural+import OpenTheory.Natural.Divides import qualified OpenTheory.Primitive.Random as Random import qualified OpenTheory.Natural.Uniform as Uniform import OpenTheory.Primitive.Test -import qualified IntegerDivides-import qualified NaturalDivides import Arithmetic.Random import Arithmetic.Prime+import qualified Arithmetic.ContinuedFraction as ContinuedFraction import qualified Arithmetic.Modular as Modular import qualified Arithmetic.Montgomery as Montgomery+import qualified Arithmetic.Smooth as Smooth+import qualified Arithmetic.SquareRoot as SquareRoot -propIntegerEgcdDivides :: Integer -> Integer -> Bool-propIntegerEgcdDivides a b =-    let (g,_) = IntegerDivides.egcd a b in-    IntegerDivides.divides g a && IntegerDivides.divides g b+propEgcdDivides :: Natural -> Natural -> Bool+propEgcdDivides a b =+    divides g a && divides g b+  where+    (g,_) = egcd a b -propIntegerEgcdEquation :: Integer -> Integer -> Bool-propIntegerEgcdEquation a b =-    let (g,(s,t)) = IntegerDivides.egcd a b in-    s * a + t * b == g+propEgcdEquation :: Natural -> Natural -> Bool+propEgcdEquation ap b =+    s * a == t * b + g+  where+    a = ap + 1+    (g,(s,t)) = egcd a b -propIntegerEgcdBound :: Integer -> Integer -> Bool-propIntegerEgcdBound a b =-    let (_,(s,t)) = IntegerDivides.egcd a b in-    abs s <= max ((abs b + 1) `div` 2) 1 &&-    abs t <= max ((abs a + 1) `div` 2) 1+propEgcdBound :: Natural -> Natural -> Bool+propEgcdBound ap b =+    s < max b 2 && t < a+  where+    a = ap + 1+    (_,(s,t)) = egcd a b -propNaturalEgcdDivides :: Natural -> Natural -> Bool-propNaturalEgcdDivides a b =-    let (g,_) = NaturalDivides.egcd a b in-    NaturalDivides.divides g a && NaturalDivides.divides g b+propSmoothInjective :: Natural -> Natural -> Bool+propSmoothInjective k np =+    Smooth.toNatural (Smooth.fromNatural k n) == n+  where+    n = np + 1 -propNaturalEgcdEquation :: Natural -> Natural -> Bool-propNaturalEgcdEquation ap b =-    let a = ap + 1 in-    let (g,(s,t)) = NaturalDivides.egcd a b in-    s * a == t * b + g+propFloorSqrt :: Natural -> Bool+propFloorSqrt n =+    sq s <= n && n < sq (s + 1)+  where+    s = SquareRoot.floor n+    sq i = i * i -propNaturalEgcdBound :: Natural -> Natural -> Bool-propNaturalEgcdBound ap b =-    let a = ap + 1 in-    let (_,(s,t)) = NaturalDivides.egcd a b in-    s < max b 2 && t < a+propCeilingSqrt :: Natural -> Bool+propCeilingSqrt n =+    (s == 0 || sq (s - 1) < n) && n <= sq s+  where+    s = SquareRoot.ceiling n+    sq i = i * i -propIntegerChineseRemainder :: Int -> Random.Random -> Bool-propIntegerChineseRemainder w r =-    n `mod` a == x && n `mod` b == y && n < a * b+propContinuedFractionSqrt :: Natural -> Bool+propContinuedFractionSqrt n =+    cf == spec   where-    (a,b) = randomCoprimeInteger w r1-    x = uniformInteger a r2-    y = uniformInteger b r3-    n = IntegerDivides.chineseRemainder a b x y-    (r1,r23) = Random.split r-    (r2,r3) = Random.split r23+    cf = ContinuedFraction.toDouble (SquareRoot.continuedFraction n)+    spec = sqrt (fromIntegral n) -propNaturalChineseRemainder :: Int -> Random.Random -> Bool-propNaturalChineseRemainder w r =+propChineseRemainder :: Int -> Random.Random -> Bool+propChineseRemainder w r =     n `mod` a == x && n `mod` b == y && n < a * b   where     (a,b) = randomCoprime w r1     x = Uniform.random a r2     y = Uniform.random b r3-    n = NaturalDivides.chineseRemainder a b x y+    n = chineseRemainder a b x y     (r1,r23) = Random.split r     (r2,r3) = Random.split r23 +propModularNegate :: Int -> Random.Random -> Bool+propModularNegate nw rnd =+    Modular.add n a b == 0 &&+    b < n+  where+    n = randomWidth nw r1+    a = Uniform.random n r2+    b = Modular.negate n a+    (r1,r2) = Random.split rnd++propModularInvert :: Int -> Random.Random -> Bool+propModularInvert nw rnd =+    case Modular.invert n a of+      Nothing -> True+      Just b -> Modular.multiply n a b == 1 && b < n+  where+    n = randomWidth nw r1+    a = Uniform.random n r2+    (r1,r2) = Random.split rnd+ randomMontgomeryParameters :: Int -> Random.Random -> Montgomery.Parameters randomMontgomeryParameters w r = Montgomery.standardParameters (randomOdd w r) @@ -281,35 +306,35 @@  checkWidthProps :: Int -> IO () checkWidthProps w =-   do checkWidthProp w "Check integer Chinese remainder properties"-        propIntegerChineseRemainder-      checkWidthProp w "Check natural Chinese remainder properties"-        propNaturalChineseRemainder-      checkWidthProp w "Check Montgomery invariant" propMontgomeryInvariant-      checkWidthProp w "Check Montgomery normalize" propMontgomeryNormalize-      checkWidthProp w "Check Montgomery reduce" propMontgomeryReduce-      checkWidthProp w "Check Montgomery reduce small" propMontgomeryReduceSmall-      checkWidthProp w "Check Montgomery toNatural" propMontgomeryToNatural-      checkWidthProp w "Check Montgomery fromNatural" propMontgomeryFromNatural-      checkWidthProp w "Check Montgomery zero" propMontgomeryZero-      checkWidthProp w "Check Montgomery one" propMontgomeryOne-      checkWidthProp w "Check Montgomery two" propMontgomeryTwo-      checkWidthProp w "Check Montgomery add" propMontgomeryAdd-      checkWidthProp w "Check Montgomery negate" propMontgomeryNegate-      checkWidthProp w "Check Montgomery multiply" propMontgomeryMultiply-      checkWidthProp w "Check Montgomery modexp" propMontgomeryModexp-      checkWidthProp w "Check Montgomery modexp2" propMontgomeryModexp2+   do checkWidthProp w "Chinese remainder" propChineseRemainder+      checkWidthProp w "Modular negate" propModularNegate+      checkWidthProp w "Modular invert" propModularInvert+      checkWidthProp w "Montgomery invariant" propMontgomeryInvariant+      checkWidthProp w "Montgomery normalize" propMontgomeryNormalize+      checkWidthProp w "Montgomery reduce" propMontgomeryReduce+      checkWidthProp w "Montgomery reduce small" propMontgomeryReduceSmall+      checkWidthProp w "Montgomery toNatural" propMontgomeryToNatural+      checkWidthProp w "Montgomery fromNatural" propMontgomeryFromNatural+      checkWidthProp w "Montgomery zero" propMontgomeryZero+      checkWidthProp w "Montgomery one" propMontgomeryOne+      checkWidthProp w "Montgomery two" propMontgomeryTwo+      checkWidthProp w "Montgomery add" propMontgomeryAdd+      checkWidthProp w "Montgomery negate" propMontgomeryNegate+      checkWidthProp w "Montgomery multiply" propMontgomeryMultiply+      checkWidthProp w "Montgomery modexp" propMontgomeryModexp+      checkWidthProp w "Montgomery modexp2" propMontgomeryModexp2       checkWidthProp w "Fermat's little theorem" propFermat       return ()  main :: IO () main =-    do check "Check integer egcd divides\n  " propIntegerEgcdDivides-       check "Check integer egcd equation\n  " propIntegerEgcdEquation-       check "Check integer egcd bound\n  " propIntegerEgcdBound-       check "Check natural egcd divides\n  " propNaturalEgcdDivides-       check "Check natural egcd equation\n  " propNaturalEgcdEquation-       check "Check natural egcd bound\n  " propNaturalEgcdBound+    do check "Check egcd divides\n  " propEgcdDivides+       check "Check egcd equation\n  " propEgcdEquation+       check "Check egcd bound\n  " propEgcdBound+       check "Check smooth injective\n  " propSmoothInjective+       check "Check floor square root\n  " propFloorSqrt+       check "Check ceiling square root\n  " propCeilingSqrt+       check "Check continued fraction square root\n  " propContinuedFractionSqrt        mapM_ checkWidthProps ws        return ()   where