# Shaurya Baranwal, Homework 21 April 7, 2024
# OK to post

#---------------------------------------------
# Important Functions Necessary for Homework |
#---------------------------------------------

#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:

#-----------------------------------------------------------------------------------------------------------------------------------------
# Part 3 - 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:

  do:
    x := Sqrt(rand(-1000..1000)(), p, 10):
  until x <> FAIL:

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

end:

#-----------------------------------------------------------
# Part 4 - AliceToBobEC(a,b,p,g) and BobToAliceEC(a,b,p,g) | 
#-----------------------------------------------------------

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:

TestABg := proc() local i, a, b, p, g:
  for i from 1 to 100 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, Ag - Bg = %a\n", a, b, p, g, AliceToBobEC(a, b, p, g) - BobToAlice(a, b, p, g)):
  od:
end:



#-------------------------------------------------------------------------------------------------------------------------------
# Part 5 - Hasse(a,b,K) 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) |
#-------------------------------------------------------------------------------------------------------------------------------

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

myConvert := 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: