1 /********************************************************************************************
2 * Supersingular Isogeny Key Encapsulation Library
3 *
4 * Abstract: elliptic curve and isogeny functions
5 *********************************************************************************************/
6
7 #include "sikep434r3.h"
8 #include "sikep434r3_fpx.h"
9 #include "sikep434r3_ec_isogeny.h"
10
11 /* Doubling of a Montgomery point in projective coordinates (X:Z).
12 * Input: projective Montgomery x-coordinates P = (X1:Z1), where x1=X1/Z1 and Montgomery curve constants A+2C and 4C.
13 * Output: projective Montgomery x-coordinates Q = 2*P = (X2:Z2). */
xDBL(const point_proj_t P,point_proj_t Q,const f2elm_t * A24plus,const f2elm_t * C24)14 void xDBL(const point_proj_t P, point_proj_t Q, const f2elm_t *A24plus, const f2elm_t *C24)
15 {
16 f2elm_t _t0, _t1;
17 f2elm_t *t0=&_t0, *t1=&_t1;
18
19 mp2_sub_p2(&P->X, &P->Z, t0); /* t0 = X1-Z1 */
20 mp2_add(&P->X, &P->Z, t1); /* t1 = X1+Z1 */
21 fp2sqr_mont(t0, t0); /* t0 = (X1-Z1)^2 */
22 fp2sqr_mont(t1, t1); /* t1 = (X1+Z1)^2 */
23 fp2mul_mont(C24, t0, &Q->Z); /* Z2 = C24*(X1-Z1)^2 */
24 fp2mul_mont(t1, &Q->Z, &Q->X); /* X2 = C24*(X1-Z1)^2*(X1+Z1)^2 */
25 mp2_sub_p2(t1, t0, t1); /* t1 = (X1+Z1)^2-(X1-Z1)^2 */
26 fp2mul_mont(A24plus, t1, t0); /* t0 = A24plus*[(X1+Z1)^2-(X1-Z1)^2] */
27 mp2_add(&Q->Z, t0, &Q->Z); /* Z2 = A24plus*[(X1+Z1)^2-(X1-Z1)^2] + C24*(X1-Z1)^2 */
28 fp2mul_mont(&Q->Z, t1, &Q->Z); /* Z2 = [A24plus*[(X1+Z1)^2-(X1-Z1)^2] + C24*(X1-Z1)^2]*[(X1+Z1)^2-(X1-Z1)^2] */
29 }
30
31 /* Computes [2^e](X:Z) on Montgomery curve with projective constant via e repeated doublings.
32 * Input: projective Montgomery x-coordinates P = (XP:ZP), such that xP=XP/ZP and Montgomery curve constants A+2C and 4C.
33 * Output: projective Montgomery x-coordinates Q <- (2^e)*P. */
xDBLe(const point_proj_t P,point_proj_t Q,const f2elm_t * A24plus,const f2elm_t * C24,const int e)34 void xDBLe(const point_proj_t P, point_proj_t Q, const f2elm_t *A24plus, const f2elm_t *C24, const int e)
35 {
36 int i;
37
38 copy_words((const digit_t*)P, (digit_t*)Q, 2*2*S2N_SIKE_P434_R3_NWORDS_FIELD);
39
40 for (i = 0; i < e; i++) {
41 xDBL(Q, Q, A24plus, C24);
42 }
43 }
44
45 /* Computes the corresponding 4-isogeny of a projective Montgomery point (X4:Z4) of order 4.
46 * Input: projective point of order four P = (X4:Z4).
47 * Output: the 4-isogenous Montgomery curve with projective coefficients A+2C/4C and the 3 coefficients
48 * that are used to evaluate the isogeny at a point in eval_4_isog(). */
get_4_isog(const point_proj_t P,f2elm_t * A24plus,f2elm_t * C24,f2elm_t * coeff)49 void get_4_isog(const point_proj_t P, f2elm_t *A24plus, f2elm_t *C24, f2elm_t *coeff)
50 {
51 mp2_sub_p2(&P->X, &P->Z, &coeff[1]); /* coeff[1] = X4-Z4 */
52 mp2_add(&P->X, &P->Z, &coeff[2]); /* coeff[2] = X4+Z4 */
53 fp2sqr_mont(&P->Z, &coeff[0]); /* coeff[0] = Z4^2 */
54 mp2_add(&coeff[0], &coeff[0], &coeff[0]); /* coeff[0] = 2*Z4^2 */
55 fp2sqr_mont(&coeff[0], C24); /* C24 = 4*Z4^4 */
56 mp2_add(&coeff[0], &coeff[0], &coeff[0]); /* coeff[0] = 4*Z4^2 */
57 fp2sqr_mont(&P->X, A24plus); /* A24plus = X4^2 */
58 mp2_add(A24plus, A24plus, A24plus); /* A24plus = 2*X4^2 */
59 fp2sqr_mont(A24plus, A24plus); /* A24plus = 4*X4^4 */
60 }
61
62 /* Evaluates the isogeny at the point (X:Z) in the domain of the isogeny, given a 4-isogeny phi defined
63 * by the 3 coefficients in coeff (computed in the function get_4_isog()).
64 * Inputs: the coefficients defining the isogeny, and the projective point P = (X:Z).
65 * Output: the projective point P = phi(P) = (X:Z) in the codomain. */
eval_4_isog(point_proj_t P,f2elm_t * coeff)66 void eval_4_isog(point_proj_t P, f2elm_t *coeff)
67 {
68 f2elm_t _t0, _t1;
69 f2elm_t *t0=&_t0, *t1=&_t1;
70
71 mp2_add(&P->X, &P->Z, t0); /* t0 = X+Z */
72 mp2_sub_p2(&P->X, &P->Z, t1); /* t1 = X-Z */
73 fp2mul_mont(t0, &coeff[1], &P->X); /* X = (X+Z)*coeff[1] */
74 fp2mul_mont(t1, &coeff[2], &P->Z); /* Z = (X-Z)*coeff[2] */
75 fp2mul_mont(t0, t1, t0); /* t0 = (X+Z)*(X-Z) */
76 fp2mul_mont(&coeff[0], t0, t0); /* t0 = coeff[0]*(X+Z)*(X-Z) */
77 mp2_add(&P->X, &P->Z, t1); /* t1 = (X-Z)*coeff[2] + (X+Z)*coeff[1] */
78 mp2_sub_p2(&P->X, &P->Z, &P->Z); /* Z = (X-Z)*coeff[2] - (X+Z)*coeff[1] */
79 fp2sqr_mont(t1, t1); /* t1 = [(X-Z)*coeff[2] + (X+Z)*coeff[1]]^2 */
80 fp2sqr_mont(&P->Z, &P->Z); /* Z = [(X-Z)*coeff[2] - (X+Z)*coeff[1]]^2 */
81 mp2_add(t1, t0, &P->X); /* X = coeff[0]*(X+Z)*(X-Z) + [(X-Z)*coeff[2] + (X+Z)*coeff[1]]^2 */
82 mp2_sub_p2(&P->Z, t0, t0); /* t0 = [(X-Z)*coeff[2] - (X+Z)*coeff[1]]^2 - coeff[0]*(X+Z)*(X-Z) */
83 fp2mul_mont(&P->X, t1, &P->X); /* Xfinal */
84 fp2mul_mont(&P->Z, t0, &P->Z); /* Zfinal */
85 }
86
87 /* Tripling of a Montgomery point in projective coordinates (X:Z).
88 * Input: projective Montgomery x-coordinates P = (X:Z), where x=X/Z and Montgomery curve constants A24plus = A+2C and A24minus = A-2C.
89 * Output: projective Montgomery x-coordinates Q = 3*P = (X3:Z3). */
xTPL(const point_proj_t P,point_proj_t Q,const f2elm_t * A24minus,const f2elm_t * A24plus)90 void xTPL(const point_proj_t P, point_proj_t Q, const f2elm_t *A24minus, const f2elm_t *A24plus)
91 {
92 f2elm_t _t0, _t1, _t2, _t3, _t4, _t5, _t6;
93 f2elm_t *t0=&_t0, *t1=&_t1, *t2=&_t2, *t3=&_t3, *t4=&_t4, *t5=&_t5, *t6=&_t6;
94
95 mp2_sub_p2(&P->X, &P->Z, t0); /* t0 = X-Z */
96 fp2sqr_mont(t0, t2); /* t2 = (X-Z)^2 */
97 mp2_add(&P->X, &P->Z, t1); /* t1 = X+Z */
98 fp2sqr_mont(t1, t3); /* t3 = (X+Z)^2 */
99 mp2_add(&P->X, &P->X, t4); /* t4 = 2*X */
100 mp2_add(&P->Z, &P->Z, t0); /* t0 = 2*Z */
101 fp2sqr_mont(t4, t1); /* t1 = 4*X^2 */
102 mp2_sub_p2(t1, t3, t1); /* t1 = 4*X^2 - (X+Z)^2 */
103 mp2_sub_p2(t1, t2, t1); /* t1 = 4*X^2 - (X+Z)^2 - (X-Z)^2 */
104 fp2mul_mont(A24plus, t3, t5); /* t5 = A24plus*(X+Z)^2 */
105 fp2mul_mont(t3, t5, t3); /* t3 = A24plus*(X+Z)^4 */
106 fp2mul_mont(A24minus, t2, t6); /* t6 = A24minus*(X-Z)^2 */
107 fp2mul_mont(t2, t6, t2); /* t2 = A24minus*(X-Z)^4 */
108 mp2_sub_p2(t2, t3, t3); /* t3 = A24minus*(X-Z)^4 - A24plus*(X+Z)^4 */
109 mp2_sub_p2(t5, t6, t2); /* t2 = A24plus*(X+Z)^2 - A24minus*(X-Z)^2 */
110 fp2mul_mont(t1, t2, t1); /* t1 = [4*X^2 - (X+Z)^2 - (X-Z)^2]*[A24plus*(X+Z)^2 - A24minus*(X-Z)^2] */
111 fp2add(t3, t1, t2); /* t2 = [4*X^2 - (X+Z)^2 - (X-Z)^2]*[A24plus*(X+Z)^2 - A24minus*(X-Z)^2] + A24minus*(X-Z)^4 - A24plus*(X+Z)^4 */
112 fp2sqr_mont(t2, t2); /* t2 = t2^2 */
113 fp2mul_mont(t4, t2, &Q->X); /* X3 = 2*X*t2 */
114 fp2sub(t3, t1, t1); /* t1 = A24minus*(X-Z)^4 - A24plus*(X+Z)^4 - [4*X^2 - (X+Z)^2 - (X-Z)^2]*[A24plus*(X+Z)^2 - A24minus*(X-Z)^2] */
115 fp2sqr_mont(t1, t1); /* t1 = t1^2 */
116 fp2mul_mont(t0, t1, &Q->Z); /* Z3 = 2*Z*t1 */
117 }
118
119 /* Computes [3^e](X:Z) on Montgomery curve with projective constant via e repeated triplings.
120 * Input: projective Montgomery x-coordinates P = (XP:ZP), such that xP=XP/ZP and Montgomery curve constants A24plus = A+2C and A24minus = A-2C.
121 * Output: projective Montgomery x-coordinates Q <- (3^e)*P. */
xTPLe(const point_proj_t P,point_proj_t Q,const f2elm_t * A24minus,const f2elm_t * A24plus,const int e)122 void xTPLe(const point_proj_t P, point_proj_t Q, const f2elm_t *A24minus, const f2elm_t *A24plus, const int e)
123 {
124 int i;
125
126 copy_words((const digit_t*)P, (digit_t*)Q, 2*2*S2N_SIKE_P434_R3_NWORDS_FIELD);
127
128 for (i = 0; i < e; i++) {
129 xTPL(Q, Q, A24minus, A24plus);
130 }
131 }
132
133 /* Computes the corresponding 3-isogeny of a projective Montgomery point (X3:Z3) of order 3.
134 * Input: projective point of order three P = (X3:Z3).
135 * Output: the 3-isogenous Montgomery curve with projective coefficient A/C. */
get_3_isog(const point_proj_t P,f2elm_t * A24minus,f2elm_t * A24plus,f2elm_t * coeff)136 void get_3_isog(const point_proj_t P, f2elm_t *A24minus, f2elm_t *A24plus, f2elm_t *coeff)
137 {
138 f2elm_t _t0, _t1, _t2, _t3, _t4;
139 f2elm_t *t0=&_t0, *t1=&_t1, *t2=&_t2, *t3=&_t3, *t4=&_t4;
140
141 mp2_sub_p2(&P->X, &P->Z, &coeff[0]); /* coeff0 = X-Z */
142 fp2sqr_mont(&coeff[0], t0); /* t0 = (X-Z)^2 */
143 mp2_add(&P->X, &P->Z, &coeff[1]); /* coeff1 = X+Z */
144 fp2sqr_mont(&coeff[1], t1); /* t1 = (X+Z)^2 */
145 mp2_add(&P->X, &P->X, t3); /* t3 = 2*X */
146 fp2sqr_mont(t3, t3); /* t3 = 4*X^2 */
147 fp2sub(t3, t0, t2); /* t2 = 4*X^2 - (X-Z)^2 */
148 fp2sub(t3, t1, t3); /* t3 = 4*X^2 - (X+Z)^2 */
149 mp2_add(t0, t3, t4); /* t4 = 4*X^2 - (X+Z)^2 + (X-Z)^2 */
150 mp2_add(t4, t4, t4); /* t4 = 2(4*X^2 - (X+Z)^2 + (X-Z)^2) */
151 mp2_add(t1, t4, t4); /* t4 = 8*X^2 - (X+Z)^2 + 2*(X-Z)^2 */
152 fp2mul_mont(t2, t4, A24minus); /* A24minus = [4*X^2 - (X-Z)^2]*[8*X^2 - (X+Z)^2 + 2*(X-Z)^2] */
153 mp2_add(t1, t2, t4); /* t4 = 4*X^2 + (X+Z)^2 - (X-Z)^2 */
154 mp2_add(t4, t4, t4); /* t4 = 2(4*X^2 + (X+Z)^2 - (X-Z)^2) */
155 mp2_add(t0, t4, t4); /* t4 = 8*X^2 + 2*(X+Z)^2 - (X-Z)^2 */
156 fp2mul_mont(t3, t4, A24plus); /* A24plus = [4*X^2 - (X+Z)^2]*[8*X^2 + 2*(X+Z)^2 - (X-Z)^2] */
157 }
158
159 /* Computes the 3-isogeny R=phi(X:Z), given projective point (X3:Z3) of order 3 on a Montgomery curve and
160 * a point P with 2 coefficients in coeff (computed in the function get_3_isog()).
161 * Inputs: projective points P = (X3:Z3) and Q = (X:Z).
162 * Output: the projective point Q <- phi(Q) = (X3:Z3). */
eval_3_isog(point_proj_t Q,const f2elm_t * coeff)163 void eval_3_isog(point_proj_t Q, const f2elm_t *coeff)
164 {
165 f2elm_t _t0, _t1, _t2;
166 f2elm_t *t0=&_t0, *t1=&_t1, *t2=&_t2;
167
168 mp2_add(&Q->X, &Q->Z, t0); /* t0 = X+Z */
169 mp2_sub_p2(&Q->X, &Q->Z, t1); /* t1 = X-Z */
170 fp2mul_mont(&coeff[0], t0, t0); /* t0 = coeff0*(X+Z) */
171 fp2mul_mont(&coeff[1], t1, t1); /* t1 = coeff1*(X-Z) */
172 mp2_add(t0, t1, t2); /* t2 = coeff0*(X+Z) + coeff1*(X-Z) */
173 mp2_sub_p2(t1, t0, t0); /* t0 = coeff1*(X-Z) - coeff0*(X+Z) */
174 fp2sqr_mont(t2, t2); /* t2 = [coeff0*(X+Z) + coeff1*(X-Z)]^2 */
175 fp2sqr_mont(t0, t0); /* t0 = [coeff1*(X-Z) - coeff0*(X+Z)]^2 */
176 fp2mul_mont(&Q->X, t2, &Q->X); /* X3final = X*[coeff0*(X+Z) + coeff1*(X-Z)]^2 */
177 fp2mul_mont(&Q->Z, t0, &Q->Z); /* Z3final = Z*[coeff1*(X-Z) - coeff0*(X+Z)]^2 */
178 }
179
180 /* 3-way simultaneous inversion
181 * Input: z1,z2,z3
182 * Output: 1/z1,1/z2,1/z3 (override inputs). */
inv_3_way(f2elm_t * z1,f2elm_t * z2,f2elm_t * z3)183 void inv_3_way(f2elm_t *z1, f2elm_t *z2, f2elm_t *z3)
184 {
185 f2elm_t _t0, _t1, _t2, _t3;
186 f2elm_t *t0=&_t0, *t1=&_t1, *t2=&_t2, *t3=&_t3;
187
188 fp2mul_mont(z1, z2, t0); /* t0 = z1*z2 */
189 fp2mul_mont(z3, t0, t1); /* t1 = z1*z2*z3 */
190 fp2inv_mont(t1); /* t1 = 1/(z1*z2*z3) */
191 fp2mul_mont(z3, t1, t2); /* t2 = 1/(z1*z2) */
192 fp2mul_mont(t2, z2, t3); /* t3 = 1/z1 */
193 fp2mul_mont(t2, z1, z2); /* z2 = 1/z2 */
194 fp2mul_mont(t0, t1, z3); /* z3 = 1/z3 */
195 fp2copy(t3, z1); /* z1 = 1/z1 */
196 }
197
198 /* Given the x-coordinates of P, Q, and R, returns the value A corresponding to the
199 * Montgomery curve E_A: y^2=x^3+A*x^2+x such that R=Q-P on E_A.
200 * Input: the x-coordinates xP, xQ, and xR of the points P, Q and R.
201 * Output: the coefficient A corresponding to the curve E_A: y^2=x^3+A*x^2+x. */
get_A(const f2elm_t * xP,const f2elm_t * xQ,const f2elm_t * xR,f2elm_t * A)202 void get_A(const f2elm_t *xP, const f2elm_t *xQ, const f2elm_t *xR, f2elm_t *A)
203 {
204 f2elm_t _t0, _t1, one = {0};
205 f2elm_t *t0=&_t0, *t1=&_t1;
206
207
208 fpcopy((const digit_t*)&Montgomery_one,one.e[0]);
209 fp2add(xP, xQ, t1); /* t1 = xP+xQ */
210 fp2mul_mont(xP, xQ, t0); /* t0 = xP*xQ */
211 fp2mul_mont(xR, t1, A); /* A = xR*t1 */
212 fp2add(t0, A, A); /* A = A+t0 */
213 fp2mul_mont(t0, xR, t0); /* t0 = t0*xR */
214 fp2sub(A, &one, A); /* A = A-1 */
215 fp2add(t0, t0, t0); /* t0 = t0+t0 */
216 fp2add(t1, xR, t1); /* t1 = t1+xR */
217 fp2add(t0, t0, t0); /* t0 = t0+t0 */
218 fp2sqr_mont(A, A); /* A = A^2 */
219 fp2inv_mont(t0); /* t0 = 1/t0 */
220 fp2mul_mont(A, t0, A); /* A = A*t0 */
221 fp2sub(A, t1, A); /* Afinal = A-t1 */
222 }
223
224 /* Computes the j-invariant of a Montgomery curve with projective constant.
225 * Input: A,C in GF(p^2).
226 * Output: j=256*(A^2-3*C^2)^3/(C^4*(A^2-4*C^2)), which is the j-invariant of the Montgomery curve
227 * B*y^2=x^3+(A/C)*x^2+x or (equivalently) j-invariant of B'*y^2=C*x^3+A*x^2+C*x. */
j_inv(const f2elm_t * A,const f2elm_t * C,f2elm_t * jinv)228 void j_inv(const f2elm_t *A, const f2elm_t *C, f2elm_t *jinv)
229 {
230 f2elm_t _t0, _t1;
231 f2elm_t *t0=&_t0, *t1=&_t1;
232
233 fp2sqr_mont(A, jinv); /* jinv = A^2 */
234 fp2sqr_mont(C, t1); /* t1 = C^2 */
235 fp2add(t1, t1, t0); /* t0 = t1+t1 */
236 fp2sub(jinv, t0, t0); /* t0 = jinv-t0 */
237 fp2sub(t0, t1, t0); /* t0 = t0-t1 */
238 fp2sub(t0, t1, jinv); /* jinv = t0-t1 */
239 fp2sqr_mont(t1, t1); /* t1 = t1^2 */
240 fp2mul_mont(jinv, t1, jinv); /* jinv = jinv*t1 */
241 fp2add(t0, t0, t0); /* t0 = t0+t0 */
242 fp2add(t0, t0, t0); /* t0 = t0+t0 */
243 fp2sqr_mont(t0, t1); /* t1 = t0^2 */
244 fp2mul_mont(t0, t1, t0); /* t0 = t0*t1 */
245 fp2add(t0, t0, t0); /* t0 = t0+t0 */
246 fp2add(t0, t0, t0); /* t0 = t0+t0 */
247 fp2inv_mont(jinv); /* jinv = 1/jinv */
248 fp2mul_mont(jinv, t0, jinv); /* jinv = t0*jinv */
249 }
250
251 /* Simultaneous doubling and differential addition.
252 * Input: projective Montgomery points P=(XP:ZP) and Q=(XQ:ZQ) such that xP=XP/ZP and xQ=XQ/ZQ,
253 * affine difference xPQ=x(P-Q) and Montgomery curve constant A24=(A+2)/4.
254 * Output: projective Montgomery points P <- 2*P = (X2P:Z2P) such that x(2P)=X2P/Z2P,
255 * and Q <- P+Q = (XQP:ZQP) such that = x(Q+P)=XQP/ZQP. */
xDBLADD(point_proj_t P,point_proj_t Q,const f2elm_t * xPQ,const f2elm_t * A24)256 static void xDBLADD(point_proj_t P, point_proj_t Q, const f2elm_t *xPQ, const f2elm_t *A24)
257 {
258 f2elm_t _t0, _t1, _t2;
259 f2elm_t *t0=&_t0, *t1=&_t1, *t2=&_t2;
260
261 mp2_add(&P->X, &P->Z, t0); /* t0 = XP+ZP */
262 mp2_sub_p2(&P->X, &P->Z, t1); /* t1 = XP-ZP */
263 fp2sqr_mont(t0, &P->X); /* XP = (XP+ZP)^2 */
264 mp2_sub_p2(&Q->X, &Q->Z, t2); /* t2 = XQ-ZQ */
265 mp2_add(&Q->X, &Q->Z, &Q->X); /* XQ = XQ+ZQ */
266 fp2mul_mont(t0, t2, t0); /* t0 = (XP+ZP)*(XQ-ZQ) */
267 fp2sqr_mont(t1, &P->Z); /* ZP = (XP-ZP)^2 */
268 fp2mul_mont(t1, &Q->X, t1); /* t1 = (XP-ZP)*(XQ+ZQ) */
269 mp2_sub_p2(&P->X, &P->Z, t2); /* t2 = (XP+ZP)^2-(XP-ZP)^2 */
270 fp2mul_mont(&P->X, &P->Z, &P->X); /* XP = (XP+ZP)^2*(XP-ZP)^2 */
271 fp2mul_mont(A24, t2, &Q->X); /* XQ = A24*[(XP+ZP)^2-(XP-ZP)^2] */
272 mp2_sub_p2(t0, t1, &Q->Z); /* ZQ = (XP+ZP)*(XQ-ZQ)-(XP-ZP)*(XQ+ZQ) */
273 mp2_add(&Q->X, &P->Z, &P->Z); /* ZP = A24*[(XP+ZP)^2-(XP-ZP)^2]+(XP-ZP)^2 */
274 mp2_add(t0, t1, &Q->X); /* XQ = (XP+ZP)*(XQ-ZQ)+(XP-ZP)*(XQ+ZQ) */
275 fp2mul_mont(&P->Z, t2, &P->Z); /* ZP = [A24*[(XP+ZP)^2-(XP-ZP)^2]+(XP-ZP)^2]*[(XP+ZP)^2-(XP-ZP)^2] */
276 fp2sqr_mont(&Q->Z, &Q->Z); /* ZQ = [(XP+ZP)*(XQ-ZQ)-(XP-ZP)*(XQ+ZQ)]^2 */
277 fp2sqr_mont(&Q->X, &Q->X); /* XQ = [(XP+ZP)*(XQ-ZQ)+(XP-ZP)*(XQ+ZQ)]^2 */
278 fp2mul_mont(&Q->Z, xPQ, &Q->Z); /* ZQ = xPQ*[(XP+ZP)*(XQ-ZQ)-(XP-ZP)*(XQ+ZQ)]^2 */
279 }
280
281 /* Swap points.
282 * If option = 0 then P <- P and Q <- Q, else if option = 0xFF...FF then P <- Q and Q <- P */
swap_points(point_proj_t P,point_proj_t Q,const digit_t option)283 static void swap_points(point_proj_t P, point_proj_t Q, const digit_t option)
284 {
285 unsigned int i;
286
287 for (i = 0; i < S2N_SIKE_P434_R3_NWORDS_FIELD; i++) {
288 digit_t temp = option & (P->X.e[0][i] ^ Q->X.e[0][i]);
289 P->X.e[0][i] = temp ^ P->X.e[0][i];
290 Q->X.e[0][i] = temp ^ Q->X.e[0][i];
291 temp = option & (P->X.e[1][i] ^ Q->X.e[1][i]);
292 P->X.e[1][i] = temp ^ P->X.e[1][i];
293 Q->X.e[1][i] = temp ^ Q->X.e[1][i];
294 temp = option & (P->Z.e[0][i] ^ Q->Z.e[0][i]);
295 P->Z.e[0][i] = temp ^ P->Z.e[0][i];
296 Q->Z.e[0][i] = temp ^ Q->Z.e[0][i];
297 temp = option & (P->Z.e[1][i] ^ Q->Z.e[1][i]);
298 P->Z.e[1][i] = temp ^ P->Z.e[1][i];
299 Q->Z.e[1][i] = temp ^ Q->Z.e[1][i];
300 }
301 }
302
LADDER3PT(const f2elm_t * xP,const f2elm_t * xQ,const f2elm_t * xPQ,const digit_t * m,const unsigned int AliceOrBob,point_proj_t R,const f2elm_t * A)303 void LADDER3PT(const f2elm_t *xP, const f2elm_t *xQ, const f2elm_t *xPQ, const digit_t* m,
304 const unsigned int AliceOrBob, point_proj_t R, const f2elm_t *A)
305 {
306 point_proj_t R0 = {0}, R2 = {0};
307 f2elm_t _A24 = {0};
308 f2elm_t *A24 = &_A24;
309 digit_t mask;
310 int i, nbits, swap, prevbit = 0;
311
312 if (AliceOrBob == S2N_SIKE_P434_R3_ALICE) {
313 nbits = S2N_SIKE_P434_R3_OALICE_BITS;
314 } else {
315 nbits = S2N_SIKE_P434_R3_OBOB_BITS - 1;
316 }
317
318 /* Initializing constant */
319 fpcopy((const digit_t*)&Montgomery_one, A24->e[0]);
320 mp2_add(A24, A24, A24);
321 mp2_add(A, A24, A24);
322 fp2div2(A24, A24);
323 fp2div2(A24, A24); /* A24 = (A+2)/4 */
324
325 /* Initializing points */
326 fp2copy(xQ, &R0->X);
327 fpcopy((const digit_t*)&Montgomery_one, (digit_t*)&R0->Z);
328 fp2copy(xPQ, &R2->X);
329 fpcopy((const digit_t*)&Montgomery_one, (digit_t*)&R2->Z);
330 fp2copy(xP, &R->X);
331 fpcopy((const digit_t*)&Montgomery_one, (digit_t*)&R->Z);
332 fpzero((digit_t*)(R->Z.e)[1]);
333
334 /* Main loop */
335 for (i = 0; i < nbits; i++) {
336 int bit = (m[i >> S2N_SIKE_P434_R3_LOG2RADIX] >> (i & (S2N_SIKE_P434_R3_RADIX-1))) & 1;
337 swap = bit ^ prevbit;
338 prevbit = bit;
339 mask = 0 - (digit_t)swap;
340
341 swap_points(R, R2, mask);
342 xDBLADD(R0, R2, &R->X, A24);
343 fp2mul_mont(&R2->X, &R->Z, &R2->X);
344 }
345 swap = 0 ^ prevbit;
346 mask = 0 - (digit_t)swap;
347 swap_points(R, R2, mask);
348 }
349