///Pairings on curves with k=9/////
(Pairings on curves k=15 and 27 below)

///parameters///
x:= 2^43+2^37+2^7+1;

p:=((x+1)^2+((x-1)^2*(2*x^3+1)^2)/3)/4;
p:=13522728075112739330542203653103566043035377929488752721485580527333685820130434460157787867627959369857;
r:=169440772426819274154778701260757423597114376011671712180533212939494849872257;
xp3:=x^3;////will be used for the new power
xp3:=712967262418449212250952106513994400129;
d:=(p^6+p^3+1)/r; ////to obtain the exponent x^3*d of the hard part of the final exponentiation
d:=36088428228002850386409293373098737304427590630970197659521680979527943311035364605465436062462975600109201460097188539174934296969288080217746500367150652002\
6253036975917714120482157250365252474729050203562164297004745136412100476484473\
3707398600109853810423030482453382659516090038079884577965949015874619575942567\
8687765066461827712381405748546584147130702015339331643500284545370139279058170\
9496437217512261293903885483630393669763997640642715452240264626054888817053601\
49841862416718843774760913150829027065128029353352413177741301661699;

//////Decomposition of the exponent x^3*d=xp3*d//////////
k2:=-(x-1)^2; k1:=x*k2; k0:=x*k1;  k5:=-x*k0; k4:=x*k5; k3:=x*k4+3;

/////Construction of finite fields
Fp:=FiniteField(p);
Fp3<s3>:=ExtensionField<Fp,x|x^3-7>;  
Fp9<w>:=ExtensionField<Fp3,y|y^3-s3>;

/////Verification of the decomposition of the exponent
A:=Random(Fp9);
IsZero(A^(xp3*d)-A^k0*A^(k1*p)*A^(k2*p^2)*A^(k3*p^3)*A^(k4*p^4)*A^(k5*p^5));

//////Construction of the curves
Ep:=EllipticCurve([Fp!(0),Fp!1]);
Ept:=EllipticCurve([Fp3!(0),Fp3!s3^2]); ///////The twist of Ep defined over Fp3////
Ep9:=EllipticCurve([Fp9!(0),Fp9!1]); /////The curve Ep defined over Fp9//////////
#Ep mod r;  ////The order of Ep is exactly divisible by r////

/////Looking for a point P of order r with Frobenius equal to 1
u:=#Ep/r;
u:=79807993562783988297252864;
P:=u*Random(Ep);

/////Looking for a point Q such that Frobenius(Q)=p*Q///////
ut:=#Ept/r;
ut:=1459402286881724979982121469312044562502140191528389579676827933003857909643480\
3037756338338439091512769469042749826318086088344472391971502866697442412700221\
902391558527653170453306215227600594312860608604181721806207906898339954691;

RR:=Random(Ept);
QRR:=ut*RR;
//////QRR is a point of order r and we send it to Ep9 on Fp9 to apply the Frobenius ///
x1:=QRR[1]*w^(-2);
y1:=QRR[2]*w^(-3);
xu:=x1^p;
yu:=y1^p;                                      
D1:=Ep9![xu,yu];
E1:=Ep9![x1,y1];
Q:=D1-E1;  /////// Q satisfies \pi(Q)=p*Q
Ep9![Q[1]^p,Q[2]^p]-p*Q;

///////The points P and Q are of order r
r*Q;
r*P;

/////The doubling step in the Miller Algorithm///////

Doubling:=function (S, K); 
x1:=S[1]*w^2; y1:=S[2]*w^3;
A:=3*x1^2; B:=(2*y1)^(-1); C:=A*B; x3:=C^2-2*x1; y3:=C*(x1-x3)-y1; f1:=K[2]*w^5+C*(x1*w^2-K[1]*w^4)-y1*w^2;
f2:=K[1]^2*w^4+x3*K[1]*w^2+x3^2;
f:=f1*f2;
return [x3*w^(-2),y3*w^(-3),f];
end function;

