# Shaurya Baranwal, March 31, 2024, Homework 19
# OK to post homework

Help := proc(): print(``): end:

#-------------------------------------
# Functions necessary to do homework |
#-------------------------------------

#C19.txt, March 28, 2024; Public Key Cryptography; Diffie-Hellman
Help:=proc(): print(`FindRP(d), IsPrimi(g,P), FindPrimi(P), FindSGP(d,K)`):
print(`FindPrimiSGP(P), AliceToBob(P,g), BobToALice(P,g) `):
end:

#FindRD(d): finds a random prime with d digits
FindRP:=proc(d):
nextprime(rand(10^(d-1)..10^d-1)()):
end:

#IsPrim(g,P): Is g primitive mod P? {g,g^2,..., g^(p-2)} are all different
#if g^i mod p=1 doesn't happen until i=p-1
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(P): given a prime P, finds a primitive g by brute force
FindPrimi:=proc(P) local g:
for g from 2 to P-2 do
  if IsPrim(g,P) then
    RETURN(g):
   fi:
od:

FAIL: 
end:

#FindSGP(d,K): finds a random Sophie Germain prime such that (P-1)/2 has d digits
#trying K times or return FAIL
FindSGP:=proc(d,K) local Q,i,ra:
ra:=rand(10^(d-1)..10^d-1):

for i from 1 to K do
 Q:=nextprime(ra()):
 if isprime(2*Q+1) then
    RETURN(2*Q+1):
 fi:
od:
FAIL:
end:

#FindPrimiSGP(P): finds a random primitive element mod P if P is a Sophie Germain prime
FindPrimiSGP:=proc(P) local Q,a:
Q:=(P-1)/2:
if not (isprime(P) and isprime(Q)) then
  RETURN(FAIL):
fi:

a:=rand(2..P-2)():
if a&^Q mod P=P-1 then
   RETURN(a):
 else
  RETURN(P-a):
fi:
end:

#AliceToBob(P,g): alice picks her secret integer a from 2 to P-2 does not
#tell anyone (not even Bob) and sends Bob g^a mod P
AliceToBob:=proc(P,g) local a:
a:=rand(3..P-3)():
[a,g&^a mod P]:
end:


#BobToAice(P,g): Bob picks his secret integer b from 2 to P-2 does not
#tell anyone (not even Alice) and sends Alice g^b mod P
BobToAlice:=proc(P,g) local b:
b:=rand(3..P-3)():
[b,g&^b mod P]:
end:


#---------------------------------------------------------------------------------
# Part 3: p(n), t(n), &Pi(n), implementing these sequences in Delahaye's article |
#---------------------------------------------------------------------------------

p := proc(n) local m, j:

 return 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:

 return 2 + n * trunc(1 / (1 + sum(trunc((n + 2) / p) - trunc((n + 1) / p), p=2..n+1))):

end:

myPi := proc(n) local k:

 return trunc(sum(sin((Pi * GAMMA(k)) / (2 * k))^2, k=1..n))

end:

#---------------------------------------------------------------------------
# Part 4: Verifying AliceToBob(P,g) and BobToAlice(P,g) share the same key |
#---------------------------------------------------------------------------

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);

end:

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

#-----------------------------------------------------------------------------------------------------------------------
# Part 5: DL(x,P,g) - inputs an integer x between 1 and P-1 and outputs an integer a between 1 and P-1 such that g^a=x |
#-----------------------------------------------------------------------------------------------------------------------

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: