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

problem1aN := 1084055561
problem1bN := 5686741440097
problem1cN := 1032899106233
problem1dN := 12038459
problem1eN := 1626641013131
problem1fN := 12038459

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




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

MillerRabinList[2, 561]

(* This is the first Miller-Rabin example from the notes: it should return the list {263, 166, 67, 1, 1}. *)




(* ******* Here is code for the Pollard p-1 factorization algorithm ******* *)

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[Transpose[{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 *)



(* ******* Here is code for the Pollard rho factorization algorithm ******* *)

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[Transpose[{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 *)




(* *******  Here are some useful functions for computing quotients, remainders, and GCDs in Z[i] ******* *)

(* Note that you ARE expected to be able to do these calculations by hand yourself: these functions are tabulated here in order to make it easier for you to check your calculations as you are learning how to do them *)


QuotientRemainder[7+3I, 2+2I]

(* This computes the quotient and remainder when dividing 7+3i by 2+2i in Z[i]. The QuotientRemainder[] command automatically recognizes Gaussian integers when you use them.  Note that the imaginary unit i must be capitalized as "I" for Mathematica to recognize it as the square root of -1. *)

GCD[7+3I, 2+2I]

(* This computes the GCD of 7+3i and 2+2i in Z[i].  The GCD[] command automatically recognizes Gaussian integers when you use them.  Mathematica will rescale the GCD by unit factors to make it unique: specifically, the GCD will always have positive real part and nonnegative imaginary part. *)

ExtendedGCD[7+3I, 2+2I]

(* This computes the GCD d of a,b along with Gaussian integers {x,y} such that d = xa + yb.  The ExtendedGCD[] command automatically recognizes Gaussian integers when you use them.  Mathematica will rescale the GCD by unit factors to make it unique: specifically, the GCD will always have positive real part and nonnegative imaginary part, and the coefficients x and y are accordingly rescaled also. *)




(* *******  Here are some useful functions for computing quotients, remainders, and GCDs in F[x] ******* *)

(* Note that you ARE expected to be able to do these calculations by hand yourself: these functions are tabulated here in order to make it easier for you to check your calculations as you are learning how to do them *)



PolynomialQuotientRemainder[x^3 + 2x + 7, x^2 + x + 1, x]

(* This computes the quotient and remainder when dividing x^3 + 2x + 7 by x^2 + x + 1, where the variable is x.  Note that the third argument x is needed to tell Mathematica what variable is being used.  By default, Mathematica will return quotients and remainders with complex coefficients (though of course here, the coefficients are all rational numbers). *)

PolynomialQuotientRemainder[x^3 + 2x + 7, x^2 + x + 1, x, Modulus -> 3]

(* This computes the quotient and remainder when dividing x^3 + 2x + 7 by x^2 + x + 1, where the variable is x, and the coefficients are considered modulo 3 -- in other words, it performs the computations in the polynomial ring F_3[x].  You can set the modulus to be any prime number. *)



PolynomialGCD[x^3 + 2x + 7, x^2 + x + 1]

(* This computes the polymomial GCD of x^3 + 2x + 7 and x^2 + x + 1.  Note that you do *not* pass PolynomialGCD[] the variable name x.  Mathematica automatically rescales the GCD to have leading coefficient 1. *)

PolynomialGCD[x^3 + 2x + 7, x^2 + x + 1, Modulus -> 13]

(* This computes the polymomial GCD of x^3 + 2x + 7 and x^2 + x + 1 with coefficients modulo 13.  Note that you do *not* pass PolynomialGCD[] the variable name x.  Mathematica automatically rescales the GCD to have leading coefficient 1. *)

PolynomialExtendedGCD[x^3 + 2x + 7, x^2 + x + 1, x]

(* This computes the polymomial GCD of x^3 + 2x + 7 by x^2 + x + 1 along with polynomials s and t such that sa + tb = d.  Note that you *do* pass PolynomialExtendedGCD[] the variable name x. *)