///////The Addition step in the Miller algorithm///////
Addition:=function(S,R, K);
x1:=S[1]*w^2; y1:=S[2]*w^3;
x2:=R[1]*w^2; y2:=R[2]*w^3;
A:=(y2-y1); B:=(x2-x1)^(-1); C:=A*B; x3:=C^2-x1-x2; 
y3:= C*(x1-x3)-y1; 
f1:=K[2]*w^5+C*(x1*w^2-K[1]*w^4)-y1*w^2;
f2:=K[1]^2*w^4+x3*K[1]*w^2+x3^2;
f:=f1*f2;
return [x3*w^(-2),y3*w^(-3),f];
end function;

///// (addition and doubling of point) are correct/////
l:=Addition(Q,2*Q,P);
IsZero(l[2]^2-l[1]^3-1);
l1:=Doubling(Q,P);
IsZero(l1[2]^2-l1[1]^3-1);

///The Main function to compute the optimal Ate pairing described in the article//
AteOpt:=function(Q ,S);
f:=1;
bin:=Intseq(Abs(x),2);
 T:=Q;
for i:=#bin-2 to 0 by -1 do
h:=Doubling(T ,S);
T:=[h[1],h[2]];
f:=h[3]*f^2;
if bin[i+1] eq 1 then 
h:=Addition(T,Q ,S); 
T:=[h[1],h[2]];
f:=h[3]*f;
end if;
end for;
//// instead of f:=f^(Integers()!((p^9-1)/r)), we compute a power;
A:=f^(p^3-1);  ////Easy part of the final exponentiation
f:=A^((xp3)*d); ////Hard part with the new exponent, formula already verified above///
return(f);
end function;
e1:=AteOpt(Q,P);

 //// The map is not trivial/////
e1;


//////Verification of the bilinearity 

eP25:=AteOpt(Q,25*P);
eQ25:=AteOpt(25*Q,P);
IsZero(eP25-e1^25);
IsZero(eQ25-e1^25);



///////PAIRING FOR CURVES k=15///
(See below for pairing with k=27)

////parameters
x:= 2^48+2^41+2^9+2^8+1;
p:=(x^12-2*x^11+x^10+x^7-2*x^6+x^5+x^2+x+1)/3;
r:=x^8-x^7+x^5-x^4+x^3-x+1;
p:=90513453509027837323015498966524467238615948271053123433025004142597079434172244878153242162500117528739179878419894291934572639098823619146067881891097102019540945709105921;
r:=41933031483931929173872138167306456628506916051359945068680984192147211979589804423702853969731915081303589283302401;
/////decomposition of the exponent of the hard part///////
xp3:=3*x^3;    ////will be used for the new power
xp3:=68482538809601545850589033778651877419850499;
d:=(p^10+p^5+1)/r; ////to obtain the exponent 3*x^3*d of the hard part of the final exponentiation
d:=880187294945139572654489462551036237678580244656016079451515968250171679053201813503642675099138703933904205252345733940911051610660916499102950174121607885151306649503349817719557134666706016903747729439038326336456738141887112827814101279404443852334788044278371858903527345454212580673250264180010017564774608041829926009481571580717058300978807438563990684211019860700632237932006141088944023350891154898234214088873296766288970796777852722598095974330559368086378910486985547304629065211601534933966114730252225307218545902631077624173802367054418396062484473230999374421036035552162154642125532702994251312595222329481924045503850227949090354735722379357048392054451328106561136158736571930239438454030720181341215393150681162473872036196887335433833332652673302378845352253067068741564131175543057264566104643365144912396421920216852009053623557986930901388219824180957603434539474388765298058944428200297967764876490925514341887789619979083477849993410549416088911416369757932514537202624321457757436185121112735364440213021833677198880688777650445356358185389103395751729681475535954854786018210892275662173433905100462264484087554151463991343435332985143720996084464929573096700456358773095729650704547141734065731099062306871808729615493651232135075230277518205526941953060515501136958100934598708638283522205431750985757346085130796859428603780396810473604773689019257878184254942102647448423309723112304936191035536415765943638184660563655421409143355460279154281459702001320862951005913197058888775193250606104356135271701768409074283369116179479843115468520492944081469726606466836213263272479033603;
//////Decomposition of the exponent 3*x^3*d=3*xp3*d//////////
k2:=-(((x-1)^2)*(x^2+x+1)); k1:=x*k2; k0:=x*k1;  k9:=-x*k0; k8:=x*k9; k7:=x*k8;
k6:=x*k7;  k5:=x*k6+3; M:=k2+ k5+k8;  k4:=M-k1-k7;   k3:=M-(k0+ k6+k9);

