(* You should be able to copy and paste this text file directly into Mathematica.  These notes will be interpreted as commented code. *)


(* ********* Here are functions to compute the discriminant and j-invariant of an elliptic curve ******** *)

ellcurvediscriminant[A_,B_] := Simplify[-16(4 A^3 + 27 B^2)]

ellcurvejinvariant[A_,B_] := Simplify[-1728 (4 A^3)/(4 A^3 + 27 B^2)]



(* ********* Here are commands for adding points on elliptic curves ********* *)

(*Performs the addition pt1+pt2 on the curve y^2=x^3+Ax+B.Enter a \
point as \[Infinity] or as {x,y}*)

ellcurveadd::pterr = "The point `1` does not lie on the curve \
y^2=x^3+`2`x+`3`.";(*custom error message for point addition*)
ellcurveadd[{A_, B_}][pt1_, pt2_] := 
 Module[{pt, x1, x2, y1, y2, m, x3, y3}, 
  Catch[If[pt1 == \[Infinity], Throw[pt2]]; 
   If[pt2 == \[Infinity], 
    Throw[pt1]];(*special cases for dealing with infinity as summand*)
   x1 = pt1[[1]]; y1 = pt1[[2]]; x2 = pt2[[1]]; y2 = pt2[[2]];
   If[y1^2 - (x1^3 + A x1 + B) != 0, 
    Message[ellcurveadd::pterr, {x1, y1}, A, B]; 
    Abort[]];(*check whether each point lies on the curve*)
   If[y2^2 - (x2^3 + A x2 + B) != 0, 
    Message[ellcurveadd::pterr, {x2, y2}, A, B]; Abort[]];
   If[(x1 == x2 && y1 != y2) || (x1 == x2 && y1 == y2 && y1 == 0), 
    Throw[\[Infinity]]];(*special cases for dealing with infinity as \
output*)m = 
    If[x1 == x2 && 
      y1 == y2, (3 x1^2 + A)/(2 y1), (y2 - y1)/(x2 - 
        x1)];(*compute the slope of the line*)x3 = m^2 - x1 - x2; 
   y3 = -m (x3 - x1) - y1;
   pt = {x3, y3};(*coordinates of the sum*)pt]]


ellcurveaddmodp::moderr = "Modulus `1` is not prime"; 
(*custom error message for nonprime modulus*)

ellcurveaddmodp[{A_, B_}, p_][pt1_, pt2_] :=
 If[PrimeQ[p] == False, Message[ellcurveaddmodp::moderr, p]; Abort[],
  Module[{pt, x1, x2, y1, y2, m, x3, y3}, 
   Catch[If[pt1 == \[Infinity], Throw[pt2]]; 
    If[pt2 == \[Infinity], 
     Throw[pt1]];(*special cases for dealing with infinity as summand*)
    x1 = Mod[pt1[[1]], p]; y1 = Mod[pt1[[2]], p]; 
    x2 = Mod[pt2[[1]], p]; y2 = Mod[pt2[[2]], p];
    If[Mod[y1^2 - (x1^3 + A x1 + B), p] != 0, 
     Message[ellcurveadd::pterr, {x1, y1}, A, B]; 
     Abort[]];(*check whether each point lies on the curve*)
    If[Mod[y2^2 - (x2^3 + A x2 + B), p] != 0, 
     Message[ellcurveadd::pterr, {x2, y2}, A, B]; Abort[]];
    If[(x1 == x2 && y1 != y2) || (x1 == x2 && y1 == y2 && y1 == 0), 
     Throw[\[Infinity]]];(*special cases for dealing with infinity as \
output*)m = 
     If[x1 == x2 && y1 == y2, 
      Expand[(3 x1^2 + A)/(2 y1), Modulus -> p], 
      Expand[(y2 - y1)/(x2 - x1), 
       Modulus -> p]];(*compute the slope of the line*)
    x3 = Mod[m^2 - x1 - x2, p]; y3 = Mod[-m (x3 - x1) - y1, p];
    pt = {x3, y3};(*coordinates of the sum*)pt]]]

(*Performs the addition pt1+pt2 on the curve y^2=x^3+Ax+B,modulo \
p.Enter a point as \[Infinity] or as {x,y}*)



(* ********* Here is a command for computing the number of points on an elliptic curve mod p ********* *)

ellcurvenumpoints[{A_,B_}, p_] := p + 1 + Sum[JacobiSymbol[x^3 + A x + B, p], {x,0,p-1}]

(* Computes the number of points on y^2 = x^3 + Ax + B modulo p *)



(* ********* Here are some useful other commands ******** *)

JacobiSymbol[11, 3]

(* This computes the Jacobi symbol (11 / 3). *)



Reduce[y^2 == x^3 + 4x + 4, Modulus -> 13]

(* This will find all solutions {x,y} to this equation, which is the same as computing all the points {x,y} on the curve y^2 = x^3 + 4x + 4 modulo 13 *)



PrimeQ[201]

(* Tests whether a given integer is prime *)


FactorInteger[143]

(* This will factor an integer *)
(* Mathematica's factorization algorithm knows how to do trial division, Pollard p-1, Pollard rho, quadratic sieving, and elliptic curve factorization *)



PowerMod[2, 10, 99]

(* This will compute 2^10 modulo 99.  PowerMod is intelligent and will use successive squaring *)
(* Mod will often also recognize if you give it a power and internally use PowerMod, but it is better to use PowerMod directly *)


PowerMod[2, -1, 99]

(* This will compute 2^(-1) modulo 99. *)


PowerMod[4, 1/2, 99]

(* This will compute 4^(1/2) modulo 99, which is to say, a square root of 4 modulo 99 *)



N[Pi, 40]

(* This will compute pi to 40 digits of precision *)