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


(* ********* Here are the numbers from the problems ********* *)

Nval5c := 57689254024036183888718503966884177523086419886438161001073372632511969905080297;
eval5c := 25250848999164909164249104247600348919876678272890811307070973522692563710935769;





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

(* The output for each of the elliptic curve operations modulo n is a pair {P,d} where P={x,y} is the desired point and d is a "flag" to capture a divisor of n if one was found while attempting to perform the operation *)

(* If the point addition works then the result will have d=1 and give the result {P,1}.  If there is a failure due to a required division by a quantity that is nonzero but not a unit modulo n, the result will be {\[Infinity], d} where d is the gcd of the non-invertible quantity and n *)


ellcurvedoublemodnwithflag[{A_, B_}, n_][pt1_] := 
 Module[{x1, y1, m, x3, y3}, 
  Catch[If[pt1 == \[Infinity], Throw[{\[Infinity], 1}]];
   x1 = pt1[[1]]; y1 = pt1[[2]];
   If[Mod[y1^2 - (x1^3 + A x1 + B), n] != 0, Message[ellcurveadd::pterr, {x1, y1}, A, B]; Abort[]]; (*check whether point lies on the curve*)
   If[1 < GCD[n, 2 y1] < n, Throw[{\[Infinity], GCD[n, 2 y1]}]]; (* if denominator has proper divisor of n, change flag to divisor's value and set point = infinity *)
   m = Mod[(3 x1^2 + A)*PowerMod[2 y1, -1, n], n];
   x3 = Mod[m^2 - 2 x1, n]; y3 = Mod[-m (x3 - x1) - y1, n];
   {{x3, y3}, 1} (*coordinates of the result along with flag = 1 to indicate no divisor found *) ]]

(* This doubles a point pt1 on the curve y^2 = x^3 + Ax + B modulo n if possible and returns {2pt1, 1} if so, and {\[Infinity], d} if not, where d is a proper divisor of n that appeared in a denominator *)


ellcurveaddmodnwithflag[{A_, B_}, n_][pt1_, pt2_] := 
 Module[{pt, x1, x2, y1, y2, m, x3, y3}, 
  Catch[If[pt1 == \[Infinity], Throw[{pt2, 1}]]; 
   If[pt2 == \[Infinity], Throw[{pt1, 1}]];(*special cases for dealing with infinity as summand*)
   x1 = Mod[pt1[[1]], n]; y1 = Mod[pt1[[2]], n]; 
   x2 = Mod[pt2[[1]], n]; y2 = Mod[pt2[[2]], n];
   If[Mod[y1^2 - (x1^3 + A x1 + B), n] != 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), n] != 0, Message[ellcurveadd::pterr, {x2, y2}, A, B]; Abort[]];
   If[(x1 == x2 && y1 == Mod[-y2, n]) || (x1 == x2 && y1 == y2 && y1 == 0), Throw[\[Infinity]]];(*special cases for dealing with infinity as output*)
   If[1 < GCD[n, 2 y1] < n, Throw[{\[Infinity], GCD[n, 2 y1]}]]; 
   If[1 < GCD[n, x2 - x1] < n, Throw[{\[Infinity], GCD[n, x2 - x1]}]]; 
   If[1 < GCD[n, y2 - y1] < n, Throw[{\[Infinity], GCD[n, y2 - y1]}]]; (*  if denominator has proper divisor of n, change flag to divisor's value and set point = infinity *)
   m = If[x1 == x2 && y1 == y2, 
     Mod[(3 x1^2 + A)*PowerMod[2 y1, -1, n], n], 
     Mod[(y2 - y1)*PowerMod[(x2 - x1), -1, n], 
      n]];(*compute the slope of the line*) 
   x3 = Mod[m^2 - x1 - x2, n]; y3 = Mod[-m (x3 - x1) - y1, n];
   {{x3, y3}, 1} (*coordinates of the result along with flag = 1 to indicate no divisor found *)]]

(* This adds pt1+pt2 on the curve y^2 = x^3 + Ax + B modulo n if possible and returns {pt1+pt2, 1} if so, and {\[Infinity], d} if not, where d is a proper divisor of n that appeared in a denominator *)