/////Construction of finite fields
Fp:=FiniteField(p);
Fp5<s3>:=ExtensionField<Fp,x|x^5-7>;  
Fp15<w>:=ExtensionField<Fp5,y|y^3-s3>;
/////Verification of the decomposition of the exponent
A:=Random(Fp15);
IsZero(A^(xp3*d)-A^k0*A^(k1*p)*A^(k2*p^2)*A^(k3*p^3)*A^(k4*p^4)*A^(k5*p^5)*A^(k6*p^6)*A^(k7*p^7)*A^(k8*p^8)*A^(k9*p^9));
//////Construction of the curves
Ep:=EllipticCurve([Fp!(0),Fp!1]);//Elliptic curve in Weierstrass form defined over Fp//
Ept:=EllipticCurve([Fp5!(0),Fp5!s3^2]); ///////The twist of Ep defined over Fp5////
Ep15:=EllipticCurve([Fp15!(0),Fp15!1]); /////The curve Ep defined over Fp15//////////
#Ep mod r;  ////The order of Ep is exactly divisible by r////

/////Looking for a point P of order r with Frobenius equal to 1
u:=#Ep/r;
u:=2158523968001482395107968013106160047076941111376761913344;
P:=u*Random(Ep);

/////Looking for a point Q such that Frobenius(Q)=p*Q///////
ut:=#Ept/r;
ut:=144880327035595369687124514092608601320347421444919827171019139801233855339400453769461260866261920797830681757217248509501145335175294461836198494600987633521172607071457191947132692379981865696460994477661256950884281101941271580655246040471466889007297439292888642885595170891943103096245403064763933016718187034151049084728878050960405167678168543856315735746317720716848904537596071324111093867972311624774127158424057675990480982314958400481595498404178902356213405389021224433480443476259800289516322287308958404534882220140415210248711322437602611330829237711856795189443788503535610257305155875300353336395464889330584329045357639727790691183456033659975748414378604885792200272787279460094757441590786027591950852917870234765846340947413251;
RR:=Random(Ept);
QRR:=ut*RR;
//////QRR is a point of order r and we send it to Ep15 on Fp15 to apply the Frobenius ///
x1:=QRR[1]*w^(-2);
y1:=QRR[2]*w^(-3);
// (x1,y1) is a point of Ep15 defined on Fp15 and we take the Frobenius////
xu:=x1^p;
yu:=y1^p;                                      
D1:=Ep15![xu,yu];
E1:=Ep15![x1,y1];
Q:=D1-E1;  /////// Q satisfies \pi(Q)=p*Q
Ep15![Q[1]^p,Q[2]^p]-p*Q;

///////The points P and Q are of order r
r*Q;
r*P;
///////////////////////////////////////////////The doubling step 
Doubling:=function (S, K); 
x1:=S[1]*w^2; y1:=S[2]*w^3;
A:=3*x1^2; B:=(2*y1)^(-1); C:=A*B; x3:=C^2-2*x1; y3:=C*(x1-x3)-y1; f1:=K[2]*w^5+C*(x1*w^2-K[1]*w^4)-y1*w^2;
f2:=K[1]^2*w^4+x3*K[1]*w^2+x3^2;
f:=f1*f2;
return [x3*w^(-2),y3*w^(-3),f];
end function;

