# OK to post

############################################################## OLD CODE ##############################################################

with(plottools):
with(plots):
with(plots, implicitplot):

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
 printf("Something bad happened\n"):
 return FAIL:
fi:
sol := sol minus {P[1],Q[1]}:
xR := sol[1]:
return [xR, -sqrt(xR^3 + a * xR + b)]:
end:

############################################################## END OLD CODE ##############################################################

# PART 1: Diag() PROCEDURE:
# -----------------------

toPlot:=proc(point, name):
	return op([plot(<<point[1]> | <point[2]>>, style = "point", color = "Black", symbol = "solidcircle"), textplot([op(point), name], `align`={`above`,`right`})]);
end;

Diag:=proc(a,b) local A,B,C,AB,graph:
	do:
		A:=RP(a,b);
		B:=RP(a,b);
	until A[1] <> B[1];

	AB:=AddE(a,b,A,B);
	C:=[AB[1],-AB[2]];
	graph:=implicitplot({y^2=x^3+a*x+b});
	display({graph,toPlot(A, "A"),toPlot(B, "B"),toPlot(C, "C"),toPlot(AB, "A+B"),line(A,C,color="blue",linestyle="solid"),line(C,AB,color="blue",linestyle="solid")});
end;

# PART 2: IntSol() PROCEDURE:
# -----------------------

IntSol:=proc(a,b,K) local S,x,y:
	S:={};
	for x from -K to K do
		y:=sqrt(x^3+a*x+b);
		if x^3+a*x+b>=0 and y=floor(y) then
			S:=S union {[x,y],[x,-y]};
		fi;
	od;

	return S;
end;

# PART 3: SolChain() PROCEDURE:
# -----------------------

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;

	return L;
end;