# Shaurya Baranwal, Homework 20 April 7, 2024
# OK to post

#---------------------------------------------
# Important Functions Necessary for Homework |
#---------------------------------------------

#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:

#---------------------------------------------------------------------------------------------
# Part 1 - Diag(a,b) that draws a picture general a, and b, and A and B picked using RP(a,b) |
#---------------------------------------------------------------------------------------------

with(plottools):
with(plots):
with(plots, implicitplot):
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:

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

Diag := proc(a, b) local A, B, C, AB, ip:

  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]]:
  ip := implicitplot({y^2 = x^3 + a * x + b}):

  display({ip, 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(a,b,K) that finds all integer solutions [x,y] with -K<=x<=K of the diopphantine equation |
#-----------------------------------------------------------------------------------------------------------

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 = trunc(y) then: 
      S := S union {[x, y], [x, -y]}: 
    fi:
  od:

  return S:

end:

#-------------------------------------------------------------------------------------------------------------------------------
# Part 3 - SolChain(a,b,S1,S2,K) starts out and repeatedly uses AddE(a,b,L[-2],L[-1]) to get a list of length K+2 of solutions |
#-------------------------------------------------------------------------------------------------------------------------------

SolChain := proc(a, b, S1, S2, K) local L:

  L := [S1, S2]:

  for i from 1 to K do:
    L := [op(L), AddE(a, b, L[-2], L[-1])]:
  od:

  return L:

end: