# Okay to post homework
# Ryan Badi, January 28, 2024, Homework #3

Help() := proc():
  printf("Available processes: EA(m, n) EAnew(m, n) EstimateProbGCD(N, K) MakeRSAkey(D1) MakeRSAkeyG(D1, G)\n"):
end:

with(NumberTheory):

# m > n m = n * q + r (r = m mod n)
# gcd(m, n) = gcd(n, r)  
# gcd(n, 0) = n
# EA(m,n): inputs m > n and outputs the gcd(m, n)
EA := proc(m, n) local q, r:
  if n = 0 then
    return m:
  fi:
  r := m mod n:
  EA(n, r):
end:

# EAnew(m, n): inputs m > n and outputs the gcd(m, n) using iteration rather than recursion
EAnew := proc(m, n) local m1, n1, m1new:
  m1 := m:
  n1 := n:
  while n1 > 0 do:
    m1new := n1:
    n1 := m1 mod n1:
    m1 := m1new:
  od:
end:

EstimateProbGCD := proc(N, K) local successes, i, a, b:
  successes := 0:
  for i from 1 to K do:
    a := rand(1..N)():
	b := rand(1..N)():
	if AreCoprime(a, b) then successes := successes + 1 fi:
  od:
  return successes / K:
end:

# MakeRSAkey: The key [n, e, d]: a -> a^e mod n
# n must be a product of two primes
# inputs D1 and outputs an RSA key
# [n, e, d], where n is a product of two primes
# with D1 digits. Try:
# MakeRSAkey(100);
MakeRSAkey := proc(D1) local n, d, S, m, p, q, e:
  p := nextprime(rand(10^(D1 - 1)..10^D1 - 1)()):
  q := nextprime(rand(10^(D1 - 1)..10^D1 - 1)()):
  if p = q then:
    return FAIL:
  fi:
  n := p * q:
  m := (p - 1) * (q - 1):
  S := rand(m / 2 .. m - 1)():
  for e from S to m - 1 while gcd(e, m) > 1 do od:
  d := e&^(-1) mod m:
  return [n, e, d]:
end:

MakeRSAkeyG := proc(D1, G) local n, d, S, m, e, g, p, P:
  n := 1:
  m := 1:
  for g from 1 to G do:
    p := do nextprime(rand(10^(D1 - 1)..10^D1 - 1)()) until not has(p, P):
    P[g] := p:
	n := n * p:
	m := m * (p - 1):
  od:
  S := rand(m / 2 .. m - 1)():
  for e from S to m - 1 while gcd(e, m) > 1 do od:
  d := e&^(-1) mod m:
  return [n, e, d]:
end:

Help():

printf("Task 2\n EAnew is significantly better than EA because it cannot overwhelm the machine's stack with recursive calls.\n"):
a := rand(10^50000..10^50001)():
b := rand(10^50000..10^50001)():
printf("time(EAnew(a, b)) - time(gcd(a, b))\n%a\n", time(EAnew(a, b)) - time(gcd(a, b))):
printf("Maple's implementation is significantly faster than our own.\n"):

smallest := 1:
largest := 0:
average := 0:
for i from 1 to 100 do:
  trial := EstimateProbGCD(10^10,10^4):
  if trial < smallest then smallest := trial fi:
  if trial > largest then largest := trial fi:
  average := average + trial:
od:
average := average / 100:
printf("\nTask 3\n Smallest Probability %f\nLargest Probability %f\n Average Probability %f\n", smallest, largest, average):

printf("\nTask 5\n"):
M := MakeRSAkeyG(100, 2):
printf("MakeRSAkeyG(100, 2)\nn %a\n e %a\nd %a\n", M[1], M[2], M[3]):