ellcurvemultbykmodnwithflag[{A_, B_}, n_][pt1_, k_] := 
 Module[{ctr = k, ptbwithflag = {\[Infinity], 1}, 
   ptcwithflag = {pt1, 1}, ptbfactored = 0, ptcfactored = 0}, 
  Catch[If[pt1 == \[Infinity], Throw[{\[Infinity],1}]];(*quick exit if the point is infinity*)
   If[Mod[pt1[[2]]^2 - (pt1[[1]]^3 + A pt1[[1]] + B), n] != 0, 
    Message[ellcurveadd::pterr, pt1, A, B]; Abort[]];(*check whether the point lies on the curve*)
   If[k < 0, ctr = -k; ptcwithflag = {{pt1[[1]], -pt1[[2]]}, 1}];(*deal with when k is negative*) 
   While[ctr > 0, If[EvenQ[ctr], ctr = ctr/2;
     ptcwithflag = ellcurvedoublemodnwithflag[{A, B}, n][ptcwithflag[[1]]],
     ctr = ctr - 1; ptbwithflag = ellcurveaddmodnwithflag[{A, B}, n][ptbwithflag[[1]], ptcwithflag[[1]]]];
    If[ptcwithflag[[2]] > 1, ptcfactored = 1; Break[]];
    If[ptbwithflag[[2]] > 1, ptbfactored = 1; Break[]]];
   (*computes k.pt1 recursively with successive squaring and outputs result, including flag *)
   If[ptcfactored == 1, ptcwithflag, ptbwithflag](* if ptc found a factor, output it, otherwise output ptb *)]]

(* This computes k*pt1 on the curve y^2 = x^3 + Ax + B modulo n if possible and returns {k*pt1, 1} if so, and {\[Infinity], d} if not, where d is a proper divisor of n that appeared in a denominator *)





(* ****** Here are commands for doing elliptic curve factorization ****** *)


ellcurvefactorptlist[n_, bound_][pt_, A_] := 
 Module[{ptlist = {pt}, currptwithflag = {pt, 1}, ctr = 1, B},
  B := Mod[pt[[2]]^2 - pt[[1]]^3 - A*pt[[1]], n];
  While[ctr < bound,
   If[currptwithflag[[2]] > 1, Break[]]; (* if a factor is already found, stop computing *)
   ctr = ctr + 1;
   currptwithflag = 
    ellcurvemultbykmodnwithflag[{A, B}, n][currptwithflag[[1]], ctr]; (* compute the next point *)
   AppendTo[ptlist, currptwithflag[[1]]]]; (* add the next point to the list of points *)
  ptlist]  (* return the list of points *)

(* This computes the points (1!)*pt, (2!)*pt, ... , (bound!)*pt on the curve y^2 = x^3 + Ax + B modulo n where B is defined so that pt lies on the curve.  If the computation cannot be completed due to a required division during the computation, it will return the list of points until that illegal division occurred. *)


ellcurvefactorfind[n_, bound_][pt_, A_] := 
 Module[{currptwithflag = {pt, 1}, ctr = 1, B},
  B := Mod[pt[[2]]^2 - pt[[1]]^3 - A*pt[[1]], n];
  While[ctr < bound,
   If[currptwithflag[[2]] > 1, Break[]]; (* if a factor is already found, stop computing *)
   ctr = ctr + 1;
   currptwithflag = 
    ellcurvemultbykmodnwithflag[{A, B}, n][currptwithflag[[1]], ctr]];
  {StringJoin["divisor = ", ToString[currptwithflag[[2]]]], 
   StringJoin["iterations = ", ToString[ctr]]}]

(* This uses Lenstra's elliptic curve factorization algorithm to try to factor n using at most bound iterations, by computing (1!)*pt, (2!)*pt, ... , (bound!)*pt on the curve y^2 = x^3 + Ax + B modulo n where B is defined so that pt lies on the curve, and checking whether a divisor of n shows up in a denominator.  It will print the factor it finds along with the number of steps it needed. *)

ellcurvefactorfindBvalue[n_, bound_][pt_, A_] := pt[[2]]^2 - pt[[1]]^3 - A*pt[[1]]

(* This will give the value of B used by the algorithm above, to make it easy to write down the curve as needed *)


ellcurvefactorptlist[170999, 20][{1, 4}, 4]

(* This will return a list of 10 points, the last being infinity, corresponding to the values (1!)*{1,4}, (2!)*{1,4}, ... , (10!)*{1,4} on the curve y^2 = x^3 + 4x + 11 modulo 170999.  There are only 10 values because the computation required a division by a nonzero nonunit modulo n, indicating a factorization was found *)

ellcurvefactorfind[170999, 100][{1, 4}, 4]

(* This will return the divisor 557 found by computing the list above, and indicate that 10 steps were needed. *)



(*  ****** Here is some code for problem 5c ****** *)

CFConvergent[alpha_, n_] := FromContinuedFraction[ContinuedFraction[alpha, n]]

(* computes the nth continued fraction convergent to alpha *)


CFConvergentList[alpha_, max_] := Table[CFConvergent[alpha, n], {n, 1, max}]

(* computes the first n continued fraction convergents to alpha *)


kedcalculate[eeeval_, list_] := Cases[(Denominator[#]*eeeval - 1)/Numerator[#] & @ list, _Integer]

(* computes the expression (k*e-1)/d on a list of rational numbers of the form d/k and keeps only the integer-valued ones *)




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

JacobiSymbol[19,7]

(* This will compute the Jacobi symbol 19 on 7 *)



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


N[Pi, 40]

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