packages feed

crypto-numbers (empty) → 0.1.0

raw patch · 12 files changed

+918/−0 lines, 12 filesdep +HUnitdep +QuickCheckdep +basesetup-changed

Dependencies added: HUnit, QuickCheck, base, bytestring, criterion, crypto-numbers, crypto-random-api, mtl, test-framework, test-framework-hunit, test-framework-quickcheck2, vector

Files

+ Benchmarks/Benchmarks.hs view
@@ -0,0 +1,68 @@+module Main where++import Criterion.Main++import Crypto.Number.Serialize+import Crypto.Number.Generate+import qualified Data.ByteString as B+import Crypto.Number.ModArithmetic+import Data.Bits++primes = [3, 5, 7, 29, 31, 211, 2309, 2311, 30029, 200560490131, 304250263527209]+carmichaelNumbers = [41041, 62745, 63973, 75361, 101101, 126217, 172081, 188461, 278545, 340561]++lg1, lg2 :: Integer+lg1 = 21389083291083903845902381390285907190274907230982112390820985903825329874812973821790321904790217490217409721904832974210974921740972109481490128430982190472109874802174907490271904124908210958093285098309582093850918902581290859012850829105809128590218590281905812905810928590128509128940821903829018390849839578967358920127598901248259797158249684571948075896458741905823982671490352896791052386357019528367902+lg2 = 21392813098390824190840192812389082390812940821904891028439028490128904829104891208940835932882910839218309812093118249089871209347472901874902407219740921840928149087284397490128903843789289014374839281492038091283923091809734832974180398210938901284839274091749021709++bitsAndShift8 n i = (n `shiftR` i, n .&. 0xff)+modAndShift8 n i = (n `shiftR` i, n `mod` 0x100)++main = defaultMain+    [ bgroup "std ops"+        [ bench "mod" $ nf (mod lg1) lg2+        , bench "rem" $ nf (rem lg1) lg2+        , bench "div" $ nf (div lg1) lg2+        , bench "quot" $ nf (quot lg1) lg2+        , bench "divmod" $ nf (divMod lg1) lg2+        , bench "quotRem" $ nf (quotRem lg1) lg2+        ]+    , bgroup "divMod by 256"+        [ bench "divmod 256" $ nf (divMod lg1) 256+        , bench "quotRem 256" $ nf (quotRem lg1) 256+        , bench "modAndShift 8" $ nf (modAndShift8 lg1) 8+        , bench "bitsAndShift 8" $ nf (bitsAndShift8 lg1) 8+        ]+    , bgroup "serialization bs->i"+        [ bench "8"    $ nf os2ip b8+        , bench "32"   $ nf os2ip b32+        , bench "64"   $ nf os2ip b64+        , bench "256"  $ nf os2ip b256+        , bench "1024" $ nf os2ip b1024+        ]+    , bgroup "serialization i->bs"+        [ bench "10"     $ nf i2osp (2^10)+        , bench "100"    $ nf i2osp (2^100)+        , bench "1000"   $ nf i2osp (2^1000)+        , bench "10000"  $ nf i2osp (2^10000)+        , bench "100000" $ nf i2osp (2^100000)+        ]+    , bgroup "serialization i->bs of size"+        [ bench "10"     $ nf (i2ospOf_ 4) (2^10)+        , bench "100"    $ nf (i2ospOf_ 16) (2^100)+        , bench "1000"   $ nf (i2ospOf_ 128) (2^1000)+        , bench "10000"  $ nf (i2ospOf_ 1560) (2^10000)+        , bench "100000" $ nf (i2ospOf_ 12502) (2^100000)+        ]+    , bgroup "exponentiation"+        [ bench "2^1234 mod 2^999" $ nf (exponantiation 2 1234) (2^999)+        , bench "130^5432 mod 100^9990" $ nf (exponantiation 130 5432) (100^9999)+        , bench "2^1234 mod 2^999" $ nf (exponantiation_rtl_binary 2 1234) (2^999)+        , bench "130^5432 mod 100^9990" $ nf (exponantiation_rtl_binary 130 5432) (100^9999)+        ]+    ]+    where b8    = B.replicate 8 0xf7+          b32   = B.replicate 32 0xf7+          b64   = B.replicate 64 0x7f+          b256  = B.replicate 256 0x7f+          b1024 = B.replicate 1024 0x7f
+ Crypto/Number/Basic.hs view
@@ -0,0 +1,88 @@+{-# LANGUAGE BangPatterns #-}+-- |+-- Module      : Crypto.Number.Basic+-- License     : BSD-style+-- Maintainer  : Vincent Hanquez <vincent@snarc.org>+-- Stability   : experimental+-- Portability : Good++module Crypto.Number.Basic+    ( sqrti+    , gcde+    , gcde_binary+    , areEven+    ) where++import Data.Bits++-- | sqrti returns two integer (l,b) so that l <= sqrt i <= b+-- the implementation is quite naive, use an approximation for the first number+-- and use a dichotomy algorithm to compute the bound relatively efficiently.+sqrti :: Integer -> (Integer, Integer)+sqrti i+    | i < 0     = error "cannot compute negative square root"+    | i == 0    = (0,0)+    | i == 1    = (1,1)+    | i == 2    = (1,2)+    | otherwise = loop x0+        where+            nbdigits = length $ show i+            x0n = (if even nbdigits then nbdigits - 2 else nbdigits - 1) `div` 2+            x0  = if even nbdigits then 2 * 10 ^ x0n else 6 * 10 ^ x0n+            loop x = case compare (sq x) i of+                LT -> iterUp x+                EQ -> (x, x)+                GT -> iterDown x+            iterUp lb = if sq ub >= i then iter lb ub else iterUp ub+                where ub = lb * 2+            iterDown ub = if sq lb >= i then iterDown lb else iter lb ub+                where lb = ub `div` 2+            iter lb ub+                | lb == ub   = (lb, ub)+                | lb+1 == ub = (lb, ub)+                | otherwise  =+                    let d = (ub - lb) `div` 2 in+                    if sq (lb + d) >= i+                        then iter lb (ub-d)+                        else iter (lb+d) ub+            sq a = a * a++-- | get the extended GCD of two integer using integer divMod+gcde :: Integer -> Integer -> (Integer, Integer, Integer)+gcde a b = if d < 0 then (-x,-y,-d) else (x,y,d) where+    (d, x, y)                     = f (a,1,0) (b,0,1)+    f t              (0, _, _)    = t+    f (a', sa, ta) t@(b', sb, tb) =+        let (q, r) = a' `divMod` b' in+        f t (r, sa - (q * sb), ta - (q * tb))++-- | get the extended GCD of two integer using the extended binary algorithm (HAC 14.61)+-- get (x,y,d) where d = gcd(a,b) and x,y satisfying ax + by = d+gcde_binary :: Integer -> Integer -> (Integer, Integer, Integer)+gcde_binary a' b'+    | b' == 0   = (1,0,a')+    | a' >= b'  = compute a' b'+    | otherwise = (\(x,y,d) -> (y,x,d)) $ compute b' a'+    where+        getEvenMultiplier !g !x !y+            | areEven [x,y] = getEvenMultiplier (g `shiftL` 1) (x `shiftR` 1) (y `shiftR` 1)+            | otherwise     = (x,y,g)+        halfLoop !x !y !u !i !j+            | areEven [u,i,j] = halfLoop x y (u `shiftR` 1) (i `shiftR` 1) (j `shiftR` 1)+            | even u          = halfLoop x y (u `shiftR` 1) ((i + y) `shiftR` 1) ((j - x) `shiftR` 1)+            | otherwise       = (u, i, j)+        compute a b =+            let (x,y,g) = getEvenMultiplier 1 a b in+            loop g x y x y 1 0 0 1++        loop g _ _ 0  !v _  _  !c !d = (c, d, g * v)+        loop g x y !u !v !a !b !c !d =+            let (u2,a2,b2) = halfLoop x y u a b+                (v2,c2,d2) = halfLoop x y v c d+             in if u2 >= v2+                then loop g x y (u2 - v2) v2 (a2 - c2) (b2 - d2) c2 d2+                else loop g x y u2 (v2 - u2) a2 b2 (c2 - a2) (d2 - b2)++-- | check if a list of integer are all even+areEven :: [Integer] -> Bool+areEven = and . map even
+ Crypto/Number/Generate.hs view
@@ -0,0 +1,35 @@+-- |+-- Module      : Crypto.Number.Generate+-- License     : BSD-style+-- Maintainer  : Vincent Hanquez <vincent@snarc.org>+-- Stability   : experimental+-- Portability : Good++module Crypto.Number.Generate+    ( generateMax+    , generateBetween+    , generateOfSize+    ) where++import Crypto.Number.Serialize+import Crypto.Random.API+import qualified Data.ByteString as B+import Data.Bits ((.|.))++-- | generate a positive integer between 0 and m.+-- using as many bytes as necessary to the same size as m, that are converted to integer.+generateMax :: CPRG g => g -> Integer -> (Integer, g)+generateMax rng m = withRandomBytes rng (lengthBytes m) $ \bs ->+    os2ip bs `mod` m++-- | generate a number between the inclusive bound [low,high].+generateBetween :: CPRG g => g -> Integer -> Integer -> (Integer, g)+generateBetween rng low high = (low + v, rng')+    where (v, rng') = generateMax rng (high - low + 1)++-- | generate a positive integer of a specific size in bits.+-- the number of bits need to be multiple of 8. It will always returns+-- an integer that is close to 2^(1+bits/8) by setting the 2 highest bits to 1.+generateOfSize :: CPRG g => g -> Int -> (Integer, g)+generateOfSize rng bits = withRandomBytes rng (bits `div` 8) $ \bs ->+    os2ip $ snd $ B.mapAccumL (\acc w -> (0, w .|. acc)) 0xc0 bs
+ Crypto/Number/ModArithmetic.hs view
@@ -0,0 +1,67 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DeriveDataTypeable #-}+-- |+-- Module      : Crypto.Number.ModArithmetic+-- License     : BSD-style+-- Maintainer  : Vincent Hanquez <vincent@snarc.org>+-- Stability   : experimental+-- Portability : Good++module Crypto.Number.ModArithmetic+    ( exponantiation_rtl_binary+    , exponantiation+    , inverse+    , inverseCoprimes+    ) where++import Control.Exception (throw, Exception)+import Crypto.Number.Basic (gcde_binary)+import Data.Bits+import Data.Typeable++-- | Raised when two numbers are supposed to be coprimes but are not.+data CoprimesAssertionError = CoprimesAssertionError+    deriving (Show,Typeable)++instance Exception CoprimesAssertionError++-- note on exponantiation: 0^0 is treated as 1 for mimicking the standard library;+-- the mathematic debate is still open on whether or not this is true, but pratically+-- in computer science it shouldn't be useful for anything anyway.++-- | exponantiation_rtl_binary computes modular exponantiation as b^e mod m+-- using the right-to-left binary exponentiation algorithm (HAC 14.79)+exponantiation_rtl_binary :: Integer -> Integer -> Integer -> Integer+exponantiation_rtl_binary 0 0 m = 1 `mod` m+exponantiation_rtl_binary b e m = loop e b 1+    where sq x          = (x * x) `mod` m+          loop !0 _  !a = a `mod` m+          loop !i !s !a = loop (i `shiftR` 1) (sq s) (if odd i then a * s else a)++-- | exponantiation computes modular exponantiation as b^e mod m+-- using repetitive squaring.+exponantiation :: Integer -> Integer -> Integer -> Integer+exponantiation b e m+    | b == 1    = b+    | e == 0    = 1+    | e == 1    = b `mod` m+    | even e    = let p = (exponantiation b (e `div` 2) m) `mod` m+                   in (p^(2::Integer)) `mod` m+    | otherwise = (b * exponantiation b (e-1) m) `mod` m++-- | inverse computes the modular inverse as in g^(-1) mod m+inverse :: Integer -> Integer -> Maybe Integer+inverse g m = if d > 1 then Nothing else Just (x `mod` m)+    where (x,_,d) = gcde_binary g m++-- | Compute the modular inverse of 2 coprime numbers.+-- This is equivalent to inverse except that the result+-- is known to exists.+--+-- if the numbers are not defined as coprime, this function+-- will raise a CoprimesAssertionError.+inverseCoprimes :: Integer -> Integer -> Integer+inverseCoprimes g m =+    case inverse g m of+        Nothing -> throw CoprimesAssertionError+        Just i  -> i
+ Crypto/Number/Polynomial.hs view
@@ -0,0 +1,133 @@+{-# LANGUAGE BangPatterns #-}+-- |+-- Module      : Crypto.Number.Polynomial+-- License     : BSD-style+-- Maintainer  : Vincent Hanquez <vincent@snarc.org>+-- Stability   : experimental+-- Portability : Good++module Crypto.Number.Polynomial+    ( Monomial(..)+    -- * polynomial operations+    , Polynomial+    , toList+    , fromList+    , addPoly+    , subPoly+    , mulPoly+    , squarePoly+    , expPoly+    , divPoly+    , negPoly+    ) where++import Data.List (intercalate, sort)+import Data.Vector ((!), Vector)+import qualified Data.Vector as V+import Control.Arrow (first)++data Monomial = Monomial {-# UNPACK #-} !Int !Integer+    deriving (Eq)++data Polynomial = Polynomial (Vector Monomial)+    deriving (Eq)++instance Ord Monomial where+    compare (Monomial w1 v1) (Monomial w2 v2) =+        case compare w1 w2 of+            EQ -> compare v1 v2+            r  -> r++instance Show Monomial where+    show (Monomial w v) = show v ++ "x^" ++ show w++instance Show Polynomial where+    show (Polynomial p) = intercalate "+" $ map show $ V.toList p++toList :: Polynomial -> [Monomial]+toList (Polynomial p) = V.toList p++fromList :: [Monomial] -> Polynomial+fromList = Polynomial . V.fromList . reverse . sort . filterZero+    where+        filterZero = filter (\(Monomial _ v) -> v /= 0)++getWeight :: Polynomial -> Int -> Maybe Integer+getWeight (Polynomial p) n = look 0+    where+        plen = V.length p+        look !i+            | i >= plen = Nothing+            | otherwise =+                let (Monomial w v) = p ! i in+                case compare w n of+                    LT -> Nothing+                    EQ -> Just v+                    GT -> look (i+1)+        ++mergePoly :: (Integer -> Integer -> Integer) -> Polynomial -> Polynomial -> Polynomial+mergePoly f (Polynomial p1) (Polynomial p2) = fromList $ loop 0 0+    where+        l1 = V.length p1+        l2 = V.length p2+        loop !i1 !i2+            | i1 == l1 && i2 == l2 = []+            | i1 == l1             = (p2 ! i2) : loop i1 (i2+1)+            | i2 == l2             = (p1 ! i1) : loop (i1+1) i2+            | otherwise            =+                let (coef, i1inc, i2inc) = addCoef (p1 ! i1) (p2 ! i2) in+                coef : loop (i1+i1inc) (i2+i2inc)+        addCoef m1@(Monomial w1 v1) (Monomial w2 v2) =+            case compare w1 w2 of+                LT -> (Monomial w2 (f 0 v2), 0, 1)+                EQ -> (Monomial w1 (f v1 v2), 1, 1)+                GT -> (m1, 1, 0)++addPoly :: Polynomial -> Polynomial -> Polynomial+addPoly = mergePoly (+)++subPoly :: Polynomial -> Polynomial -> Polynomial+subPoly = mergePoly (-)++negPoly :: Polynomial -> Polynomial+negPoly (Polynomial p) = Polynomial $ V.map negateMonomial p+    where negateMonomial (Monomial w v) = Monomial w (-v)++mulPoly :: Polynomial -> Polynomial -> Polynomial+mulPoly p1@(Polynomial v1) p2@(Polynomial v2) =+    fromList $ filter (\(Monomial _ v) -> v /= 0) $ map (\i -> Monomial i (c i)) $ reverse [0..(m+n)]+    where+        (Monomial m _) = v1 ! 0+        (Monomial n _) = v2 ! 0+        c r = foldl (\acc i -> (b $ r-i) * (a $ i) + acc) 0 [0..r]+            where+                a = maybe 0 id . getWeight p1+                b = maybe 0 id . getWeight p2++squarePoly :: Polynomial -> Polynomial+squarePoly p = p `mulPoly` p++expPoly :: Polynomial -> Integer -> Polynomial+expPoly p e = loop p e+    where+        loop t 0 = t+        loop t n = loop (squarePoly t) (n-1)++divPoly :: Polynomial -> Polynomial -> (Polynomial, Polynomial)+divPoly p1 p2@(Polynomial pp2) = first fromList $ divLoop p1+    where divLoop d1@(Polynomial pp1)+            | V.null pp1 = ([], d1)+            | otherwise  =+                let (Monomial w1 v1) = pp1 ! 0 in+                let (Monomial w2 v2) = pp2 ! 0 in+                let w = w1 - w2 in+                let (v,r) = v1 `divMod` v2 in+                if w >= 0 && r == 0+                    then+                        let mono = (Monomial w v) in+                        let remain = d1 `subPoly` (p2 `mulPoly` (fromList [mono])) in+                        let (l, finalRem) = divLoop remain in+                        (mono : l, finalRem)+                    else+                        ([], d1)
+ Crypto/Number/Prime.hs view
@@ -0,0 +1,197 @@+{-# LANGUAGE BangPatterns #-}+-- |+-- Module      : Crypto.Number.Prime+-- License     : BSD-style+-- Maintainer  : Vincent Hanquez <vincent@snarc.org>+-- Stability   : experimental+-- Portability : Good++module Crypto.Number.Prime+    ( generatePrime+    , generateSafePrime+    , isProbablyPrime+    , findPrimeFrom+    , findPrimeFromWith+    , primalityTestNaive+    , primalityTestMillerRabin+    , primalityTestFermat+    , isCoprime+    ) where++import Crypto.Random.API+import Data.Bits+import Crypto.Number.Generate+import Crypto.Number.Basic (sqrti, gcde_binary)+import Crypto.Number.ModArithmetic (exponantiation)++-- | returns if the number is probably prime.+-- first a list of small primes are implicitely tested for divisibility,+-- then a fermat primality test is used with arbitrary numbers and+-- then the Miller Rabin algorithm is used with an accuracy of 30 recursions+isProbablyPrime :: CPRG g => g -> Integer -> (Bool, g)+isProbablyPrime rng !n+    | any (\p -> p `divides` n) (filter (< n) firstPrimes) = (False, rng)+    | primalityTestFermat 50 (n`div`2) n                   = primalityTestMillerRabin rng 30 n+    | otherwise                                            = (False, rng)++-- | generate a prime number of the required bitsize+generatePrime :: CPRG g => g -> Int -> (Integer, g)+generatePrime rng bits =+    let (sp, rng') = generateOfSize rng bits+     in findPrimeFrom rng' sp++-- | generate a prime number of the form 2p+1 where p is also prime.+-- it is also knowed as a Sophie Germaine prime or safe prime.+--+-- The number of safe prime is significantly smaller to the number of prime,+-- as such it shouldn't be used if this number is supposed to be kept safe.+generateSafePrime :: CPRG g => g -> Int -> (Integer, g)+generateSafePrime rng bits =+    let (sp, rng') = generateOfSize rng bits+        (p, rng'') = findPrimeFromWith rng' (\g i -> isProbablyPrime g (2*i+1)) (sp `div` 2)+     in (2*p+1, rng'')++-- | find a prime from a starting point where the property hold.+findPrimeFromWith :: CPRG g => g -> (g -> Integer -> (Bool,g)) -> Integer -> (Integer, g)+findPrimeFromWith rng prop !n+    | even n        = findPrimeFromWith rng prop (n+1)+    | otherwise     = case isProbablyPrime rng n of+        (False, rng')    -> findPrimeFromWith rng' prop (n+2)+        (True, rng')     ->+            case prop rng' n of+                (False, rng'') -> findPrimeFromWith rng'' prop (n+2)+                (True, rng'')  -> (n, rng'')++-- | find a prime from a starting point with no specific property.+findPrimeFrom :: CPRG g => g -> Integer -> (Integer, g)+findPrimeFrom rng n = findPrimeFromWith rng (\g _ -> (True, g)) n++-- | Miller Rabin algorithm return if the number is probably prime or composite.+-- the tries parameter is the number of recursion, that determines the accuracy of the test.+primalityTestMillerRabin :: CPRG g => g -> Int -> Integer -> (Bool, g)+primalityTestMillerRabin rng tries !n+    | n <= 3     = error "Miller-Rabin requires tested value to be > 3"+    | even n     = (False, rng)+    | tries <= 0 = error "Miller-Rabin tries need to be > 0"+    | otherwise  = let (witnesses, rng') = generateTries tries rng+                    in (loop witnesses, rng')+        where !nm1 = n-1+              !nm2 = n-2++              (!s,!d) = (factorise 0 nm1)++              generateTries 0 g = ([], g)+              generateTries t g = let (v,g')   = generateBetween g 2 nm2+                                      (vs,g'') = generateTries (t-1) g'+                                   in (v:vs, g'')++              -- factorise n-1 into the form 2^s*d+              factorise :: Integer -> Integer -> (Integer, Integer)+              factorise !si !vi+                  | vi `testBit` 0 = (si, vi)+                  | otherwise     = factorise (si+1) (vi `shiftR` 1) -- probably faster to not shift v continously, but just once.+              expmod = exponantiation++              -- when iteration reach zero, we have a probable prime+              loop []     = True+              loop (w:ws) = let x = expmod w d n+                             in if x == (1 :: Integer) || x == nm1+                                   then loop ws+                                   else loop' ws ((x*x) `mod` n) 1++              -- loop from 1 to s-1. if we reach the end then it's composite+              loop' ws !x2 !r+                  | r == s    = False+                  | x2 == 1   = False+                  | x2 /= nm1 = loop' ws ((x2*x2) `mod` n) (r+1)+                  | otherwise = loop ws++{-+    n < z -> witness to test+              1373653 [2,3]+              9080191 [31,73]+              4759123141 [2,7,61]+              2152302898747 [2,3,5,7,11]+              3474749660383 [2,3,5,7,11,13]+              341550071728321 [2,3,5,7,11,13,17]+-}++-- | Probabilitic Test using Fermat primility test.+-- Beware of Carmichael numbers that are Fermat liars, i.e. this test+-- is useless for them. always combines with some other test.+primalityTestFermat :: Int -- ^ number of iterations of the algorithm+                    -> Integer -- ^ starting a+                    -> Integer -- ^ number to test for primality+                    -> Bool+primalityTestFermat n a p = and $ map expTest [a..(a+fromIntegral n)]+    where !pm1 = p-1+          expTest i = exponantiation i pm1 p == 1++-- | Test naively is integer is prime.+-- while naive, we skip even number and stop iteration at i > sqrt(n)+primalityTestNaive :: Integer -> Bool+primalityTestNaive n+    | n <= 1    = False+    | n == 2    = True+    | even n    = False+    | otherwise = search 3+        where !ubound = snd $ sqrti n+              search !i+                  | i > ubound    = True+                  | i `divides` n = False+                  | otherwise     = search (i+2)++-- | Test is two integer are coprime to each other+isCoprime :: Integer -> Integer -> Bool+isCoprime m n = case gcde_binary m n of (_,_,d) -> d == 1++-- | list of the first primes till 2903..+firstPrimes :: [Integer]+firstPrimes =+    [ 2    , 3    , 5    , 7    , 11   , 13   , 17   , 19   , 23   , 29+    , 31   , 37   , 41   , 43   , 47   , 53   , 59   , 61   , 67   , 71+    , 73   , 79   , 83   , 89   , 97   , 101  , 103  , 107  , 109  , 113+    , 127  , 131  , 137  , 139  , 149  , 151  , 157  , 163  , 167  , 173+    , 179  , 181  , 191  , 193  , 197  , 199  , 211  , 223  , 227  , 229+    , 233  , 239  , 241  , 251  , 257  , 263  , 269  , 271  , 277  , 281+    , 283  , 293  , 307  , 311  , 313  , 317  , 331  , 337  , 347  , 349+    , 353  , 359  , 367  , 373  , 379  , 383  , 389  , 397  , 401  , 409+    , 419  , 421  , 431  , 433  , 439  , 443  , 449  , 457  , 461  , 463+    , 467  , 479  , 487  , 491  , 499  , 503  , 509  , 521  , 523  , 541+    , 547  , 557  , 563  , 569  , 571  , 577  , 587  , 593  , 599  , 601+    , 607  , 613  , 617  , 619  , 631  , 641  , 643  , 647  , 653  , 659+    , 661  , 673  , 677  , 683  , 691  , 701  , 709  , 719  , 727  , 733+    , 739  , 743  , 751  , 757  , 761  , 769  , 773  , 787  , 797  , 809+    , 811  , 821  , 823  , 827  , 829  , 839  , 853  , 857  , 859  , 863+    , 877  , 881  , 883  , 887  , 907  , 911  , 919  , 929  , 937  , 941+    , 947  , 953  , 967  , 971  , 977  , 983  , 991  , 997  , 1009 , 1013+    , 1019 , 1021 , 1031 , 1033 , 1039 , 1049 , 1051 , 1061 , 1063 , 1069+    , 1087 , 1091 , 1093 , 1097 , 1103 , 1109 , 1117 , 1123 , 1129 , 1151+    , 1153 , 1163 , 1171 , 1181 , 1187 , 1193 , 1201 , 1213 , 1217 , 1223+    , 1229 , 1231 , 1237 , 1249 , 1259 , 1277 , 1279 , 1283 , 1289 , 1291+    , 1297 , 1301 , 1303 , 1307 , 1319 , 1321 , 1327 , 1361 , 1367 , 1373+    , 1381 , 1399 , 1409 , 1423 , 1427 , 1429 , 1433 , 1439 , 1447 , 1451+    , 1453 , 1459 , 1471 , 1481 , 1483 , 1487 , 1489 , 1493 , 1499 , 1511+    , 1523 , 1531 , 1543 , 1549 , 1553 , 1559 , 1567 , 1571 , 1579 , 1583+    , 1597 , 1601 , 1607 , 1609 , 1613 , 1619 , 1621 , 1627 , 1637 , 1657+    , 1663 , 1667 , 1669 , 1693 , 1697 , 1699 , 1709 , 1721 , 1723 , 1733+    , 1741 , 1747 , 1753 , 1759 , 1777 , 1783 , 1787 , 1789 , 1801 , 1811+    , 1823 , 1831 , 1847 , 1861 , 1867 , 1871 , 1873 , 1877 , 1879 , 1889+    , 1901 , 1907 , 1913 , 1931 , 1933 , 1949 , 1951 , 1973 , 1979 , 1987+    , 1993 , 1997 , 1999 , 2003 , 2011 , 2017 , 2027 , 2029 , 2039 , 2053+    , 2063 , 2069 , 2081 , 2083 , 2087 , 2089 , 2099 , 2111 , 2113 , 2129+    , 2131 , 2137 , 2141 , 2143 , 2153 , 2161 , 2179 , 2203 , 2207 , 2213+    , 2221 , 2237 , 2239 , 2243 , 2251 , 2267 , 2269 , 2273 , 2281 , 2287+    , 2293 , 2297 , 2309 , 2311 , 2333 , 2339 , 2341 , 2347 , 2351 , 2357+    , 2371 , 2377 , 2381 , 2383 , 2389 , 2393 , 2399 , 2411 , 2417 , 2423+    , 2437 , 2441 , 2447 , 2459 , 2467 , 2473 , 2477 , 2503 , 2521 , 2531+    , 2539 , 2543 , 2549 , 2551 , 2557 , 2579 , 2591 , 2593 , 2609 , 2617+    , 2621 , 2633 , 2647 , 2657 , 2659 , 2663 , 2671 , 2677 , 2683 , 2687+    , 2689 , 2693 , 2699 , 2707 , 2711 , 2713 , 2719 , 2729 , 2731 , 2741+    , 2749 , 2753 , 2767 , 2777 , 2789 , 2791 , 2797 , 2801 , 2803 , 2819+    , 2833 , 2837 , 2843 , 2851 , 2857 , 2861 , 2879 , 2887 , 2897 , 2903+    ]++{-# INLINE divides #-}+divides :: Integer -> Integer -> Bool+divides i n = n `mod` i == 0
+ Crypto/Number/Serialize.hs view
@@ -0,0 +1,73 @@+module Crypto.Number.Serialize+    ( i2osp+    , os2ip+    , i2ospOf+    , i2ospOf_+    , lengthBytes+    ) where++import Data.ByteString (ByteString)+import qualified Data.ByteString as B+import qualified Data.ByteString.Internal as B+import Data.Bits+import Foreign.Storable+import Foreign.Ptr++{-# INLINE divMod256 #-}+divMod256 :: Integer -> (Integer, Integer)+divMod256 n = (n `shiftR` 8, n .&. 0xff)++-- | os2ip converts a byte string into a positive integer+{-# INLINE os2ip #-}+os2ip :: ByteString -> Integer+os2ip = B.foldl' (\a b -> (256 * a) .|. (fromIntegral b)) 0++-- | i2osp converts a positive integer into a byte string+i2osp :: Integer -> ByteString+i2osp m+    | m < 0     = error "i2osp: cannot convert a negative integer to a bytestring"+    | otherwise = B.reverse $ B.unfoldr fdivMod256 m+    where fdivMod256 0 = Nothing+          fdivMod256 n = Just (fromIntegral a,b) where (b,a) = divMod256 n+++-- | just like i2osp, but take an extra parameter for size.+-- if the number is too big to fit in @len bytes, nothing is returned+-- otherwise the number is padded with 0 to fit the @len required.+--+-- FIXME: use unsafeCreate to fill the bytestring+i2ospOf :: Int -> Integer -> Maybe ByteString+i2ospOf len m+    | lenbytes < len  = Just $ B.replicate (len - lenbytes) 0 `B.append` bytes+    | lenbytes == len = Just bytes+    | otherwise       = Nothing+    where+        lenbytes = B.length bytes+        bytes    = i2osp m++-- | just like i2ospOf except that it doesn't expect a failure.+-- for example if you just took a modulo of the number that represent+-- the size (example the RSA modulo n).+{-# INLINE i2ospOf_ #-}+i2ospOf_ :: Int -> Integer -> ByteString+i2ospOf_ len m = B.unsafeCreate len fillPtr+    where fillPtr srcPtr = loop m (srcPtr `plusPtr` (len-1))+            where loop n ptr = do+                      let (nn,a) = divMod256 n+                      poke ptr (fromIntegral a)+                      if ptr == srcPtr+                          then return ()+                          else (if nn == 0 then fillerLoop else loop nn) (ptr `plusPtr` (-1))+                  fillerLoop ptr = do+                      poke ptr 0+                      if ptr == srcPtr+                          then return ()+                          else fillerLoop (ptr `plusPtr` (-1))++-- | returns the number of bytes to store an integer with i2osp+--+-- FIXME: really slow implementation. use log or bigger shifts.+lengthBytes :: Integer -> Int+lengthBytes n+    | n < 256   = 1+    | otherwise = 1 + lengthBytes (n `shiftR` 8)
+ LICENSE view
@@ -0,0 +1,24 @@+Copyright (c) 2010-2012 Vincent Hanquez <vincent@snarc.org>++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:+1. Redistributions of source code must retain the above copyright+   notice, this list of conditions and the following disclaimer.+2. Redistributions in binary form must reproduce the above copyright+   notice, this list of conditions and the following disclaimer in the+   documentation and/or other materials provided with the distribution.++THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ Tests/RNG.hs view
@@ -0,0 +1,38 @@+module RNG where++import Data.Word+import Data.List (foldl')+import qualified Data.ByteString as B+import Crypto.Random.API+import Control.Arrow (first)++{- this is a just test rng. this is absolutely not a serious RNG. DO NOT use elsewhere -}+data Rng = Rng (Int, Int)++getByte :: Rng -> (Word8, Rng)+getByte (Rng (mz, mw)) = (r, g)+    where mz2 = 36969 * (mz `mod` 65536)+          mw2 = 18070 * (mw `mod` 65536)+          r   = fromIntegral (mz2 + mw2)+          g   = Rng (mz2, mw2)++getBytes :: Int -> Rng -> ([Word8], Rng)+getBytes 0 g = ([], g)+getBytes n g =+    let (b, g')  = getByte g+        (l, g'') = getBytes (n-1) g'+     in (b:l, g'')++instance CPRG Rng where+    cprgGenBytes g len    = first B.pack $ getBytes len g+    cprgSupplyEntropy g e = reseed e g+    cprgNeedReseed _      = maxBound++reseed :: B.ByteString -> Rng -> Rng+reseed bs (Rng (a,b)) = Rng (fromIntegral a', b')+        where a' = foldl' (\v i -> ((fromIntegral v) + (fromIntegral i) * 36969) `mod` 65536) a l+              b' = foldl' (\v i -> ((fromIntegral v) + (fromIntegral i) * 18070) `mod` 65536) b l+              l  = B.unpack bs++rng :: Rng+rng = Rng (1,2) 
+ Tests/Tests.hs view
@@ -0,0 +1,135 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE ViewPatterns #-}++import Test.Framework (defaultMain, testGroup)+import Test.Framework.Providers.QuickCheck2 (testProperty)+import Test.Framework.Providers.HUnit (testCase)++import Test.QuickCheck+import Test.HUnit+--import Test.QuickCheck.Test++import Control.Applicative ((<$>))++import qualified Data.ByteString as B++import Crypto.Number.ModArithmetic+import Crypto.Number.Basic+import Crypto.Number.Generate+import Crypto.Number.Prime+import Crypto.Number.Serialize++import RNG++prop_gcde_binary_valid :: (Positive Integer, Positive Integer) -> Bool+prop_gcde_binary_valid (Positive a, Positive b) =+    and [v==v', a*x' + b*y' == v', a*x + b*y == v, gcd a b == v]+    where (x,y,v)    = gcde_binary a b+          (x',y',v') = gcde a b++prop_modexp_rtl_valid :: (NonNegative Integer,+                          NonNegative Integer,+                          Positive Integer)+                      -> Bool+prop_modexp_rtl_valid (NonNegative a, NonNegative b, Positive m) =+    exponantiation_rtl_binary a b m == ((a ^ b) `mod` m)++prop_modinv_valid :: (Positive Integer, Positive Integer) -> Bool+prop_modinv_valid (Positive a, Positive m)+    | m > 1           = case inverse a m of+                             Just ainv -> (ainv * a) `mod` m == 1+                             Nothing   -> True+    | otherwise       = True++prop_sqrti_valid :: Positive Integer -> Bool+prop_sqrti_valid (Positive i) = l*l <= i && i <= u*u where (l, u) = sqrti i++prop_generate_prime_valid :: Seed -> Bool+prop_generate_prime_valid i =+    -- because of the next naive test, we can't generate easily number above 32 bits+    -- otherwise it slows down the tests to uselessness. when AKS or ECPP is implemented+    -- we can revisit the number here+    primalityTestNaive $ withRNG i (\g -> generatePrime g 32)++prop_miller_rabin_valid :: (Seed, PositiveSmall) -> Bool+prop_miller_rabin_valid (seed, PositiveSmall i)+    | i <= 3    = True+    | otherwise =+        -- miller rabin only returns with certitude that the integer is composite.+        let b = withRNG seed (\g -> isProbablyPrime g i)+         in (b == False && primalityTestNaive i == False) || b == True++prop_generate_valid :: (Seed, Positive Integer) -> Bool+prop_generate_valid (seed, Positive h) =+    let v = withRNG seed (\g -> generateMax g h)+     in (v >= 0 && v < h)++withAleasInteger :: Rng -> Seed -> (Rng -> (a,Rng)) -> a+withAleasInteger g (Seed i) f = fst $ f $ reseed (i2osp $ fromIntegral i) g++withRNG :: Seed -> (Rng -> (a,Rng)) -> a+withRNG seed f = withAleasInteger rng seed f++newtype PositiveSmall = PositiveSmall Integer+                      deriving (Show,Eq)++instance Arbitrary PositiveSmall where+    arbitrary = PositiveSmall . fromIntegral <$> (resize (2^(20 :: Int)) (arbitrary :: Gen Int))++data Range = Range Integer Integer+           deriving (Show,Eq)++instance Arbitrary Range where+    arbitrary = do (Positive x) <- arbitrary :: Gen (Positive Int)+                   (Positive r) <- arbitrary :: Gen (Positive Int)+                   return $ Range (fromIntegral x) (fromIntegral r)++newtype Seed = Seed Integer+             deriving (Eq)++instance Show Seed where+    show s = "Seed " ++ show s++instance Arbitrary Seed where+    arbitrary = arbitrary >>= \(Positive i) -> return (Seed i)++serializationKATTests = concatMap f vectors+    where f (v, bs) = [ testCase ("i2osp " ++ show v) (i2osp v  @=? bs)+                      , testCase ("os2ip " ++ show v) (os2ip bs @=? v)+                      ]+          vectors =+            [ (0x10000, "\SOH\NUL\NUL")+            , (0x1234, "\DC24")+            , (0xf123456, "\SI\DC24V")+            , (0xf21908421feabd21490, "\SI!\144\132!\254\171\210\DC4\144")+            , (0x7521908421feabd21490, "u!\144\132!\254\171\210\DC4\144")+            ]++main :: IO ()+main = defaultMain+    [ testGroup "serialization"+        [ testProperty "unbinary.binary==id" (\(Positive i) -> os2ip (i2osp i) == i)+        , testProperty "length integer" (\(Positive i) -> B.length (i2osp i) == lengthBytes i)+        , testGroup "KAT" serializationKATTests+        ]+    , testGroup "gcde binary"+        [ testProperty "gcde" prop_gcde_binary_valid+        ]+    , testGroup "exponantiation"+        [ testProperty "right-to-left" prop_modexp_rtl_valid+        ]+    , testGroup "inverse"+        [ testProperty "inverse" prop_modinv_valid+        ]+    , testGroup "sqrt integer"+        [ testProperty "sqrt" prop_sqrti_valid+        ]+    , testGroup "generation"+        [ testProperty "max" prop_generate_valid+        --, testProperty "between" (\seed (Range l h) -> let generated = withRNG seed (\rng -> generateBetween rng l (l+h))+        --                                                in (generated > l && generated < h))+        ]+    , testGroup "primality test"+        [ testProperty "miller-rabin" prop_miller_rabin_valid+        ]+    ]
+ crypto-numbers.cabal view
@@ -0,0 +1,58 @@+Name:                crypto-numbers+Version:             0.1.0+Description:         Cryptographic numbers: functions and algorithms+License:             BSD3+License-file:        LICENSE+Copyright:           Vincent Hanquez <vincent@snarc.org>+Author:              Vincent Hanquez <vincent@snarc.org>+Maintainer:          Vincent Hanquez <vincent@snarc.org>+Synopsis:            Cryptographic numbers: functions and algorithms+Category:            Cryptography+Build-Type:          Simple+Homepage:            http://github.com/vincenthz/hs-crypto-numbers+Cabal-Version:       >=1.8+Extra-Source-Files:  Tests/*.hs++Library+  Build-Depends:     base >= 4 && < 5+                   , bytestring+                   , vector+                   , crypto-random-api+  Exposed-modules:   Crypto.Number.ModArithmetic+                     Crypto.Number.Serialize+                     Crypto.Number.Generate+                     Crypto.Number.Basic+                     Crypto.Number.Polynomial+                     Crypto.Number.Prime+  ghc-options:       -Wall++Test-Suite test-crypto-numbers+  type:              exitcode-stdio-1.0+  hs-source-dirs:    Tests+  Main-Is:           Tests.hs+  Build-depends:     base >= 4 && < 5+                   , crypto-random-api+                   , crypto-numbers+                   , bytestring+                   , vector+                   , QuickCheck >= 2+                   , HUnit+                   , test-framework >= 0.3.3+                   , test-framework-quickcheck2 >= 0.2.9+                   , test-framework-hunit+  ghc-options:       -Wall -O2++Benchmark bench-crypto-numbers+  hs-source-dirs:    Benchmarks+  Main-Is:           Benchmarks.hs+  type:              exitcode-stdio-1.0+  Build-depends:     base >= 4 && < 5+                   , bytestring+                   , crypto-random-api+                   , crypto-numbers+                   , criterion+                   , mtl++source-repository head+  type:     git+  location: git://github.com/vincenthz/hs-crypto-numbers