/////ATTENTION//Sage code for addition formula and algorithms below
///////Magma Code to obtain points that the sum does not belong to the curve. i.e the addition is not complete. Please change the value of the prime p to obtain other points///////////////

AddProjThetaFirst := function(P,Q,c0,c2)
    X0 := P[1]; X1 := P[2]; X2 := P[3]; X3 := P[4];
    Y0 := Q[1]; Y1 := Q[2]; Y2 := Q[3]; Y3 := Q[4];
 
    Z0 := (X0^2*Y0^2 + X2^2*Y2^2) - 4*(c2/c0)*X1*X3*Y1*Y3;
    Z1 := c0*(X0*X1*Y0*Y1 + X2*X3*Y2*Y3) - 2*c2*(X2*X3*Y0*Y1 + X0*X1*Y2*Y3);
    Z2 := (X1^2*Y1^2 + X3^2*Y3^2) - 4*(c2/c0)*X0*X2*Y0*Y2;
    Z3 := c0*(X0*X3*Y0*Y3 + X1*X2*Y1*Y2) - 2*c2*(X0*X3*Y1*Y2 + X1*X2*Y0*Y3);
    return [Z0,Z1,Z2,Z3];
end function;
 
AddProjThetaSecond := function(P,Q,c0,c2)
    X0 := P[1]; X1 := P[2]; X2 := P[3]; X3 := P[4];
    Y0 := Q[1]; Y1 := Q[2]; Y2 := Q[3]; Y3 := Q[4];
   
   Z0 := c0*(X0*X1*Y0*Y3 + X2*X3*Y1*Y2) - 2*c2*(X0*X1*Y1*Y2 + X2*X3*Y0*Y3);
   Z1 := X1^2*Y0^2 + X3^2*Y2^2 - 4*(c2/c0)*X0*X2*Y1*Y3;
   Z2 := c0*(X0*X3*Y2*Y3 + X1*X2*Y0*Y1) - 2*c2*(X0*X3*Y0*Y1 + X1*X2*Y2*Y3);
   Z3 := X0^2*Y3^2 + X2^2*Y1^2 - 4*(c2/c0)*X1*X3*Y0*Y2;
    return [Z0,Z1,Z2,Z3];
end function;
 
 
AddProjThetaThird := function(P,Q,c0,c2)
    X0 := P[1]; X1 := P[2]; X2 := P[3]; X3 := P[4];
    Y0 := Q[1]; Y1 := Q[2]; Y2 := Q[3]; Y3 := Q[4];
   
   Z0 := X1^2*Y3^2 + X3^2*Y1^2 - 4*(c2/c0)*X0*X2*Y0*Y2;
   Z1 := c0*(X0*X3*Y1*Y2 + X1*X2*Y0*Y3) - 2*c2*(X0*X3*Y0*Y3 + X1*X2*Y1*Y2);
   Z2 := X0^2*Y2^2 + X2^2*Y0^2 - 4*(c2/c0)*X1*X3*Y1*Y3;
   Z3 := c0*(X0*X1*Y2*Y3 + X2*X3*Y1*Y0) - 2*c2*(X0*X1*Y0*Y1 + X2*X3*Y2*Y3);
    return [Z0,Z1,Z2,Z3];
end function;
 
 
AddProjThetaForth := function(P,Q,c0,c2)
    X0 := P[1]; X1 := P[2]; X2 := P[3]; X3 := P[4];
    Y0 := Q[1]; Y1 := Q[2]; Y2 := Q[3]; Y3 := Q[4];
   
   Z0 := c0*(X0*X3*Y0*Y1 + X1*X2*Y2*Y3) - 2*c2*(X0*X3*Y2*Y3 + X1*X2*Y0*Y1);
   Z1 := X0^2*Y1^2 + X2^2*Y3^2 - 4*(c2/c0)*X1*X3*Y0*Y2;
   Z2 := c0*(X0*X1*Y1*Y2 + X2*X3*Y0*Y3) - 2*c2*(X0*X1*Y0*Y3 + X2*X3*Y1*Y2);
   Z3 := X3^2*Y0^2 + X1^2*Y2^2 - 4*(c2/c0)*X0*X2*Y1*Y3;
    return [Z0,Z1,Z2,Z3];
end function;
 
 
AddThetaModel := function(P,Q,c0,c2)
  E := Curve(Parent(P));
    X0 := P[1]; X1 := P[2]; X2 := P[3]; X3 := P[4];
    Y0 := Q[1]; Y1 := Q[2]; Y2 := Q[3]; Y3 := Q[4];
 
   R1 := AddProjThetaFirst(P,Q,c0,c2);
   R2 := AddProjThetaSecond(P,Q,c0,c2);
   R3 := AddProjThetaThird(P,Q,c0,c2);
   R4 := AddProjThetaForth(P,Q,c0,c2);
 
   if R1 in E then
      //print "Formule 1 - One";
      return E!R1;
   elif R2 in E then
      print "Formule 2 - Two";
      return E!R2;
   elif R3 in E then
      print "Formule 3 - Three";
      return E!R3;
   elif R4 in E then
      print "Formule 4 - Forth";
      return E!R4;
   else
      return [0,0,0,0];
   end if;
end function;
 
 // Choice of the curve :
p := NextPrime(Random(2^3));
Fp := GF(p);
c0c2 := []; /// list of possible [c0,c2]
/// generation of  c0c2
for c0 in Fp do
 for c in Fp do
  if c0*c*(c0^2 + 4*c^2) eq 1 then c0c2 := Append(c0c2,[c0,c]); end if;
end for; end for;
c02 :=Random(c0c2); 
c0 := c02[1]; c2 := c02[2];
lbd := 1/(c0*c2);
 Prj<x0,x1,x2,x3> := ProjectiveSpace(Fp,3);
E1 := x0^2 + x2^2 - lbd*x1*x3;
E2 := x1^2 + x3^2 - lbd*x0*x2;
 E := Curve(Prj,[E1,E2]);
 for P in Points(E) do
   for Q in Points(E) do
      PQ := AddThetaModel(P,Q,c0,c2);
      if PQ notin E then print false, P,Q,PQ,E,c0,c2; end if;
   end for;
end for;





////////Formula 3 Verification of addition formula over any field

R.<c0,c2,X0,X1,X2,X3,Y0,Y1,Y2,Y3> = QQ[]
lbd = c0^2 + 4*c2^2; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
F1 = numerator(Y0^2 + Y2^2 - lbd*Y1*Y3); 
F2 = numerator(Y1^2 + Y3^2 - lbd*Y0*Y2)
S = R.quo([E1,E2,F1,F2,LB])
Z0 = c0*(X0*X1*Y0*Y3 + X2*X3*Y1*Y2) - 2*c2*(X0*X1*Y1*Y2 + X2*X3*Y0*Y3);
Z1 = X1^2*Y0^2 + X3^2*Y2^2 - 4*(c2/c0)*X0*X2*Y1*Y3;
Z2 = c0*(X0*X3*Y2*Y3 + X1*X2*Y0*Y1) - 2*c2*(X0*X3*Y0*Y1 + X1*X2*Y2*Y3);
Z3 =X0^2*Y3^2 + X2^2*Y1^2 - 4*(c2/c0)*X1*X3*Y0*Y2;
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0

////////Formula 4 Verification of addition formula over any field
R.<c0,c2,X0,X1,X2,X3,Y0,Y1,Y2,Y3> = QQ[]
lbd = c0^2 + 4*c2^2; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
F1 = numerator(Y0^2 + Y2^2 - lbd*Y1*Y3); 
F2 = numerator(Y1^2 + Y3^2 - lbd*Y0*Y2)
S = R.quo([E1,E2,F1,F2,LB])
Z0 = X1^2*Y3^2 + X3^2*Y1^2 - 4*(c2/c0)*X0*X2*Y0*Y2;
Z1 = c0*(X0*X3*Y1*Y2 + X1*X2*Y0*Y3) - 2*c2*(X0*X3*Y0*Y3 + X1*X2*Y1*Y2);
Z2 = X0^2*Y2^2 + X2^2*Y0^2 - 4*(c2/c0)*X1*X3*Y1*Y3;
Z3 = c0*(X0*X1*Y2*Y3 + X2*X3*Y1*Y0) - 2*c2*(X0*X1*Y0*Y1 + X2*X3*Y2*Y3);
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0

