Elliptic curves   are sometimes used in   cryptography   as a way to perform   digital signatures.

The purpose of this task is to implement a simplified (without modular arithmetic) version of the elliptic curve arithmetic which is required by the   elliptic curve DSA   protocol.

In a nutshell, an elliptic curve is a bi-dimensional curve defined by the following relation between the x and y coordinates of any point on the curve:

   




y

2


=

x

3


+
a
x
+
b


{\displaystyle y^{2}=x^{3}+ax+b}


a and b are arbitrary parameters that define the specific curve which is used.

For this particular task, we'll use the following parameters:

    a=0,   b=7 
The most interesting thing about elliptic curves is the fact that it is possible to define a   group   structure on it.

To do so we define an   internal composition   rule with an additive notation +,   such that for any three distinct points P, Q and R on the curve, whenever these points are aligned, we have:

    P + Q + R = 0  
Here   0   (zero)   is the infinity point,   for which the x and y values are not defined.   It's basically the same kind of point which defines the horizon in   projective geometry.

We'll also assume here that this infinity point is unique and defines the   neutral element   of the addition.

This was not the definition of the addition, but only its desired property.   For a more accurate definition, we proceed as such:

Given any three aligned points P, Q and R,   we define the sum   S = P + Q   as the point (possibly the infinity point) such that   S, R   and the infinity point are aligned.

Considering the symmetry of the curve around the x-axis, it's easy to convince oneself that two points S and R can be aligned with the infinity point if and only if S and R are symmetric of one another towards the x-axis   (because in that case there is no other candidate than the infinity point to complete the alignment triplet).

S is thus defined as the symmetric of R towards the x axis.

The task consists in defining the addition which, for any two points of the curve, returns the sum of these two points.   You will pick two random points on the curve, compute their sum and show that the symmetric of the sum is aligned with the two initial points.

You will use the a and b parameters of secp256k1, i.e. respectively zero and seven.

Hint:   You might need to define a "doubling" function, that returns P+P for any given point P.

Extra credit:   define the full elliptic curve arithmetic (still not modular, though) by defining a "multiply" function that returns,

for any point P and integer n,   the point P + P + ... + P     (n times).




 
(* 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 ()
 