1{-# LANGUAGE BangPatterns #-} 2-- | 3-- Module : Crypto.Number.ModArithmetic 4-- License : BSD-style 5-- Maintainer : Vincent Hanquez <vincent@snarc.org> 6-- Stability : experimental 7-- Portability : Good 8 9module Crypto.Number.ModArithmetic 10 ( 11 -- * Exponentiation 12 expSafe 13 , expFast 14 -- * Inverse computing 15 , inverse 16 , inverseCoprimes 17 , inverseFermat 18 -- * Squares 19 , jacobi 20 , squareRoot 21 ) where 22 23import Control.Exception (throw, Exception) 24import Crypto.Number.Basic 25import Crypto.Number.Compat 26 27-- | Raised when two numbers are supposed to be coprimes but are not. 28data CoprimesAssertionError = CoprimesAssertionError 29 deriving (Show) 30 31instance Exception CoprimesAssertionError 32 33-- | Compute the modular exponentiation of base^exponent using 34-- algorithms design to avoid side channels and timing measurement 35-- 36-- Modulo need to be odd otherwise the normal fast modular exponentiation 37-- is used. 38-- 39-- When used with integer-simple, this function is not different 40-- from expFast, and thus provide the same unstudied and dubious 41-- timing and side channels claims. 42-- 43-- Before GHC 8.4.2, powModSecInteger is missing from integer-gmp, 44-- so expSafe has the same security as expFast. 45expSafe :: Integer -- ^ base 46 -> Integer -- ^ exponent 47 -> Integer -- ^ modulo 48 -> Integer -- ^ result 49expSafe b e m 50 | odd m = gmpPowModSecInteger b e m `onGmpUnsupported` 51 (gmpPowModInteger b e m `onGmpUnsupported` 52 exponentiation b e m) 53 | otherwise = gmpPowModInteger b e m `onGmpUnsupported` 54 exponentiation b e m 55 56-- | Compute the modular exponentiation of base^exponent using 57-- the fastest algorithm without any consideration for 58-- hiding parameters. 59-- 60-- Use this function when all the parameters are public, 61-- otherwise 'expSafe' should be preferred. 62expFast :: Integer -- ^ base 63 -> Integer -- ^ exponent 64 -> Integer -- ^ modulo 65 -> Integer -- ^ result 66expFast b e m = gmpPowModInteger b e m `onGmpUnsupported` exponentiation b e m 67 68-- | @exponentiation@ computes modular exponentiation as /b^e mod m/ 69-- using repetitive squaring. 70exponentiation :: Integer -> Integer -> Integer -> Integer 71exponentiation b e m 72 | b == 1 = b 73 | e == 0 = 1 74 | e == 1 = b `mod` m 75 | even e = let p = exponentiation b (e `div` 2) m `mod` m 76 in (p^(2::Integer)) `mod` m 77 | otherwise = (b * exponentiation b (e-1) m) `mod` m 78 79-- | @inverse@ computes the modular inverse as in /g^(-1) mod m/. 80inverse :: Integer -> Integer -> Maybe Integer 81inverse g m = gmpInverse g m `onGmpUnsupported` v 82 where 83 v 84 | d > 1 = Nothing 85 | otherwise = Just (x `mod` m) 86 (x,_,d) = gcde g m 87 88-- | Compute the modular inverse of two coprime numbers. 89-- This is equivalent to inverse except that the result 90-- is known to exists. 91-- 92-- If the numbers are not defined as coprime, this function 93-- will raise a 'CoprimesAssertionError'. 94inverseCoprimes :: Integer -> Integer -> Integer 95inverseCoprimes g m = 96 case inverse g m of 97 Nothing -> throw CoprimesAssertionError 98 Just i -> i 99 100-- | Computes the Jacobi symbol (a/n). 101-- 0 ≤ a < n; n ≥ 3 and odd. 102-- 103-- The Legendre and Jacobi symbols are indistinguishable exactly when the 104-- lower argument is an odd prime, in which case they have the same value. 105-- 106-- See algorithm 2.149 in "Handbook of Applied Cryptography" by Alfred J. Menezes et al. 107jacobi :: Integer -> Integer -> Maybe Integer 108jacobi a n 109 | n < 3 || even n = Nothing 110 | a == 0 || a == 1 = Just a 111 | n <= a = jacobi (a `mod` n) n 112 | a < 0 = 113 let b = if n `mod` 4 == 1 then 1 else -1 114 in fmap (*b) (jacobi (-a) n) 115 | otherwise = 116 let (e, a1) = asPowerOf2AndOdd a 117 nMod8 = n `mod` 8 118 nMod4 = n `mod` 4 119 a1Mod4 = a1 `mod` 4 120 s' = if even e || nMod8 == 1 || nMod8 == 7 then 1 else -1 121 s = if nMod4 == 3 && a1Mod4 == 3 then -s' else s' 122 n1 = n `mod` a1 123 in if a1 == 1 then Just s 124 else fmap (*s) (jacobi n1 a1) 125 126-- | Modular inverse using Fermat's little theorem. This works only when 127-- the modulus is prime but avoids side channels like in 'expSafe'. 128inverseFermat :: Integer -> Integer -> Integer 129inverseFermat g p = expSafe g (p - 2) p 130 131-- | Raised when the assumption about the modulus is invalid. 132data ModulusAssertionError = ModulusAssertionError 133 deriving (Show) 134 135instance Exception ModulusAssertionError 136 137-- | Modular square root of @g@ modulo a prime @p@. 138-- 139-- If the modulus is found not to be prime, the function will raise a 140-- 'ModulusAssertionError'. 141-- 142-- This implementation is variable time and should be used with public 143-- parameters only. 144squareRoot :: Integer -> Integer -> Maybe Integer 145squareRoot p 146 | p < 2 = throw ModulusAssertionError 147 | otherwise = 148 case p `divMod` 8 of 149 (v, 3) -> method1 (2 * v + 1) 150 (v, 7) -> method1 (2 * v + 2) 151 (u, 5) -> method2 u 152 (_, 1) -> tonelliShanks p 153 (0, 2) -> \a -> Just (if even a then 0 else 1) 154 _ -> throw ModulusAssertionError 155 156 where 157 x `eqMod` y = (x - y) `mod` p == 0 158 159 validate g y | (y * y) `eqMod` g = Just y 160 | otherwise = Nothing 161 162 -- p == 4u + 3 and u' == u + 1 163 method1 u' g = 164 let y = expFast g u' p 165 in validate g y 166 167 -- p == 8u + 5 168 method2 u g = 169 let gamma = expFast (2 * g) u p 170 g_gamma = g * gamma 171 i = (2 * g_gamma * gamma) `mod` p 172 y = (g_gamma * (i - 1)) `mod` p 173 in validate g y 174 175tonelliShanks :: Integer -> Integer -> Maybe Integer 176tonelliShanks p a 177 | aa == 0 = Just 0 178 | otherwise = 179 case expFast aa p2 p of 180 b | b == p1 -> Nothing 181 | b == 1 -> Just $ go (expFast aa ((s + 1) `div` 2) p) 182 (expFast aa s p) 183 (expFast n s p) 184 e 185 | otherwise -> throw ModulusAssertionError 186 where 187 aa = a `mod` p 188 p1 = p - 1 189 p2 = p1 `div` 2 190 n = findN 2 191 192 x `mul` y = (x * y) `mod` p 193 194 pow2m 0 x = x 195 pow2m i x = pow2m (i - 1) (x `mul` x) 196 197 (e, s) = asPowerOf2AndOdd p1 198 199 -- find a quadratic non-residue 200 findN i 201 | expFast i p2 p == p1 = i 202 | otherwise = findN (i + 1) 203 204 -- find m such that b^(2^m) == 1 (mod p) 205 findM b i 206 | b == 1 = i 207 | otherwise = findM (b `mul` b) (i + 1) 208 209 go !x b g !r 210 | b == 1 = x 211 | otherwise = 212 let r' = findM b 0 213 z = pow2m (r - r' - 1) g 214 x' = x `mul` z 215 b' = b `mul` g' 216 g' = z `mul` z 217 in go x' b' g' r' 218