\\\\\\\Formula 5 verification of addition over any field
R.<c0,c2,X0,X1,X2,X3,Y0,Y1,Y2,Y3> = QQ[]
lbd = c0^2 + 4*c2^2; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
F1 = numerator(Y0^2 + Y2^2 - lbd*Y1*Y3); 
F2 = numerator(Y1^2 + Y3^2 - lbd*Y0*Y2)
S = R.quo([E1,E2,F1,F2,LB])
Z0 =c0*(X0*X3*Y0*Y1 + X1*X2*Y2*Y3) - 2*c2*(X0*X3*Y2*Y3 + X1*X2*Y0*Y1);
Z1 =X0^2*Y1^2 + X2^2*Y3^2 - 4*(c2/c0)*X1*X3*Y0*Y2;
Z2 = c0*(X0*X1*Y1*Y2 + X2*X3*Y0*Y3) - 2*c2*(X0*X1*Y0*Y3 + X2*X3*Y1*Y2);
Z3 = X3^2*Y0^2 + X1^2*Y2^2 - 4*(c2/c0)*X0*X2*Y1*Y3;
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0

////////Formula 19 Verification of doubling formula over any field
R.<c0,c2,X0,X1,X2,X3> = QQ[]
lbd = c0^2 + 4*c2^2; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
S = R.quo([E1,E2,LB])
Z0 = c0*X1*X3*(X0^2+X2^2) - 2*c2*X0*X2*(X1^2+X3^2);
Z1 = X1^2*X0^2 + X3^2*X2^2 - 4*(c2/c0)*X0*X2*X1*X3;
Z2 = c0*X0*X2*(X3^2 + X1^2) - 2*c2*X1*X3*(X0^2 + X2^2);
Z3 = X0^2*X3^2 + X2^2*X1^2 - 4*(c2/c0)*X1*X3*X0*X2;
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0

////////Formula 20 Verification of doubling formula over any field
R.<c0,c2,X0,X1,X2,X3> = QQ[]
lbd = c0^2 + 4*c2^2; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
S = R.quo([E1,E2,LB])
Z0 = 2*X1^2*X3^2  - 4*(c2/c0)*X0^2*X2^2;
Z1 = 2*c0*X0*X3*X1*X2  - 2*c2*(X0^2*X3^2 + X1^2*X2^2);
Z2 = 2*X0^2*X2^2  - 4*(c2/c0)*X1^2*X3^2;
Z3 = 2*c0*X0*X1*X2*X3 - 2*c2*(X0^2*X1^2 + X2^2*X3^2);
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0

////////Formula 21 Verification of addition formula over binary fields
R.<c0,c2,X0,X1,X2,X3,Y0,Y1,Y2,Y3> = GF(2)[]
lbd = c0^2 ; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
F1 = numerator(Y0^2 + Y2^2 - lbd*Y1*Y3); 
F2 = numerator(Y1^2 + Y3^2 - lbd*Y0*Y2)
S = R.quo([E1,E2,F1,F2,LB])
Z0 = c0*(X0*X1*Y0*Y3 + X2*X3*Y1*Y2) ;
Z1 = (X1*Y0 + X3*Y2)^2 ;
Z2 = c0*(X0*X3*Y2*Y3 + X1*X2*Y0*Y1) ;
Z3 =(X0*Y3+ X2*Y1)^2 ;
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0







////////Formula 22 Verification of addition formula over binary fields
R.<c0,c2,X0,X1,X2,X3,Y0,Y1,Y2,Y3> = GF(2)[]
lbd = c0^2 ; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
F1 = numerator(Y0^2 + Y2^2 - lbd*Y1*Y3); 
F2 = numerator(Y1^2 + Y3^2 - lbd*Y0*Y2)
S = R.quo([E1,E2,F1,F2,LB])
Z0 = (X1*Y3 + X3*Y1)^2 ;
Z1 = c0*(X0*X3*Y1*Y2 + X1*X2*Y0*Y3) ;
Z2 = (X0*Y2 + X2*Y0)^2 ;
Z3 = c0*(X0*X1*Y2*Y3 + X2*X3*Y1*Y0) ;
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0

\\\\\\\Formula 23 verification of addition formula over binary fields
R.<c0,c2,X0,X1,X2,X3,Y0,Y1,Y2,Y3> = GF(2)[]
lbd = c0^2 ; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
F1 = numerator(Y0^2 + Y2^2 - lbd*Y1*Y3); 
F2 = numerator(Y1^2 + Y3^2 - lbd*Y0*Y2)
S = R.quo([E1,E2,F1,F2,LB])
Z0 =c0*(X0*X3*Y0*Y1 + X1*X2*Y2*Y3) ;
Z1 =(X0*Y1 + X2*Y3)^2 ;
Z2 = c0*(X0*X1*Y1*Y2 + X2*X3*Y0*Y3) ;
Z3 = (X3*Y0 + X1*Y2)^2 ;
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0

