# OK to post homework
# Alex Varjabedian, 6-Apr-2024, Homework 21

Help := proc() print(`RandPt(a, b, p)\nAliceToBobEC(a, b, p, g)\nBobToAliceEC(a, b, p, g)\nHasse(a, b, K)`) end:

DoesLie:=proc(a,b,P): P[2]^2-P[1]^3-a*P[1]-b:end:

Sqrt:=proc(n,p,K) local a,i,ra,x:
ra:=rand(1..p-1):

if n&^((p-1)/2) mod p=p-1 then
   RETURN(FAIL):
fi:

for i from 1 to K do
 a:=ra():
if (a^2-n mod p)&^((p-1)/2) mod p=p-1 then
   x:=expand((a+sqrt(a^2-n))^((p+1)/2)):
   
  RETURN(op(1,x) mod p):
  fi:
od:
FAIL:
end:

Mul:=proc(g,A,a,b,p) local B:
if A=1 then
  RETURN(g):
elif A mod 2=0 then
 B:=Mul(g,A/2,a,b,p):
 RETURN(AEp(a,b,B,B,p)):
else
 RETURN(AEp(a,b,Mul(g,A-1,a,b,p),g,p)):
fi:

end:

AEp:=proc(a,b,P,Q,p) local x,y,s,R,xR,yR:
if P=infinity then
  RETURN(Q):
fi:
if Q=infinity then
  RETURN(P):
fi:

if P[1]<>Q[1] then
 s:=(P[2]-Q[2])/(P[1]-Q[1]) mod p:
 xR:=s^2-P[1]-Q[1] mod p:
 yR:=s*(P[1]-xR)-P[2] mod p:
 RETURN([xR,yR]):
 else
  if P[1]=Q[1] and P[2]<>Q[2] then
   RETURN(infinity):
   elif P[1]=Q[1] and P[2]=0 and Q[2]=0 then
    RETURN(infinity):
    elif P=Q then
      s:=(3*P[1]^2+a)/(2*P[2]) mod p:
      xR:=(s^2-2*P[1]) mod p:
      yR:=s*(P[1]-xR)-P[2] mod p:
      RETURN([xR,yR]):
     else
      RETURN(FAIL):
    fi:
  fi:
end: 

############################
# ------------------------ #
# PART 1 - RandPt(a, b, p) #
# ------------------------ #
############################

RandPt := proc(a, b, p) local x:
do:
 x := Sqrt(rand(-1000..1000)(), p, 10):
until x <> FAIL:
[x, sqrt((x^3 + a*x + b) mod p)]:
end:

###############################################################
# ----------------------------------------------------------- #
# PART 2 - AliceToBobEC(a, b, p, g), BobToAliceEC(a, b, p, g) #
# ----------------------------------------------------------- #
###############################################################

print(`PART 2`);

AliceToBobEC := proc(a, b, p, g) Mul(g, rand(1..1000)(), a, b, p) end:
BobToAliceEC := proc(a, b, p, g) Mul(g, rand(1..1000)(), a, b, p) end:

# Testing
for i from 1 to 10 do:
 a := rand(1..100)():
 b := rand(1..100)():
 p := nextprime(rand(1..100)()):
 g := RandPt(a, b, p):
 printf("a: %a, b: %a, p: %a, g: %a, AB(g) = BA(g): %a\n", a, b, p, g, evalb(AliceToBobEC(a, b, p, BobToAliceEC(a, b, p, g)) = BobToAliceEC(a, b, p, AliceToBobEC(a, b, p, g)))):
od:

###########################
# ----------------------- #
# PART 3 - Hasse(a, b, K) #
# ----------------------- #
###########################

print(`PART 3`);

Hasse := proc(a, b, K) local S, N, p, i, min, max, dl, mdl:
N := [[0$100]$5]:
p := 3:
while p <= K do:
 S := {}:
 for i from 0 to p do:
  dl := (i^3 + a*i + b) mod p:
  mdl := (-(i^3 + a*i + b)) mod p:
  if evalb(DoesLie(a, b, [i, dl]) = 0) then S := S union {[i, dl]}: fi:
  if evalb(DoesLie(a, b, [i, mdl]) = 0) then S := S union {[i, mdl]}: fi:
 od:
 N[conv(p)[1]][conv(p)[2]] := nops(S):
 p := nextprime(p):
od:
for i from 3 to K do:
 if N[conv(i)[1]][conv(i)[2]] <> 0 then N[conv(i)[1]][conv(i)[2]] := (N[conv(i)[1]][conv(i)[2]]-i-1)/sqrt(i) fi:
od:
min := 3:
max := 3:
for i from 3 to K do:
 if evalb(evalf(N[conv(i)[1]][conv(i)[2]]) > evalf(N[conv(max)[1]][conv(max)[2]])) and i = nextprime(i-1) then max := i: fi:
 if evalb(evalf(N[conv(i)[1]][conv(i)[2]]) < evalf(N[conv(min)[1]][conv(min)[2]])) and i = nextprime(i-1) then min := i: fi:
od:
[min, max]:
end:

conv := proc(p) local L:
  L := [0, 0]:
  L[1] := ceil(p/100):
  L[2] := p mod 100:
  if L[2] = 0 then L[2] := 100: fi:
  return L:
end:

# Testing
for a from 1 to 10 do:
 for b from 1 to 10 do:
  T[a, b] := Hasse(a, b, 500):
  printf("a: %a b: %a [minp, maxp]: %a\n", a, b, T[a, b]):
 od:
od:
op(T):
printf("T = %a\n", op(T)); # Final table