-----------------------------------------------------------------------------
-- |
-- Module      :  Codec.Crypto.ECC.Base
-- Copyright   :  (c) Marcel Fourné 20[09..13]
-- License     :  BSD3
-- Maintainer  :  Marcel Fourné (mail@marcelfourne.de)
-- Stability   :  experimental
-- Portability :  Good
--
-- ECC Base algorithms & point formats
-- 
-----------------------------------------------------------------------------

{-# OPTIONS_GHC -O2 -fllvm -optlo-O3 -feager-blackholing #-}
{-# LANGUAGE GADTs, PatternGuards, FlexibleInstances, BangPatterns #-}

module Codec.Crypto.ECC.Base (EC(..),
                              getBitLength,
                              geta,
                              getb,
                              getp,
                              getr,
                              ECPF(..),
                              getCurve,
                              getx,
                              gety,
                              getz,
                              getaz4,
                              getxA,
                              getyA,
                              padd,
                              pdouble,
                              modinv, 
                              pmul, 
                              ison) 
       where 

import qualified Data.F2 as F2
import qualified Data.Bits as B (testBit)
import qualified Data.List as L (length)
import qualified Numeric as N (showIntAtBase)
import qualified Data.Char as DC (intToDigit)
import qualified Crypto.Types as CT (BitLength)
import Data.Serialize as S (Serialize,put,get)
import Control.Applicative ((<$>),(<*>))

-- |all Elliptic Curves, the parameters being the BitLength L, A, B and P
data EC a where
     -- the Integer Curves, having the form y^2=x^3+A*x+B mod P
     ECi :: CT.BitLength -> Integer -> Integer -> Integer -> Integer -> EC Integer
     -- the Curves on F2, having the form  y^2+x*y=x^3+a*x^2+b mod P; relevant for "ison"
     ECb :: CT.BitLength -> F2.F2 -> F2.F2 -> F2.F2 -> F2.F2 -> EC F2.F2
instance Eq (EC a) where
  (ECi l a b p r) == (ECi l' a' b' p' r') = l==l' && a==a' && b==b' && p==p' && r==r'
  (ECb l a b p r) == (ECb l' a' b' p' r') = l==l' && a==a' && b==b' && p==p' && r==r'
  _ == _ = False
instance Show (EC a) where
  show (ECi l a b p r) = "Curve with length" ++ show l ++", y^2=x^3+" ++ show a ++ "*x+" ++ show b ++ " mod " ++ show p ++ " and group order " ++ show r
  show (ECb l a b p r) = "Curve with length" ++ show l ++", y^2=x^3+" ++ show a ++ "*x+" ++ show b ++ " mod " ++ show p ++ " and group order " ++ show r
-- for now only an EC Integer instance, since F2 is not instance of Serialize; also: a very simple one
instance Serialize (EC Integer) where
  put (ECi l a b p r) = S.put l >> S.put a >> S.put b >> S.put p >> S.put r
  get = (ECi) <$> S.get <*> S.get <*> S.get <*> S.get <*> S.get
instance Serialize (EC F2.F2) where
  put (ECb l a b p r) = S.put l >> S.put a >> S.put b >> S.put p >> S.put r
  get = (ECb) <$> S.get <*> S.get <*> S.get <*> S.get <*> S.get

-- |get bitlength          
getBitLength :: EC a -> Int
getBitLength (ECi l _ _ _ _) = l
getBitLength (ECb l _ _ _ _) = l

-- |get Curve parameter A
geta :: EC a -> a
geta (ECi _ a _ _ _) = a
geta (ECb _ a _ _ _) = a

-- |get Curve parameter B
getb :: EC a -> a
getb (ECi _ _ b _ _) = b
getb (ECb _ _ b _ _) = b

-- |get Curve parameter P
getp :: EC a -> a
getp (ECi _ _ _ p _) = p
getp (ECb _ _ _ p _) = p

-- |get Curve order r
getr :: EC a -> a
getr (ECi _ _ _ _ r) = r
getr (ECb _ _ _ _ r) = r

-- every point has a curve on which it is valid (has to be tested manually), plus possibly some coordinates
-- parametrised by the kind of numbers one which it may be computed
-- point formats may be translated through functions
-- |data of all Elliptic Curve Points
data ECPF a where 
  -- Elliptic Curve Point Affine coordinates, two parameters x and y
  ECPa :: EC Integer -> Integer -> Integer -> ECPF Integer
  -- Elliptic Curve Point Projective coordinates, three parameters x, y and z, like affine (x/z,y/z)
  ECPp ::EC Integer -> Integer -> Integer -> Integer -> ECPF Integer
  -- Elliptic Curve Point Jacobian coordinates, three parameter x, y and z, like affine (x/z^2,y/z^3)
  ECPj :: EC Integer -> Integer -> Integer -> Integer -> ECPF Integer
  -- Elliptic Curve 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)
  ECPmj :: EC Integer -> Integer -> Integer -> Integer -> Integer -> ECPF Integer
  -- Elliptic Curve Point Affine coordinates in F2, two parameters x and y
  ECPaF2 :: EC F2.F2 -> F2.F2 -> F2.F2 -> ECPF F2.F2
  -- Elliptic Curve Point Projective coordinates in F2, three parameters x, y and z, like affine (x/z,y/z)
  ECPpF2 :: EC F2.F2 -> F2.F2 -> F2.F2 -> F2.F2 -> ECPF F2.F2
  -- conserve the elliptic curve, but the point at infinity does not need coordinates
  -- Elliptic Curve Point at Infinity on an Integer Curve
  ECPInfI :: EC Integer -> ECPF Integer
  -- Elliptic Curve Point at Infinity on an F2 Curve
  ECPInfF2 :: EC F2.F2 -> ECPF F2.F2
instance Eq (ECPF a) where
  (ECPa curve x y) == (ECPa curve' x' y') = curve==curve' && x==x' && y==y'
  (ECPp curve x y z) == (ECPp curve' x' y' z') = curve==curve' && x==x' && y==y' && z==z'
  (ECPj curve x y z) == (ECPj curve' x' y' z') = curve==curve' && x==x' && y==y' && z==z'
  (ECPmj curve x y z az4) == (ECPmj curve' x' y' z' az4') = curve==curve' && x==x' && y==y' && z==z' && az4==az4'
  (ECPaF2 curve x y) == (ECPaF2 curve' x' y') = curve==curve' && x==x' && y==y'
  (ECPpF2 curve x y z) == (ECPpF2 curve' x' y' z') = curve==curve' && x==x' && y==y' && z==z'
  (ECPInfI curve) == (ECPInfI curve') = curve==curve'
  (ECPInfF2 curve) == (ECPInfF2 curve') = curve==curve'
  _ == _ = False
instance Show (ECPF a) where
  show (ECPa curve x y) = show (curve,x,y)
  show (ECPp curve x y z) = show (curve,x,y,z)
  show (ECPj curve x y z) = show (curve,x,y,z)
  show (ECPmj curve x y z az4) = show (curve,x,y,z,az4)
  show (ECPaF2 curve x y) = show (curve,x,y) 
  show (ECPpF2 curve x y z) = show (curve,x,y,z)
  show (ECPInfI curve) = show "Point at Infinity on the " ++ show curve
  show (ECPInfF2 curve) = show "Point at Infinity on the " ++ show curve
-- for now only an ECPF Integer instance, since F2 is not instance of Serialize; also: a very simple one
instance Serialize (ECPF Integer) where
  -- not using getxA,getzA for a single put, because "decode . encode = id" and ECPInfI!
  -- the first char is a simple tag
  put pt@(ECPa _ _ _) = S.put "a" >> (S.put $ getCurve pt) >> (S.put $ getx pt) >> (S.put $ gety pt)
  put pt@(ECPp _ _ _ _) = S.put "p" >> (S.put $ getCurve pt) >> (S.put $ getx pt) >> (S.put $ gety pt) >> (S.put $ getz pt)
  put pt@(ECPj _ _ _ _) = S.put "j" >> (S.put $ getCurve pt) >> (S.put $ getx pt) >> (S.put $ gety pt) >> (S.put $ getz pt)
  put pt@(ECPmj _ _ _ _ _) = S.put "j" >> (S.put $ getCurve pt) >> (S.put $ getx pt) >> (S.put $ gety pt) >> (S.put $ getz pt) >> (S.put $ getaz4 pt)
  put pt@(ECPInfI _) = S.put "i" >> (S.put $ getCurve pt)
  -- in the following part a monad is needed, because the tag t implicates the output type and how many get are done
  get = do
    t <- S.get
    case t of
      "a" -> (ECPa) <$> S.get <*> S.get <*> S.get
      "p" -> (ECPp) <$> S.get <*> S.get <*> S.get <*> S.get
      "j" -> (ECPj) <$> S.get <*> S.get <*> S.get <*> S.get
      "m" -> (ECPmj) <$> S.get <*> S.get <*> S.get <*> S.get <*> S.get
      "i" -> (ECPInfI) <$> S.get
      _ -> fail "Wrong format!"
instance Serialize (ECPF F2.F2) where
  put pt@(ECPaF2 _ _ _) = S.put "a2" >> (S.put $ getCurve pt) >> (S.put $ getx pt) >> (S.put $ gety pt)
  put pt@(ECPpF2 _ _ _ _) = S.put "p2" >> (S.put $ getCurve pt) >> (S.put $ getx pt) >> (S.put $ gety pt) >> (S.put $ getz pt)
  put pt@(ECPInfF2 _) = S.put "i2" >> (S.put $ getCurve pt)
  get = do
    t <- S.get
    case t of
      "a2" -> (ECPaF2) <$> S.get <*> S.get <*> S.get
      "p2" -> (ECPpF2) <$> S.get <*> S.get <*> S.get <*> S.get
      "i2" -> (ECPInfF2) <$> S.get
      _ -> fail "Wrong format!"
  
-- |get contents of the curve
getCurve :: ECPF a -> EC a
getCurve (ECPa curve _ _) = curve
getCurve (ECPp curve _ _ _) = curve
getCurve (ECPj curve _ _ _) = curve
getCurve (ECPmj curve _ _ _ _) = curve
getCurve (ECPaF2 curve _ _) = curve
getCurve (ECPpF2 curve _ _ _) = curve
getCurve (ECPInfI c) = c
getCurve (ECPInfF2 c) = c

-- |generic getter, returning the x-value
getx :: ECPF a -> a
getx (ECPa _ x _) = x
getx (ECPp _ x _ _) = x
getx (ECPj _ x _ _) = x
getx (ECPmj _ x _ _ _) = x
getx (ECPaF2 _ x _) = x
getx (ECPpF2 _ x _ _) = x
getx (ECPInfI _) = undefined
getx (ECPInfF2 _) = undefined

-- |generic getter, returning the y-value
gety :: ECPF a -> a
gety (ECPa _ _ y) = y
gety (ECPp _ _ y _) = y
gety (ECPj _ _ y _) = y
gety (ECPmj _ _ y _ _) = y
gety (ECPaF2 _ _ y) = y
gety (ECPpF2 _ _ y _) = y
gety (ECPInfI _) = undefined
gety (ECPInfF2 _) = undefined

-- |generic getter, returning the z-value for points having them
getz :: ECPF a -> a
getz (ECPa _ _ _) = undefined
getz (ECPp _ _ _ z) = z
getz (ECPj _ _ _ z) = z
getz (ECPmj _ _ _ z _) = z
getz (ECPaF2 _ _ _) = undefined
getz (ECPpF2 _ _ _ z) = z
getz (ECPInfI _) = undefined
getz (ECPInfF2 _) = undefined

-- |generic getter, returning the a*z^4-value for points having them
getaz4 :: ECPF a -> a
getaz4 (ECPa _ _ _) = undefined
getaz4 (ECPp _ _ _ _) = undefined
getaz4 (ECPj _ _ _ _) = undefined
getaz4 (ECPmj _ _ _ _ az4) = az4
getaz4 (ECPaF2 _ _ _) = undefined
getaz4 (ECPpF2 _ _ _ _) = undefined
getaz4 (ECPInfI _) = undefined
getaz4 (ECPInfF2 _) = undefined

-- |generic getter, returning the affine x-value
getxA :: ECPF a -> a
getxA pt@(ECPa _ _ _) = getx pt
getxA pt@(ECPp _ _ _ _) = 
  let p = getp $ getCurve pt
      x = getx pt
      z = getz pt 
  in (x * (modinv z p)) `mod` p
getxA pt@(ECPj _ _ _ _) = 
  let p = getp $ getCurve pt
      x = getx pt
      z = getz pt 
  in (x * (modinv (z^(2::Int)) p)) `mod` p
getxA pt@(ECPmj _ _ _ _ _) = 
  let p = getp $ getCurve pt
      x = getx pt
      z = getz pt 
  in (x * (modinv (z^(2::Int)) p)) `mod` p
getxA pt@(ECPaF2 _ _ _) = getx pt
getxA pt@(ECPpF2 _ _ _ _) = 
  let p = getp $ getCurve pt
      x = getx pt
      z = getz pt 
  in (x * (F2.bininv z p)) `F2.mod` p
getxA (ECPInfI _) = undefined
getxA (ECPInfF2 _) = undefined

-- |generic getter, returning the affine y-value
getyA :: ECPF a -> a
getyA pt@(ECPa _ _ _) = gety pt
getyA pt@(ECPp _ _ _ _) = 
  let p = getp $ getCurve pt
      y = gety pt
      z = getz pt 
  in (y * (modinv z p)) `mod` p
getyA pt@(ECPj _ _ _ _) = 
  let p = getp $ getCurve pt
      y = gety pt
      z = getz pt 
  in (y * (modinv (z^(3::Int)) p)) `mod` p
getyA pt@(ECPmj _ _ _ _ _) = 
  let p = getp $ getCurve pt
      y = gety pt
      z = getz pt 
  in (y * (modinv (z^(3::Int)) p)) `mod` p
getyA pt@(ECPaF2 _ _ _) = gety pt
getyA pt@(ECPpF2 _ _ _ _) = 
  let p = getp $ getCurve pt
      y = gety pt
      z = getz pt 
  in (y * (F2.bininv z p)) `F2.mod` p
getyA (ECPInfI _) = undefined
getyA (ECPInfF2 _) = undefined

-- |add an elliptic point onto itself, base for padd a a
pdouble :: (ECPF a) -> (ECPF a)
pdouble pt@(ECPInfI _) = pt
pdouble pt@(ECPInfF2 _) = pt
pdouble pt@(ECPa _ _ _) = 
  let curve = getCurve pt
      alpha = geta curve 
      p = getp curve
      x1 = getx pt
      y1 = gety pt
      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 ECPa curve x3 y3
pdouble pt@(ECPp _ _ _ _) = 
  let curve = getCurve pt
      alpha = geta curve 
      p = getp curve
      x1 = getx pt
      y1 = gety pt
      z1 = getz pt
      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 ECPp curve x3 y3 z3
pdouble pt@(ECPj _ _ _ _) = 
  let curve = getCurve pt
      alpha = geta curve 
      p = getp curve
      x1 = getx pt
      y1 = gety pt
      z1 = getz pt
      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 ECPj curve x3 y3 z3
pdouble pt@(ECPmj _ _ _ _ _) = 
  let curve = getCurve pt
      p = getp curve
      x1 = getx pt
      y1 = gety pt
      z1 = getz pt
      z1' = getaz4 pt
      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 ECPmj curve x3 y3 z3 z3'
pdouble pt@(ECPaF2 _ _ _) = 
  let curve = getCurve pt
      alpha = geta curve 
      p = getp curve
      x1 = getx pt
      y1 = gety pt
      lambda = (x1 + (y1 * (F2.bininv x1 p)))
      x3 = (lambda `F2.pow` (fromInteger 2)) + lambda + alpha `F2.mod` p
      y3 = (lambda * (x1 + x3)) + x3 + y1 `F2.mod` p
  in ECPaF2 curve x3 y3 
pdouble pt@(ECPpF2 _ _ _ _) = 
  let curve = getCurve pt
      alpha = geta curve 
      p = getp curve
      x1 = getx pt
      y1 = gety pt
      z1 = getz pt
      a = (x1 `F2.pow` (fromInteger 2)) `F2.mod` p
      b = (a + (y1 * z1)) `F2.mod` p
      c = (x1 * z1) `F2.mod` p
      d = (c `F2.pow` (fromInteger 2)) `F2.mod` p
      e = ((b `F2.pow` (fromInteger 2)) + (b * c) + (alpha * d)) `F2.mod` p
      x3 = (c * e) `F2.mod` p
      y3 = (((b + c) * e) + ((a `F2.pow` (fromInteger 2)) * c)) `F2.mod` p
      z3 = (c * d) `F2.mod` p
  in ECPpF2 curve x3 y3 z3

-- |"generic" verify, if generic ECP is on EC via getxA and getyA
ison :: ECPF a -> Bool
ison pt@(ECPa _ _ _) = 
  let curve = getCurve pt
      alpha = geta curve
      beta = getb curve
      p = getp curve
      x = getxA pt
      y = getyA pt
  in (y^(2::Int)) `mod` p == (x^(3::Int)+alpha*x+beta) `mod` p
ison pt@(ECPp _ _ _ _) = 
  let curve = getCurve pt
      alpha = geta curve
      beta = getb curve
      p = getp curve
      x = getxA pt
      y = getyA pt
  in (y^(2::Int)) `mod` p == (x^(3::Int)+alpha*x+beta) `mod` p
ison pt@(ECPj _ _ _ _) = 
  let curve = getCurve pt
      alpha = geta curve
      beta = getb curve
      p = getp curve
      x = getxA pt
      y = getyA pt
  in (y^(2::Int)) `mod` p == (x^(3::Int)+alpha*x+beta) `mod` p
ison pt@(ECPmj _ _ _ _ _) = 
  let curve = getCurve pt
      alpha = geta curve
      beta = getb curve
      p = getp curve
      x = getxA pt
      y = getyA pt
  in (y^(2::Int)) `mod` p == (x^(3::Int)+alpha*x+beta) `mod` p
ison pt@(ECPaF2 _ _ _) = 
  let curve = getCurve pt
      alpha = geta curve
      beta = getb curve
      p = getp curve
      x = getxA pt
      y = getyA pt
  in ((y `F2.pow` (fromInteger 2)) + (x * y)) `F2.mod` p == ((x `F2.pow` (fromInteger 3)) + (alpha * (x `F2.pow` (fromInteger 2))) + beta) `F2.mod` p
ison pt@(ECPpF2 _ _ _ _) = 
  let curve = getCurve pt
      alpha = geta curve
      beta = getb curve
      p = getp curve
      x = getxA pt
      y = getyA pt
  in ((y `F2.pow` (fromInteger 2)) + (x * y)) `F2.mod` p == ((x `F2.pow` (fromInteger 3)) + (alpha * (x `F2.pow` (fromInteger 2))) + beta) `F2.mod` p
ison (ECPInfI _) = True
ison (ECPInfF2 _) = True

-- |extended euclidean algorithm, recursive variant
eeukl :: (Integral a ) => a -> a -> (a, a, a)
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 :: (Integral a) => a -- ^the number to invert
       -> a -- ^the modulus
       -> a -- ^the inverted value
modinv a m = let (x,y,_) = eeukl a m
             in if x == 1 
                then mod y m
                else undefined

-- |add 2 elliptic points
padd :: (ECPF a) -> (ECPF a) -> (ECPF a)
padd pt@(ECPInfI _) _ = pt
padd _ pt@(ECPInfI _) = pt
padd pt@(ECPInfF2 _) _ = pt
padd _ pt@(ECPInfF2 _) = pt
padd a@(ECPa _ _ _) b@(ECPa _ _ _) 
        | x1==x2,y1==(-y2),curve==curve' = ECPInfI curve
        | a==b = pdouble a
        | 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 if curve==curve' then ECPa curve x3 y3
               else undefined
  where curve = getCurve a
        p = getp curve
        x1 = getx a
        y1 = gety a
        curve' = getCurve b
        x2 = getx b
        y2 = gety b
padd p1@(ECPp _ _ _ _) p2@(ECPp _ _ _ _)
        | x1==x2,y1==(-y2),curve==curve' = ECPInfI curve
        | p1==p2 = pdouble p1
        | 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 if curve==curve' then ECPp curve x3 y3 z3
               else undefined
  where curve = getCurve p1
        p = getp curve
        x1 = getx p1
        y1 = gety p1
        z1 = getz p1
        curve' = getCurve p2
        x2 = getx p2
        y2 = gety p2                 
        z2 = getz p2
padd p1@(ECPj _ _ _ _) p2@(ECPj _ _ _ _)
        | x1==x2,y1==(-y2),curve==curve' = ECPInfI curve
        | p1==p2 = pdouble p1
        | 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 if curve==curve' then ECPj curve x3 y3 z3
               else undefined
  where curve = getCurve p1
        p = getp curve
        x1 = getx p1
        y1 = gety p1
        z1 = getz p1
        curve' = getCurve p2
        x2 = getx p2
        y2 = gety p2                 
        z2 = getz p2
padd p1@(ECPmj _ _ _ _ _) p2@(ECPmj _ _ _ _ _)
        | x1==x2,y1==(-y2),curve==curve' = ECPInfI curve
        | p1==p2 = pdouble p1
        | 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 if curve==curve' then ECPmj curve x3 y3 z3 z3'
               else undefined
  where curve = getCurve p1
        alpha = geta curve
        p = getp curve
        x1 = getx p1
        y1 = gety p1
        z1 = getz p1
        curve' = getCurve p2
        x2 = getx p2
        y2 = gety p2                 
        z2 = getz p2
padd a@(ECPaF2 _ _ _) b@(ECPaF2 _ _ _) 
        | x1==x2, y1==(x2 + y2), curve==curve' = ECPInfF2 curve
        | a==b = pdouble a
        | otherwise = 
            let lambda = ((y1 + y2) * (F2.bininv (x1 + x2) p)) `F2.mod` p
                x3 = ((lambda `F2.pow` (fromInteger 2))  + lambda +  x1  +  x2 + alpha) `F2.mod` p
                y3 = ((lambda * (x1 + x3)) + x3 + y1) `F2.mod` p
            in if curve==curve' then ECPaF2 curve x3 y3
               else undefined
  where curve = getCurve a
        alpha = geta curve
        p = getp curve
        x1 = getx a
        y1 = gety a
        curve' = getCurve b
        x2 = getx b
        y2 = gety b
padd p1@(ECPpF2 _ _ _ _) p2@(ECPpF2 _ _ _ _)
        | x1==x2,y1==(x2 + y2) = ECPInfF2 curve
        | p1==p2 = pdouble p1
        | otherwise = 
            let a = ((y1 * z2) + (z1 * y2)) `F2.mod` p
                b = ((x1 * z2)  +  (z1 * x2)) `F2.mod` p
                c = (x1 * z1) `F2.mod` p
                d = (c `F2.pow` (fromInteger 2)) `F2.mod` p
                e = ((((a `F2.pow` (fromInteger 2)) + (a * b) + (alpha * c)) * d) + (b * c)) `F2.mod` p
                x3 = (b * e) `F2.mod` p
                y3 = (((c * ((a * x1) + (y1 * b))) * z2) + ((a + b) * e)) `F2.mod` p
                z3 = ((b `F2.pow` (fromInteger 3)) * d) `F2.mod` p
            in if curve==curve' then ECPpF2 curve x3 y3 z3
               else undefined
  where curve = getCurve p1
        alpha = geta curve
        p = getp curve
        x1 = getx p1
        y1 = gety p1
        z1 = getz p1
        curve' = getCurve p2
        x2 = getx p2
        y2 = gety p2                 
        z2 = getz p2
padd _ _ = undefined

-- |this is a generic handle for Point Multiplication. The implementation may change.
pmul :: (ECPF a) -> Integer -> (ECPF a)
pmul = montgladder

-- montgomery ladder, timing-attack-resistant (except for caches...)
montgladder :: (ECPF a) -> Integer -> (ECPF a)
montgladder b@(ECPa _ _ _) k'  = 
  let p = getp $ getCurve b
      k = k' `mod` (p - 1)
      ex !p1 !p2 i
        | i < 0 = p1
        | not (B.testBit k i) = ex (pdouble p1) (padd p1 p2) (i - 1)
        | otherwise = ex (padd p1 p2) (pdouble p2) (i - 1)
  in ex b (pdouble b) ((L.length (binary k)) - 2)
montgladder b@(ECPp _ _ _ _) k'  = 
  let p = getp $ getCurve b
      k = k' `mod` (p - 1)
      ex !p1 !p2 i
        | i < 0 = p1
        | not (B.testBit k i) = ex (pdouble p1) (padd p1 p2) (i - 1)
        | otherwise = ex (padd p1 p2) (pdouble p2) (i - 1)
  in ex b (pdouble b) ((L.length (binary k)) - 2)
montgladder b@(ECPj _ _ _ _) k'  = 
  let p = getp $ getCurve b
      k = k' `mod` (p - 1)
      ex !p1 !p2 i
        | i < 0 = p1
        | not (B.testBit k i) = ex (pdouble p1) (padd p1 p2) (i - 1)
        | otherwise = ex (padd p1 p2) (pdouble p2) (i - 1)
  in ex b (pdouble b) ((L.length (binary k)) - 2)
montgladder b@(ECPmj _ _ _ _ _) k'  = 
  let p = getp $ getCurve b
      k = k' `mod` (p - 1)
      ex !p1 !p2 i
        | i < 0 = p1
        | not (B.testBit k i) = ex (pdouble p1) (padd p1 p2) (i - 1)
        | otherwise = ex (padd p1 p2) (pdouble p2) (i - 1)
  in ex b (pdouble b) ((L.length (binary k)) - 2)
montgladder b@(ECPaF2 _ _ _) k'  = 
  let p = getp $ getCurve b
      k = k' `mod` ((F2.toInteger p) - 1)
      ex !p1 !p2 i
        | i < 0 = p1
        | not (B.testBit k i) = ex (pdouble p1) (padd p1 p2) (i - 1)
        | otherwise = ex (padd p1 p2) (pdouble p2) (i - 1)
  in ex b (pdouble b) ((L.length (binary k)) - 2)
montgladder b@(ECPpF2 _ _ _ _) k'  = 
  let p = getp $ getCurve b
      k = k' `mod` ((F2.toInteger p) - 1)
      ex !p1 !p2 i
        | i < 0 = p1
        | not (B.testBit k i) = ex (pdouble p1) (padd p1 p2) (i - 1)
        | otherwise = ex (padd p1 p2) (pdouble p2) (i - 1)
  in ex b (pdouble b) ((L.length (binary k)) - 2)
montgladder b@(ECPInfI _) _ = b
montgladder b@(ECPInfF2 _) _ = b
           
-- |binary representation of an integer
-- |taken from http://haskell.org/haskellwiki/Fibonacci_primes_in_parallel
binary :: Integer -> String
binary = flip (N.showIntAtBase 2 DC.intToDigit) []