#OK to post homework
#GLORIA LIU, Apr 6 2024, Assignment 21

read `C21.txt`:

#=====================================#
#1. Read and understand Ewa Syta's notes
#Done

#=====================================#
#2. Read and understand the wikipedia article Cippola's algorithm. Then read the Maple procedure (that I 
#finished after class) Sqrt(n,p,K). Note that it does use the Legendre symbol directly but rather the so-
#called Euler criterion (n is a square mod p if (n^((p-1)/2) mod p=1)
#Done

#=====================================#
#3. Write procedure
#RandPt(a,b,p) that inputs integers a and b and a prime p and outputs a random point on the elliptic curve 
#y^2=x^3+a*x+b mod p

RandPT:=proc(a,b,p) local x,y,n,ra:
ra:=rand(1..p-1):
x:=ra():
n:=x^3+a*x+b mod p:
y:=Sqrt(n,p,10):
while y=FAIL do
 x:=ra():
 n:=x^3+a*x+b mod p:
 y:=Sqrt(n,p,10):
od:
[x,y]:
end:

#=====================================#
#4. Write procedures
#AliceToBobEC(a,b,p,g) and BobToAliceEC(a,b,p,g)
#where after both agree on a,b, (the elliptic curve), on p (the prime) and some random point g on 
#EC(a,b,p) where
#Alice picks a random (large) A (and keeps it to herself) and sends Bob Ag (i.e. g+....+g repeated A times)
#Bob picks a random (large) B (and keeps it to himself) and sends Alice Bg (i.e. g+....+g repeated B times)
#Experiment it with a few random cases and verify that the shared key (AB(g)=BA(g)) is indeed the same.

AliceToBobEC:=proc(a,b,p,g) local A,Ag:
A:=rand(10^3..10^4)():
Ag:=Mul(g,A,a,b,p):
[A,Ag]:
end:

BobToAliceEC:=proc(a,b,p,g) local B,Bg:
B:=rand(10^3..10^4)():
Bg:=Mul(g,B,a,b,p):
[B,Bg]:
end:

result:={}:
for i from 1 to 10 do
 p:=nextprime(rand(10^5..10^6)()):
 a:=rand(1..p-1)():
 b:=rand(1..p-1)():
 g:=RandPT(a,b,p):
 Alice:=AliceToBobEC(a,b,p,g):
 Bob:=BobToAliceEC(a,b,p,g):
 aAnswer:=Mul(Bob[2],Alice[1],a,b,p):
 bAnswer:=Mul(Alice[2],Bob[1],a,b,p):
 result:=result union {evalb(aAnswer=bAnswer)}:
od:
#should be true
print(result);

#=====================================#
#5. Write a procedure
#Hasse(a,b,K)
#that finds the primes for which the number of points in y^2=x^3+a*x+b mod p minus (p+1), divided by 
#sqrt(p), for all primes between 3 and K is maximal and minimal. Tabulate
#Hasse(a,b,500) for 1<=a,b<=10

Hasse:=proc(a,b,K) local p,max,min,maxp,minp,c:
p:=3:
max:=0:
min:=K:
while p<=K do
 c:=evalf(nops(EC(a,b,p))-(p+1))/(sqrt(p)):
 c:=evalf(c):
 if c > max then
  max:=c:
  maxp:=p:
 fi:
 if c < min then
  min:=c:
  minp:=p:
 fi:
 p:=nextprime(p):
od:
[maxp, minp]:
end:

T:=table():
for a from 1 to 10 do
 for b from 1 to 10 do
  T[a,b]:=Hasse(a,b,500):
  print(T[a,b]);
 od:
od:

#=====================================#
#6. Write a procedure
#SqrtF(n,p,K)
#whose input and output are the same as Sqrt(n,P,K) but does it much faster (for large p).
#Skipped