fastPower :: (Integral a) => a -> a -> a fastPower a n | n < 0 = error "Fast power with negative exponent" | otherwise = fastPower' 1 a n -- Invariant: acc * a^n = const fastPower' :: (Integral a) => a -> a -> a -> a fastPower' acc a 0 = acc fastPower' acc a n | even n = fastPower' acc (a*a) (n `div` 2) | otherwise = fastPower' (acc*a) a (n - 1) powerMod :: (Integral a) => a -> a -> a -> a powerMod a n m | n < 0 = error "Fast power with negative exponent" | otherwise = powerMod' 1 a n m -- Invariant: acc * a^n (mod m) = const powerMod' :: (Integral a) => a -> a -> a -> a -> a powerMod' acc _ 0 _ = acc powerMod' acc a n m | even n = powerMod' acc (a*a `mod` m) (n `div` 2) m | otherwise = powerMod' (acc*a `mod` m) a (n - 1) m -- Maximal prime number lesser than 2^64 pMax64 = 18446744073709551557 rootPMax = 9223372036854775778 randomNumber :: (Integral a) => a -> a randomNumber n = powerMod (fromIntegral rootPMax) n (fromIntegral pMax64) pollardPFactor :: (Integral a) => a -> a pollardPFactor m = pollardPFactor' b 1 n m where b = 97 n = 10000000 pollardPFactor' :: (Integral a) => a -> a -> a -> a -> a pollardPFactor' p k n m | k >= n = d | d /= 1 = d | otherwise = pollardPFactor' (powerMod p (k+1) m) (k+1) n m where d = gcd (p-1) m extEucl :: (Integral a) => a -> a -> (a, a, a) extEucl 0 0 = error "gcd(0, 0) is undefined" extEucl m n = extEucl' (m, n, 1, 0, 0, 1) -- Invariant -- gcd(m, n) == gcd(a, b) -- a == u1*m + v1*n -- b == u2*m + v2*n -------------------------- -- (a, b) <- (b, r) -- a = q*b + r, r = a - q*b -- u1' = u2, v1' = v2 -- r = a - q*b = (u1*m + v1*n) - q*(u2*m + v2*n) = -- = (u1 - q*u2)*m + (v1 - q*v2)*n -- u2' = u1 - q*u2, v2' = v1 - q*v2 ------------------------------------ extEucl' :: (Integral n) => (n, n, n, n, n, n) -> (n, n, n) extEucl' (a, 0, u1, v1, _, _) | a > 0 = (a, u1, v1) | otherwise = ((-a), (-u1), (-v1)) extEucl' (a, b, u1, v1, u2, v2) = let q = a `div` b r = a `mod` b u1' = u2 v1' = v2 u2' = u1 - q*u2 v2' = v1 - q*v2 in extEucl' (b, r, u1', v1', u2', v2') invmod :: (Integral a) => a -> a -> a invmod x 0 = error "invmod: zero module" invmod x m | d /= 1 = 0 -- element x is not invertible modulo m | otherwise = u where (d, u, _) = extEucl x m