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