%% Elliptic curve messing about

%% y^2 = x^3 + a x^2 + b

numeric u, v;
numeric a, b;
numeric minx, maxx, miny, maxy;
transform gt;
pair O;
%% prologues := 1;
a = -2; b = 1;
scale = 2 cm;
labeloffset := 6 pt;

ahlength := 2 mm;
ahangle := 30;

O = (0, 0);

vardef ecy(expr x) = sqrt (x*x*x + a*x + b) enddef;
vardef ecadd(expr u, v) =
  numeric ux, uy, vx, vy;
  numeric dx, dy;
  numeric m;
  ux := xpart u; uy = ypart u; 
  vx := xpart v; vy = ypart v;
  if ux = vx:
    m := (3*ux*ux + a)/(2*uy);
  else:
    m := (vy - uy)/(vx - ux);
  fi
  dx := m**2 - ux - vx;
  dy := m * (ux - dx) - uy;
  (dx, dy)
enddef;
 
minx = -2; maxx = 3;
maxy = ecy(maxx); miny = -maxy;
u = 7/5 cm;
v = u * (maxx - minx)/(maxy - miny);
gt = identity xscaled u yscaled v shifted (10 cm, 10 cm);

vardef eclabel@#(expr l, z) =
  pickup pencircle scaled 1/4 pt;
  draw fullcircle scaled 1 mm shifted (z transformed gt);
  draw thelabel@#(l, z transformed gt);
enddef;

vardef drawcurve =
  boolean ok;
  path p;
  numeric t;
  pickup pencircle scaled 1/4 pt;
  drawarrow ((minx, 0) -- (maxx, 0)) transformed gt;
  drawarrow ((0, miny) -- (0, maxy)) transformed gt;
  ok := false;
  pickup pencircle scaled 2/3 pt;
  for x = minx step 0.01 until maxx:
    t := x*x*x + a*x + b;
    if t >= 0:
      y := sqrt t;
      if ok:
	p := p .. (x, y);
      else:
	p := (x, y);
      fi
      ok := true;
    else:
      if ok:
	draw (reverse p reflectedabout (O, O + right) .. p .. cycle)
	  transformed gt;
      fi
      ok := false;
    fi
  endfor
  if ok:
    draw (reverse p reflectedabout (O, O + right) .. p) transformed gt;
  fi
enddef;

vardef setbackground(expr c) =
  save bboxmargin;
  picture pic;
  bboxmargin := 1/2 cm;
  pic := currentpicture;
  currentpicture := nullpicture;
  fill (bbox pic) withcolor c;
  draw pic;
enddef;

beginfig(0);
  drawcurve;
  x0 := -5/4; x1 := 1/4;
  y0 := -ecy(x0); y1 := ecy(x1);
  z2 = ecadd(z0, z1);
  x3 := x2; y3 := -y2;
  pickup pencircle scaled 1/4 pt;
  draw ((-0.1)[z0, z3] -- 1.1[z0, z3]) transformed gt dashed evenly;
  draw ((-0.1)[z2, z3] -- 1.1[z2, z3]) transformed gt dashed evenly;
  eclabel.bot(btex $P_0$ etex, z0);
  eclabel.top(btex $P_1$ etex, z1);
  eclabel.rt(btex $P'$ etex, z3);
  eclabel.lft(btex $P_2 = -P' = P_0 + P_1$ etex, z2);
  setbackground(white);
endfig;

beginfig(1);
  drawcurve;
  x0 := -18/16;
  y0 = ecy(x0);
  z1 = ecadd(z0, z0);
  x2 := x1; y2 := -y1;
  pickup pencircle scaled 1/4 pt;
  draw ((-0.1)[z0, z2] -- 1.1[z0, z2]) transformed gt dashed evenly;
  draw ((-0.1)[z2, z1] -- 1.1[z2, z1]) transformed gt dashed evenly;
  eclabel.top(btex $P_0$ etex, z0);
  eclabel.rt(btex $P'$ etex, z2);
  eclabel.lft(btex $P_1 = -P' = 2 P_0$ etex, z1);
  setbackground(white);
endfig;

beginfig(2);
  drawcurve;
  x0 := 5/2;
  y0 = ecy(x0);
  x1 := x0; y1 := -y0;
  pickup pencircle scaled 1/4 pt;
  draw ((x0, miny) -- (x0, maxy)) transformed gt dashed evenly;
  eclabel.rt(btex $P_0$ etex, z0);
  eclabel.rt(btex $P_1$ etex, z1);
  setbackground(white);
endfig;

end;