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

problem1c  := 09925133775403023161145765436168132558312424526975951883114357491259354928694566
problem1N  := 25721846750453992857920211151404185773973809170109756788979587593106010382198809
problem1e  := 257
problem12ec:= 23134950268396883542330416160947543385253909800862777028533354726015417063152123
problem1w  := 02280248284034080022224826080806382820080830084208281010482840340428243040380834

problem2aN := 1084055561
problem2bN := 5686741440097
problem2cN := 1032899106233
problem2dN := 12038459
problem2eN := 1626641013131
problem2fN := 12038459

problem3N1 := 147451228887363586625323456966525905720989842312760509775958662775459536677624741
problem3N2 := 181724486732607374235034401344439931270145141565372874381350646276632766328969281
problem3N3 := 258424126740178352128100370736889906817607518086806632752038758788555704304604649
problem3N4 := 324234657928347051123113232023409710234012389751239847120398471917665655581200339
problem3N5 := 408869971164328247524265450583823930434406844303142816841351879439544818685702841
problem3N6 := 542408184634943257672698834917404611542248228873337459368210624910406937582942097

problem5N1 := 2^11 - 1
problem5N2 := 2^13 - 1
problem5N3 := 2^17 - 1
problem5N4 := 2^19 - 1
problem5N5 := 2^23 - 1
problem5N6 := 2^29 - 1

problem7dN := 843156784620274963828079044664499378320177127026840734436833335222593049312927235387489615873



(* ******* Here is code for Miller-Rabin ******* *)

v2[n_] := If[IntegerQ[n/2], 1 + v2[n/2], 0]
MillerRabinList[a_, n_] := Table[Mod[1 + PowerMod[a, (n - 1)/2^(v2[n - 1] - k), n], n] - 1, {k, 0, v2[n - 1]}]

(* This function with the helper function above generates the Miller-Rabin list for n with base a.  It also marks -1 as -1 rather than n-1 mod n, to make it easier to identify when -1 shows up on the list. *)


(* ******* Here is code for factorization algorithms ******* *)

PollardPminus1FactorFind[B_, a_, n_] := 
 Module[{x = a, q = 1, i = 1}, While[(q < 2) && (i < B),
   i = i + 1;
   x = PowerMod[x, i, n];
   q = GCD[x-1, n]];
  q]

(* This is a direct implementation of the Pollard p-1 algorithm starting with x=a that runs until it either finds a divisor of n that is >1 or iterates B times *)

PollardPminus1FactorFind[20, 2, 4913429]

(* This is the Pollard p-1 example from the notes, and should give the nontrivial divisor 2521 *)




PollardPminus1IntegerList[B_, a_, n_] := Table[PowerMod[a, k!, n], {k, 1, B}]
PollardPminus1IntegerGCDList[B_, a_, n_] := GCD[n, PollardPminus1IntegerList[B, a, n] - 1]
PollardPminus1Table[B_, a_, n_] := 
 MatrixForm[{Join[{"i"}, Table[i, {i, 1, B}]], 
   Join[{"a^{i!} mod n"}, PollardPminus1IntegerList[B, a, n]], 
   Join[{"gcd(a^i!-1, n)"}, PollardPminus1IntegerGCDList[B, a, n]]}]

(* This function with the two helper functions above generate a table of the values obtained while running the Pollard p-1 algorithm: the first row is the step number i, the second row is the value of a^(i!) mod n, and the third row is the gcd *)

PollardPminus1Table[8, 2, 4913429]

(* This generates the table of computations for the Pollard p-1 example from the notes *)




PollardRhoFactorFind[B_, a_, n_] := 
 Module[{x = a, y = a, q = 1, i = 0}, While[(q < 2) && (i < B),
   i = i + 1;
   x = Mod[x^2 + 1, n];
   y = Mod[(y^2 + 1)^2 + 1, n];
   q = GCD[n, x - y]];
  q]

(* This is a direct implementation of the Pollard rho algorithm starting with x=a that runs until it either finds a divisor of n that is >1 or iterates B times *)

PollardRhoFactorFind[11, 2, 1242301]

(* This is the Pollard rho example from the notes, and should give the nontrivial divisor 281 *)



PollardRhoXList[B_, a_, n_] := PollardRhoXList[B, a, n] = NestList[Mod[#^2 + 1, n] &, a, B]
PollardRhoYList[B_, a_, n_] := PollardRhoYList[B, a, n] = NestList[Mod[(#^2 + 1)^2 + 1, n] &, a, B]
PollardRhoGCDList[B_, a_, n_] := GCD[PollardRhoXList[B, a, n] - PollardRhoYList[B, a, n], n]
PollardRhoTable[B_, a_, n_] := 
 MatrixForm[{Join[{"i"}, Table[i, {i, 0, B}]], 
   Join[{"x_i mod n"}, PollardRhoXList[B, a, n]], 
   Join[{"y_i mod n)"}, PollardRhoYList[B, a, n]], 
   Join[{"gcd(x_i-y_i, n)"}, PollardRhoGCDList[B, a, n]]}]

(* This generates a table of the values obtained while running the Pollard rho  algorithm: the first row is the step number i, the second row is the value of x_i = (x_(i-1)^2) + 1 mod n, the third row is the value of y_i = ((y_(i-1)^2)+1)^2+1) mod n, and the last row is the gcd *)

PollardRhoTable[11, 2, 1242301]

(* This generates the table of computations for the Pollard rho example from the notes *)