////////Formula 24 Verification of doubling formula over binary fields
R.<c0,c2,X0,X1,X2,X3> = GF(2)[]
lbd = c0^2 ; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
S = R.quo([E1,E2,LB])
Z0 = c0*X1*X3*(X0+X2)^2 ;
Z1 = (X1*X0 + X3*X2)^2 ;
Z2 = c0*X0*X2*(X3 + X1)^2 ;
Z3 = (X0*X3 + X2*X1)^2 ;
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0




////////Table 1 : Verification of Algorithm in Table 1 addition formula 3 over any field
R.<c0,c2,X0,X1,X2,X3,Y0,Y1,Y2,Y3> = QQ[]
lbd = c0^2 + 4*c2^2; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
F1 = numerator(Y0^2 + Y2^2 - lbd*Y1*Y3); 
F2 = numerator(Y1^2 + Y3^2 - lbd*Y0*Y2)
S = R.quo([E1,E2,F1,F2,LB])
A=X0*Y3; B=X3*Y2; C=X2*Y1; D=X1*Y0; E=A^2; F=B^2;
G=C^2; H=D^2; I=(1/2)*((A+B)^2-E-F); J=(1/2)*((A+D)^2-E-H);
K=(1/2)*((C+D)^2-G-H); L=(1/2)*((C+B)^2-G-F);
U1=X0*X1; V1=X2*X3; U2=Y0*Y1; V2=Y2*Y3
N=(X0+X2)*(X3+X1)-U1-V1;
O=(Y0+Y2)*(Y3+Y1)-U2-V2;
P=U2+V2; 
Q1=U1+V1;
R1=N*P-I-K;
S1=O*Q1-J-L; 
T=(1/2)*((A+C)^2-E-G)
U=(1/2)*((B+D)^2-F-H)
Z0=c0*(J+L)-2*c2*S1; Z1=H+F-4*(c2/c0)*T;
Z2=c0*(I+K)-2*c2*R1; Z3=E+G-4*(c2/c0)*U;
U3=Z0*Z1; V3=Z2*Z3
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0



////////Table 2 : Verification of Algorithm in Table 2 addition formula 4 over any field
R.<c0,c2,X0,X1,X2,X3,Y0,Y1,Y2,Y3> = QQ[]
lbd = c0^2 + 4*c2^2; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
F1 = numerator(Y0^2 + Y2^2 - lbd*Y1*Y3); 
F2 = numerator(Y1^2 + Y3^2 - lbd*Y0*Y2)
S = R.quo([E1,E2,F1,F2,LB])
A=X1*Y3; B=X3*Y1; C=X0*Y2; D=X2*Y0; E=A^2; F=B^2;
G=C^2; H=D^2; I=(1/2)*((A+C)^2-E-G); J=(1/2)*((B+D)^2-F-H);
K=(1/2)*((A+B)^2-E-F); L=(1/2)*((A+D)^2-E-H);
U1=X0*X1; V1=X2*X3; U2=Y0*Y1; V2=Y2*Y3
N=(X0+X2)*(X3+X1)-U1-V1;
O=(Y0+Y2)*(Y3+Y1)-U2-V2;
P=N*O; 
T=(1/2)*((B+C)^2-F-G)
Q1=P-L-T;
R1=(U1+V1)*(U2+V2);
S1=R1-I-J
U=(1/2)*((D+C)^2-G-H)
Z0=E+F-4*(c2/c0)*U;   Z1=c0*(T+L)-2*c2*Q1;
Z2=G+H-4*(c2/c0)*K; Z3=c0*(I+J)-2*c2*S1;
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0




