{-# LANGUAGE PatternGuards #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Codec.Encryption.ECC.Base
-- Copyright   :  (c) Marcel Fourné 2009
-- License     :  BSD3
-- Maintainer  :  Marcel Fourné (hecc@bitrot.dyndns.org
--
-- ECC Base algorithms & point formats
--
-----------------------------------------------------------------------------

module Codec.Encryption.ECC.Base (ECP(..),
                                  EC(..),
                                  modinv, 
                                  pmul, 
                                  ison,
                                  genkey,
                                  EPa(..), 
                                  EPp(..), 
                                  EPj(..), 
                                  EPmj(..))
    where 

import Control.Monad.Random
import Data.Bits
import Numeric
import Char

-- |extended euclidean algorithm, recursive variant
eeukl :: Integer -> Integer -> (Integer, Integer, Integer)
eeukl a 0 = (a,1,0)
eeukl a b = let (d,s,t) = eeukl b (a `mod` b)
            in (d,t,s-(div a b)*t)

-- |computing the modular inverse of @a@ `mod` @m@
modinv :: Integer -- ^the number to invert
       -> Integer -- ^the modulus
       -> Integer -- ^the inverted value
modinv a m = let (x,y,_) = eeukl a m
             in if x == 1 
                then mod y m
                else undefined

-- |class of all Elliptic Curves, has the form y^2=x^3+A*x+B mod P, the parameters being A, B and P
data EC = EC (Integer, Integer, Integer)
        deriving (Eq)
instance Show EC where show (EC (a,b,p)) = "y^2=x^3+" ++ show a ++ "*x+" ++ show b ++ " mod " ++ show p

-- |class of all Elliptic Curve Points
class ECP a where
    -- |function returning the appropriate INF in the specific ECP-Format, for generic higher-level-algorithms
    inf :: a
    -- |generic getter, returning the affine x-value
    getx :: a -> EC -> Integer
    -- |generic getters, returning the affine y-value
    gety :: a -> EC -> Integer
    -- |add an elliptic point onto itself, base for padd a a c
    pdouble :: a -> EC -> a
    -- |add 2 elliptic points
    padd :: a -> a -> EC -> a

-- |Elliptic Point Affine coordinates, two parameters x and y
data EPa = EPa (Integer, Integer) 
         | Infa
           deriving (Eq)
instance Show EPa where show (EPa (a,b)) = show (a,b)
                        show Infa = "Null"
instance ECP EPa where 
    inf = Infa
    getx (EPa (x,_)) _ = x
    getx Infa _ = undefined
    gety (EPa (_,y)) _ = y
    gety Infa _ = undefined
    pdouble (EPa (x1,y1)) (EC (alpha,_,p)) = 
        let lambda = ((3*x1^(2::Int)+alpha)*(modinv (2*y1) p)) `mod` p
            x3 = (lambda^(2::Int) - 2*x1) `mod` p
            y3 = (lambda*(x1-x3)-y1) `mod` p
        in EPa (x3,y3)
    pdouble Infa _ = Infa
    padd Infa a _ = a
    padd a Infa _ = a
    padd a@(EPa (x1,y1)) b@(EPa (x2,y2)) c@(EC (_,_,p)) 
        | x1==x2,y1==(-y2) = Infa
        | a==b = pdouble a c
        | otherwise = 
            let lambda = ((y2-y1)*(modinv (x2-x1) p)) `mod` p
                x3 = (lambda^(2::Int) - x1 - x2) `mod` p
                y3 = (lambda*(x1-x3)-y1) `mod` p
            in EPa (x3,y3)

-- |Elliptic Point Projective coordinates, three parameters x, y and z, like affine (x/z,y/z)
data EPp = EPp (Integer, Integer, Integer) 
         | Infp
           deriving (Eq)
instance Show EPp where show (EPp (a,b,c)) = show (a,b,c)
                        show Infp = "Null"
instance ECP EPp where
    inf = Infp
    getx (EPp (x,_,z)) (EC (_,_,p)) = (x * (modinv z p)) `mod` p
    getx Infp _ = undefined
    gety (EPp (_,y,z)) (EC (_,_,p)) = (y * (modinv z p)) `mod` p
    gety Infp _ = undefined
    pdouble (EPp (x1,y1,z1)) (EC (alpha,_,p)) = 
        let a = (alpha*z1^(2::Int)+3*x1^(2::Int)) `mod` p
            b = (y1*z1) `mod` p
            c = (x1*y1*b) `mod` p
            d = (a^(2::Int)-8*c) `mod` p
            x3 = (2*b*d) `mod` p
            y3 = (a*(4*c-d)-8*y1^(2::Int)*b^(2::Int)) `mod` p
            z3 = (8*b^(3::Int)) `mod` p
        in EPp (x3,y3,z3)
    pdouble Infp _ = Infp
    padd Infp a _ = a
    padd a Infp _ = a
    padd p1@(EPp (x1,y1,z1)) p2@(EPp (x2,y2,z2)) curve@(EC (_,_,p)) 
        | x1==x2,y1==(-y2) = Infp
        | p1==p2 = pdouble p1 curve
        | otherwise = 
            let a = (y2*z1 - y1*z2) `mod` p
                b = (x2*z1 - x1*z2) `mod` p
                c = (a^(2::Int)*z1*z2 - b^(3::Int) - 2*b^(2::Int)*x1*z2) `mod` p
                x3 = (b*c) `mod` p
                y3 = (a*(b^(2::Int)*x1*z2-c)-b^(3::Int)*y1*z2) `mod` p
                z3 = (b^(3::Int)*z1*z2) `mod` p
            in EPp (x3,y3,z3)
    
-- |Elliptic Point Jacobian coordinates, three parameter x, y and z, like affine (x/z^2,y/z^3)
data EPj = EPj (Integer, Integer, Integer) 
         | Infj
           deriving (Eq)
instance Show EPj where show (EPj (a,b,c)) = show (a,b,c)
                        show Infj = "Null"
instance ECP EPj where
    inf = Infj
    getx (EPj (x,_,z)) (EC (_,_,p)) = (x * (modinv (z^(2::Int)) p)) `mod` p
    getx Infj _ = undefined
    gety (EPj (_,y,z)) (EC (_,_,p)) = (y * (modinv (z^(3::Int)) p)) `mod` p
    gety Infj _ = undefined
    pdouble (EPj (x1,y1,z1)) (EC (alpha,_,p)) = 
        let a = 4*x1*y1^(2::Int) `mod` p
            b = (3*x1^(2::Int) + alpha*z1^(4::Int)) `mod` p
            x3 = (-2*a + b^(2::Int)) `mod` p
            y3 = (-8*y1^(4::Int) + b*(a-x3)) `mod` p
            z3 = 2*y1*z1 `mod` p
        in EPj (x3,y3,z3)
    pdouble Infj _ = Infj
    padd Infj a _ = a
    padd a Infj _ = a 
    padd p1@(EPj (x1,y1,z1)) p2@(EPj (x2,y2,z2)) curve@(EC (_,_,p)) 
        | x1==x2,y1==(-y2) = Infj
        | p1==p2 = pdouble p1 curve
        | otherwise = 
            let a = (x1*z2^(2::Int)) `mod` p
                b = (x2*z1^(2::Int)) `mod` p
                c = (y1*z2^(3::Int)) `mod` p
                d = (y2*z1^(3::Int)) `mod` p
                e = (b - a) `mod` p
                f = (d - c) `mod` p
                x3 = (-e^(3::Int) - 2*a*e^(2::Int) + f^(2::Int)) `mod` p
                y3 = (-c*e^(3::Int) + f*(a*e^(2::Int) - x3)) `mod` p
                z3 = (z1*z2*e) `mod` p
            in EPj (x3,y3,z3)

-- |Elliptic Point Modified Jacobian coordinates, four parameters x,y,z and A*z^4 (A being the first curve-parameter), like affine coordinates (x/z^2,y/z^3)
data EPmj = EPmj (Integer, Integer, Integer, Integer) 
         | Infmj
           deriving (Eq)
instance Show EPmj where show (EPmj (a,b,c,d)) = show (a,b,c,d)
                         show Infmj = "Null"
instance ECP EPmj where
    inf = Infmj
    getx (EPmj (x,_,z,_)) (EC (_,_,p)) = (x * (modinv (z^(2::Int)) p)) `mod` p
    getx Infmj _ = undefined
    gety (EPmj (_,y,z,_)) (EC (_,_,p)) = (y * (modinv (z^(3::Int)) p)) `mod` p
    gety Infmj _ = undefined
    pdouble (EPmj (x1,y1,z1,z1')) (EC (_,_,p)) = 
        let s = 4*x1*y1^(2::Int) `mod` p
            u = 8*y1^(4::Int) `mod` p
            m = (3*x1^(2::Int) + z1') `mod` p
            t = (-2*s + m^(2::Int)) `mod` p
            x3 = t
            y3 = (m*(s - t) - u) `mod` p
            z3 = 2*y1*z1 `mod` p
            z3' = 2*u*z1' `mod` p
        in EPmj (x3,y3,z3,z3')
    pdouble Infmj _ = Infmj
    padd Infmj a _ = a
    padd a Infmj _ = a 
    padd p1@(EPmj (x1,y1,z1,_)) p2@(EPmj (x2,y2,z2,_)) curve@(EC (alpha,_,p)) 
        | x1==x2,y1==(-y2) = Infmj
        | p1==p2 = pdouble p1 curve
        | otherwise = 
            let u1 = (x1*z2^(2::Int)) `mod` p
                u2 = (x2*z1^(2::Int)) `mod` p
                s1 = (y1*z2^(3::Int)) `mod` p
                s2 = (y2*z1^(3::Int)) `mod` p
                h = (u2 - u1) `mod` p
                r = (s2 - s1) `mod` p
                x3 = (-h^(3::Int) - 2*u1*h^(2::Int) + r^(2::Int)) `mod` p
                y3 = (-s1*h^(3::Int) + r*(u1*h^(2::Int) - x3)) `mod` p
                z3 = (z1*z2*h) `mod` p
                z3' = (alpha*z3^(4::Int)) `mod` p
            in EPmj (x3,y3,z3,z3')

-- |this is a generic handle for Point Multiplication. The implementation may change.
pmul :: (ECP a) => a -- ^the point to multiply
     -> Integer -- ^times to multiply the point
     -> EC -- ^the curve to operate on
     -> a -- ^the result-point
pmul = montgladder
-- pmul = dnadd

-- |double and add for generic ECP
dnadd :: (ECP a) => a -> Integer -> EC -> a
dnadd b k' c@(EC (_,_,p)) = 
        let k = k' `mod` (p - 1)
            ex a i
                | i < 0 = a
                | not (testBit k i) = ex (pdouble a c) (i - 1)
                | otherwise = ex (padd (pdouble a c) b c) (i - 1)
        in ex inf (length (binary k) - 1)

-- montgomery ladder, timing-attack-resistant (except for caches...)
montgladder :: (ECP a) => a -> Integer -> EC -> a
montgladder b k' c@(EC (_,_,p)) =
        let k = k' `mod` (p - 1)
            ex p1 p2 i
                | i < 0 = p1
                | not (testBit k i) = ex (pdouble p1 c) (padd p1 p2 c) (i - 1)
                | otherwise = ex (padd p1 p2 c) (pdouble p2 c) (i - 1)
        in ex b (pdouble b c) ((length (binary k)) - 2)

-- binary representation of an integer
-- taken from http://haskell.org/haskellwiki/Fibonacci_primes_in_parallel
binary :: Integer -> String 
binary = flip (showIntAtBase 2 intToDigit) []

-- |generic verify, if generic ECP is on EC via getx and gety
ison :: (ECP a) => a -- ^ the elliptic curve point which we check
     -> EC -- ^the curve to test on
     -> Bool -- ^is the point on the curve?
ison pt curve@(EC (alpha,beta,p)) = let x = getx pt curve
                                        y = gety pt curve
                                    in (y^(2::Int)) `mod` p == (x^(3::Int)+alpha*x+beta) `mod` p

-- | given a generator and a curve, generate a point randomly
genkey :: (ECP a) => a -- ^a generator (a point on the curve which multiplied gets to be every other point on the curve)
       -> EC -- ^the curve
       -> IO a -- ^the random point which will be the key
genkey a c@(EC (_,_,p)) = do
  n <- evalRandIO $ getRandomR (1,p)
  return $ pmul a n c