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

with(plots):
with(plottools):

# 1.
Diag:=proc(a,b) local P,Q,C,v3,P1,P2,P3,P4,P5,P6,P7,L1,s,x,y:
 P:=RP(a,b):
 Q:=RP(a,b):
 v3:=AddE(a,b,P,Q):
 C:=[v3[1],-v3[2]]:


 P4:= plot(Vector([P[1],Q[1],C[1],v3[1]]), Vector([P[2], Q[2],C[2],v3[2]]), style = point, symbol = solidcircle, color = "Blue"):
 P6:= implicitplot(y^2=x^3+a*x+b):
 s:=(P[2]-Q[2])/(P[1]-Q[1]):
 L1:= plot(s*(x-C[1])+C[2]): # why are straight lines so hard to plot well !!!
 display([P4,P6,L1], color=[blue,red,blue]):

end:

# 2.
IntSol := proc(a,b,K) local x,y, S:
 S:={}:
 for x from -K to K do:
  if x^3+a*x+b>=0 then
    y:=sqrt(x^3+a*x+b): 
    if type(y,nonnegint) then
     S:=S union {[x,y]}:
    fi:
  fi:
 od:
 S:
end:

# 3. The solutions appear to repeat for sufficently large K.

SolChain:=proc(a,b,S1,S2,K) local L:
 L:=[S1,S2]:
 while nops(L)<K+2 do:
  L:=[op(L), AddE(a,b,L[-2],L[-1])]:
 od:
end:

######## C20 Code below:
#C20.txt: April 1, 2024
Help:=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: