CoCalc Public Fileswww / talks / circles / handout.m
Author: William A. Stein
1TEENIE_WEENIE := 10^(-8);
2
3function other_half(points)
4   return Reverse([<a[1],-a[2]> : a in points]);
5end function;
6
7intrinsic PointString(component::SeqEnum, ymin::FldPrElt, ymax::FldPrElt)
8                            -> MonStgElt
9{}
10   AssertAttribute(FldPr,"OutputPrecision", 5);
11   ans := "";
12   for p in [q : q in component | q[2] ge ymin and q[2] le ymax] do
13      x, y := Explode(p);
14      if Abs(x) lt TEENIE_WEENIE then
15         x := 0;
16      end if;
17      if Abs(y) lt TEENIE_WEENIE then
18         y := 0;
19      end if;
20      ans := ans cat "(" cat Sprint(x) cat "," cat Sprint(y) cat ")\n";
21   end for;
22   return ans;
23end intrinsic;
24
25function RealRoots(E)
26   a := aInvariants(E);
27   assert a[1] eq 0 and a[3] eq 0;
28   R<x> := PolynomialRing(RealField());
29   f := x^3 + a[2]*x^2 + a[4]*x + a[5];
30   real_roots := [r[1] : r in Roots(f) | IsReal(r[1])];
31   return real_roots;
32end function;
33
34function LocalMax(E)
35   a := aInvariants(E);
36   assert a[1] eq 0 and a[3] eq 0;
37   a2 := a[2]*1.0; a4 := a[4]*1.0; a6 := a[5]*1.0;
38   if 4*a2^2 - 12*a4 ge 0 then
39      x0 := (-2*a2 + Sqrt(4*a2^2 - 12*a4))/6.0;
40      x1 := (-2*a2 - Sqrt(4*a2^2 - 12*a4))/6.0;
41      "x = ",x0,x1;
42      return Max(Sqrt(Abs(x0^3+a2*x0^2 + a4*x0+a6)),Sqrt(Abs(x1^3+a2*x1^2 + a4*x1+a6)));
43   end if;
44   return 0;
45end function;
46
47intrinsic PlotEllipticCurve(E::CrvEll,
48                         xmin::FldPrElt, xmax::FldPrElt, ybound::FldPrElt,
49                          res::RngIntElt)
50   -> SeqEnum, SeqEnum
51{}
52   /* For each x between xmin and xmax at res discrete intervals, compute
53      the two roots y of the Wierstrass equation for E. Do this using the
55   */
56
57   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
58   require a1 eq 0 and a3 eq 0:
59        "Currently a1 and a3 of the a-invariants of argument 1 must be 0.";
60
61   component_1 := [];
62   component_2 := [];
63   phase := 1;
64
65   interval := (xmax - xmin)*1.0/res;
66   x := xmax;
67   while x ge xmin and phase lt 4 do
68      y := Sqrt(x^3+a2*x^2+a4*x+a6);
69      if IsReal(y) then
70         if phase eq 1 then
71//            if Abs(y) le ybound then
72               Append(~component_1, <x,y>);
73//            end if;
74         elif phase eq 3 then
75//            if Abs(y) le ybound then
76               Append(~component_2, <x,y>);
77//            end if;
78         elif phase eq 2 then
79//            if Abs(y) le ybound then
80
81//               Append(~component_2, <RealRoots(E)[2],0.0>);
82
83               Append(~component_2, <x,y>);
84//            end if;
85            phase := 3;
86         end if;
87      else   // not real
88         if phase eq 1 then
89            component_1 := component_1 cat other_half(component_1);
90            phase := 2;
91         elif phase eq 3 then
92            component_2 := component_2 cat other_half(component_2);
93            phase := 4;
94         end if;
95      end if;
96      x -:= interval;
97   end while;
98
99   return component_1, component_2;
100
101end intrinsic;
102
103function CurveModel(E)
104   R<x,y> := PolynomialRing(Rationals(),2);
105   a:=aInvariants(E);
106   s := Sprintf("%o = %o", y^2 + a[1]*x*y + a[3]*y,
107		             x^3 + a[2]*x^2 + a[4]*x+a[5]);
108   return &cat[s[i] : i in [1..#s] | s[i] ne "*"];
109end function;
110
111function String(x)
112   AssertAttribute(FldPr,"OutputPrecision", 5);
113   return Sprint(x);
114end function;
115
116function Point(x,y)
117   return Sprintf("(%o,%o)", String(x), String(y));
118end function;
119
120intrinsic MakePage(E::CrvEll, P1::PtEll, P2::PtEll,
121		 res::RngIntElt, fname::MonStgElt)
122{}
123   file := Open(fname cat ".tex","w");
124   fprintf file, "%o", MakePage(E,P1,P2,res);
125   Flush(file);
126   System("latex " cat fname cat ".tex; dvips " cat fname cat ".dvi -o " cat fname cat ".ps" );
127end intrinsic;
128
129intrinsic MakePage(E::CrvEll, P1::PtEll, P2::PtEll,
130		 res::RngIntElt) -> MonStgElt
131{A LaTeX file that gives a colorful diagram of E along with a choice
132 of two points on E and the suggestion to add together these two points.
133 The two points and their sum all fit in the diagram, and the diagram
134 contains a grid.}
135
136   /* Compute the following:
137      unit
138      x0, y0    lower-left corner
139      x1, y1    upper-right corner
140      P3        sum of P1 and P2
141    */
142
143   a := aInvariants(E);
144   real_roots := RealRoots(E);
145   local_max := LocalMax(E);
146"local_max = ",local_max;
147
148   P3 := P1 + P2;
149   mP3 := -P3;
150   x0 := Min([-1,P1[1], P2[1], P3[1], mP3[1]] cat real_roots)*1.3 - 0.1;
151   x1 := Max([1,P1[1], P2[1], P3[1], mP3[1]] cat real_roots)*1.3 + 0.1;
152
153   y0 := Min([P1[2], P2[2], P3[2], mP3[2], -local_max, -1])*1.3 - 0.1;
154   y1 := Max([P1[2], P2[2], P3[2], mP3[2], local_max, 1])*1.3 + 0.1;
155
156//   y1 := Max([Abs(P1[2]), Abs(P2[2]), Abs(P3[2]), Abs(mP3[2])])*1.0 + 1;
157//   y0 := -y1;
158   unit := Min(11.0/(x1-x0), 18.0/(y1-y0));
159   bit  := Max(y1-y0,x1-x0)/25.0;
160   c1, c2 := PlotEllipticCurve(E, x0, x1, y1, res);
161   disk := Sprint(bit/5.0);
162
163   page := "\\documentclass[12pt]{article}\n" cat
164           "\\usepackage{pstricks}\n"         cat
165           "\\pagestyle{empty}\n"             cat
166           "\\textwidth=1.2\\textwidth\n"     cat
167           "\\hoffset=-.7in\n"                cat
168           "\\textheight=1.25\\textheight\n"   cat
169           "\\voffset=-.8in\n"                cat
170           "\\begin{document}\n"              cat
171           "Add $" cat Point(P1[1],P1[2]) cat "$ and $" cat 172 Point(P2[1],P2[2]) cat "$ together!\\vspace{2ex}\n\n" cat
173           "\n{\\LARGE$$" cat CurveModel(E) cat "\\vspace{2ex}$$}\n\n" cat
174           "\\begin{center}\n" cat
175           "\\psset{unit=" cat String(unit) cat "}\n" cat
176           "\\pspicture" cat Point(x0,y0) cat Point(x1,y1) cat "\n" cat
177           "\\psgrid[gridcolor=gray, subgriddiv=1]\n\n" cat
178           "\\pscircle*[linecolor=red]" cat Point(P1[1],P1[2]) cat
179                     "{" cat disk cat "}\\rput" cat Point(P1[1]+bit,P1[2]+bit/5.0) cat
180                     "{$" cat Point(P1[1],P1[2]) cat "$}\n" cat
181           "\\pscircle*[linecolor=red]" cat Point(P2[1],P2[2]) cat
182                     "{" cat disk cat "}\\rput" cat Point(P2[1]+bit,P2[2]+bit/5.0) cat
183                     "{$" cat Point(P2[1],P2[2]) cat "$}\n\n" cat
184           "\\psline[linewidth=0.03]{->}" cat Point(x0,0) cat Point(x1,0) cat
185           "    \\rput" cat Point(x1+bit/2,0) cat "{$x$}\n" cat
186           "\\psline[linewidth=0.03]{->}" cat Point(0,y0) cat Point(0,y1) cat
187           "    \\rput" cat Point(0,y1+bit/2) cat "{$y$}\n\n" cat
188           "\\pscurve[linecolor=blue]\n" cat PointString(c1, y0, y1) cat "\n\n" cat
189           (#c2 gt 0 select
190              "\\psccurve[linecolor=blue]\n" cat PointString(c2, y0, y1) cat "\n\n"
191           else "") cat
192           "\\endpspicture\n" cat
193           "\\end{center}\n\n" cat
194           "\\end{document}\n";
195
196   return page;
197
198end intrinsic;
199
200
201intrinsic FindGoodCurves(arange::SeqEnum,
202                         brange::SeqEnum,
203                         hbound::FldPrElt) -> SeqEnum
204{Sequence of triples <E, P1, P2>, where E is an elliptic curve
205defined by an equation of the form y^2=x^3+a*X+b with a in arange
206and b in brange, and P1, P2, and P1+P2 are points on E of small height.}
207
208   ans := [* *];
209   analyzed := [];
210   for a in arange do
211      for b in brange do
212         print a,b;
213         if -4*a^3 - 27*b^2 ne 0 then
214            E := EllipticCurve([a,b]);
215            min := MinimalModel(E);
216            if Index(analyzed,min) eq 0 then
217               Append(~analyzed, min);
218               points := IntegralPoints(E : Bound := hbound);
219               if #points ge 3 then
220                  P1 := points[1];
221                  P2 := points[2];
222                  P3 := points[3];
223                  if Height(P3) le hbound then
224                     found := <aInvariants(E), P1, P2>;
225                     Append(~ans, found);
226                     print "found ", found;
227                  end if;
228               end if;
229            else
230               print "skipping -- already considered equivalent curve";
231            end if;
232         else
233            print "skipping -- singular";
234         end if;
235      end for;
236   end for;
237   return ans;
238end intrinsic;
239
240intrinsic MakeHandouts(good::List, prefix::MonStgElt)
241{}
242
243   for i in [1..#good] do
244      fname := Sprintf("%o%o", prefix, i);
245      print "Making: ",fname;
246      g := good[i];
247      MakePage(g[1], g[2], g[3], 300, fname);
248   end for;
249end intrinsic;
250