#OK to post homework
#GLORIA LIU, Mar 29 2024, Assignment 19

read `C19.txt`:

#=====================================#
#1. If you haven't done it yet email me your choice of project, and whether you are willing to be the 
#team leader.
#I will do this

#=====================================#
#2. Read and understand Jean-Paul Delahaye's fascinating column ( Pour La Science, April 2024). [en 
#Francais, but you can click on "English"] (here is the .pdf version (in French)
#Done

#=====================================#
#3. Write Maple procedures pn, t(n), &Pi(n), implementing these sequences in Delahaye's article. #Verify that they give what they are supposed to give.

pnhelper:=proc(j)
floor(((1+(j-1)!)/j)-floor((j-1)!/j)):
end:

pn:=proc(n) local j,m:
1+add(floor(floor((n/(1+add(pnhelper(j),j=2..m))))^(1/n)), m=1..(2^n)):
end:

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

Pin:=proc(n) local k:
floor(add(sin(((Pi*GAMMA(k))/(2*k)))^2, k=1..n)):
end:

with(NumberTheory):
#Checking if functions work
#pn becomes slow even for small numbers
#this set should only have true
{seq(evalb(pn(i)=ithprime(i)) ,i=1..5)};

#tn works as expected. Returns a prime number that is n+2 from n if n+2 is prime. Else returns 2.
#these sets should only have true
#{seq(evalb(tn(ithprime(i) - 2)=ithprime(i)), i=1..100)};
#2 more than an even number is never prime so this is just a simple way to test if the function
#returns true for some cases where it's supposed to
{seq(evalb(tn(i*2)=2), i=2..100)};

#Pin matches Maple's PrimeCounting function.
#this set should only have true for i>2
{seq(evalb(Pin(i)=PrimeCounting(i)),i=2..100)};

#=====================================#
#4. By using AliceToBob(P,g) and BobToAlice(P,g) verify that they share the same key (do it several times)
result:={}:
for i from 1 to 10 do
 #Pick Sophie Germain prime with between 5 and 50 digits
 P:=FindSGP(rand(5..50)(),5000):
 g:=FindPrimiSGP(P):
 aKeys:=AliceToBob(P,g):
 bKeys:=BobToAlice(P,g):
 aAnswer:=bKeys[2]&^aKeys[1] mod P;
 bAnswer:=aKeys[2]&^bKeys[1] mod P;
 result:=result union {evalb(aAnswer=bAnswer)}:
od:
#This should only have true
print(result);

#=====================================#
#5. Write a procedure
#DL(x,P,g)
#that inputs an integer x between 1 and P-1 and outputs an integer a between 1 and P-1 such that ga=x 
#(or returns FAIL if none exits). Convince yourself that for large P it is very slow.

#Brute force way, there are faster ways 
DL:=proc(x,P,g) local a:
if not type(x,integer) or x < 1 or 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:

#Yes it is very slow for large P.

#=====================================#
#6. [Optional challenge, 5 dolalrs for the best implementation] Implement (in Maple) the Rabin-Miller (pseudo-) primality test
#MillerRabin(p,k): Miller Rabin primality test
#p is an odd number we are checking for primality
#It will print the witness if there is one
#k is the number of t imes we want to check for a non primality witness for p
#Returns false for composite and true for (probably) prime
MillerRabin:=proc(p,k) local s,d,i,b,f,a,j,b1:
for i from 1 to k do
 a:=rand(2..p-1)():
 if MRHelper(p,a)=false then
  print(a);
  RETURN(false):
 fi:
od:
true:
end:

MRHelper:=proc(p,a) local s,d,i,b,f,j,b1:
d:=p-1:
s:=0:
while d mod 2=0 do
 d:=d/2:
 s:=s+1:
od:
 b:=a&^d mod p:
 if b = 1 or b = p-1 then RETURN(true): fi:
 for j from 1 to s do
  b1:=b&^2 mod p:
  if b1 = -1 mod p then
   #Not all of result is -1 so a is not a witness
   RETURN(true):
   f:=false:
   break:
  fi:
  b:=b1:
 od:
false:
end:

MillerRabin(121, 10); #false
MillerRabin(19, 10); #true
MillerRabin(547, 10); #true
MillerRabin(911, 10); #true
MillerRabin(1000, 10); #false
MillerRabin(129, 10); #false