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

problem1N1 := 147451228887363586625323456966525905720989842312760509775958662775459536677624741
problem1N2 := 181724486732607374235034401344439931270145141565372874381350646276632766328969281
problem1N3 := 258424126740178352128100370736889906817607518086806632752038758788555704304604649
problem1N4 := 324234657928347051123113232023409710234012389751239847120398471917665655581200339
problem1N5 := 408869971164328247524265450583823930434406844303142816841351879439544818685702841
problem1N6 := 542408184634943257672698834917404611542248228873337459368210624910406937582942097


problem2Na := 1084055561
problem2Nb := 5686741440097
problem2Nc := 1032899106233
problem2Nd := 12038459
problem2Ne := 1626641013131
problem2Nf := 12038459


problem4N  := 1359692821
problem4m  := 117285016
problem4c  := 823845737


problem5N  := 1085444233
problem5m  := 277891194
problem5c  := 878460400


problem6a1 := 2^11 - 1
problem6a2 := 2^13 - 1
problem6a3 := 2^17 - 1
problem6a4 := 2^19 - 1
problem6a5 := 2^23 - 1
problem6a6 := 2^29 - 1


problem8bN := 843156784620274963828079044664499378320177127026840734436833335222593049312927235387489615873



(* ******* Here is code implementing Miller-Rabin, for problems 1, 3, and 6 ******* *)

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 constructing primes and products of primes with a given number of digits, for problem 3 ******* *)

RandomPrime[{500, 999}]

(* This function generates a random prime between 500 and 999 inclusive *)

Ranprimed[d_] := RandomPrime[{Ceiling[10^(d-1)], Floor[10^d - 1]}]

(* This function generates a random d-digit prime *)

Rancomposite2d[d_] := Ranprimed[d/2] * Ranprimed[d/2]

(* This function generates an approximately d-digit product of two primes *)

Rancomposite3d[d_] := Ranprimed[d/3] * Ranprimed[d/3] * Ranprimed[d/3]

(* This function generates an approximately d-digit product of three primes *)




(* ******* Here is code illustrating the Timing function, which measures how long a computation lasts ******* *)

largeval := 2^128 + 1

(* This defines a large integer for later use.  Note that you can apply functions directly to named variables. *)

Timing[PrimeQ[largeval]]

(* This measures how long Mathematica takes to test whether the given large number is prime *)

Timing[FactorInteger[largeval]]

(* This measures how long Mathematica takes to factor the given large integer *)




(* ******* Here is code for the Pollard p-1 factorization algorithm, for problems 2 and 3 ******* *)

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, for problems 2, 3, and 6 ******* *)

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 commands for powers mod m ******** *)

Mod[100,6]

(* This will compute 100 modulo 6, yielding 2 *)



PrimeQ[201]

(* This will test whether 201 is prime, which it is not *)



FactorInteger[143]

(* This will factor the integer 143, which is 11 times 13 *)
(* Mathematica's factorization algorithm knows how to do trial division, Pollard p-1, Pollard rho, quadratic sieving, and elliptic curve factorization, and is typically fairly fast on numbers with 50 digits or fewer *)



PowerMod[2, 10, 99]

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



PowerMod[2, -1, 99]

(* This will compute 2^(-1) modulo 99, which is to say, the multiplicative inverse of 2 modulo 99. *)



EulerPhi[40]

(* This will compute phi(40). *)



ChineseRemainder[{1,2,3},{4,5,9}]

(* This will compute the smallest positive integer solution to x = 1 mod 4, x = 2 mod 5, x = 3 mod 9 *)




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

GCD[180, 250]

(* This computes the GCD of 180 and 250 *)



ExtendedGCD[180, 250]

(* This performs the extended Euclidean algorithm, and computes the GCD = 10 along with integers {x,y} such that 180x + 250y = 10 *)



PrimitiveRoot[25]

(* This calculates the smallest primitive root modulo 25 *)



MultiplicativeOrder[7, 25]

(* This calculates the order of 7 modulo 25 *)



Solve[x^3 == 1, x, Modulus -> 49]

(* This calculates all of the numbers x satisfying the equation x^3 = 1 modulo 49 *)


Reduce[x^3 == 1, Modulus -> 49]

(* This simplifies the equation x^3 = 1 modulo 49, which here will give a list of all the solutions *)