# OK to post

############################################################## OLD CODE ##############################################################

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: 

############################################################## END OLD CODE ##############################################################

# PART 1: RandPt() PROCEDURE:
# -----------------------

RandPt:=proc(a,b,p) local x:
	do
		x:=Sqrt(rand(-1000..1000)(),p,10);
	until x <> FAIL;

	return [x,sqrt((x^3+a*x+b) mod p)];
end:

# PART 2: AliceToBobEC() and BobToAliceEC PROCEDURES:
# -----------------------

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

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

# PART 1: Hasse() PROCEDURE:
# -----------------------

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;

	return [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;