1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators #-}
import Data.Bits
import Data.Proxy
import Data.Ratio (denominator, numerator)
import GHC.TypeLits
import Text.Printf (PrintfArg, printf)
-- FiniteFields
--https://stackoverflow.com/questions/39823408/prime-finite-field-z-pz-in-haskell-with-operator-overloading
newtype FieldElement (p :: Nat) = FieldElement Integer deriving Eq
instance KnownNat n => Num (FieldElement n) where
FieldElement x + FieldElement y = fromInteger $ x + y
FieldElement x * FieldElement y = fromInteger $ x * y
abs x = x
signum _ = 1
negate (FieldElement x) = fromInteger $ negate x
fromInteger a = FieldElement (mod a n) where n = natVal (Proxy :: Proxy n)
instance KnownNat n => Fractional (FieldElement n) where
recip a = a ^ (n - 2) where n = natVal (Proxy :: Proxy n)
fromRational r = fromInteger (numerator r) / fromInteger (denominator r)
instance KnownNat n => Show (FieldElement n) where
show (FieldElement a) | n == (2 ^ 256 - 2 ^ 32 - 977) = printf "0x%064x" a
| otherwise = "FieldElement_" ++ show n ++ " " ++ show a
where n = natVal (Proxy :: Proxy n)
assert :: Bool -> Bool
assert False = error "WRONG"
assert x = x
aa =
let a = FieldElement 2 :: FieldElement 31
b = FieldElement 15
in (a + b == FieldElement 17, a /= b, a - b == FieldElement 18)
bb =
let a = FieldElement 19 :: FieldElement 31
b = FieldElement 24
in a * b
-- Elliptic curve
data ECPoint a
= Infinity
| ECPoint
{ x :: a
, y :: a
, a :: a
, b :: a
}
deriving (Eq)
instance KnownNat n => Show (ECPoint (FieldElement n)) where
show Infinity = "ECPoint(Infinity)"
show p
| n == (2 ^ 256 - 2 ^ 32 - 977) = "S256Point" ++ points
| otherwise = "ECPoint_" ++ show n ++ points ++ params
where
n = natVal (Proxy :: Proxy n)
points = "(" ++ si (x p) ++ ", " ++ si (y p) ++ ")"
params = "a_" ++ si (a p) ++ "|b_" ++ si (b p)
si (FieldElement r) | n == (2 ^ 256 - 2 ^ 32 - 977) = printf "0x%064x" r
| otherwise = show r
validECPoint :: (Eq a, Num a) => ECPoint a -> Bool
validECPoint Infinity = True
validECPoint (ECPoint x y a b) = y ^ 2 == x ^ 3 + a * x + b
add :: (Eq a, Fractional a) => ECPoint a -> ECPoint a -> ECPoint a
add Infinity p = p
add p Infinity = p
add p q
| a p /= a q || b p /= b q = error "point not on same curve"
| x p == x q && y p /= y q = Infinity
| x p /= x q = new_point $ (y q - y p) / (x q - x p)
| x p == x q && y p == 0 = Infinity
| p == q = new_point $ (3 * x p ^ 2 + a p) / (2 * y p)
| otherwise = error "Unexpected case of points"
where
new_point slope =
let new_x = slope ^ 2 - x p - x q
new_y = slope * (x p - new_x) - y p
in ECPoint new_x new_y (a p) (b p)
binex :: (Eq a, Fractional a) => Integer -> ECPoint a -> ECPoint a -> ECPoint a
binex m value result | m == 0 = result
| m .&. 1 == 1 = loop (add result value)
| otherwise = loop result
where loop = binex (m `shiftR` 1) (add value value)
crmul :: (Eq a, Fractional a) => Integer -> ECPoint a -> ECPoint a
crmul m ec = binex m ec Infinity
tre = FieldElement 3 :: FieldElement 31
cc =
let a = ECPoint tre (-7) 5 7
b = ECPoint 18 77 5 7
c = ECPoint (-1) (-1) 5 7
in ( validECPoint a
, validECPoint b
, validECPoint c
, a /= b
, a == a
, add Infinity a
, add a (ECPoint 3 7 5 7)
, add (ECPoint 3 7 5 7) c
, add c c
)
dd =
let prime = 223
a = FieldElement 0 :: FieldElement prime
b = FieldElement 7
x = FieldElement 192
y = FieldElement 105
point = ECPoint x y a b
in point
type S256Field = FieldElement (2 ^ 256- 2^ 32 - 977)
type S256Point = ECPoint S256Field
s256point :: S256Field -> S256Field -> S256Point
s256point x y = ECPoint x y 0 7
li :: S256Field
li = 12
ri= ECPoint 3 7 5 7 :: S256Point
ncons = 0xfffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141
gcons = s256point 0x79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f81798 0x483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8
|