e(P):=P[2]^2-P[1]^3-a*P[1]^2-b*P[1]-c;  /* 楕円曲線の方程式 */
Y:[y1,y2,y3]; E:[e([x1,y1]),e([x2,y2]),e([x3,y3])];
elladd(P,Q):=                           /* 点の加法 P+Q を返す */
  block([L:0],[N:0],[R1:0],[R2:0],      /* 局所変数の定義 */
        L:ratsimp(Q[2]-P[2])/(Q[1]-P[1]),
        N:ratsimp(P[2]-L*P[1]),
        R1:ratsimp(L^2-P[1]-Q[1]-a),
        R2:-ratsimp(L*R1+N),
        return([R1,R2]));
srem(X,Y,Z):=divide(X,Y,Z)[2];   /* 剰余を与える関数を自作 */
ellred(X,i):=                    /* mod E[i] で yi の冪を下げる */
  block([S:0],[T:0],
        S:num(X), T:denom(X),        /* 分子，分母に別々に適用 */
     S:srem(S,E[i],Y[i]),T:srem(T,E[i],Y[i]),
        return(S/T));                    /* 分数に直して返す */
S:elladd([x1,y1],[x2,y2])$
S:elladd(S,[x3,y3])$              /* (P1+P2)+P3 */       
T:elladd([x2,y2],[x3,y3])$
T:elladd([x1,y1],T)$              /* P1+(P2+P3) */
D:ratsimp(S-T)$          /* 上２つの座標の差, 成分毎に計算してくれる */
    (for i:1 thru 3 do      /* E[i] で mod しないと一致せず */
         (D[j]:ratsimp(ellred(D[j],i))));
D;                                 /* D=[0,0] なら証明成功 */