/////////////////////////////////The Addition step in the Miller algorithm///////
Addition:=function(S,R, K);
x1:=S[1]*w^2; y1:=S[2]*w^3;
x2:=R[1]*w^2; y2:=R[2]*w^3;
A:=(y2-y1); B:=(x2-x1)^(-1); C:=A*B; x3:=C^2-x1-x2; 
y3:= C*(x1-x3)-y1; 
f1:=K[2]*w^5+C*(x1*w^2-K[1]*w^4)-y1*w^2;
f2:=K[1]^2*w^4+x3*K[1]*w^2+x3^2;
f:=f1*f2;
return [x3*w^(-2),y3*w^(-3),f];
end function;

///// ( addition and doubling of point) are correct/////
l:=Addition(Q,2*Q,P);
IsZero(l[2]^2-l[1]^3-1);
l1:=Doubling(Q,P);
IsZero(l1[2]^2-l1[1]^3-1);

///The Main function to compute the optimal Ate pairing described in the article//
AteOpt:=function(Q ,S);
f:=1;
bin:=Intseq(Abs(x),2);
 T:=Q;
for i:=#bin-2 to 0 by -1 do
h:=Doubling(T ,S);
T:=[h[1],h[2]];
f:=h[3]*f^2;
if bin[i+1] eq 1 then 
h:=Addition(T,Q ,S); 
T:=[h[1],h[2]];
f:=h[3]*f;
end if;
end for;
//// instead of f:=f^(Integers()!((p^15-1)/r)), we compute a power;
A:=f^(p^5-1);  ////Easy part of the final exponentiation
f:=A^((xp3)*d); ////Hard part with the new exponent, formula already verified above///
return(f);
end function;
e1:=AteOpt(Q,P);
  //// The map  is not trivial/////
e1;

//////Verification of the bilinearity 
eP25:=AteOpt(Q,25*P);
eQ25:=AteOpt(25*Q,P);
IsZero(eP25-e1^25);
IsZero(eQ25-e1^25);





////////PARINGS FOR CURVES WITH k=27

//////Parameters

x:= 2^29+2^19+2^17+2^14;
p:=(1/3)*(x-1)^2*(x^18+x^9+1)+x;
r:=(1/3)*(x^18+x^9+1);
p:=1352500393920932776288342881517739001698216127874508012163729668076196545366650832718124682950296167601825262362589324645313089939923437585202457165226294516117144846798645931;
r:=4680707861603025179856470723242682683313968691227380097658414641202569036266960847786881049178731378039261425821772294435544682009706547447370448789224204971;

/////Construction of finite fields
Fp:=FiniteField(p);
Fp3<s31>:=ExtensionField<Fp,x|x^3-7>;  
Fp9<s3>:=ExtensionField<Fp3,y|y^3-s31>;
Fp27<w>:=ExtensionField<Fp9,z|z^3-s3>;

//////Construction of the curves
Ep:=EllipticCurve([Fp!(0),Fp!(-2)]);
Ept:=EllipticCurve([Fp9!(0),Fp9!(-2*s3^2)]);
Ep27:=EllipticCurve([Fp27!(0),Fp27!(-2)]);
#Ep mod r;  ////The order of Ep is exactly divisible by r////

////Looking for a point P of order r with Frobenius equal to 1
u:=#Ep/r;
#Ep-p+x;
u:=288952105944449025;
P:=u*Random(Ep);

