The following commands for the Magma computer-algebra system check various group properties for the usual doubling and addition formulas for affine coordinates.
Addition is on the curve. The following commands start from affine points (x1,y1) and (x2,y2) on the curve; feed (x1,y1) and (x2,y2) to the affine addition formulas to obtain (x3,y3); and check that (x3,y3) is on the curve:
K<a,b,x1,x2>:=FieldOfFractions(PolynomialRing(Rationals(),4));
R<y1,y2>:=PolynomialRing(K,2);
S:=quo<R|y1^2-x1^3-a*x1-b,y2^2-x2^3-a*x2-b>;
S!(y1^2-x1^3-a*x1-b);
S!(y2^2-x2^3-a*x2-b);
lambda:=(y2-y1)/(x2-x1);
x3:=lambda^2-x1-x2; y3:=lambda*(x1-x3)-y1;
S!(y3^2-x3^3-a*x3-b);
Doubling is on the curve. The following commands start from an affine point (x1,y1) on the curve; feed (x1,y1) to the affine doubling formulas to obtain (x3,y3); and check that (x3,y3) is on the curve:
K<a,b,x1>:=FieldOfFractions(PolynomialRing(Rationals(),3));
R<y1>:=PolynomialRing(K,1);
S:=quo<R|y1^2-x1^3-a*x1-b>;
x2:=x1; y2:=y1;
S!(y1^2-x1^3-a*x1-b);
S!(y2^2-x2^3-a*x2-b);
lambda:=(3*x1^2+a)/(2*y1);
x3:=lambda^2-x1-x2; y3:=lambda*(x1-x3)-y1;
S!(y3^2-x3^3-a*x3-b);
Addition associates with addition.
K<a,b,x1,x2,x4>:=FieldOfFractions(PolynomialRing(Rationals(),5));
R<y1,y2,y4>:=PolynomialRing(K,3);
S:=quo<R|y1^2-x1^3-a*x1-b,y2^2-x2^3-a*x2-b,y4^2-x4^3-a*x4-b>;
lambda:=(y2-y1)/(x2-x1); x3:=lambda^2-x1-x2; y3:=lambda*(x1-x3)-y1;
lambda:=(y4-y3)/(x4-x3); x7:=lambda^2-x3-x4; y7:=lambda*(x3-x7)-y3;
lambda:=(y4-y2)/(x4-x2); x6:=lambda^2-x2-x4; y6:=lambda*(x2-x6)-y2;
lambda:=(y6-y1)/(x6-x1); x9:=lambda^2-x1-x6; y9:=lambda*(x1-x9)-y1;
S!(y1^2-x1^3-a*x1-b);
S!(y2^2-x2^3-a*x2-b);
S!(y3^2-x3^3-a*x3-b);
S!(y4^2-x4^3-a*x4-b);
S!(y6^2-x6^3-a*x6-b);
S!(y7^2-x7^3-a*x7-b);
S!(y9^2-x9^3-a*x9-b);
S!(x9-x7); S!(y9-y7);
Addition associates with doubling.
K<a,b,x1,x2>:=FieldOfFractions(PolynomialRing(Rationals(),4));
x4:=x2;
R<y1,y2>:=PolynomialRing(K,2);
y4:=y2;
S:=quo<R|y1^2-x1^3-a*x1-b,y2^2-x2^3-a*x2-b>;
lambda:=(y2-y1)/(x2-x1); x3:=lambda^2-x1-x2; y3:=lambda*(x1-x3)-y1;
lambda:=(y4-y3)/(x4-x3); x7:=lambda^2-x3-x4; y7:=lambda*(x3-x7)-y3;
lambda:=(3*x2^2+a)/(2*y2); x6:=lambda^2-x2-x4; y6:=lambda*(x2-x6)-y2;
lambda:=(y6-y1)/(x6-x1); x9:=lambda^2-x1-x6; y9:=lambda*(x1-x9)-y1;
S!(y1^2-x1^3-a*x1-b);
S!(y2^2-x2^3-a*x2-b);
S!(y3^2-x3^3-a*x3-b);
S!(y4^2-x4^3-a*x4-b);
S!(y6^2-x6^3-a*x6-b);
S!(y7^2-x7^3-a*x7-b);
S!(y9^2-x9^3-a*x9-b);
S!(x9-x7); S!(y9-y7);