#OK to post homework
#GLORIA LIU, Apr 4 2024, Assignment 20

read `C20.txt`:

#=====================================#
#1. Using Maple's graphing command, write a procedure
#Diag(a,b)
#that draws a picture like Fig. 3 (with the points A, B,C, A+B labeled and indicated on the diagram) 
#for a general elliptic curve y^2=x^3+a*x+b , for general a, and b, and A and B picked using RP(a,b)

with(plots):
with(plottools):

Diag3:=proc(a,b) local A,B,Btemp,C,AB,slope,ALabel,BLabel,CLabel,ABLabel,bottomPlot,linePlot,topPlot, verticalLine:
#Estimating the values from figure 3
A:=[-6/5,sqrt((-6/5)^3+a*(-6/5)+b)]:
B:=[-1/3,sqrt((-1/3)^3+a*(-1/3)+b)]:
AB:=AddE(a,b,A,B):
C:=[AB[1], -AB[2]]:
slope:=(C[2]-A[2])/(C[1]-A[1]):
topPlot:=plot(sqrt(x^3 + a*x + b), x=-2..2):
bottomPlot:=plot(-sqrt(x^3 + a*x + b), x=-2..2):
linePlot:=plot(slope*x+A[2]-slope*A[1], x=A[1]..C[1]):
verticalLine:=line(C,AB):
#points:=pointplot([A,B,C,AB]):
ALabel:=textplot([A[1], A[2], "A"]):
BLabel:=textplot([B[1], B[2], "B"]):
CLabel:=textplot([C[1], C[2], "C"]):
ABLabel:=textplot([AB[1], AB[2], "A+B"]):
display(topPlot, bottomPlot, linePlot, ALabel, BLabel, CLabel, ABLabel, verticalLine);
end:

#Shows something similar to figure 3 in the paper
Diag3(-1,1);

#General version
Diag:=proc(a,b) local A,B,Btemp,C,AB,slope,ALabel,BLabel,CLabel,ABLabel,bottomPlot,linePlot,topPlot,verticalLine,Atemp:
#Estimating the values from figure 3
A:=RP(a,b):
B:=RP(a,b):
if A[1] > B[1] then
 Atemp:=A:
 A:=B:
 B:=Atemp:
fi:
AB:=AddE(a,b,A,B):
C:=[AB[1], -AB[2]]:
slope:=(C[2]-A[2])/(C[1]-A[1]):
topPlot:=plot(sqrt(x^3 + a*x + b), x=-2..2):
bottomPlot:=plot(-sqrt(x^3 + a*x + b), x=-2..2):
linePlot:=plot(slope*x+A[2]-slope*A[1], x=max(A[1],B[1])..C[1]):
verticalLine:=line(C,AB):
ALabel:=textplot([A[1], A[2], "A"]):
BLabel:=textplot([B[1], B[2], "B"]):
CLabel:=textplot([C[1], C[2], "C"]):
ABLabel:=textplot([AB[1], AB[2], "A+B"]):
display(topPlot, bottomPlot, linePlot, ALabel, BLabel, CLabel, ABLabel, verticalLine);
end:

Diag(3,5);

#=====================================#
#2. Write a procedure
#IntSol(a,b,K)
#that finds all integer solutions [x,y] with -K<=x<=K of the diopphantine equation
#y^2=x^3+a*x+b

IntSol:=proc(a,b,K) local x,y,sol,y1,y2:
sol:={}:
for x from -K to K do
y1:=sqrt(x^3+a*x+b):
y2:=-sqrt(x^3+a*x+b):
if type(y1,integer) then
 sol:=sol union {[x,y1],[x,y2]}:
fi:
od:
sol:
end:

IntSol(3,5,20);
IntSol(-1,1,20);


#=====================================#
#3. Write a procedure
#SolChain(a,b,S1,S2,K)
#That starts out with two solutions S1 and S2 of y^2=x^3+a*x+b and repeatedly uses 
#AddE(a,b,L[-2],L[-1]) to get a list of length K+2 of solutions in rational numbers. Can you get 
#infinitely solutions this way?

SolChain:=proc(a,b,S1,S2,K) local k,L,ng:
L:=[S1,S2]:
for k from 1 to K+2 do
 ng:=AddE(a,b,L[-2],L[-1]):
 L:=[op(L), ng]:
od:
L:
end:

#SolChain(3,5,RP(3,5),RP(3,5),10);
#I think it should be possible to find infinitely many solutions this way, unless the solutions repeat