struct curve { real a=0; real b=8; real y2(real x) { return x^3+a*x+b; } real disc() { return -16*(4*a*a*a+27*b*b); } real lowx () { return sqrt(-a/3); } int comps() { if (a < 0) { real x=sqrt(-a/3); return y2(x) < 0 ? 2 : 1; } return 1; } void locus(picture pic=currentpicture, real m, real M, int n=100, pen p=currentpen) { path flip(path p, bool close) { path pp=reverse(yscale(-1)*p)..p; return close ? pp..cycle : pp; } path section(real m, real M, int n) { guide g; real width=(M-m)/n; for(int i=0; i <= n; ++i) { real x=m+width*i; real yy=y2(x); if (yy > 0) g=g..(x,sqrt(yy)); } return g; } if (comps() == 1) { draw(pic,flip(section(m,M,n),false),p); } else { real x=lowx(); // The minimum on x^3+ax+b if (m < x) draw(pic,flip(section(m,min(x,M),n),true),p); if (x < M) draw(pic,flip(section(max(x,m),M,n),false),p); } } pair neg(pair P) { return finite(P.y) ? yscale(-1)*P : P; } pair add(pair P, pair Q) { if (P.x == Q.x && P.x != Q.x) return (0,infinity); else { real lambda=P == Q ? (3*P.x^2+a)/(2*P.y) : (Q.y-P.y)/(Q.x-P.x); real Rx=lambda^2-P.x-Q.x; return (Rx,(P.x-Rx)*lambda-P.y); } } } import graph; import math; size(0,200); curve c; c.a=-1; c.b=4; pair oncurve(real x) { return (x,sqrt(c.y2(x))); } picture output; axes(); c.locus(-4,3,.3red+.7blue); pair P=oncurve(-1),Q=oncurve(1.2); pair PP=c.add(P,P),sum=c.add(P,Q); save(); drawline(P,Q,dashed); drawline(c.neg(sum),sum,dashed); dot("\$P\$", P, NW); dot("\$Q\$", Q, SSE); dot(c.neg(sum)); dot("\$P+Q\$", sum, 2SW); add(output,currentpicture.fit(),(-0.5cm,0),W); restore(); save(); drawline(P,c.neg(PP),dashed); drawline(c.neg(PP),PP,dashed); dot("\$P\$", P, NW); dot(c.neg(PP)); dot("\$2P\$", PP, SW); add(output,currentpicture.fit(),(0.5cm,0),E); shipout(output); restore();