-----------------------------------------------------------------------------
-- |
-- Module      :  Codec.Crypto.ECC.Base
-- Copyright   :  (c) Marcel Fourné 20[09..13]
-- License     :  BSD3
-- Maintainer  :  Marcel Fourné (hecc@bitrot.dyndns.org)
--
-- ECC Base algorithms & point formats
-- 
-----------------------------------------------------------------------------

{-# LANGUAGE GADTs, PatternGuards, FlexibleInstances #-}

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

import qualified Data.F2 as F2
import Data.Bits (testBit)
import Data.List as L (length)
import Numeric (showIntAtBase)
import Data.Char (intToDigit)
import Crypto.Types (BitLength)
import Data.Serialize (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 :: (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 :: (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)) = put l >> put a >> put b >> put p >> put r
  get = (ECi) <$> ((,,,,) <$> get <*> get <*> get <*> get <*> 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 _) = put 'a' >> (put $ getCurve pt) >> (put $ getx pt) >> (put $ gety pt)
  put pt@(ECPp _) = put 'p' >> (put $ getCurve pt) >> (put $ getx pt) >> (put $ gety pt) >> (put $ getz pt)
  put pt@(ECPj _) = put 'j' >> (put $ getCurve pt) >> (put $ getx pt) >> (put $ gety pt) >> (put $ getz pt)
  put pt@(ECPmj _) = put 'j' >> (put $ getCurve pt) >> (put $ getx pt) >> (put $ gety pt) >> (put $ getz pt) >> (put $ getaz4 pt)
  put pt@(ECPInfI _) = put 'i' >> (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 <- get
    case t of
      'a' -> (ECPa) <$> ((,,) <$> get <*> get <*> get)
      'p' -> (ECPp) <$> ((,,,) <$> get <*> get <*> get <*> get)
      'j' -> (ECPj) <$> ((,,,) <$> get <*> get <*> get <*> get)
      'm' -> (ECPmj) <$> ((,,,,) <$> get <*> get <*> get <*> get <*> get)
      'i' -> (ECPInfI) <$> 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.mul` (F2.bininv z p)) `F2.reduceBy` 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.mul` (F2.bininv z p)) `F2.reduceBy` 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 `F2.add` (y1 `F2.mul` (F2.bininv x1 p)))
                            x3 = (lambda `F2.pow` (F2.fromInteger 2)) `F2.add` lambda `F2.add` alpha `F2.reduceBy` p
                            y3 = (lambda `F2.mul` (x1 `F2.add` x3)) `F2.add` x3 `F2.add` y1 `F2.reduceBy` 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` (F2.fromInteger 2)) `F2.reduceBy` p
                            b = (a `F2.add` (y1 `F2.mul` z1)) `F2.reduceBy` p
                            c = (x1 `F2.mul` z1) `F2.reduceBy` p
                            d = (c `F2.pow` (F2.fromInteger 2)) `F2.reduceBy` p
                            e = ((b `F2.pow` (F2.fromInteger 2)) `F2.add` (b `F2.mul` c) `F2.add` (alpha `F2.mul` d)) `F2.reduceBy` p
                            x3 = (c `F2.mul` e) `F2.reduceBy` p
                            y3 = (((b `F2.add` c) `F2.mul` e) `F2.add` ((a `F2.pow` (F2.fromInteger 2)) `F2.mul` c)) `F2.reduceBy` p
                            z3 = (c `F2.mul` d) `F2.reduceBy` 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` (F2.fromInteger 2)) `F2.add` (x `F2.mul` y)) `F2.reduceBy` p == ((x `F2.pow` (F2.fromInteger 3)) `F2.add` (alpha `F2.mul` (x `F2.pow` (F2.fromInteger 2))) `F2.add` beta) `F2.reduceBy` 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` (F2.fromInteger 2)) `F2.add` (x `F2.mul` y)) `F2.reduceBy` p == ((x `F2.pow` (F2.fromInteger 3)) `F2.add` (alpha `F2.mul` (x `F2.pow` (F2.fromInteger 2))) `F2.add` beta) `F2.reduceBy` 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 _) 
        | ((F2.length x1 == F2.length x2) && (x1==x2)), (F2.length y1 == F2.length y2 && F2.length x2 == F2.length y2) && (y1==(x2 `F2.add` y2)), curve==curve' = ECPInfF2 curve
        | (F2.length x1 == F2.length x2) && (F2.length y1 == F2.length y2) && a==b = pdouble a
        | otherwise = 
            let lambda = ((y1 `F2.add` y2) `F2.mul` (F2.bininv (x1 `F2.add` x2) p)) `F2.reduceBy` p
                x3 = ((lambda `F2.pow` (F2.fromInteger 2))  `F2.add` lambda `F2.add`  x1  `F2.add`  x2 `F2.add` alpha) `F2.reduceBy` p
                y3 = ((lambda `F2.mul` (x1 `F2.add` x3)) `F2.add` x3 `F2.add` y1) `F2.reduceBy` 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 _)
        | ((F2.length x1 == F2.length x2) && (x1==x2)),((F2.length y1 == F2.length y2 && F2.length x2 == F2.length y2) && y1==(x2 `F2.add` y2)) = ECPInfF2 curve
        | (F2.length x1 == F2.length x2) && (F2.length y1 == F2.length y2) && p1==p2 = pdouble p1
        | otherwise = 
            let a = ((y1 `F2.mul` z2) `F2.add` (z1 `F2.mul` y2)) `F2.reduceBy` p
                b = ((x1 `F2.mul` z2)  `F2.add`  (z1 `F2.mul` x2)) `F2.reduceBy` p
                c = (x1 `F2.mul` z1) `F2.reduceBy` p
                d = (c `F2.pow` (F2.fromInteger 2)) `F2.reduceBy` p
                e = ((((a `F2.pow` (F2.fromInteger 2)) `F2.add` (a `F2.mul` b) `F2.add` (alpha `F2.mul` c)) `F2.mul` d) `F2.add` (b `F2.mul` c)) `F2.reduceBy` p
                x3 = (b `F2.mul` e) `F2.reduceBy` p
                y3 = (((c `F2.mul` ((a `F2.mul` x1) `F2.add` (y1 `F2.mul` b))) `F2.mul` z2) `F2.add` ((a `F2.add` b) `F2.mul` e)) `F2.reduceBy` p
                z3 = ((b `F2.pow` (F2.fromInteger 3)) `F2.mul` d) `F2.reduceBy` 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 (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 (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 (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 (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 (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 (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 (showIntAtBase 2 intToDigit) []