#OK to post homework
#Ramesh Balaji,4/7/2024,hw20

# Part 1

with(plottools):
with(plots):
RP:=proc(a,b) local x: x:=rand(0..200)()/100: [x,sqrt(x^3+a*x+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:

Diag := proc(a,b):
    J:=plot([-sqrt(x^3 + a*x + b), sqrt(x^3 + a*x + b)]);

    A:=RP(a,b);
    B:=RP(a,b);
    AplusB:=AddE(a,b,A,B);



    display({J,point(A),point(B), point(AplusB), textplot([op(A), "A"]), textplot([op(B), "B"]), textplot([op(AplusB), "A + B"])});

end:

Diag(1,1);


# Part 2: brute force
IntSol := proc(a,b,K):
    S := {};
    for i from -K to K do:
        for j from -K to K do:
            if i^2 = j^3 + a*j + b then:
                S := S union {[i,j]};
            fi:
        od:
    od:
    return S;
end:

s:=IntSol(1,1,2);

# Part 3: SolChain

SolChain := proc(a,b,S1,S2,K):
    L := [S1, S2];
    print(L[-2]);
    print(L[-1]);
    for i from 1 to K do:
        L := [op(L), AddE(a,b,L[-2],L[-1])];
    od:
    return L;
end:

SolChain(1,1,RP(1,1),RP(1,1),3);