# Okay to post homework
# Ryan Badi, Homework 21,  April 7, 2024

DoesLie := proc(a, b, P): P[2]^2 - P[1]^3 - a * P[1] - b: 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:

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:

# Returns a random point [x, sqrt((+-(x^3 + a * x + b)) mod p)] where x is any random integer from the start of the domain to 100
# The upper bound 100 was chosen arbitrarily to reduce computation speeds. Any integer in the domain will suffice as an upper bound.
RandPt := proc(a, b, p) local k:
  k := trunc(rand(evalf(lhs(solve(evalc(Im(sqrt(x^3 + a*x + b))), {x})[1][1]))..100)()):
  return [k, sqrt([k^3 + a * k + b, -k^3 - a * k + b][rand(1..2)()] mod p)]:
end:

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, evalb(ABg = BAg) = %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:
end:

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[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 evalb(evalf(N[myConvert(i)[1]][myConvert(i)[2]]) > evalf(N[myConvert(max)[1]][myConvert(max)[2]])) and i = nextprime(i-1) then max := i: fi:
	if evalb(evalf(N[myConvert(i)[1]][myConvert(i)[2]]) < evalf(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:

TestHasse := proc() local a, b, T:
  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:
  return T:
end: