#Homework20
#Nuray Kutlu
#OK TO POST

with(plottools):
with(plots):
#You might need to run a couple of times, 
#the next point is often randomly chosen outside of the screen
Diag:=proc(a,b) local A,B,C,ApB,y,x,sol,s, seg1,seg2,elp,D,seg3 :
A:= RP(a,b):
B:= RP(a,b):
s:=(A[2]-B[2])/(A[1]-B[1]):
y:=A[2]+(x-A[1])*s:
sol:={solve(y^2-x^3-a*x-b,x)};
sol:=sol minus {A[1],B[1]}:
C:=[sol[1],sqrt(sol[1]^3+a*sol[1]+b)] :
seg1:=plot([A,B], color = yellow):
seg2:=plot([B,C], color = green):
elp:=implicitplot(Y^2 = X^3 + a*X+b):
D:=AddE(a,b,A,B):
seg3:=plot([C,D], color = purple):
display(seg1,seg2,seg3,elp):
end:


IntSol:= proc(a,b,K) local solL, i:
solL:= []:
for i from -1*K to K do:
if type(sqrt(i^3+a*i+b), integer) then
solL:= [op(solL), [i,sqrt(i^3+a*i+b)]]:
fi:
od:
solL:
end:

#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 L,i:
L:=[S1,S2]:
for i from 1 to K do:
 L:=[op(L), AddE(a,b,L[-2],L[-1])]:
od:
L:
end:


#OLD CODE

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: