diff --git a/Benchmarks/Benchmarks.hs b/Benchmarks/Benchmarks.hs --- a/Benchmarks/Benchmarks.hs +++ b/Benchmarks/Benchmarks.hs @@ -3,9 +3,10 @@ import Criterion.Main import Crypto.Number.Serialize -import Crypto.Number.Generate +-- import Crypto.Number.Generate import qualified Data.ByteString as B import Crypto.Number.ModArithmetic +import Crypto.Number.F2m import Data.Bits primes = [3, 5, 7, 29, 31, 211, 2309, 2311, 30029, 200560490131, 304250263527209] @@ -14,6 +15,7 @@ lg1, lg2 :: Integer lg1 = 21389083291083903845902381390285907190274907230982112390820985903825329874812973821790321904790217490217409721904832974210974921740972109481490128430982190472109874802174907490271904124908210958093285098309582093850918902581290859012850829105809128590218590281905812905810928590128509128940821903829018390849839578967358920127598901248259797158249684571948075896458741905823982671490352896791052386357019528367902 lg2 = 21392813098390824190840192812389082390812940821904891028439028490128904829104891208940835932882910839218309812093118249089871209347472901874902407219740921840928149087284397490128903843789289014374839281492038091283923091809734832974180398210938901284839274091749021709 +fx = 11692013098647223345629478661730264157247460344009 -- x^163+x^7+x^6+x^3+1 bitsAndShift8 n i = (n `shiftR` i, n .&. 0xff) modAndShift8 n i = (n `shiftR` i, n `mod` 0x100) @@ -59,6 +61,14 @@ , 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) + ] + , bgroup "F2m" + [ bench "addition" $ nf (addF2m lg1) lg2 + , bench "multiplication" $ nf (mulF2m fx lg1) lg2 + , bench "square" $ nf (squareF2m fx) lg1 + , bench "square multiplication" $ nf (mulF2m fx lg1) lg1 + , bench "reduction" $ nf (modF2m fx) lg1 + , bench "inversion" $ nf (invF2m fx) lg1 ] ] where b8 = B.replicate 8 0xf7 diff --git a/Crypto/Number/Basic.hs b/Crypto/Number/Basic.hs --- a/Crypto/Number/Basic.hs +++ b/Crypto/Number/Basic.hs @@ -1,4 +1,8 @@ {-# LANGUAGE BangPatterns #-} +{-# LANGUAGE CPP #-} +#if MIN_VERSION_integer_gmp(0,5,1) +{-# LANGUAGE UnboxedTuples #-} +#endif -- | -- Module : Crypto.Number.Basic -- License : BSD-style @@ -13,7 +17,11 @@ , areEven ) where +#if MIN_VERSION_integer_gmp(0,5,1) +import GHC.Integer.GMP.Internals +#else import Data.Bits +#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 @@ -48,17 +56,29 @@ 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' @@ -82,6 +102,7 @@ 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 diff --git a/Crypto/Number/F2m.hs b/Crypto/Number/F2m.hs new file mode 100644 --- /dev/null +++ b/Crypto/Number/F2m.hs @@ -0,0 +1,122 @@ +{-# LANGUAGE CPP #-} +#ifdef VERSION_integer_gmp +{-# LANGUAGE MagicHash #-} +#endif +-- | +-- Module : Crypto.Number.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) + +#ifdef VERSION_integer_gmp +import GHC.Exts +import GHC.Integer.Logarithms (integerLog2#) +#endif + +-- | 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. +mulF2m :: Integer -- ^ 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 :: Integer -- ^ 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 over F₂m using extended Euclidean algorithm. +invF2m :: Integer -- ^ Irreducible binary polynomial + -> Integer -> Maybe Integer +invF2m _ 0 = Nothing +invF2m fx n = go n fx 1 0 + where + go u v g1 g2 + | u == 1 = Just $ modF2m fx g1 + | otherwise = if j < 0 + then go u (v `xor` shift u (-j)) + g1 (g2 `xor` shift g1 (-j)) + else 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'. +divF2m :: Integer -- ^ Irreducible binary polynomial + -> Integer -- ^ Dividend + -> Integer -- ^ Quotient + -> Maybe Integer +divF2m fx n1 n2 = mulF2m fx n1 <$> invF2m fx n2 +{-# INLINE divF2m #-} + +log2 :: Integer -> Int +#if defined(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/ModArithmetic.hs b/Crypto/Number/ModArithmetic.hs --- a/Crypto/Number/ModArithmetic.hs +++ b/Crypto/Number/ModArithmetic.hs @@ -1,5 +1,6 @@ {-# LANGUAGE BangPatterns #-} {-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE CPP #-} -- | -- Module : Crypto.Number.ModArithmetic -- License : BSD-style @@ -8,16 +9,29 @@ -- Portability : Good module Crypto.Number.ModArithmetic - ( exponantiation_rtl_binary + ( + -- * 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 -import Data.Typeable +#endif -- | Raised when two numbers are supposed to be coprimes but are not. data CoprimesAssertionError = CoprimesAssertionError @@ -25,34 +39,97 @@ instance Exception CoprimesAssertionError --- note on exponantiation: 0^0 is treated as 1 for mimicking the standard library; +-- | 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. +expSafe :: Integer -- ^ base + -> Integer -- ^ exponant + -> Integer -- ^ modulo + -> Integer -- ^ result +#if MIN_VERSION_integer_gmp(0,5,1) +expSafe b e m + | odd m = powModSecInteger b e m + | 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. --- | exponantiation_rtl_binary computes modular exponantiation as b^e mod m +-- | exponentiation_rtl_binary computes modular exponentiation 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 +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 --- | exponantiation computes modular exponantiation as b^e mod m +-- | exponentiation computes modular exponentiation as b^e mod m -- using repetitive squaring. -exponantiation :: Integer -> Integer -> Integer -> Integer -exponantiation b e m +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 = (exponantiation b (e `div` 2) m) `mod` m + | even e = let p = (exponentiation b (e `div` 2) m) `mod` m in (p^(2::Integer)) `mod` m - | otherwise = (b * exponantiation b (e-1) m) `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 -inverse g m = if d > 1 then Nothing else Just (x `mod` m) - where (x,_,d) = gcde_binary g m +#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 diff --git a/Crypto/Number/Prime.hs b/Crypto/Number/Prime.hs --- a/Crypto/Number/Prime.hs +++ b/Crypto/Number/Prime.hs @@ -1,4 +1,8 @@ +{-# LANGUAGE CPP #-} {-# LANGUAGE BangPatterns #-} +#if MIN_VERSION_integer_gmp(0,5,1) +{-# LANGUAGE MagicHash #-} +#endif -- | -- Module : Crypto.Number.Prime -- License : BSD-style @@ -19,11 +23,17 @@ ) where import Crypto.Random.API -import Data.Bits 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 @@ -64,11 +74,22 @@ -- | 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 +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) @@ -105,6 +126,7 @@ | x2 == 1 = False | x2 /= nm1 = loop' ws ((x2*x2) `mod` n) (r+1) | otherwise = loop ws +#endif {- n < z -> witness to test diff --git a/Tests/Tests.hs b/Tests/Tests.hs --- a/Tests/Tests.hs +++ b/Tests/Tests.hs @@ -1,5 +1,4 @@ {-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE ViewPatterns #-} import Test.Framework (defaultMain, testGroup) import Test.Framework.Providers.QuickCheck2 (testProperty) @@ -19,6 +18,7 @@ import Crypto.Number.Generate import Crypto.Number.Prime import Crypto.Number.Serialize +import Crypto.Number.F2m import RNG @@ -65,6 +65,12 @@ let v = withRNG seed (\g -> generateMax g h) in (v >= 0 && v < h) +prop_invF2m_valid :: Fx -> PositiveLarge -> Bool +prop_invF2m_valid (Fx fx) (PositiveLarge a) = maybe True ((1 ==) . mulF2m fx a) (invF2m fx a) + +prop_squareF2m_valid :: Fx -> PositiveLarge -> Bool +prop_squareF2m_valid (Fx fx) (PositiveLarge a) = mulF2m fx a a == squareF2m fx a + withAleasInteger :: Rng -> Seed -> (Rng -> (a,Rng)) -> a withAleasInteger g (Seed i) f = fst $ f $ reseed (i2osp $ fromIntegral i) g @@ -77,6 +83,27 @@ instance Arbitrary PositiveSmall where arbitrary = PositiveSmall . fromIntegral <$> (resize (2^(20 :: Int)) (arbitrary :: Gen Int)) +newtype PositiveLarge = PositiveLarge Integer + deriving (Show,Eq) + +instance Arbitrary PositiveLarge where + arbitrary = PositiveLarge <$> sized (\n -> choose (1, fromIntegral n^(100::Int))) + +newtype Fx = Fx Integer deriving (Show,Eq) + +instance Arbitrary Fx where + arbitrary = elements $ map Fx + [ 283 -- [8,4,3,1,0] Rijndael + -- SEC2 polynomials + , 11692013098647223345629478661730264157247460344009 -- [163,7,6,3,0] + , 13803492693581127574869511724554050904902217944359662576256527028453377 -- [233,74,0] + , 883423532389192164791648750371459257913741948437809479060803169365786625 -- [239,36,0] + , 883423532389192164791649115746868590639471499359017658131558014629445633 -- [239,158,0] + , 15541351137805832567355695254588151253139254712417116170014499277911234281641667989665 -- [283,12,7,5,0] + , 1322111937580497197903830616065542079656809365928562438569297590548811582472622691650378420879430724437687334722581078999041 -- [409,87,0] + , 7729075046034516689390703781863974688597854659412869997314470502903038284579120849072387533163845155924927232063004354354730157322085975311485817346934161497393961629647909 -- [571,10,5,2,0] + ] + data Range = Range Integer Integer deriving (Show,Eq) @@ -132,5 +159,9 @@ ] , testGroup "primality test" [ testProperty "miller-rabin" prop_miller_rabin_valid + ] + , testGroup "F2m" + [ testProperty "invF2m" prop_invF2m_valid + , testProperty "squareF2m" prop_squareF2m_valid ] ] diff --git a/crypto-numbers.cabal b/crypto-numbers.cabal --- a/crypto-numbers.cabal +++ b/crypto-numbers.cabal @@ -1,5 +1,5 @@ Name: crypto-numbers -Version: 0.2.1 +Version: 0.2.2 Description: Cryptographic numbers: functions and algorithms License: BSD3 License-file: LICENSE @@ -13,6 +13,10 @@ Cabal-Version: >=1.8 Extra-Source-Files: Tests/*.hs +Flag integer-gmp + Description: Are we using integer-gmp? + Default: True + Library Build-Depends: base >= 4 && < 5 , bytestring @@ -23,7 +27,11 @@ Crypto.Number.Generate Crypto.Number.Basic Crypto.Number.Polynomial + Crypto.Number.F2m Crypto.Number.Prime + if impl(ghc) && flag(integer-gmp) + Build-depends: integer-gmp + , ghc-prim ghc-options: -Wall Test-Suite test-crypto-numbers @@ -49,7 +57,6 @@ type: exitcode-stdio-1.0 Build-depends: base >= 4 && < 5 , bytestring - , crypto-random , crypto-numbers , criterion , mtl