# OK to post homework
# Alex Varjabedian, 30-Mar-2024, Homework 19

Help := proc() print(`p(n)\nt(n)\nprimePi(n)\nVerify(P)\nDL(x, P, g)`) end:

IsPrim:=proc(g,P) local i:
for i from 1 to P-2 do
 if g&^i mod P=1 then
    RETURN(false):
 fi: 
od:
true:
end:

FindPrimi:=proc(P) local g:
for g from 2 to P-2 do
  if IsPrim(g,P) then
    RETURN(g):
   fi:
od:
FAIL: 
end:

AliceToBob:=proc(P,g) local a:
a:=rand(3..P-3)():
[a,g&^a mod P]:
end:

BobToAlice:=proc(P,g) local b:
b:=rand(3..P-3)():
[b,g&^b mod P]:
end:

###################################
# ------------------------------- #
# PART 1 - p(n), t(n), primePi(n) #
# ------------------------------- #
###################################

print(`PART 3`);

p := proc(n) local m, j:
1 + sum(floor(floor((n/(1 + sum(floor((1 + (j-1)!)/j - floor((j-1)!/j)), j = 2..m))))^(1/n)), m = 1..2^n):
end:

t := proc(n) local p:
2 + n*floor(1/(1 + sum(floor((n+2)/p - floor((n+1)/p)), p = 2..n+1))):
end:

primePi := proc(n):
 floor(sum(sin((Pi * GAMMA(k)) / (2 * k))^2, k=1..n))
end:

# Testing
printf("First 6 primes with p(n): %a\n", {seq(p(n), n = 1..6)});
printf("First 8 primes with s(n): %a\n", {seq(`if`(n mod 2 = 1, t(n), NULL), n = 1..20)});
printf("Number of primes less than 20: %a\n", primePi(20));
printf("Number of primes less than 30: %a\n", primePi(30));

#########################
# --------------------- #
# PART 4 - Verification #
# --------------------- #
#########################

print(`PART 4`);

Verify := proc(P) local g, AtoB, BtoA, shared, shared_a, shared_b:
g := FindPrimi(P):
AtoB := AliceToBob(P, g):
BtoA := BobToAlice(P, g):
shared_a := BtoA[2]&^AtoB[1] mod P:
shared_b := AtoB[2]&^BtoA[1] mod P:
shared := g&^(AtoB[1]*BtoA[1]) mod P;
printf("%a %a %a\n", shared_a, shared_b, shared); # all variables should have the same value
end:

for i from 1 to 5 do
 P := nextprime(rand(10^4..10^6)()):
 Verify(P):
od:

########################
# -------------------- #
# PART 5 - DL(x, P, g) #
# -------------------- #
########################

DL := proc(x,P,g) local a:
if not (x >= 1 and x <= P-1) then
 return FAIL:
fi:
for a from 1 to P-1 do
 if g^a=x then
  return a:
 fi:
od:
return FAIL:
end:

# After testing with P = 1000003, I am definitely convinced that this is a slow process for large P (Maple crashed and my computer almost crashed)