////////Table 3: Verification of Algorithm in Table 3 addition formula 5 over any field
R.<c0,c2,X0,X1,X2,X3,Y0,Y1,Y2,Y3> = QQ[]
lbd = c0^2 + 4*c2^2; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
F1 = numerator(Y0^2 + Y2^2 - lbd*Y1*Y3); 
F2 = numerator(Y1^2 + Y3^2 - lbd*Y0*Y2)
S = R.quo([E1,E2,F1,F2,LB])
A=X0*Y1; B=X2*Y3; C=X3*Y0; D=X1*Y2; E=A^2; F=B^2;
G=C^2; H=D^2; I=(1/2)*((A+B)^2-E-F); J=(1/2)*((A+C)^2-E-G);
K=(1/2)*((A+D)^2-E-H); L=(1/2)*((C+B)^2-G-F);
U1=X0*X1; V1=X2*X3; U2=Y0*Y1; V2=Y2*Y3
N=(X0+X2)*(X3+X1)-U1-V1;
O=(Y0+Y2)*(Y3+Y1)-U2-V2;
P=U2+V2; 
Q1=U1+V1;
T=(1/2)*((B+D)^2-F-H)
R1=N*P-J-T;
S1=O*Q1-K-L; 
U=(1/2)*((D+C)^2-G-H)
Z0=c0*(J+T)-2*c2*R1; Z1=E+F-4*(c2/c0)*U;
Z2=c0*(K+L)-2*c2*S1; Z3=G+H-4*(c2/c0)*I;
U3=Z0*Z1; V3=Z2*Z3
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0

///////Table 4 Verification of algorithm for formula 21 Addition over binary fields
R.<c0, c2, X0, X1, X2, X3, Y0, Y1, Y2, Y3> = GF(2)[]
 lbd = c0^2 + 4*c2^2;
 LB = numerator(lbd - 1/(c0*c2))
 E1 = numerator(X0^2+X2^2 - lbd*X1*X3);
 E2 = numerator(X1^2+X3^2 - lbd*X0*X2);
 F1 = numerator(Y0^2+Y2^2-lbd*Y1*Y3);
 F2 = numerator(Y1^2+Y3^2-lbd*Y0*Y2);
 S = R. quo([E1, E2, F1, F2, LB])
A = X0*Y3; B=X2*Y1; C=X3*Y2; D=X1*Y0; Z0=c0*(A*D+B*C);
Z2=c0*(A*C+B*D)
Z1=(D+C)^2; Z3=(A+B)^2
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2+Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0;  S(numerator(G2)) == 0


///////Table 4 Verification of algorithm for formula 22 Addition over binary fields
R.<c0, c2, X0, X1, X2, X3, Y0, Y1, Y2, Y3> = GF(2)[]
 lbd = c0^2 + 4*c2^2;
 LB = numerator(lbd - 1/(c0*c2))
 E1 = numerator(X0^2+X2^2 - lbd*X1*X3);
 E2 = numerator(X1^2+X3^2 - lbd*X0*X2);
 F1 = numerator(Y0^2+Y2^2-lbd*Y1*Y3);
 F2 = numerator(Y1^2+Y3^2-lbd*Y0*Y2);
 S = R. quo([E1, E2, F1, F2, LB])
A = X1*Y3; B=X3*Y1; C=X0*Y2; D=X2*Y0; 
Z1=c0*(B*C+A*D);
Z2=(C+D)^2
Z0=(A+B)^2; Z3=c0*(C*A+D*B)
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2+Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0;  S(numerator(G2)) == 0

///////Table 4 Verification of algorithm for formula 23 Addition over binary fields

R.<c0, c2, X0, X1, X2, X3, Y0, Y1, Y2, Y3> = GF(2)[]
 lbd = c0^2 + 4*c2^2;
 LB = numerator(lbd - 1/(c0*c2))
 E1 = numerator(X0^2+X2^2 - lbd*X1*X3);
 E2 = numerator(X1^2+X3^2 - lbd*X0*X2);
 F1 = numerator(Y0^2+Y2^2-lbd*Y1*Y3);
 F2 = numerator(Y1^2+Y3^2-lbd*Y0*Y2);
 S = R. quo([E1, E2, F1, F2, LB])
A = X0*Y1; B=X2*Y3; C=X3*Y0; D=X1*Y2; Z0=c0*(A*C+B*D);
Z2=c0*(A*D+B*C)
Z1=(A+B)^2; Z3=(C+D)^2
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2+Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0;  S(numerator(G2)) == 0



