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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
|
-- |
-- Module: Math.NumberTheory.Curves.Montgomery
-- Copyright: (c) 2017 Andrew Lelechenko
-- Licence: MIT
-- Maintainer: Andrew Lelechenko <andrew.lelechenko@gmail.com>
--
-- Arithmetic on Montgomery elliptic curves.
-- This is an internal module, exposed only for purposes of testing.
--
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -fno-warn-type-defaults #-}
{-# OPTIONS_HADDOCK hide #-}
module Math.NumberTheory.Curves.Montgomery
( Point
, pointX
, pointZ
, pointN
, pointA24
, SomePoint(..)
, newPoint
, add
, double
, multiply
) where
import Data.Proxy
import GHC.Exts
import GHC.Integer.Logarithms
import GHC.TypeNats (KnownNat, SomeNat(..), Nat, natVal, someNatVal)
import Math.NumberTheory.Utils (recipMod)
-- | We use the Montgomery form of elliptic curve:
-- b Y² = X³ + a X² + X (mod n).
-- See Eq. (10.3.1.1) at p. 260 of <http://www.ams.org/journals/mcom/1987-48-177/S0025-5718-1987-0866113-7/S0025-5718-1987-0866113-7.pdf Speeding the Pollard and Elliptic Curve Methods of Factorization> by P. L. Montgomery.
--
-- Switching to projective space by substitutions Y = y \/ z, X = x \/ z,
-- we get b y² z = x³ + a x² z + x z² (mod n).
-- The point on projective elliptic curve is characterized by three coordinates,
-- but it appears that only x- and z-components matter for computations.
-- By the same reason there is no need to store coefficient b.
--
-- That said, the chosen curve is represented by a24 = (a + 2) \/ 4
-- and modulo n at type level, making points on different curves
-- incompatible.
data Point (a24 :: Nat) (n :: Nat) = Point
{ pointX :: !Integer -- ^ Extract x-coordinate.
, pointZ :: !Integer -- ^ Extract z-coordinate.
}
-- | Extract (a + 2) \/ 4, where a is a coefficient in curve's equation.
pointA24 :: forall a24 n. KnownNat a24 => Point a24 n -> Integer
pointA24 _ = toInteger $ natVal (Proxy :: Proxy a24)
-- | Extract modulo of the curve.
pointN :: forall a24 n. KnownNat n => Point a24 n -> Integer
pointN _ = toInteger $ natVal (Proxy :: Proxy n)
-- | In projective space 'Point's are equal, if they are both at infinity
-- or if respective ratios 'pointX' \/ 'pointZ' are equal.
instance KnownNat n => Eq (Point a24 n) where
Point _ 0 == Point _ 0 = True
Point _ 0 == _ = False
_ == Point _ 0 = False
p@(Point x1 z1) == Point x2 z2 = let n = pointN p in (x1 * z2 - x2 * z1) `rem` n == 0
-- | For debugging.
instance (KnownNat a24, KnownNat n) => Show (Point a24 n) where
show p = "(" ++ show (pointX p) ++ ", " ++ show (pointZ p) ++ ") (a24 "
++ show (pointA24 p) ++ ", mod "
++ show (pointN p) ++ ")"
-- | Point on unknown curve.
data SomePoint where
SomePoint :: (KnownNat a24, KnownNat n) => Point a24 n -> SomePoint
instance Show SomePoint where
show (SomePoint p) = show p
-- | 'newPoint' @s@ @n@ creates a point on an elliptic curve modulo @n@, uniquely determined by seed @s@.
-- Some choices of @s@ and @n@ produce ill-parametrized curves, which is reflected by return value 'Nothing'.
--
-- We choose a curve by Suyama's parametrization. See Eq. (3)-(4) at p. 4
-- of <http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdf Implementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware>
-- by K. Gaj, S. Kwon et al.
newPoint :: Integer -> Integer -> Maybe SomePoint
newPoint s n = do
a24denRecip <- recipMod a24den n
a24 <- case a24num * a24denRecip `rem` n of
-- (a+2)/4 = 0 corresponds to singular curve with A = -2
0 -> Nothing
-- (a+2)/4 = 1 corresponds to singular curve with A = 2
1 -> Nothing
t -> Just t
SomeNat (_ :: Proxy a24Ty) <- if a24 < 0
then Nothing
else Just $ someNatVal $ fromInteger a24
SomeNat (_ :: Proxy nTy) <- if n < 0
then Nothing
else Just $ someNatVal $ fromInteger n
return $ SomePoint (Point x z :: Point a24Ty nTy)
where
u = s * s `rem` n - 5
v = 4 * s
d = v - u
x = u * u * u `mod` n
z = v * v * v `mod` n
a24num = d * d * d * (3 * u + v) `mod` n
a24den = 16 * x * v `rem` n
-- | If @p0@ + @p1@ = @p2@, then 'add' @p0@ @p1@ @p2@ equals to @p1@ + @p2@.
-- It is also required that z-coordinates of @p0@, @p1@ and @p2@ are coprime with modulo
-- of elliptic curve; and x-coordinate of @p0@ is non-zero.
-- If preconditions do not hold, return value is undefined.
--
-- Remarkably such addition does not require 'KnownNat' @a24@ constraint.
--
-- Computations follow Algorithm 3 at p. 4
-- of <http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdf Implementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware>
-- by K. Gaj, S. Kwon et al.
add :: KnownNat n => Point a24 n -> Point a24 n -> Point a24 n -> Point a24 n
add p0@(Point x0 z0) (Point x1 z1) (Point x2 z2) = Point x3 z3
where
n = pointN p0
a = (x1 - z1) * (x2 + z2) `rem` n
b = (x1 + z1) * (x2 - z2) `rem` n
apb = a + b
amb = a - b
c = apb * apb `rem` n
d = amb * amb `rem` n
x3 = c * z0 `mod` n
z3 = d * x0 `mod` n
-- | Multiply by 2.
--
-- Computations follow Algorithm 3 at p. 4
-- of <http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdf Implementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware>
-- by K. Gaj, S. Kwon et al.
double :: (KnownNat a24, KnownNat n) => Point a24 n -> Point a24 n
double p@(Point x z) = Point x' z'
where
n = pointN p
a24 = pointA24 p
r = x + z
s = x - z
u = r * r `rem` n
v = s * s `rem` n
t = u - v
x' = u * v `mod` n
z' = (v + a24 * t `rem` n) * t `mod` n
-- | Multiply by given number, using binary algorithm.
multiply :: (KnownNat a24, KnownNat n) => Word -> Point a24 n -> Point a24 n
multiply 0 _ = Point 0 0
multiply 1 p = p
multiply (W# w##) p =
case wordLog2# w## of
l# -> go (l# -# 1#) p (double p)
where
go 0# !p0 !p1 = case w## `and#` 1## of
0## -> double p0
_ -> add p p0 p1
go i# p0 p1 = case uncheckedShiftRL# w## i# `and#` 1## of
0## -> go (i# -# 1#) (double p0) (add p p0 p1)
_ -> go (i# -# 1#) (add p p0 p1) (double p1)
|
