Elliptic curve arithmetic: Difference between revisions
Content added Content deleted
(Updated C++ entry, add subtraction fixed negative integer multiplication) |
(OCaml implementation) |
||
Line 1,270: | Line 1,270: | ||
a + b + d = Zero |
a + b + d = Zero |
||
a * 12345 = (10.759, 35.388)</pre> |
a * 12345 = (10.759, 35.388)</pre> |
||
=={{header|OCaml}}== |
|||
Original version by [http://rosettacode.org/wiki/User:Vanyamil User:Vanyamil]. |
|||
Some float-precision issues but overall works. |
|||
<lang OCaml> |
|||
(* Task : Elliptic_curve_arithmetic *) |
|||
(* |
|||
Using the secp256k1 elliptic curve (a=0, b=7), |
|||
define the addition operation on points on the curve. |
|||
Extra credit: define the full elliptic curve arithmetic |
|||
(still not modular, though) by defining a "multiply" function. |
|||
*) |
|||
(*** Helpers ***) |
|||
type ec_point = Point of float * float | Inf |
|||
type ec_curve = { a : float; b : float } |
|||
(* By default, cube root doesn't work for negative bases *) |
|||
let cube_root : float -> float = |
|||
let third = 1. /. 3. in |
|||
let f x = |
|||
if x > 0. |
|||
then x ** third |
|||
else ~-. (~-. x ** third) |
|||
in |
|||
f |
|||
(* Finds the left-most x on this curve *) |
|||
let ec_minx ({a; b} : ec_curve) : float = |
|||
let factor = ~-. b *. 0.5 in |
|||
let discr = (factor ** 2.) +. (a ** 3. /. 27.) in |
|||
if discr <= 0. |
|||
then failwith "Not a simple curve" |
|||
else |
|||
let root = sqrt discr in |
|||
cube_root (factor +. root) +. cube_root (factor -. root) |
|||
(* Negates the point by negating y coord *) |
|||
let ec_neg : ec_point -> ec_point = function |
|||
| Inf -> Inf |
|||
| Point (x, y) -> Point (x, ~-. y) |
|||
(*** Actual task at hand ***) |
|||
(* Generates a random point in the vicinity of x=0 *) |
|||
let ec_random ({a; b} as c : ec_curve) : ec_point = |
|||
let minx = ec_minx c in |
|||
let x = Random.float (~-. minx *. 2.) +. minx in |
|||
let rhs = x ** 3. +. a *. x +. b in |
|||
Point (x, sqrt rhs) |
|||
(* Verifies that the point is on curve. |
|||
Due to rounding errors, sometimes these calculations aren't perfect. |
|||
*) |
|||
let on_curve ?(debug : bool = false) ({a; b} : ec_curve) : ec_point -> bool = function |
|||
| Inf -> true |
|||
| Point (x, y) -> |
|||
let lhs = y *. y in |
|||
let rhs = x ** 3. +. a *. x +. b in |
|||
let delta = abs_float (lhs -. rhs) in |
|||
( |
|||
if debug then Printf.printf "Delta = %.8f" delta; |
|||
delta < 0.000001 |
|||
) |
|||
(* Doubles a point on the curve (adds a point to itself) *) |
|||
let ec_double ({a; b} as c : ec_curve) : ec_point -> ec_point = function |
|||
| Inf -> Inf |
|||
| Point (x, y) as p -> |
|||
if not (on_curve c p) |
|||
then failwith "Point not on this curve." |
|||
else if y = 0. |
|||
then Inf |
|||
else |
|||
let s = (3. *. x *. x +. a) /. (2. *. y) in |
|||
let x' = s *. s -. 2. *. x in |
|||
let y' = y +. s *. (x' -. x) in |
|||
Point (x', -. y') |
|||
(* Adds any two points on the curve *) |
|||
let ec_add ({a; b} as c : ec_curve) (p : ec_point) (q : ec_point) : ec_point = |
|||
match p, q with |
|||
| Inf, x | x, Inf -> x |
|||
| Point (px, py), Point (qx, qy) -> |
|||
if not (on_curve c p) || not (on_curve c q) |
|||
then failwith "Point not on this curve." |
|||
else if abs_float (px -. qx) < 0.000001 then |
|||
begin |
|||
if abs_float (py +. qy) < 0.000001 |
|||
then Inf |
|||
else |
|||
(* py must equal qy here, otherwise something goes real bad *) |
|||
ec_double c p |> ec_neg |
|||
end |
|||
else |
|||
let s = (py -. qy) /. (px -. qx) in |
|||
let rx = s *. s -. px -. qx in |
|||
let ry = py +. s *. (rx -. px) in |
|||
Point (rx, -. ry) |
|||
(* Extra credit : multiplies a point by a scalar *) |
|||
let ec_mul ({a; b} as c : ec_curve) (p : ec_point) (n : int) : ec_point = |
|||
let rec helper n curPow acc = |
|||
if n = 0 then acc |
|||
else |
|||
let doubled = ec_double c curPow in |
|||
if n mod 2 = 0 |
|||
then helper (n / 2) doubled acc |
|||
else helper (n / 2) doubled (ec_add c acc curPow) |
|||
in |
|||
helper n p Inf |
|||
(*** Output ***) |
|||
let string_of_point : ec_point -> string = function |
|||
| Inf -> "Zero" |
|||
| Point (x, y) -> Printf.sprintf "(%.4f, %.4f)" x y |
|||
let print_output () = |
|||
let c = { a = 0.; b = 7. } in |
|||
let p = ec_random c in |
|||
let q = ec_random c in |
|||
let r = ec_add c p q in |
|||
let t = ec_neg r in |
|||
Printf.printf "p = %s\n" (string_of_point p); |
|||
Printf.printf "q = %s\n" (string_of_point q); |
|||
Printf.printf "r = p + q = %s\n" (string_of_point r); |
|||
Printf.printf "t = -r = %s\n" (string_of_point t); |
|||
Printf.printf "r + t = %s\n" (ec_add c r t |> string_of_point); |
|||
Printf.printf "p + (q + t) = %s\n" (ec_add c q t |> ec_add c p |> string_of_point); |
|||
Printf.printf "p * 12345 = %s\n" (ec_mul c p 12345 |> string_of_point) |
|||
let _ = |
|||
print_output (); |
|||
print_output () |
|||
</lang> |
|||
{{out}} |
|||
<pre> |
|||
p = (-0.0489, 2.6457) |
|||
q = (-0.0236, 2.6457) |
|||
r = p + q = (0.0726, -2.6458) |
|||
t = -r = (0.0726, 2.6458) |
|||
r + t = Zero |
|||
p + (q + t) = Zero |
|||
p * 12345 = (-1.7905, -1.1224) |
|||
p = (1.2715, 3.0092) |
|||
q = (1.6370, 3.3745) |
|||
r = p + q = (-1.9103, 0.1696) |
|||
t = -r = (-1.9103, -0.1696) |
|||
r + t = Zero |
|||
p + (q + t) = Zero |
|||
p * 12345 = (-0.3948, -2.6341) |
|||
</pre> |
|||
=={{header|PARI/GP}}== |
=={{header|PARI/GP}}== |