From 90d02607ba4a26b2870537f2be086f69a52340c7 Mon Sep 17 00:00:00 2001 From: Vincent Hanquez Date: Mon, 9 Feb 2015 05:47:11 +0000 Subject: [PATCH] merge crypto-numbers minus all the random parts --- Crypto/Math/Polynomial.hs | 133 +++++++++++++++++++ Crypto/Number/Basic.hs | 135 ++++++++++++++++++++ Crypto/Number/F2m.hs | 110 ++++++++++++++++ Crypto/Number/Generate.hs | 56 ++++++++ Crypto/Number/ModArithmetic.hs | 153 ++++++++++++++++++++++ Crypto/Number/Prime.hs | 226 +++++++++++++++++++++++++++++++++ Crypto/Number/Serialize.hs | 166 ++++++++++++++++++++++++ cryptonite.cabal | 6 + 8 files changed, 985 insertions(+) create mode 100644 Crypto/Math/Polynomial.hs create mode 100644 Crypto/Number/Basic.hs create mode 100644 Crypto/Number/F2m.hs create mode 100644 Crypto/Number/Generate.hs create mode 100644 Crypto/Number/ModArithmetic.hs create mode 100644 Crypto/Number/Prime.hs create mode 100644 Crypto/Number/Serialize.hs diff --git a/Crypto/Math/Polynomial.hs b/Crypto/Math/Polynomial.hs new file mode 100644 index 0000000..133dca3 --- /dev/null +++ b/Crypto/Math/Polynomial.hs @@ -0,0 +1,133 @@ +{-# LANGUAGE BangPatterns #-} +-- | +-- Module : Crypto.Math.Polynomial +-- License : BSD-style +-- Maintainer : Vincent Hanquez +-- Stability : experimental +-- Portability : Good + +module Crypto.Math.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) diff --git a/Crypto/Number/Basic.hs b/Crypto/Number/Basic.hs new file mode 100644 index 0000000..618e7e9 --- /dev/null +++ b/Crypto/Number/Basic.hs @@ -0,0 +1,135 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE CPP #-} +#ifndef MIN_VERSION_integer_gmp +#define MIN_VERSION_integer_gmp(a,b,c) 0 +#endif +#if MIN_VERSION_integer_gmp(0,5,1) +{-# LANGUAGE UnboxedTuples #-} +#endif +#ifdef VERSION_integer_gmp +{-# LANGUAGE MagicHash #-} +#endif +-- | +-- Module : Crypto.Number.Basic +-- License : BSD-style +-- Maintainer : Vincent Hanquez +-- Stability : experimental +-- Portability : Good + +module Crypto.Number.Basic + ( sqrti + , gcde + , gcde_binary + , areEven + , log2 + ) where + +#if MIN_VERSION_integer_gmp(0,5,1) +import GHC.Integer.GMP.Internals +#else +import Data.Bits +#endif +#ifdef VERSION_integer_gmp +import GHC.Exts +import GHC.Integer.Logarithms (integerLog2#) +#endif + +-- | 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 'a' 'b' find (x,y,gcd(a,b)) where ax + by = d +-- +gcde :: Integer -> Integer -> (Integer, Integer, Integer) +#if MIN_VERSION_integer_gmp(0,5,1) +gcde a b = (s, t, g) + where (# g, s #) = gcdExtInteger a b + t = (g - s * a) `div` b +#else +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)) +#endif + +-- | 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) +#if MIN_VERSION_integer_gmp(0,5,1) +gcde_binary = gcde +#else +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) +#endif + +-- | check if a list of integer are all even +areEven :: [Integer] -> Bool +areEven = and . map even + +log2 :: Integer -> Int +#ifdef VERSION_integer_gmp +log2 0 = 0 +log2 x = I# (integerLog2# x) +#else +-- http://www.haskell.org/pipermail/haskell-cafe/2008-February/039465.html +log2 = imLog 2 + where + imLog b x = if x < b then 0 else (x `div` b^l) `doDiv` l + where + l = 2 * imLog (b * b) x + doDiv x' l' = if x' < b then l' else (x' `div` b) `doDiv` (l' + 1) +#endif +{-# INLINE log2 #-} diff --git a/Crypto/Number/F2m.hs b/Crypto/Number/F2m.hs new file mode 100644 index 0000000..ee62c96 --- /dev/null +++ b/Crypto/Number/F2m.hs @@ -0,0 +1,110 @@ +{-# LANGUAGE CPP #-} +#ifndef MIN_VERSION_integer_gmp +#define MIN_VERSION_integer_gmp(a,b,c) 0 +#endif +-- | +-- Module : Crypto.Math.F2m +-- License : BSD-style +-- Maintainer : Danny Navarro +-- Stability : experimental +-- Portability : Good +-- +-- This module provides basic arithmetic operations over F₂m. Performance is +-- not optimal and it doesn't provide protection against timing +-- attacks. The 'm' parameter is implicitly derived from the irreducible +-- polynomial where applicable. +module Crypto.Number.F2m + ( addF2m + , mulF2m + , squareF2m + , modF2m + , invF2m + , divF2m + ) where + +import Control.Applicative ((<$>)) +import Data.Bits ((.&.),(.|.),xor,shift,testBit) +import Crypto.Number.Basic + +type BinaryPolynomial = Integer + +-- | Addition over F₂m. This is just a synonym of 'xor'. +addF2m :: Integer -> Integer -> Integer +addF2m = xor +{-# INLINE addF2m #-} + +-- | Binary polynomial reduction modulo using long division algorithm. +modF2m :: Integer -- ^ Irreducible binary polynomial + -> Integer -> Integer +modF2m fx = go + where + lfx = log2 fx + go n | s == 0 = n `xor` fx + | s < 0 = n + | otherwise = go $ n `xor` shift fx s + where + s = log2 n - lfx +{-# INLINE modF2m #-} + +-- | Multiplication over F₂m. +-- +-- n1 * n2 (in F(2^m)) +mulF2m :: BinaryPolynomial -- ^ Irreducible binary polynomial + -> Integer -> Integer -> Integer +mulF2m fx n1 n2 = modF2m fx + $ go (if n2 `mod` 2 == 1 then n1 else 0) (log2 n2) + where + go n s | s == 0 = n + | otherwise = if testBit n2 s + then go (n `xor` shift n1 s) (s - 1) + else go n (s - 1) +{-# INLINABLE mulF2m #-} + +-- | Squaring over F₂m. +-- TODO: This is still slower than @mulF2m@. + +-- Multiplication table? C? +squareF2m :: BinaryPolynomial -- ^ Irreducible binary polynomial + -> Integer -> Integer +squareF2m fx = modF2m fx . square +{-# INLINE squareF2m #-} + +square :: Integer -> Integer +square n1 = go n1 ln1 + where + ln1 = log2 n1 + go n s | s == 0 = n + | otherwise = go (x .|. y) (s - 1) + where + x = shift (shift n (2 * (s - ln1) - 1)) (2 * (ln1 - s) + 2) + y = n .&. (shift 1 (2 * (ln1 - s) + 1) - 1) +{-# INLINE square #-} + +-- | Inversion of @n over F₂m using extended Euclidean algorithm. +-- +-- If @n doesn't have an inverse, Nothing is returned. +invF2m :: BinaryPolynomial -- ^ Irreducible binary polynomial + -> Integer -> Maybe Integer +invF2m _ 0 = Nothing +invF2m fx n + | n >= fx = Nothing + | otherwise = go n fx 1 0 + where + go u v g1 g2 + | u == 1 = Just $ modF2m fx g1 + | j < 0 = go u (v `xor` shift u (-j)) g1 (g2 `xor` shift g1 (-j)) + | otherwise = go (u `xor` shift v j) v (g1 `xor` shift g2 j) g2 + where + j = log2 u - log2 v +{-# INLINABLE invF2m #-} + +-- | Division over F₂m. If the dividend doesn't have an inverse it returns +-- 'Nothing'. +-- +-- Compute n1 / n2 +divF2m :: BinaryPolynomial -- ^ Irreducible binary polynomial + -> Integer -- ^ Dividend + -> Integer -- ^ Quotient + -> Maybe Integer +divF2m fx n1 n2 = mulF2m fx n1 <$> invF2m fx n2 +{-# INLINE divF2m #-} diff --git a/Crypto/Number/Generate.hs b/Crypto/Number/Generate.hs new file mode 100644 index 0000000..5d6fea3 --- /dev/null +++ b/Crypto/Number/Generate.hs @@ -0,0 +1,56 @@ +-- | +-- Module : Crypto.Number.Generate +-- License : BSD-style +-- Maintainer : Vincent Hanquez +-- Stability : experimental +-- Portability : Good + +module Crypto.Number.Generate + ( {-generateMax + , generateBetween + , generateOfSize + , generateBits-} + ) where + +import Crypto.Number.Basic +import Crypto.Number.Serialize +import qualified Data.ByteString as B +import Data.Bits ((.|.), (.&.), shiftR) + +{- +-- | generate a positive integer x, s.t. 0 <= x < m +generateMax :: CPRG g => g -> Integer -> (Integer, g) +generateMax rng 1 = (0, rng) +generateMax rng m + | (result' >= m) = generateMax rng' m + | otherwise = (result', rng') + where + bytesLength = lengthBytes m + bitsLength = (log2 (m-1) + 1) + bitsPoppedOff = 8 - (bitsLength `mod` 8) + randomInt bytes = withRandomBytes rng bytes $ \bs -> os2ip bs + + (result, rng') = randomInt bytesLength + result' = result `shiftR` bitsPoppedOff + +-- | 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 + +-- | Generate a number with the specified number of bits +generateBits :: CPRG g => g -> Int -> (Integer, g) +generateBits rng nbBits = withRandomBytes rng nbBytes' $ \bs -> modF (os2ip bs) + where (nbBytes, strayBits) = nbBits `divMod` 8 + nbBytes' | strayBits == 0 = nbBytes + | otherwise = nbBytes + 1 + modF | strayBits == 0 = id + | otherwise = (.&.) (2^nbBits - 1) +-} diff --git a/Crypto/Number/ModArithmetic.hs b/Crypto/Number/ModArithmetic.hs new file mode 100644 index 0000000..ab0ea0f --- /dev/null +++ b/Crypto/Number/ModArithmetic.hs @@ -0,0 +1,153 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE CPP #-} +#ifndef MIN_VERSION_integer_gmp +#define MIN_VERSION_integer_gmp(a,b,c) 0 +#endif +-- | +-- Module : Crypto.Number.ModArithmetic +-- License : BSD-style +-- Maintainer : Vincent Hanquez +-- Stability : experimental +-- Portability : Good + +module Crypto.Number.ModArithmetic + ( + -- * exponentiation + expSafe + , expFast + , exponentiation_rtl_binary + , exponentiation + -- * deprecated name for exponentiation + , exponantiation_rtl_binary + , exponantiation + -- * inverse computing + , inverse + , inverseCoprimes + ) where + +import Control.Exception (throw, Exception) +import Data.Typeable + +#if MIN_VERSION_integer_gmp(0,5,1) +import GHC.Integer.GMP.Internals +#else +import Crypto.Number.Basic (gcde_binary) +import Data.Bits +#endif + +-- | Raised when two numbers are supposed to be coprimes but are not. +data CoprimesAssertionError = CoprimesAssertionError + deriving (Show,Typeable) + +instance Exception CoprimesAssertionError + +-- | Compute the modular exponentiation of base^exponant using +-- algorithms design to avoid side channels and timing measurement +-- +-- Modulo need to be odd otherwise the normal fast modular exponentiation +-- is used. +-- +-- When used with integer-simple, this function is not different +-- from expFast, and thus provide the same unstudied and dubious +-- timing and side channels claims. +-- +-- with GHC 7.10, the powModSecInteger is missing from integer-gmp +-- (which is now integer-gmp2), so is has the same security as old +-- ghc version. +expSafe :: Integer -- ^ base + -> Integer -- ^ exponant + -> Integer -- ^ modulo + -> Integer -- ^ result +#if MIN_VERSION_integer_gmp(0,5,1) +expSafe b e m +#if !(MIN_VERSION_integer_gmp(1,0,0)) + | odd m = powModSecInteger b e m +#endif + | otherwise = powModInteger b e m +#else +expSafe = exponentiation +#endif + +-- | Compute the modular exponentiation of base^exponant using +-- the fastest algorithm without any consideration for +-- hiding parameters. +-- +-- Use this function when all the parameters are public, +-- otherwise 'expSafe' should be prefered. +expFast :: Integer -- ^ base + -> Integer -- ^ exponant + -> Integer -- ^ modulo + -> Integer -- ^ result +expFast = +#if MIN_VERSION_integer_gmp(0,5,1) + powModInteger +#else + exponentiation +#endif + +-- note on exponentiation: 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. + +-- | exponentiation_rtl_binary computes modular exponentiation as b^e mod m +-- using the right-to-left binary exponentiation algorithm (HAC 14.79) +exponentiation_rtl_binary :: Integer -> Integer -> Integer -> Integer +#if MIN_VERSION_integer_gmp(0,5,1) +exponentiation_rtl_binary = expSafe +#else +exponentiation_rtl_binary 0 0 m = 1 `mod` m +exponentiation_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) +#endif + +-- | exponentiation computes modular exponentiation as b^e mod m +-- using repetitive squaring. +exponentiation :: Integer -> Integer -> Integer -> Integer +#if MIN_VERSION_integer_gmp(0,5,1) +exponentiation = expSafe +#else +exponentiation b e m + | b == 1 = b + | e == 0 = 1 + | e == 1 = b `mod` m + | even e = let p = (exponentiation b (e `div` 2) m) `mod` m + in (p^(2::Integer)) `mod` m + | otherwise = (b * exponentiation b (e-1) m) `mod` m +#endif + +--{-# DEPRECATED exponantiation_rtl_binary "typo in API name it's called exponentiation_rtl_binary #-} +exponantiation_rtl_binary :: Integer -> Integer -> Integer -> Integer +exponantiation_rtl_binary = exponentiation_rtl_binary + +--{-# DEPRECATED exponentiation "typo in API name it's called exponentiation #-} +exponantiation :: Integer -> Integer -> Integer -> Integer +exponantiation = exponentiation + +-- | inverse computes the modular inverse as in g^(-1) mod m +inverse :: Integer -> Integer -> Maybe Integer +#if MIN_VERSION_integer_gmp(0,5,1) +inverse g m + | r == 0 = Nothing + | otherwise = Just r + where r = recipModInteger g m +#else +inverse g m + | d > 1 = Nothing + | otherwise = Just (x `mod` m) + where (x,_,d) = gcde_binary g m +#endif + +-- | 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 diff --git a/Crypto/Number/Prime.hs b/Crypto/Number/Prime.hs new file mode 100644 index 0000000..8909901 --- /dev/null +++ b/Crypto/Number/Prime.hs @@ -0,0 +1,226 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE BangPatterns #-} +#ifndef MIN_VERSION_integer_gmp +#define MIN_VERSION_integer_gmp(a,b,c) 0 +#endif +#if MIN_VERSION_integer_gmp(0,5,1) +{-# LANGUAGE MagicHash #-} +#endif +-- | +-- Module : Crypto.Number.Prime +-- License : BSD-style +-- Maintainer : Vincent Hanquez +-- Stability : experimental +-- Portability : Good + +module Crypto.Number.Prime + ( +{- + generatePrime + , generateSafePrime + , isProbablyPrime + , findPrimeFrom + , findPrimeFromWith + , primalityTestMillerRabin +-} + primalityTestNaive + , primalityTestFermat + , isCoprime + ) where + +import Crypto.Number.Generate +import Crypto.Number.Basic (sqrti, gcde_binary) +import Crypto.Number.ModArithmetic (exponantiation) + +#if MIN_VERSION_integer_gmp(0,5,1) +import GHC.Integer.GMP.Internals +import GHC.Base +#else +import Data.Bits +#endif + +{- +-- | 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 = +#if MIN_VERSION_integer_gmp(0,5,1) + (nextPrimeInteger n, rng) +#else + findPrimeFromWith rng (\g _ -> (True, g)) n +#endif + +-- | 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) +#if MIN_VERSION_integer_gmp(0,5,1) +primalityTestMillerRabin rng (I# tries) !n = + case testPrimeInteger n tries of + 0# -> (False, rng) + _ -> (True, rng) +#else +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 +#endif +-} + +{- + 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 diff --git a/Crypto/Number/Serialize.hs b/Crypto/Number/Serialize.hs new file mode 100644 index 0000000..28e9124 --- /dev/null +++ b/Crypto/Number/Serialize.hs @@ -0,0 +1,166 @@ +{-# LANGUAGE CPP #-} +#ifndef MIN_VERSION_integer_gmp +#define MIN_VERSION_integer_gmp(a,b,c) 0 +#endif +#if MIN_VERSION_integer_gmp(0,5,1) +{-# LANGUAGE MagicHash, UnboxedTuples, BangPatterns #-} +#endif +-- | +-- Module : Crypto.Number.Serialize +-- License : BSD-style +-- Maintainer : Vincent Hanquez +-- Stability : experimental +-- Portability : Good +-- +-- fast serialization primitives for integer +module Crypto.Number.Serialize + ( i2osp + , os2ip + , i2ospOf + , i2ospOf_ + , lengthBytes + ) where + +import Data.ByteString (ByteString) +import qualified Data.ByteString.Internal as B +import qualified Data.ByteString as B +import Foreign.Ptr + +#if MIN_VERSION_integer_gmp(0,5,1) +#if __GLASGOW_HASKELL__ >= 710 +import Control.Monad (void) +#endif +import GHC.Integer.GMP.Internals +import GHC.Base +import GHC.Ptr +import System.IO.Unsafe +import Foreign.ForeignPtr +#else +import Foreign.Storable +import Data.Bits +#endif + +#if !MIN_VERSION_integer_gmp(0,5,1) +{-# INLINE divMod256 #-} +divMod256 :: Integer -> (Integer, Integer) +divMod256 n = (n `shiftR` 8, n .&. 0xff) +#endif + +-- | os2ip converts a byte string into a positive integer +os2ip :: ByteString -> Integer +#if MIN_VERSION_integer_gmp(0,5,1) +os2ip bs = unsafePerformIO $ withForeignPtr fptr $ \ptr -> + let !(Ptr ad) = (ptr `plusPtr` ofs) +#if __GLASGOW_HASKELL__ >= 710 + in importIntegerFromAddr ad (int2Word# n) 1# +#else + in IO $ \s -> importIntegerFromAddr ad (int2Word# n) 1# s +#endif + where !(fptr, ofs, !(I# n)) = B.toForeignPtr bs +{-# NOINLINE os2ip #-} +#else +os2ip = B.foldl' (\a b -> (256 * a) .|. (fromIntegral b)) 0 +{-# INLINE os2ip #-} +#endif + +-- | i2osp converts a positive integer into a byte string +i2osp :: Integer -> ByteString +#if MIN_VERSION_integer_gmp(0,5,1) +i2osp 0 = B.singleton 0 +i2osp m = B.unsafeCreate (I# (word2Int# sz)) fillPtr + where !sz = sizeInBaseInteger m 256# +#if __GLASGOW_HASKELL__ >= 710 + fillPtr (Ptr srcAddr) = void $ exportIntegerToAddr m srcAddr 1# +#else + fillPtr (Ptr srcAddr) = IO $ \s -> case exportIntegerToAddr m srcAddr 1# s of + (# s2, _ #) -> (# s2, () #) +#endif +{-# NOINLINE i2osp #-} +#else +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 +#endif + + +-- | 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 +#if MIN_VERSION_integer_gmp(0,5,1) +i2ospOf len m + | sz <= len = Just $ i2ospOf_ len m + | otherwise = Nothing + where !sz = I# (word2Int# (sizeInBaseInteger m 256#)) +#else +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 +#endif + +-- | just like i2ospOf except that it doesn't expect a failure: i.e. +-- an integer larger than the number of output bytes requested +-- +-- for example if you just took a modulo of the number that represent +-- the size (example the RSA modulo n). +i2ospOf_ :: Int -> Integer -> ByteString +#if MIN_VERSION_integer_gmp(0,5,1) +i2ospOf_ len m = unsafePerformIO $ B.create len fillPtr + where !sz = (sizeInBaseInteger m 256#) + isz = I# (word2Int# sz) + fillPtr ptr + | len < isz = error "cannot compute i2ospOf_ with integer larger than output bytes" + | len == isz = + let !(Ptr srcAddr) = ptr in +#if __GLASGOW_HASKELL__ >= 710 + void (exportIntegerToAddr m srcAddr 1#) +#else + IO $ \s -> case exportIntegerToAddr m srcAddr 1# s of + (# s2, _ #) -> (# s2, () #) +#endif + | otherwise = do + let z = len-isz + _ <- B.memset ptr 0 (fromIntegral len) + let !(Ptr addr) = ptr `plusPtr` z +#if __GLASGOW_HASKELL__ >= 710 + void (exportIntegerToAddr m addr 1#) +#else + IO $ \s -> case exportIntegerToAddr m addr 1# s of + (# s2, _ #) -> (# s2, () #) +#endif +{-# NOINLINE i2ospOf_ #-} +#else +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)) +{-# INLINE i2ospOf_ #-} +#endif + +-- | returns the number of bytes to store an integer with i2osp +-- +-- with integer-simple, this function is really slow. +lengthBytes :: Integer -> Int +#if MIN_VERSION_integer_gmp(0,5,1) +lengthBytes n = I# (word2Int# (sizeInBaseInteger n 256#)) +#else +lengthBytes n + | n < 256 = 1 + | otherwise = 1 + lengthBytes (n `shiftR` 8) +#endif diff --git a/cryptonite.cabal b/cryptonite.cabal index bdeef42..f912839 100644 --- a/cryptonite.cabal +++ b/cryptonite.cabal @@ -32,6 +32,12 @@ Library Crypto.Cipher.RC4 Crypto.MAC.Poly1305 Crypto.MAC.HMAC + Crypto.Number.Basic + Crypto.Number.F2m + Crypto.Number.Generate + Crypto.Number.ModArithmetic + Crypto.Number.Prime + Crypto.Number.Serialize Crypto.KDF.PBKDF2 Crypto.KDF.Scrypt Crypto.Hash