1(* Elliptic algebra functions: FEE format. 2 3 y^2 = x^3 + c x^2 + a x + b. 4 5 Montgomery: b = 0, a = 1; 6 Weierstrass: c = 0; 7 Atkin3: c = a = 0; 8 Atkin4: c = b = 0; 9 10 Parameters c, a, b, p must be global. 11 *) 12 13elleven[pt_] := Block[{x1 = pt[[1]], z1 = pt[[2]], e, f }, 14 e = Mod[(x1^2 - a z1^2)^2 - 4 b (2 x1 + c z1) z1^3, p]; 15 f = Mod[4 z1 (x1^3 + c x1^2 z1 + a x1 z1^2 + b z1^3), p]; 16 Return[{e,f}] 17]; 18 19ellodd[pt_, pu_, pv_] := Block[ 20 {x1 = pt[[1]], z1 = pt[[2]], 21 x2 = pu[[1]], z2 = pu[[2]], 22 xx = pv[[1]], zz = pv[[2]], i, j}, 23 i = Mod[zz ((x1 x2 - a z1 z2)^2 - 24 4 b(x1 z2 + x2 z1 + c z1 z2) z1 z2), p]; 25 j = Mod[xx (x1 z2 - x2 z1)^2, p]; 26 Return[{i,j}] 27]; 28 29bitList[k_] := Block[{li = {}, j = k}, 30 While[j > 0, 31 li = Append[li, Mod[j,2]]; 32 j = Floor[j/2]; 33 ]; 34 Return[Reverse[li]]; 35 ]; 36 37elliptic[pt_, k_] := Block[{porg, ps, pp, q}, 38 39 If[k ==1, Return[pt]]; 40 If[k ==2, Return[elleven[pt]]]; 41 porg = pt; 42 ps = elleven[pt]; 43 pp = pt; 44 bitlist = bitList[k]; 45 Do[ 46 If[bitlist[[q]] == 1, 47 pp = ellodd[ps, pp, porg]; 48 ps = elleven[ps], 49 ps = ellodd[pp, ps, porg]; 50 pp = elleven[pp] 51 ], 52 {q,2,Length[bitlist]} 53 ]; 54 Return[Mod[pp,p]] 55]; 56ellinv[n_] := PowerMod[n,-1,p]; 57ex[pt_] := Mod[pt[[1]] * ellinv[pt[[2]]], p]; 58