# OK to post homework
# Alex Varjabedian, 25-Jan-2024, Homework 3

Help := proc() print(`EA(m, n)\nEAnew(m, n)\nEstimateProbGCD(N, K)\nLimProb()\nMakeRSAkeyG(D1, g)\nMakeSignature(D1, message)`) end:
with(plots): with(linalg):

###########################
# ----------------------- #
# PART 1 - Maple examples #
# ----------------------- #
###########################

print(`PART 1`);

# 5.9 - Solving differential equations

f := 'f': y := f(x);
dy := diff(y, x);
ddy := diff(%, x);
dsolve(ddy + 5*dy + 6*y = sin(x)*exp(-3*x), y);

# 6.1 - 2-dimensional plotting

plot(sin(x), x = -2*Pi..2*Pi);
plot([cos(t), 0.5*sin(t), t = 0..2*Pi]);
polarplot(cos(5*t), t = 0..2*Pi);

# 6.2 - 3-dimensional plotting

x := 'x': y := 'y':
plot3d(exp(-(x^2 + y^2 - 1)^2), x = -2..2, y = -2..2);
plot3d([sqrt(1 + u^2)*cos(t), sqrt(1 + u^2)*sin(t), u], u = -1..1, t = 0..2*Pi);
spacecurve([cos(t), sin(t), t], t = 0..4*Pi, numpoints = 200);
contourplot(exp(-(x^2 + y^2 - 1)^2), x = -(1.3)..(1.3), y = -(1.3)..(1.3), filled = true, coloring = [blue, red]);

# 6.3 - Animation

animate(1/(1+x*t), x = 0..10, t = 0..1, frames = 100);

# 7.1 - Vectors, arrays, and matrices

v := vector([1, 2, 3]);
a := 'a': b := 'b': c := 'c': d := 'd': e := 'e': f := 'f':
A := matrix(2, 3, [a, b, c, d, e, f]);
print(v);
print(A);
op(A);
eval(A);

# Matrix operations

A := matrix(2, 2, [1, 2, 3, 4]):
B := matrix(2, 2, [-2, 3, -5, 1]):
A + B;
evalm(%);
multiply(A, B);

########################
# -------------------- #
# PART 2 - EAnew(m, n) #
# -------------------- #
########################

print(`PART 2`);

EA:=proc(m,n) local q,r:
if n=0 then
 RETURN(m):
fi:
r:=m mod n:
EA(n,r):
end:

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:

gcd_time := 0: EAnew_time := 0: EA_time := 0:

for i from 1 to 10 do
 a:=rand(10^50000..10^50001)(): b:=rand(10^50000..10^50001)(): gcd_time := gcd_time + time(gcd(a,b)): EAnew_time := EAnew_time + time(EAnew(a,b)): EA_time := EA_time + time(EA(a, b)):
od:

gcd_time := gcd_time; EAnew_time := EAnew_time; EA_time := EA_time;
print(`After 10 trials, we can see that the built-in gcd is the fastest, followed by EAnew, and then EA. EA is the slowest because a function programmed recursively is usually slower than one programmed iteratively, because recursion relies on pushing and popping stack frames that store the information of the current state of the function (time-consuming operations), but iteration does not. I am not quite sure why the built-in gcd is the fastest (and we cannot see the implementation because Maple is not open-source) but, if I had to guess, I would say that the built-in likely uses a more efficient algorithm (please see https://en.wikipedia.org/wiki/Binary_GCD_algorithm#Efficiency).`);

##################################
# ------------------------------ #
# PART 3 - EstimateProbGCD(N, K) #
# ------------------------------ #
##################################

print(`PART 3`);

with(NumberTheory):

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

for i from 1 to 10 do
 if i = 1 then
  printf("[%f, ", EstimateProbGCD(10^10, 10^4));
 elif i = 10 then
  printf("%f]", EstimateProbGCD(10^10, 10^4));
 else
  printf("%f, ", EstimateProbGCD(10^10, 10^4));
 fi:
od:

print(`The answers are very similar.`);

#################################
# ----------------------------- #
# PART 4 - Optional Challenge 1 #
# ----------------------------- #
#################################

print(`PART 4`);

print(`This limit may take a while to compute...`);
s := 0:
for i from 1 to 100 do s := s + EstimateProbGCD(10^1000, 1000) od:
limit_of_probability := evalf(s/100);
print(`It is difficult to say what the exact value of this limit would be...`);

###############################
# --------------------------- #
# PART 5 - MakeRSAkeyG(D1, g) #
# --------------------------- #
###############################

print(`PART 5`);

MakeRSAkeyG := proc(D1, g) local i, primes, n, d, S, m, e:
primes := []:
for i from 1 to g do primes := [op(primes), nextprime(rand(10^(D1-1)..10^D1-1)())] od:
if nops(primes) <> nops(convert(primes, set)) then return FAIL fi:
n := 1: m := 1:
for i from 1 to nops(primes) do n := n*primes[i]: m := m*(primes[i]-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:

key := MakeRSAkeyG(400, 25);
ciphertext := 12345678901234567890&^(key[2]) mod key[1];
message := ciphertext&^(key[3]) mod key[1];

# RSA is fascinating!

print(`The reason why g > 2 is less secure than g = 2 is because a brute force attack would take O(n^(1/g)) time to complete, so a greater value of g decreases the time complexity of a brute force attack. This means that smaller values of g are more secure.`);

#################################
# ----------------------------- #
# PART 6 - Optional Challenge 2 #
# ----------------------------- #
#################################

print(`PART 6`);

# Returns a list containing the signature and generated RSA key
MakeSignature := proc(D1, message) local RSAkey:
RSAkey := MakeRSAkeyG(D1, 2):
return [message&^RSAkey[3] mod RSAkey[1], RSAkey]:
end:

# Signature in action!

# First, Bob generates a signature-key pair for his message "12345678901234567890"
bob_pair := MakeSignature(40, 12345678901234567890);
bob_signature := bob_pair[1];
bob_key := bob_pair[2];

# Now, Alice also needs an RSA key, so she needs to generate one
alice_key := MakeRSAkeyG(40, 2);

# Bob needs to encrypt his signature using Alice's public key for privacy
ciphertext_signature := bob_signature&^alice_key[2] mod alice_key[1];

# After Alice receives the encrypted signature, she can decrypt it using her private key
decrypted_signature := ciphertext_signature&^alice_key[3] mod alice_key[1];

# Finally, Alice can decrypt Bob's signature to obtain his original message by using his public key
message := decrypted_signature&^bob_key[2] mod bob_key[1];