///////Table 5 Verification of algorithm for formula 19 Doubling over any field
R.<c0,c2,X0,X1,X2,X3> = QQ[]
lbd = c0^2 + 4*c2^2; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
S = R.quo([E1,E2,LB])
A=X1*X3; B=X0*X2; C=A^2; D=B^2; Z2=lbd*c0*D-2*lbd*c2*C;
G=A*B; Z0=lbd*c0*C-2*lbd*c2*D;
U1=X0*X1; V1=X2*X3;
F=V1^2; E=U1^2;
Z1=E+F-4*(c2/c0)*G;
Z3=lbd^2*G-E-F-4*(c2/c0)*U1*V1
V3=Z2*Z3; U3=Z0*Z1

G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0


///////Table 5 Verification of algorithm for formula 20 Doubling over any field

R.<c0,c2,X0,X1,X2,X3> = QQ[]
lbd = c0^2 + 4*c2^2; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
S = R.quo([E1,E2,LB])
A=X1*X3; B=X0*X2; C=A^2; D=B^2; 
U1=X0*X1; V1=X2*X3;
F=V1^2; E=U1^2; U3=Z0*Z1;
G=A*B; 
Z0=2*C-4*(c2/c0)*D;
Z2=2*D-4*(c2/c0)*C;
Z1=2*c0*U1*V1-2*c2*(((X0+X1)^2-2*U1)*((X2+X3)^2-2*V1)-D-C)
Z3=2*c0*U1*V1-2*c2*(E+F)
U3=Z0*Z1;V3=Z2*Z3
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0



//////// Table 6 : verification of Formula 24 doubling formula over binary fields
R.<c0,c2,X0,X1,X2,X3> = GF(2)[]
lbd = c0^2 ; 
LB = numerator(lbd -1/(c0*c2))
E1 = numerator(X0^2 + X2^2 - lbd*X1*X3); 
E2 = numerator(X1^2 + X3^2 - lbd*X0*X2)
S = R.quo([E1,E2,LB])
A=X1*X3; B=X0*X2; C=A^2; D=B^2; U1=X0*X1; V1=X2*X3; E=U1^2; F=V1^2
Z0=lbd*c0*C;
Z1=(U1+V1)^2;
Z2=lbd*c0*D;
Z3=lbd^2*A*B-E-F
V3=Z2*Z3; U3=Z0*Z1
G1 = Z0^2 + Z2^2 - lbd*Z1*Z3; G2 = Z1^2 + Z3^2 - lbd*Z0*Z2
S(numerator(G1)) == 0, S(numerator(G2)) == 0






////////Magma Code for the completeness of the set of Four addition formulas, one from the previous work and the three obtained in this work. Change the value of p any time to see that in all cases we have an output from one formula which gives a point on the curve.//////

AddProjThetaFirst := function(P,Q,c0,c2)
    X0 := P[1]; X1 := P[2]; X2 := P[3]; X3 := P[4];
    Y0 := Q[1]; Y1 := Q[2]; Y2 := Q[3]; Y3 := Q[4];
 
    Z0 := (X0^2*Y0^2 + X2^2*Y2^2) - 4*(c2/c0)*X1*X3*Y1*Y3;
    Z1 := c0*(X0*X1*Y0*Y1 + X2*X3*Y2*Y3) - 2*c2*(X2*X3*Y0*Y1 + X0*X1*Y2*Y3);
    Z2 := (X1^2*Y1^2 + X3^2*Y3^2) - 4*(c2/c0)*X0*X2*Y0*Y2;
    Z3 := c0*(X0*X3*Y0*Y3 + X1*X2*Y1*Y2) - 2*c2*(X0*X3*Y1*Y2 + X1*X2*Y0*Y3);
    return [Z0,Z1,Z2,Z3];
end function;
 
AddProjThetaSecond := function(P,Q,c0,c2)
    X0 := P[1]; X1 := P[2]; X2 := P[3]; X3 := P[4];
    Y0 := Q[1]; Y1 := Q[2]; Y2 := Q[3]; Y3 := Q[4];
   
   Z0 := c0*(X0*X1*Y0*Y3 + X2*X3*Y1*Y2) - 2*c2*(X0*X1*Y1*Y2 + X2*X3*Y0*Y3);
   Z1 := X1^2*Y0^2 + X3^2*Y2^2 - 4*(c2/c0)*X0*X2*Y1*Y3;
   Z2 := c0*(X0*X3*Y2*Y3 + X1*X2*Y0*Y1) - 2*c2*(X0*X3*Y0*Y1 + X1*X2*Y2*Y3);
   Z3 := X0^2*Y3^2 + X2^2*Y1^2 - 4*(c2/c0)*X1*X3*Y0*Y2;
    return [Z0,Z1,Z2,Z3];
