# Homework 21 -- Spring 2024
# Pablo Blanco. OK to post.


# 3. For some reason, this implementation has a much higher rate of failure than expected.
# 

RandPt:= proc(a,b,p) local x,y,failc:
 x:=rand(1..p-1)():
 y:=Sqrt(x^3+a*x+b,p,100): # 2^{-100} chance of failure IF the input has a square
 # given the density of quadratic residues, shouldn't x^3+a*x+b be a square mod p with probability around 1/2?
 failc:=1:
 while y = FAIL do:
  x:=rand(0..p-1)():
  y:=Sqrt(n,p,100): # 2^{-100} chance of failure
  failc:=failc+1:
  if failc>=p then:
   RETURN(FAIL):
  fi:
 od:

 
 [x,y]:
end:


# 4. Verified that this works by picking a,b,p and g in EC(a,b,p) and checking the following equality
# Mul(Q[2],P[1],a,b,p) = Mul(P[2],Q[1],a,b,p)
# where
# P:=AliceToBobEC(a,b,p,g);
# Q:=BobToAliceEC(a,b,p,g);

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

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

# 5. I tabulated the results below funktion().

with(NumberTheory):
Hasse:=proc(a,b,K) local ec,bnd,P,i,M,p,argL:
 bnd:=PrimeCounting(500):
 
 P:=[seq(ithprime(i),i=2..bnd)]: # the set of primes to check
 M:= [seq( evalf(funktion(a,b,p)) , p in P)]: # apparently maple just can't sort expressions with radicals
 argL:=sort(M, 'output = permutation'):
 
 return([ P[argL[1]], P[argL[bnd-1]] ]):
end:

# function which Hasse minimizes/maximizes
funktion := proc(a,b,p):
 ((nops(EC(a,b,p)) mod p) - p - 1 )/sqrt(p):
end:


### Table: the  T[x]:= Hasse(a,x,500)    with a=1..10
# T[1] := [[491, 173], [499, 409], [479, 283], [487, 439], [461, 313], [491, 241], [431, 353], [499, 311], [479, 461], [479, 271]]:
# T[2] := [[499, 337], [499, 277], [467, 431], [479, 109], [457, 223], [491, 457], [443, 277], [487, 211], [479, 431], [499, 421]]:
# T[3] := [[487, 193], [499, 409], [491, 463], [499, 367], [487, 307], [461, 433], [467, 397], [487, 293], [467, 431], [491, 467]]:
# T[4] := [[499, 241], [487, 277], [479, 487], [487, 499], [491, 157], [479, 449], [439, 491], [479, 443], [461, 293], [479, 457]]:
# T[5] := [[491, 353], [491, 277], [463, 307], [491, 463], [499, 401], [449, 433], [491, 137], [463, 331], [499, 379], [487, 397]]:
# T[6] := [[499, 179], [487, 379], [479, 277], [491, 449], [467, 379], [491, 487], [499, 433], [491, 167], [487, 479], [461, 487]]:
# T[7] := [[487, 109], [461, 233], [467, 277], [457, 397], [461, 487], [499, 389], [487, 223], [487, 167], [499, 167], [499, 487]]:
# T[8] := [[491, 431], [491, 397], [487, 277], [463, 229], [499, 157], [499, 277], [499, 179], [487, 409], [479, 271], [499, 317]]:
# T[9] := [[457, 211], [499, 431], [487, 431], [499, 487], [487, 389], [491, 251], [457, 191], [479, 409], [467, 307], [499, 317]]:
# T[10]:= [[487, 397], [499, 257], [491, 193], [491, 227], [457, 293], [467, 157], [491, 271], [491, 401], [491, 283], [487, 233]]

