#OK to post homework
#Joseph Koutsoutis, 01-28-2024, Assignment 3

with(NumberTheory):
Help:=proc():
print(` EstimateProbGCD(N,K), MakeRSAkeyG(D1,g), EncryptAndSign(M,na,da,nb,eb),`):
print(` DecryptAndVerifySignature(M,S,nb,db,na,ea) `):
end:

#1

#Here are some small examples based on pages 60-90
#of Frank Garvan's awesome Maple booklet

#with(plots):
#p1 := plot([sqrt(x), 3*log(x)], x=0..400, title="The Square Root and log functions", linestyle=dash):
#p2 := textplot([[360,16,"y=3log(x)"],[130,10,"y=sqrt(x)"]]):
#display(p1, p2);
#p1 := 'p1': p2 := 'p2':

#implicitplot3d(y^2 + x^2 - z^2 = 1, x=-2..2, y=-2..2, z=-1..1);

#animate3d([r*cos(t+a), r*sin(t+a), r^2*cos(2*t)+sin(a)], r=0..1, t=0..2*Pi, a=0..2*Pi, 
#frames=10, style=patch, title="The Flying Pizza");

#with(Student:-LinearAlgebra):

#f := (i,j) -> x^(i*j):
#A := Matrix(4,4,f):
#factor(Determinant(A));
#f := 'f': A := 'A':

#A := RandomMatrix(3,3);
#MatrixInverse(A);
#ReducedRowEchelonForm(A);
#A.A;
#JordanForm(A);
#A := 'A':

#2

#EAnew is better than the function EA created in class because EA used recursion while
#EAnew only uses iteration. This means that EA is both more likely to fail (due to
#recursion depth) and is probably a bit slower (due to creating extra stack frames).
 
#Overall, I found that EAnew was still slower than the built-in gcd.
#I skimmed Donald Knuth's TAOCP vol 2 which has a section on computing the gcd,
#and an algorithm that speeds up computation for bigints was presented 
#(a second, faster algorithm was left as an exercise that appears to save time by
#using bit shifts and addition/subtraction instead of division).
#I didn't have time to read the full analysis, but I believe that the algorithm
#presented produces better results by doing division with single precision
#integers whenever possible (the algorithm tries to use only the leading digits, and
#Knuth states that the step using multi-precision operations instead of single-precision 
#operations is an "extremely rare occurence"), 
#and my guess is that the built-in gcd does something similar.

#3

EstimateProbGCD := proc(N,K) local i,count:
 count := 0:
 for i from 1 to K do:
  if AreCoprime(rand(1..N)(), rand(1..N)()) then count += 1 fi:
 od: 

 RETURN(count / K):
end:

#after running EstimateProbGCD(10^10,10^4) a couple times, I found that the outputted probabilities
#were roughly in the .595-.615 range.

#4

#I couldn't solve this problem but here was my attempt:

#Suppose a,b are randomly chosen integers between 1 and N.
#Let p_i be the ith prime, and let A_i be the event that p_i is not a factor of both a and b.
#By the fundamental theorem of arithmetic, a and b have unique prime decompositions,
#so we have that P(gcd(a,b) = 1) = P(A_1,A_2,A_3,...).
#As N goes to infinity, the probability that a prime p divides a specific
#integer approaches 1/p, so by independence, P(A_i) = (1-1/p_i^2).
#Additionally, as N goes to infinity, I think that we have
#P(A_i|A_1,A_2,...,A_{i-1}) = P(A_i) (I don't know how to formalize this, but
#I think this is true since we can calculate P(A_i) based on the mutually exclusive
#events where we label each prime p_1,...,p_{i-1} as being included/excluded as a 
#factor of a and included/excluded as a factor of b. The conditions A_1,...,A_{i-1} 
#only restrict primes from being included in both, but the probability of 
#p_i dividing a (or b) should still be 1/p_i).
#With this, we have that P(gcd(a,b) = 1) = prod(1-1/p_i^2,p_i), which is the same as:
#P(gcd(a,b) = 1) = 1 - sum(1/p_i,pi) + sum(1/(p_i*p_j),p_i < p_j) - ...

#This is where I gave up, but when I was skimming the section of TAOCP on 
#computing the gcd, I found that Knuth gave this problem as an exercise with hints.
#Knuth gives the hint that for two absolutely convergent series:
#(sum(a_k/k^z, k=1,2,...))(sum(b_m/m^z, m=1,2,...)) = sum(sum(a_d b_{n/d}, d|n)/n^z, n=1,2,...)
#and the exercise says to prove:
#(sum(u(k)/k^2, k=1,2,...))(sum(1/m^2, m=1,2,...)) = 1 where u is the Moebius function.
#This equals 1 since a_d = u(d) and b_{n/d} = 1, meaning that when computing sum(a_d b_{n/d}, d|n), 
#we can think of this sum as counting how many different ways we can choose unique primes from the 
#prime factorization of n, but we are subtracting the ways that involve an odd number of primes.
#Hence this sum can be understood as the binomial expansion of (1-1)^f(n) where f(n) is the
#number of distinct primes dividing n.
#Maple is able to compute sum(1/m^2, m=1,2,...) = pi^2/6, so
#P(gcd(a,b) = 1) = 6/pi^2 which matches our results from problem 4.


#5

#MakeRSAkeyG takes as input an integer D1 and an integer g used to compute n.
#This outputs n,e,d where n is the product of g primes with D1 digits, 
#e is the encryption key, and d is the decryption key.
MakeRSAkeyG:=proc(D1,g) local n,d,S,m,p,e,i:
n := 1:
m := 1:
for i from 1 to g do:
 p := nextprime(rand(10^(D1-1)..10^D1-1)()):
 n *= p:
 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:

[n,e,d]:
end:

#I do not know why MakeRSAkey(D1,3) is weaker than MakeRSAkey(D1,2).
#My current understanding is that RSA is secure since searching for
#large prime factors is a difficult problem that needs to be solved in order to decrypt a
#message, and since the minimum length of the primes is the same, I don't know how this problem
#has become easier (or if a different approach has become available).

#6

#EncryptAndSign takes as input a message N where N<min(na,nb),
#integers na,da that correspond to the n,d outputs of MakeRSAkey for the current
#user who is encrypting/signing the message, and
#integers nb,eb that correspond to the n,e outputs of MakeRSAkey for the user who
#will be receiving the message.
#The output is a 2 entry list.
#The first entry is the message M encrypted using eb,nb, and
#the second entry is the message S signed by using da,na.
EncryptAndSign := proc(M,na,da,nb,eb):
 RETURN([M&^eb mod nb, M&^da mod na]):
end:

#DecryptAndVerifySignature takes an encrypted message M and a signature S
#as input along with
#integers nb,db that correspond to the n,d outputs of MakeRSAkey for the current
#user who is decrypting/verifying the message, and
#integers na,ea that correspond to the n,e outputs of MakeRSAkey for the user who
#sent/signed the message.
#The output is a 2 entry list.
#The first entry is the message decrypted by using db,nb, and
#the second entry is an integer obtained by taking the decrypted message and
#subtracting S after encrypting using ea,na.
#If the second entry is 0, the signature matches the message.
DecryptAndVerifySignature := proc(M,S,nb,db,na,ea):
 RETURN([M&^db mod nb, (M&^db mod nb) - (S&^ea mod na)]):
end:

#7 done