/////Looking for a point Q such that Frobenius(Q)=p*Q///////
ut:=#Ept/r;
ut:=3235378284218996041020628334555048353746687415791702110741362745607241012273289319490998100883360270140705444557817564333373970649196612709047960904504727376868836045879159893663749880485466097512959238393120957579174173722220769961965352109080756921414878821233494544942360241327904011276807478858565560475041657714859107418444304817769450537761391040286670975405474480702112256762237413877369520068658997811386234281044157320443969096913334657615996156910629987637321879968194559996451203622322855173490681315929340942845941901854762602482230345544778738233277695905149714849526040614153838451261090255301423115701413230650497436257358410673648969677575521082062660402719830010745131401857581241136771533959019632974838068321575021890239654250923239698058238532343293640583332746460549591118080473911339134214326805999615934889023975672403116797582756674727493991534767352065457123484256656059605341018889815246545798807124610587408576954994040926079696509810443325864338801147014393695521862315284030321439470300121888698228859049470995597746867551962112920666605249053365847645455641652839345532054083933046497019380937980801267829539781034425459803946933565477870059918126212656924303929710879629583716821028188063832088155868621647460338354760227739293275552535630628539912296348917701212830650737364550173593007332885398096588506291388885219345384480041310305023866955910802442221913227414055036258755703;
RR:=Random(Ept);
QRR:=ut*RR;
//////QRR is a point of order r and we send it to Ep27 on Fp27 to apply the ///Frobenius ///
x1:=QRR[1]*w^(-2);
y1:=QRR[2]*w^(-3);
xu:=x1^p;
yu:=y1^p;                                      
D1:=Ep27![xu,yu];
E1:=Ep27![x1,y1];
Q:=D1-E1;  /////// Q satisfies \pi(Q)=p*Q
Ep27![Q[1]^p,Q[2]^p]-p*Q;

///////The points P and Q are of order r
r*Q;
r*P;

/////The doubling step in the Miller Algorithm///////
Doubling:=function (S, K); 
x1:=S[1]*w^2; y1:=S[2]*w^3;
A:=3*x1^2; B:=(2*y1)^(-1); C:=A*B; x3:=C^2-2*x1; y3:=C*(x1-x3)-y1; f1:=K[2]*w^5+C*(x1*w^2-K[1]*w^4)-y1*w^2;
f2:=K[1]^2*w^4+x3*K[1]*w^2+x3^2;
f:=f1*f2;
return [x3*w^(-2),y3*w^(-3),f];
end function;

///////The Addition step in the Miller algorithm///////
Addition:=function(S,R, K);
x1:=S[1]*w^2; y1:=S[2]*w^3;
x2:=R[1]*w^2; y2:=R[2]*w^3;
A:=(y2-y1); B:=(x2-x1)^(-1); C:=A*B; x3:=C^2-x1-x2; 
y3:= C*(x1-x3)-y1; 
f1:=K[2]*w^5+C*(x1*w^2-K[1]*w^4)-y1*w^2;
f2:=K[1]^2*w^4+x3*K[1]*w^2+x3^2;
f:=f1*f2;
return [x3*w^(-2),y3*w^(-3),f];
end function;

///// ( addition and doubling of point) are correct/////
l:=Addition(Q,2*Q,P);
IsZero(l[2]^2-l[1]^3+2);
l1:=Doubling(Q,P);
IsZero(l1[2]^2-l1[1]^3+2);

///The Main function to compute the optimal Ate pairing described in the article//
AteOpt:=function(Q ,S);
f:=1;
bin:=Intseq(Abs(x),2);
 T:=Q;
for i:=#bin-2 to 0 by -1 do
h:=Doubling(T ,S);
T:=[h[1],h[2]];
f:=h[3]*f^2;
if bin[i+1] eq 1 then 
h:=Addition(T,Q ,S); 
T:=[h[1],h[2]];
f:=h[3]*f;
end if;
end for;
f:=f^(Integers()!((p^27-1)/r)); 
return(f);
end function;
e1:=AteOpt(Q,P);

//// The map  is not trivial/////
e1;

//////Verification of the bilinearity 
eP25:=AteOpt(Q,25*P);
eQ25:=AteOpt(25*Q,P);
IsZero(eP25-e1^25);
IsZero(eQ25-e1^25);