# all together in a list:
# T[b][a]:=Hasse(a,b,500): # <----
# T:= [seq(T[i],i=1..10)];
# T:= [[[491, 173], [499, 409], [479, 283], [487, 439], [461, 313], [491, 241], [431, 353], [499, 311], [479, 461], [479, 271]], [[499, 337], [499, 277], [467, 431], [479, 109], [457, 223], [491, 457], [443, 277], [487, 211], [479, 431], [499, 421]], [[487, 193], [499, 409], [491, 463], [499, 367], [487, 307], [461, 433], [467, 397], [487, 293], [467, 431], [491, 467]], [[499, 241], [487, 277], [479, 487], [487, 499], [491, 157], [479, 449], [439, 491], [479, 443], [461, 293], [479, 457]], [[491, 353], [491, 277], [463, 307], [491, 463], [499, 401], [449, 433], [491, 137], [463, 331], [499, 379], [487, 397]], [[499, 179], [487, 379], [479, 277], [491, 449], [467, 379], [491, 487], [499, 433], [491, 167], [487, 479], [461, 487]], [[487, 109], [461, 233], [467, 277], [457, 397], [461, 487], [499, 389], [487, 223], [487, 167], [499, 167], [499, 487]], [[491, 431], [491, 397], [487, 277], [463, 229], [499, 157], [499, 277], [499, 179], [487, 409], [479, 271], [499, 317]], [[457, 211], [499, 431], [487, 431], [499, 487], [487, 389], [491, 251], [457, 191], [479, 409], [467, 307], [499, 317]], [[487, 397], [499, 257], [491, 193], [491, 227], [457, 293], [467, 157], [491, 271], [491, 401], [491, 283], [487, 233]]]


###############################################
############## C21 Code below
#C21.txt, 4/4/24

Help:=proc(): print(`   AE(a,b,P,Q), EC(a,b,p), AEp(a,b,P,Q,p), Mul(g,A,a,b,p), Sqrt(n,p,K) `):end:

#old code

Help20:=proc(): print(`RP(a,b), DoesLie(a,b,P), AddE(a,b,P,Q) `):end:
#Elliptic curve: y^2=x^3+a*x+b 
RP:=proc(a,b) local x: x:=rand(0..200)()/100: [x,sqrt(x^3+a*x+b)]:end:

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

AddE:=proc(a,b,P,Q) local x,y,s,eq,sol,R,xR:
 s:=(P[2]-Q[2])/(P[1]-Q[1]):
  y:=P[2]+(x-P[1])*s:
sol:={solve(y^2-x^3-a*x-b,x)};


if not (nops(sol)=3 and member(P[1],sol) and member(Q[1],sol)) then
  print(`Something bad happened`):
  RETURN(FAIL):
fi:

sol:=sol minus {P[1],Q[1]}:
xR:=sol[1]:
[xR,-sqrt(xR^3+a*xR+b)]:

end: 

###end of C20.txt

AE:=proc(a,b,P,Q) 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]):
 xR:=s^2-P[1]-Q[1]:
 yR:=s*(P[1]-xR)-P[2]:
 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]):
      xR:=(s^2-2*P[1]):
      yR:=s*(P[1]-xR)-P[2]:
      RETURN([xR,yR]):
     else
      RETURN(FAIL):
    fi:
  fi:
end:

 #EC(a,b,p): all the points (x,y) such that y^2=x^3+a*x+b mod p 
EC:=proc(a,b,p) local x,y,S:
S:={}:
 for  x from 0 to p-1 do
  for y from 0 to p-1 do
   if (y^2-x^3-a*x-b) mod p=0 then
     S:=S union {[x,y]}:
   fi:
 od:
od: 
S union {infinity}:
end:

#AEp(a,b,P,Q,p): addiditon in the elliptic ruve y^2=x^3+a*x+b mod p
 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:    

#Mul(g,A,a,b,p): inputs a member,g, of the elliptic curve y^2=x^3+a*x+b mod p
#and outputs g+g+...+g (A times)
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:
 

#Sqrt(n,p,K): the x such that x^2=n mod p
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: