TEENIE_WEENIE := 10^(-8);

function other_half(points)
   return Reverse([<a[1],-a[2]> : a in points]);
end function;

intrinsic PointString(component::SeqEnum, ymin::FldPrElt, ymax::FldPrElt) 
                            -> MonStgElt
{}
   AssertAttribute(FldPr,"OutputPrecision", 5);
   ans := "";
   for p in [q : q in component | q[2] ge ymin and q[2] le ymax] do
      x, y := Explode(p);
      if Abs(x) lt TEENIE_WEENIE then
         x := 0;
      end if; 
      if Abs(y) lt TEENIE_WEENIE then
         y := 0;
      end if; 
      ans := ans cat "(" cat Sprint(x) cat "," cat Sprint(y) cat ")\n";
   end for;
   return ans;
end intrinsic;

function RealRoots(E)
   a := aInvariants(E);
   assert a[1] eq 0 and a[3] eq 0;
   R<x> := PolynomialRing(RealField());
   f := x^3 + a[2]*x^2 + a[4]*x + a[5];
   real_roots := [r[1] : r in Roots(f) | IsReal(r[1])];
   return real_roots;
end function;

function LocalMax(E)
   a := aInvariants(E);
   assert a[1] eq 0 and a[3] eq 0;
   a2 := a[2]*1.0; a4 := a[4]*1.0; a6 := a[5]*1.0;
   if 4*a2^2 - 12*a4 ge 0 then
      x0 := (-2*a2 + Sqrt(4*a2^2 - 12*a4))/6.0;
      x1 := (-2*a2 - Sqrt(4*a2^2 - 12*a4))/6.0;
      "x = ",x0,x1;
      return Max(Sqrt(Abs(x0^3+a2*x0^2 + a4*x0+a6)),Sqrt(Abs(x1^3+a2*x1^2 + a4*x1+a6)));
   end if;
   return 0;
end function;

intrinsic PlotEllipticCurve(E::CrvEll, 
                         xmin::FldPrElt, xmax::FldPrElt, ybound::FldPrElt, 
                          res::RngIntElt) 
   -> SeqEnum, SeqEnum
{}
   /* For each x between xmin and xmax at res discrete intervals, compute
      the two roots y of the Wierstrass equation for E. Do this using the
      quadratic formula.  
   */      

   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
   require a1 eq 0 and a3 eq 0: 
        "Currently a1 and a3 of the a-invariants of argument 1 must be 0.";

   component_1 := [];   
   component_2 := [];
   phase := 1;

   interval := (xmax - xmin)*1.0/res;
   x := xmax;
   while x ge xmin and phase lt 4 do
      y := Sqrt(x^3+a2*x^2+a4*x+a6);
      if IsReal(y) then 
         if phase eq 1 then
//            if Abs(y) le ybound then
               Append(~component_1, <x,y>);
//            end if;
         elif phase eq 3 then
//            if Abs(y) le ybound then
               Append(~component_2, <x,y>);
//            end if;
         elif phase eq 2 then
//            if Abs(y) le ybound then

//               Append(~component_2, <RealRoots(E)[2],0.0>);

               Append(~component_2, <x,y>);
//            end if;
            phase := 3;
         end if;
      else   // not real
         if phase eq 1 then
            component_1 := component_1 cat other_half(component_1);
            phase := 2;
         elif phase eq 3 then
            component_2 := component_2 cat other_half(component_2);   
            phase := 4;
         end if;
      end if;
      x -:= interval;
   end while;

   return component_1, component_2;

end intrinsic;

function CurveModel(E)
   R<x,y> := PolynomialRing(Rationals(),2);
   a:=aInvariants(E);
   s := Sprintf("%o = %o", y^2 + a[1]*x*y + a[3]*y, 
		             x^3 + a[2]*x^2 + a[4]*x+a[5]);
   return &cat[s[i] : i in [1..#s] | s[i] ne "*"];
end function;

function String(x)
   AssertAttribute(FldPr,"OutputPrecision", 5);
   return Sprint(x);
end function;

function Point(x,y)
   return Sprintf("(%o,%o)", String(x), String(y));
end function;

intrinsic MakePage(E::CrvEll, P1::PtEll, P2::PtEll, 
		 res::RngIntElt, fname::MonStgElt)
{}
   file := Open(fname cat ".tex","w");
   fprintf file, "%o", MakePage(E,P1,P2,res);
   Flush(file);
   System("latex " cat fname cat ".tex; dvips " cat fname cat ".dvi -o " cat fname cat ".ps" );
end intrinsic;

intrinsic MakePage(E::CrvEll, P1::PtEll, P2::PtEll, 
		 res::RngIntElt) -> MonStgElt
{A LaTeX file that gives a colorful diagram of E along with a choice
 of two points on E and the suggestion to add together these two points.
 The two points and their sum all fit in the diagram, and the diagram
 contains a grid.}

   /* Compute the following:
      unit      
      x0, y0    lower-left corner
      x1, y1    upper-right corner
      P3        sum of P1 and P2
    */

   a := aInvariants(E);
   real_roots := RealRoots(E);
   local_max := LocalMax(E);
"local_max = ",local_max;

   P3 := P1 + P2;
   mP3 := -P3;
   x0 := Min([-1,P1[1], P2[1], P3[1], mP3[1]] cat real_roots)*1.3 - 0.1;
   x1 := Max([1,P1[1], P2[1], P3[1], mP3[1]] cat real_roots)*1.3 + 0.1;

   y0 := Min([P1[2], P2[2], P3[2], mP3[2], -local_max, -1])*1.3 - 0.1;
   y1 := Max([P1[2], P2[2], P3[2], mP3[2], local_max, 1])*1.3 + 0.1;

//   y1 := Max([Abs(P1[2]), Abs(P2[2]), Abs(P3[2]), Abs(mP3[2])])*1.0 + 1;
//   y0 := -y1;
   unit := Min(11.0/(x1-x0), 18.0/(y1-y0));
   bit  := Max(y1-y0,x1-x0)/25.0;
   c1, c2 := PlotEllipticCurve(E, x0, x1, y1, res);
   disk := Sprint(bit/5.0);

   page := "\\documentclass[12pt]{article}\n" cat
           "\\usepackage{pstricks}\n"         cat
           "\\pagestyle{empty}\n"             cat
           "\\textwidth=1.2\\textwidth\n"     cat
           "\\hoffset=-.7in\n"                cat
           "\\textheight=1.25\\textheight\n"   cat
           "\\voffset=-.8in\n"                cat
           "\\begin{document}\n"              cat
           "Add $" cat Point(P1[1],P1[2]) cat "$ and $" cat 
                       Point(P2[1],P2[2]) cat "$ together!\\vspace{2ex}\n\n" cat
           "\n{\\LARGE$$" cat CurveModel(E) cat "\\vspace{2ex}$$}\n\n" cat
           "\\begin{center}\n" cat 
           "\\psset{unit=" cat String(unit) cat "}\n" cat
           "\\pspicture" cat Point(x0,y0) cat Point(x1,y1) cat "\n" cat
           "\\psgrid[gridcolor=gray, subgriddiv=1]\n\n" cat
           "\\pscircle*[linecolor=red]" cat Point(P1[1],P1[2]) cat 
                     "{" cat disk cat "}\\rput" cat Point(P1[1]+bit,P1[2]+bit/5.0) cat
                     "{$" cat Point(P1[1],P1[2]) cat "$}\n" cat
           "\\pscircle*[linecolor=red]" cat Point(P2[1],P2[2]) cat 
                     "{" cat disk cat "}\\rput" cat Point(P2[1]+bit,P2[2]+bit/5.0) cat
                     "{$" cat Point(P2[1],P2[2]) cat "$}\n\n" cat
           "\\psline[linewidth=0.03]{->}" cat Point(x0,0) cat Point(x1,0) cat
           "    \\rput" cat Point(x1+bit/2,0) cat "{$x$}\n" cat 
           "\\psline[linewidth=0.03]{->}" cat Point(0,y0) cat Point(0,y1) cat
           "    \\rput" cat Point(0,y1+bit/2) cat "{$y$}\n\n" cat 
           "\\pscurve[linecolor=blue]\n" cat PointString(c1, y0, y1) cat "\n\n" cat
           (#c2 gt 0 select 
              "\\psccurve[linecolor=blue]\n" cat PointString(c2, y0, y1) cat "\n\n" 
           else "") cat
           "\\endpspicture\n" cat
           "\\end{center}\n\n" cat
           "\\end{document}\n";

   return page;   

end intrinsic;


intrinsic FindGoodCurves(arange::SeqEnum, 
                         brange::SeqEnum, 
                         hbound::FldPrElt) -> SeqEnum
{Sequence of triples <E, P1, P2>, where E is an elliptic curve
defined by an equation of the form y^2=x^3+a*X+b with a in arange
and b in brange, and P1, P2, and P1+P2 are points on E of small height.}

   ans := [* *];
   analyzed := [];   
   for a in arange do
      for b in brange do
         print a,b;
         if -4*a^3 - 27*b^2 ne 0 then
            E := EllipticCurve([a,b]);
            min := MinimalModel(E);
            if Index(analyzed,min) eq 0 then
               Append(~analyzed, min);
               points := IntegralPoints(E : Bound := hbound);
               if #points ge 3 then
                  P1 := points[1];
                  P2 := points[2];
                  P3 := points[3];
                  if Height(P3) le hbound then
                     found := <aInvariants(E), P1, P2>;
                     Append(~ans, found);
                     print "found ", found;
                  end if;
               end if;
            else
               print "skipping -- already considered equivalent curve";
            end if;
         else
            print "skipping -- singular";
         end if;
      end for;
   end for;
   return ans;
end intrinsic;

intrinsic MakeHandouts(good::List, prefix::MonStgElt)
{}

   for i in [1..#good] do
      fname := Sprintf("%o%o", prefix, i);
      print "Making: ",fname;
      g := good[i];
      MakePage(g[1], g[2], g[3], 300, fname);
   end for;
end intrinsic;