module Math.Algebra.ModP where
import Data.Bits
import Data.Ratio
import Data.Int
p :: Int64
p = 2^31 1
newtype Zp = Zp Int64 deriving (Eq, Show)
fromZp :: Zp -> Int
fromZp (Zp k) = fromIntegral k
mkZp :: Integral a => a -> Zp
mkZp n = Zp (mod (fromIntegral n) p)
instance Num Zp where
(+) = addZp
() = subZp
(*) = mulZp
fromInteger = mkZp . fromInteger
abs = id
signum _ = Zp 1
instance Fractional Zp where
recip (Zp a) = mkZp $ invZp_euclid a
a / b = a * recip b
fromRational r = fromInteger (numerator r) / fromInteger (denominator r)
addZp :: Zp -> Zp -> Zp
addZp (Zp a) (Zp b)
| c < 0 = Zp (c p)
| c >= p = Zp (c p)
| otherwise = Zp c
where
c = a + b
subZp :: Zp -> Zp -> Zp
subZp (Zp a) (Zp b) = Zp (if b<=a then ab else a+pb)
mulZp :: Zp -> Zp -> Zp
mulZp (Zp a0) (Zp b0) = Zp (fromInteger c) where
a = fromIntegral a0 :: Integer
b = fromIntegral b0 :: Integer
c = mod (a * b) (fromIntegral p)
invZp_euclid :: Int64 -> Int64
invZp_euclid a
| a == 0 = 0
| otherwise = go 1 0 a p
where
modp :: Int64 -> Int64
modp n = mod n p
halfp1 = shiftR (p+1) 1
go :: Int64 -> Int64 -> Int64 -> Int64 -> Int64
go !x1 !x2 !u !v
| u==1 = x1
| v==1 = x2
| otherwise = stepU x1 x2 u v
stepU :: Int64 -> Int64 -> Int64 -> Int64 -> Int64
stepU !x1 !x2 !u !v = if even u
then let u' = shiftR u 1
x1' = if even x1 then shiftR x1 1 else shiftR x1 1 + halfp1
in stepU x1' x2 u' v
else stepV x1 x2 u v
stepV :: Int64 -> Int64 -> Int64 -> Int64 -> Int64
stepV !x1 !x2 !u !v = if even v
then let v' = shiftR v 1
x2' = if even x2 then shiftR x2 1 else shiftR x2 1 + halfp1
in stepV x1 x2' u v'
else final x1 x2 u v
final :: Int64 -> Int64 -> Int64 -> Int64 -> Int64
final !x1 !x2 !u !v = if u>=v
then let u' = uv
x1' = if x1 >= x2 then modp (x1x2) else modp (x1+px2)
in go x1' x2 u' v
else let v' = vu
x2' = if x2 >= x1 then modp (x2x1) else modp (x2+px1)
in go x1 x2' u v'