end function;
 
 
AddProjThetaThird := function(P,Q,c0,c2)
    X0 := P[1]; X1 := P[2]; X2 := P[3]; X3 := P[4];
    Y0 := Q[1]; Y1 := Q[2]; Y2 := Q[3]; Y3 := Q[4];
   
   Z0 := X1^2*Y3^2 + X3^2*Y1^2 - 4*(c2/c0)*X0*X2*Y0*Y2;
   Z1 := c0*(X0*X3*Y1*Y2 + X1*X2*Y0*Y3) - 2*c2*(X0*X3*Y0*Y3 + X1*X2*Y1*Y2);
   Z2 := X0^2*Y2^2 + X2^2*Y0^2 - 4*(c2/c0)*X1*X3*Y1*Y3;
   Z3 := c0*(X0*X1*Y2*Y3 + X2*X3*Y1*Y0) - 2*c2*(X0*X1*Y0*Y1 + X2*X3*Y2*Y3);
    return [Z0,Z1,Z2,Z3];
end function;
 
 
AddProjThetaForth := function(P,Q,c0,c2)
    X0 := P[1]; X1 := P[2]; X2 := P[3]; X3 := P[4];
    Y0 := Q[1]; Y1 := Q[2]; Y2 := Q[3]; Y3 := Q[4];
   
   Z0 := c0*(X0*X3*Y0*Y1 + X1*X2*Y2*Y3) - 2*c2*(X0*X3*Y2*Y3 + X1*X2*Y0*Y1);
   Z1 := X0^2*Y1^2 + X2^2*Y3^2 - 4*(c2/c0)*X1*X3*Y0*Y2;
   Z2 := c0*(X0*X1*Y1*Y2 + X2*X3*Y0*Y3) - 2*c2*(X0*X1*Y0*Y3 + X2*X3*Y1*Y2);
   Z3 := X3^2*Y0^2 + X1^2*Y2^2 - 4*(c2/c0)*X0*X2*Y1*Y3;
    return [Z0,Z1,Z2,Z3];
end function;
 
 
AddThetaModel := function(P,Q,c0,c2)
  E := Curve(Parent(P));
    X0 := P[1]; X1 := P[2]; X2 := P[3]; X3 := P[4];
    Y0 := Q[1]; Y1 := Q[2]; Y2 := Q[3]; Y3 := Q[4];
 
   R1 := AddProjThetaFirst(P,Q,c0,c2);
   R2 := AddProjThetaSecond(P,Q,c0,c2);
   R3 := AddProjThetaThird(P,Q,c0,c2);
   R4 := AddProjThetaForth(P,Q,c0,c2);
 
   if R1 in E then
      print "Formule 1 - One";
      return E!R1;
   elif R2 in E then
      print "Formule 2 - Two";
      return E!R2;
   elif R3 in E then
      print "Formule 3 - Three";
      return E!R3;
   elif R4 in E then
      print "Formule 4 - Forth";
      return E!R4;
   else
      return [0,0,0,0];
   end if;
end function;
 

 ////An Example
// Choice of the curve :
p := NextPrime(Random(2^3));
Fp := GF(p);
 
c0c2 := []; /// liste des [c0,c2]
 
/// generation of c0 c2
for c0 in Fp do
 for c in Fp do
  if c0*c*(c0^2 + 4*c^2) eq 1 then c0c2 := Append(c0c2,[c0,c]); end if;
end for; end for;
 
c02 :=Random(c0c2); 
c0 := c02[1]; c2 := c02[2];
lbd := 1/(c0*c2);
 
Prj<x0,x1,x2,x3> := ProjectiveSpace(Fp,3);
 
E1 := x0^2 + x2^2 - lbd*x1*x3;
E2 := x1^2 + x3^2 - lbd*x0*x2;
 
E := Curve(Prj,[E1,E2]);
 
 
for P in Points(E) do
   for Q in Points(E) do
      PQ := AddThetaModel(P,Q,c0,c2);
      if PQ notin E then print false; end if;
   end for;
end for;


////////////////////THE END///////////////////////