1
/*
2
*******************************************************************************
3
\file ecp.c
4
\brief Elliptic curves over prime fields
5
\project bee2 [cryptographic library]
6
\author (C) Sergey Agievich [agievich@{bsu.by|gmail.com}]
7
\created 2012.06.26
8
\version 2016.07.05
9
\license This program is released under the GNU General Public License 
10
version 3. See Copyright Notices in bee2/info.h.
11
*******************************************************************************
12
*/
13

14
#include "bee2/core/mem.h"
15
#include "bee2/core/stack.h"
16
#include "bee2/core/util.h"
17
#include "bee2/math/ecp.h"
18
#include "bee2/math/gfp.h"
19
#include "bee2/math/pri.h"
20
#include "bee2/math/ww.h"
21
#include "bee2/math/zz.h"
22

23
/*
24
*******************************************************************************
25
Общие положения
26

27
Ссылки на все реализованные алгоритмы имеются на сайте
28
http://www.hyperelliptic.org/efd. Там же можно найти соглашения по
29
обозначению сложности алгоритмов. В этих обозначениях фигурируют
30
следующие формальные выражения:
31
	add -- сложение или вычитание \mod p,
32
	c -- умножение на малую константу c (2,3,...) \mod p,
33
	half -- умножение на 2^{-1} \mod p,
34
	*A -- умножение на коэффициент A \mod p,
35
	S -- возведение в квадрат \mod p,
36
	M -- умножение \mod p,
37
	D -- деление \mod p.
38

39
При общей оценке сложности считается, что 1D = 100M и 1S = 1M.
40
Аддитивные операции игнорируются.
41

42
Используются обозначения:
43
	A <- A + A -- сложение аффинных точек,
44
	A <- 2A -- удвоение аффинной точки;
45
	P <- P + P -- сложение проективных точек;
46
	P <- P + A -- добавление к проективной точке аффинной;
47
	P <- 2P -- удвоение проективной точки;
48
	P <- 2A -- удвоение аффинной точки с переходом к проективным координатам.
49
*******************************************************************************
50
*/
51

52
#define ecpSeemsOnA(a, ec)\
53
	(zmIsIn(ecX(a), (ec)->f) && zmIsIn(ecY(a, (ec)->f->n), (ec)->f))
54

55
#define ecpSeemsOn3(a, ec)\
56
	(ecpSeemsOnA(a, ec) && zmIsIn(ecZ(a, (ec)->f->n), (ec)->f))
57

58
/*
59
*******************************************************************************
60
Якобиановы координаты:
61
	x = X / Z^2, y = Y / Z^3,
62
	-(X : Y : Z) = (X : -Y : Z).
63

64
В функции ecpDblJ() выполняется удвоение P <- 2P. Реализован
65
алгоритм dbl-1998-hnm [Hasegawa T., Nakajima J., Matsui M. A Practical
66
Implementation of Elliptic Curve Cryptosystems over GF(p) on a
67
16-bit Microcomputer. Public Key Cryptography, 1998].
68
Сложность алгоритма:
69
	3M + 6S + 1*A + 1half + 6add + 3*2 \approx 9M.
70

71
\todo Сравнить dbl-1998-hnm с dbl-2007-bl, сложность которого
72
	1M + 8S + 1A + 10add + 1*8 + 2*2 + 1*3.
73

74
В функции ecpDblJA3() выполняется удвоение P <- 2P для случая A = -3.
75
Реализован алгоритм dbl-1998-hnm2 [Hasegawa T., Nakajima J., Matsui M.
76
A Practical Implementation of Elliptic Curve Cryptosystems over GF(p)
77
on a 16-bit Microcomputer. Public Key Cryptography, 1998].
78
Сложность алгоритма:
79
	4M + 4S + 1*half + 7add + 3*2 \approx 8M.
80

81
\todo Сравнить dbl-1998-hnm2 с dbl-2001-b, сложность которого
82
	3M +5S + 8add + 1*4 + 2*8 + 1*3.
83

84
В функции ecpDblAJ() выполняется удвоение P <- 2A. Реализован алгоритм 
85
mdbl-2007-bl [Bernstein-Lange, 2007]. Сложность алгоритма:
86
	1M + 5S + 7add + 1*8 + 3*2 + 1*3 \approx 6M.
87

88
В функции ecpAddJ() выполняется сложение P <- P + P. Реализован алгоритм 
89
add-2007-bl [Bernstein-Lange, 2007]. Сложность алгоритма:
90
	11M + 5S + 9add + 4*2 \approx 16M.
91

92
\todo Сравнить add-2007-bl с add-1998-cmo-2, сложность которого
93
	12M + 4S + 6add + 1*2.
94

95
В функции ecpAddAJ() выполняется сложение P <- P + A.
96
Реализован алгоритм madd-2004-hmv [Hankerson D., Menezes A., Vanstone S.
97
Guide to Elliptic Curve Cryptography, Springer, 2004].
98
Сложность алгоритма:
99
	8M + 3S + 6add + 1*2 \approx 11M.
100

101
В функции ecpTplJ() выполняется утроение P <- 3P. Реализован алгоритм 
102
tpl-2007-bl [Bernstein-Lange, 2007]. Сложность алгоритма
103
	5M + 10S + 1*A + 15add + 2*4 + 1*6 + 1*8 + 1*16 + 1*3 \approx 15M.
104

105
В функции ecpTplJA3() выполняется утроение P <- 3P для случая A = -3. 
106
Реализован алгоритм tpl-2007-bl-2 [Bernstein-Lange, 2007]. 
107
Сложность алгоритма
108
	7M + 7S + 13add + 2*4 + 1*8 + 1*12 + 1*16 + 1*3 \approx 14M.
109

110
Целевые функции ci(l), определенные в описании реализации ecMul() в ec.c, 
111
принимают следующий вид:
112
	c1(l) = l/3 11;
113
	c2(l, w) = 103 + (2^{w-2} - 2)102 + l/(w + 1) 11;
114
	c3(l, w) = 6 + (2^{w-2} - 2)16 + l/(w + 1) 16.
115

116
Расчеты показывают, что
117
	с1(l) <= min_w c3(l), l <= 81,
118
	min_w c3(l, w) <= m	in_w c2(l, w), l <= 899.
119
Поэтому для практически используемых размерностей l (192 <= l <= 571)
120
первые две стратегии являются проигрышными. Реализована только стратегия 3.
121

122
\todo Сравнить madd-2004-hmv с madd-2007-bl, сложность которого
123
	7M + 4S + 9add + 1*4 + 3*2.
124

125
\todo При выполнении функции ecpAddAJ() может выясниться, что
126
складываются одинаковые точки. В этом случае вызывается функция
127
ecpDblAJ(). Подумать, как "помочь" этой функции, передав уже
128
рассчитанные значения, например, za^2. Аналогичное замечание касается
129
функции ecpAddJ().
130
*******************************************************************************
131
*/
132

133
// [3n]b <- [2n]a (P <- A)
134 1
static bool_t ecpFromAJ(word b[], const word a[], const ec_o* ec, void* stack)
135
{
136 1
	const size_t n = ec->f->n;
137
	// pre
138 1
	ASSERT(ecIsOperable(ec) && ec->d == 3);
139 1
	ASSERT(ecpSeemsOnA(a, ec));
140 1
	ASSERT(a == b || wwIsDisjoint2(a, 2 * n, b, 3 * n));
141
	// xb <- xa
142 1
	qrCopy(ecX(b), ecX(a), ec->f);
143
	// yb <- ya
144 1
	qrCopy(ecY(b, n), ecY(a, n), ec->f);
145
	// zb <- unity
146 1
	qrSetUnity(ecZ(b, n), ec->f);
147 1
	return TRUE;
148
}
149

150
// [2n]b <- [3n]a (A <- P)
151 1
static bool_t ecpToAJ(word b[], const word a[], const ec_o* ec, void* stack)
152
{
153 1
	const size_t n = ec->f->n;
154
	// переменные в stack
155 1
	word* t1 = (word*)stack;
156 1
	word* t2 = t1 + n;
157 1
	stack = t2 + n;
158
	// pre
159 1
	ASSERT(ecIsOperable(ec) && ec->d == 3);
160 1
	ASSERT(ecpSeemsOn3(a, ec));
161 1
	ASSERT(a == b || wwIsDisjoint2(a, 3 * n, b, 2 * n));
162
	// a == O => b <- O
163 1
	if (qrIsZero(ecZ(a, n), ec->f))
164 1
		return FALSE;
165
	// t1 <- za^{-1}
166 1
	qrInv(t1, ecZ(a, n), ec->f, stack);
167
	// t2 <- t1^2
168 1
	qrSqr(t2, t1, ec->f, stack);
169
	// xb <- xa t2
170 1
	qrMul(ecX(b), ecX(a), t2, ec->f, stack);
171
	// t2 <- t1 t2
172 1
	qrMul(t2, t1, t2, ec->f, stack);
173
	// yb <- ya t2
174 1
	qrMul(ecY(b, n), ecY(a, n), t2, ec->f, stack);
175
	// b != O
176 1
	return TRUE;
177
}
178

179 1
static size_t ecpToAJ_deep(size_t n, size_t f_deep)
180
{
181 1
	return O_OF_W(2 * n) + f_deep;
182
}
183

184
// [3n]b <- -[3n]a (P <- -P)
185 0
static void ecpNegJ(word b[], const word a[], const ec_o* ec, void* stack)
186
{
187 0
	const size_t n = ec->f->n;
188
	// pre
189 0
	ASSERT(ecIsOperable(ec) && ec->d == 3);
190 0
	ASSERT(ecpSeemsOn3(a, ec));
191 0
	ASSERT(wwIsSameOrDisjoint(a, b, 3 * n));
192
	// xb <- xa
193 0
	qrCopy(ecX(b), ecX(a), ec->f);
194
	// yb <- -ya
195 0
	zmNeg(ecY(b, n), ecY(a, n), ec->f);
196
	// zb <- za
197 0
	qrCopy(ecZ(b, n), ecZ(a, n), ec->f);
198
}
199

200
// [3n]b <- 2[3n]a (P <- 2P)
201 1
static void ecpDblJ(word b[], const word a[], const ec_o* ec, void* stack)
202
{
203 1
	const size_t n = ec->f->n;
204
	// переменные в stack
205 1
	word* t1 = (word*)stack;
206 1
	word* t2 = t1 + n;
207 1
	stack = t2 + n;
208
	// pre
209 1
	ASSERT(ecIsOperable(ec) && ec->d == 3);
210 1
	ASSERT(ecpSeemsOn3(a, ec));
211 1
	ASSERT(wwIsSameOrDisjoint(a, b, 3 * n));
212
	// za == 0 или ya == 0? => b <- O
213 1
	if (qrIsZero(ecZ(a, n), ec->f) || qrIsZero(ecY(a, n), ec->f))
214
	{
215 1
		qrSetZero(ecZ(b, n), ec->f);
216 1
		return;
217
	}
218
	// t1 <- za^2
219 1
	qrSqr(t1, ecZ(a, n), ec->f, stack);
220
	// zb <- ya za
221 1
	qrMul(ecZ(b, n), ecY(a, n), ecZ(a, n), ec->f, stack);
222
	// zb <- 2 zb
223 1
	gfpDouble(ecZ(b, n), ecZ(b, n), ec->f);
224
	// t1 <- t1^2
225 1
	qrSqr(t1, t1, ec->f, stack);
226
	// t1 <- A t1
227 1
	qrMul(t1, ec->A, t1, ec->f, stack);
228
	// t2 <- xa^2
229 1
	qrSqr(t2, ecX(a), ec->f, stack);
230
	// t1 <- t1 + t2
231 1
	zmAdd(t1, t1, t2, ec->f);
232
	// t2 <- 2 t2
233 1
	gfpDouble(t2, t2, ec->f);
234
	// t1 <- t1 + t2
235 1
	zmAdd(t1, t1, t2, ec->f);
236
	// yb <- 2 ya
237 1
	gfpDouble(ecY(b, n), ecY(a, n), ec->f);
238
	// yb <- yb^2
239 1
	qrSqr(ecY(b, n), ecY(b, n), ec->f, stack);
240
	// t2 <- yb^2
241 1
	qrSqr(t2, ecY(b, n), ec->f, stack);
242
	// t2 <- t2 / 2
243 1
	gfpHalf(t2, t2, ec->f);
244
	// yb <- yb xa
245 1
	qrMul(ecY(b, n), ecY(b, n), ecX(a), ec->f, stack);
246
	// xb <- t1^2
247 1
	qrSqr(ecX(b), t1, ec->f, stack);
248
	// xb <- xb - yb
249 1
	zmSub(ecX(b), ecX(b), ecY(b, n), ec->f);
250
	// xb <- xb - yb
251 1
	zmSub(ecX(b), ecX(b), ecY(b, n), ec->f);
252
	// yb <- yb - xb
253 1
	zmSub(ecY(b, n), ecY(b, n), ecX(b), ec->f);
254
	// yb <- yb t1
255 1
	qrMul(ecY(b, n), ecY(b, n), t1, ec->f, stack);
256
	// yb <- yb - t2
257 1
	zmSub(ecY(b, n), ecY(b, n), t2, ec->f);
258
}
259

260 1
static size_t ecpDblJ_deep(size_t n, size_t f_deep)
261
{
262 1
	return O_OF_W(2 * n) + f_deep;
263
}
264

265
// [3n]b <- 2[3n]a (P <- 2P, A = -3)
266 1
static void ecpDblJA3(word b[], const word a[], const ec_o* ec, void* stack)
267
{
268 1
	const size_t n = ec->f->n;
269
	// переменные в stack
270 1
	word* t1 = (word*)stack;
271 1
	word* t2 = t1 + n;
272 1
	stack = t2 + n;
273
	// pre
274 1
	ASSERT(ecIsOperable(ec) && ec->d == 3);
275 1
	ASSERT(ecpSeemsOn3(a, ec));
276 1
	ASSERT(wwIsSameOrDisjoint(a, b, 3 * n));
277
	// za == 0 или ya == 0? => b <- O
278 1
	if (qrIsZero(ecZ(a, n), ec->f) || qrIsZero(ecY(a, n), ec->f))
279
	{
280 1
		qrSetZero(ecZ(b, n), ec->f);
281 1
		return;
282
	}
283
	// t1 <- za^2
284 1
	qrSqr(t1, ecZ(a, n), ec->f, stack);
285
	// zb <- ya za
286 1
	qrMul(ecZ(b, n), ecY(a, n), ecZ(a, n), ec->f, stack);
287
	// zb <- 2 zb
288 1
	gfpDouble(ecZ(b, n), ecZ(b, n), ec->f);
289
	// t2 <- xa - t1
290 1
	zmSub(t2, ecX(a), t1, ec->f);
291
	// t1 <- xa + t1
292 1
	zmAdd(t1, ecX(a), t1, ec->f);
293
	// t2 <- t1 t2
294 1
	qrMul(t2, t1, t2, ec->f, stack);
295
	// t1 <- 2 t2
296 1
	gfpDouble(t1, t2, ec->f);
297
	// t1 <- t1 + t2
298 1
	zmAdd(t1, t1, t2, ec->f);
299
	// yb <- 2 ya
300 1
	gfpDouble(ecY(b, n), ecY(a, n), ec->f);
301
	// yb <- yb^2
302 1
	qrSqr(ecY(b, n), ecY(b, n), ec->f, stack);
303
	// t2 <- yb^2
304 1
	qrSqr(t2, ecY(b, n), ec->f, stack);
305
	// t2 <- t2 / 2
306 1
	gfpHalf(t2, t2, ec->f);
307
	// yb <- yb xa
308 1
	qrMul(ecY(b, n), ecY(b, n), ecX(a), ec->f, stack);
309
	// xb <- t1^2
310 1
	qrSqr(ecX(b), t1, ec->f, stack);
311
	// xb <- xb - yb
312 1
	zmSub(ecX(b), ecX(b), ecY(b, n), ec->f);
313
	// xb <- xb - yb
314 1
	zmSub(ecX(b), ecX(b), ecY(b, n), ec->f);
315
	// yb <- yb - xb
316 1
	zmSub(ecY(b, n), ecY(b, n), ecX(b), ec->f);
317
	// yb <- yb t1
318 1
	qrMul(ecY(b, n), ecY(b, n), t1, ec->f, stack);
319
	// yb <- yb - t2
320 1
	zmSub(ecY(b, n), ecY(b, n), t2, ec->f);
321
}
322

323 1
static size_t ecpDblJA3_deep(size_t n, size_t f_deep)
324
{
325 1
	return O_OF_W(2 * n) + f_deep;
326
}
327

328
// [3n]b <- 2[2n]a (P <- 2A)
329 1
static void ecpDblAJ(word b[], const word a[], const ec_o* ec, void* stack)
330
{
331 1
	const size_t n = ec->f->n;
332
	// переменные в stack
333 1
	word* t1 = (word*)stack;
334 1
	word* t2 = t1 + n;
335 1
	word* t3 = t2 + n;
336 1
	word* t4 = t3 + n;
337 1
	stack = t4 + n;
338
	// pre
339 1
	ASSERT(ecIsOperable(ec) && ec->d == 3);
340 1
	ASSERT(ecpSeemsOnA(a, ec));
341 1
	ASSERT(a == b || wwIsDisjoint2(a, 2 * n, b, 3 * n));
342
	// ya == 0? => b <- O
343 1
	if (qrIsZero(ecY(a, n), ec->f))
344
	{
345 0
		qrSetZero(ecZ(b, n), ec->f);
346 0
		return;
347
	}
348
	// t1 <- xa^2 [X1^2 = XX]
349 1
	qrSqr(t1, ecX(a), ec->f, stack);
350
	// t2 <- ya^2 [Y1^2 = YY]
351 1
	qrSqr(t2, ecY(a, n), ec->f, stack);
352
	// t3 <- t2^2 [YY^2 = YYYY]
353 1
	qrSqr(t3, t2, ec->f, stack);
354
	// t2 <- t2 + xa [X1 + YY]
355 1
	zmAdd(t2, t2, ecX(a), ec->f);
356
	// t2 <- t2^2 [(X1 + YY)^2]
357 1
	qrSqr(t2, t2, ec->f, stack);
358
	// t2 <- t2 - t1 [(X1 + YY)^2 - XX]
359 1
	zmSub(t2, t2, t1, ec->f);
360
	// t2 <- t2 - t3 [(X1 + YY)^2 - XX - YYYY]
361 1
	zmSub(t2, t2, t3, ec->f);
362
	// t2 <- 2 t2 [2((X1 + YY)^2 - XX - YYYY) = S]
363 1
	gfpDouble(t2, t2, ec->f);
364
	// t4 <- 2 t1 [2 XX]
365 1
	gfpDouble(t4, t1, ec->f);
366
	// t4 <- t4 + t1 [3 XX]
367 1
	zmAdd(t4, t4, t1, ec->f);
368
	// t4 <- t4 + A [3 XX + A = M]
369 1
	zmAdd(t4, t4, ec->A, ec->f);
370
	// t1 <- 2 t2 [2S]
371 1
	gfpDouble(t1, t2, ec->f);
372
	// xb <- t4^2 [M^2]
373 1
	qrSqr(ecX(b), t4, ec->f, stack);
374
	// xb <- xb - t1 [M^2 - 2S = T]
375 1
	zmSub(ecX(b), ecX(b), t1, ec->f);
376
	// zb <- 2 ya [2Y1]
377 1
	gfpDouble(ecZ(b, n), ecY(a, n), ec->f);
378
	// t2 <- t2 - xb [S - T]
379 1
	zmSub(t2, t2, ecX(b), ec->f);
380
	// yb <- t4 t2 [M(S - T)]
381 1
	qrMul(ecY(b, n), t4, t2, ec->f, stack);
382
	// t3 <- 2 t3 [2 YYYY]
383 1
	gfpDouble(t3, t3, ec->f);
384
	// t3 <- 2 t3 [4 YYYY]
385 1
	gfpDouble(t3, t3, ec->f);
386
	// t3 <- 2 t3 [8 YYYY]
387 1
	gfpDouble(t3, t3, ec->f);
388
	// yb <- yb - t3 [M(S - T) - 8 YYYY]
389 1
	zmSub(ecY(b, n), ecY(b, n), t3, ec->f);
390
}
391

392 1
static size_t ecpDblAJ_deep(size_t n, size_t f_deep)
393
{
394 1
	return O_OF_W(4 * n) + f_deep;
395
}
396

397
// [3n]c <- [3n]a + [3n]b (P <- P + P)
398 1
static void ecpAddJ(word c[], const word a[], const word b[], const ec_o* ec,
399
	void* stack)
400
{
401 1
	const size_t n = ec->f->n;
402
	// переменные в stack
403 1
	word* t1 = (word*)stack;
404 1
	word* t2 = t1 + n;
405 1
	word* t3 = t2 + n;
406 1
	word* t4 = t3 + n;
407 1
	stack = t4 + n;
408
	// pre
409 1
	ASSERT(ecIsOperable(ec) && ec->d == 3);
410 1
	ASSERT(ecpSeemsOn3(a, ec));
411 1
	ASSERT(ecpSeemsOn3(b, ec));
412 1
	ASSERT(wwIsSameOrDisjoint(a, c, 3 * n));
413 1
	ASSERT(wwIsSameOrDisjoint(b, c, 3 * n));
414
	// a == O => c <- b
415 1
	if (qrIsZero(ecZ(a, n), ec->f))
416
	{
417 1
		wwCopy(c, b, 3 * n);
418 1
		return;
419
	}
420
	// b == O => c <- a
421 1
	if (qrIsZero(ecZ(b, n), ec->f))
422
	{
423 0
		wwCopy(c, a, 3 * n);
424 0
		return;
425
	}
426
	// t1 <- za^2 [Z1Z1]
427 1
	qrSqr(t1, ecZ(a, n), ec->f, stack);
428
	// t2 <- zb^2 [Z2Z2]
429 1
	qrSqr(t2, ecZ(b, n), ec->f, stack);
430
	// t3 <- zb t2 [Z2 Z2Z2]
431 1
	qrMul(t3, ecZ(b, n), t2, ec->f, stack);
432
	// t3 <- ya t3 [Y1 Z2 Z2Z2 = S1]
433 1
	qrMul(t3, ecY(a, n), t3, ec->f, stack);
434
	// t4 <- za t1 [Z1 Z1Z1]
435 1
	qrMul(t4, ecZ(a, n), t1, ec->f, stack);
436
	// t4 <- yb t4 [Y2 Z1 Z1Z1 = S2]
437 1
	qrMul(t4, ecY(b, n), t4, ec->f, stack);
438
	// zc <- za + zb [Z1 + Z2]
439 1
	zmAdd(ecZ(c, n), ecZ(a, n), ecZ(b, n), ec->f);
440
	// zc <- zc^2 [(Z1 + Z2)^2]
441 1
	qrSqr(ecZ(c, n), ecZ(c, n), ec->f, stack);
442
	// zc <- zc - t1 [(Z1 + Z2)^2 - Z1Z1]
443 1
	zmSub(ecZ(c, n), ecZ(c, n), t1, ec->f);
444
	// zc <- zc - t2 [(Z1 + Z2)^2 - Z1Z1 - Z2Z2]
445 1
	zmSub(ecZ(c, n), ecZ(c, n), t2, ec->f);
446
	// t1 <- xb t1 [X1 Z2Z2 = U2]
447 1
	qrMul(t1, ecX(b), t1, ec->f, stack);
448
	// t2 <- xa t2 [X2 Z1Z1 = U1]
449 1
	qrMul(t2, ecX(a), t2, ec->f, stack);
450
	// t1 <- t1 - t2 [U2 - U1 = H]
451 1
	zmSub(t1, t1, t2, ec->f);
452
	// t1 == 0 => xa zb^2 == xb za^2
453 1
	if (qrIsZero(t1, ec->f))
454
	{
455
		// t3 == t4 => ya zb^3 == yb za^3 => a == b => c <- 2a
456 1
		if (qrCmp(t3, t4, ec->f) == 0)
457 0
			ecpDblJ(c, c == a ? b : a, ec, stack);
458
		// t3 != t4 => a == -b => c <- O
459
		else
460 1
			qrSetZero(ecZ(c, n), ec->f);
461 1
		return;
462
	}
463
	// zc <- zc t1 [((Z1 + Z2)^2 - Z1Z1 - Z2Z2)H = Z3]
464 1
	qrMul(ecZ(c, n), ecZ(c, n), t1, ec->f, stack);
465
	// t4 <- t4 - t3 [S2 - S1]
466 1
	zmSub(t4, t4, t3, ec->f);
467
	// t4 <- 2 t4 [2(S2 - S1) = r]
468 1
	gfpDouble(t4, t4, ec->f);
469
	// yc <- 2 t1 [2H]
470 1
	gfpDouble(ecY(c, n), t1, ec->f);
471
	// yc <- yc^2 [(2H)^2 = I]
472 1
	qrSqr(ecY(c, n), ecY(c, n), ec->f, stack);
473
	// t1 <- t1 yc [H I = J]
474 1
	qrMul(t1, t1, ecY(c, n), ec->f, stack);
475
	// yc <- t2 yc [U1 I = V]
476 1
	qrMul(ecY(c, n), t2, ecY(c, n), ec->f, stack);
477
	// t2 <- 2 yc [2 V]
478 1
	gfpDouble(t2, ecY(c, n), ec->f);
479
	// xc <- t4^2 [r^2]
480 1
	qrSqr(ecX(c), t4, ec->f, stack);
481
	// xc <- xc - t1 [r^2 - J]
482 1
	zmSub(ecX(c), ecX(c), t1, ec->f);
483
	// xc <- xc - t2 [r^2 - J - 2V = X3]
484 1
	zmSub(ecX(c), ecX(c), t2, ec->f);
485
	// yc <- yc - xc [V - X3]
486 1
	zmSub(ecY(c, n), ecY(c, n), ecX(c), ec->f);
487
	// yc <- t4 yc [r(V - X3)]
488 1
	qrMul(ecY(c, n), t4, ecY(c, n), ec->f, stack);
489
	// t3 <- 2 t3 [2S1]
490 1
	gfpDouble(t3, t3, ec->f);
491
	// t3 <- t3 t1 [2S1 J]
492 1
	qrMul(t3, t3, t1, ec->f, stack);
493
	// yc <- yc - t3 [r(V - X3) - 2 S1 J]
494 1
	zmSub(ecY(c, n), ecY(c, n), t3, ec->f);
495
}
496

497 1
static size_t ecpAddJ_deep(size_t n, size_t f_deep)
498
{
499 1
	return O_OF_W(4 * n) + 
500 1
		utilMax(2,
501
			f_deep,
502
			ecpDblJ_deep(n, f_deep));
503
}
504

505
// [3n]c <- [3n]a + [2n]b (P <- P + A)
506 1
static void ecpAddAJ(word c[], const word a[], const word b[], const ec_o* ec,
507
	void* stack)
508
{
509 1
	const size_t n = ec->f->n;
510
	// переменные в stack
511 1
	word* t1 = (word*)stack;
512 1
	word* t2 = t1 + n;
513 1
	word* t3 = t2 + n;
514 1
	word* t4 = t3 + n;
515 1
	stack = t4 + n;
516
	// pre
517 1
	ASSERT(ecIsOperable(ec) && ec->d == 3);
518 1
	ASSERT(ecpSeemsOn3(a, ec));
519 1
	ASSERT(ecpSeemsOnA(b, ec));
520 1
	ASSERT(wwIsSameOrDisjoint(a,  c, 3 * n));
521 1
	ASSERT(b == c || wwIsDisjoint2(b, 2 * n, c, 3 * n));
522
	// a == O => c <- (xb : yb : 1)
523 1
	if (qrIsZero(ecZ(a, n), ec->f))
524
	{
525 1
		qrCopy(ecX(c), ecX(b), ec->f);
526 1
		qrCopy(ecY(c, n), ecY(b, n), ec->f);
527 1
		qrSetUnity(ecZ(c, n), ec->f);
528 1
		return;
529
	}
530
	// t1 <- za^2
531 1
	qrSqr(t1, ecZ(a, n), ec->f, stack);
532
	// t2 <- t1 za
533 1
	qrMul(t2, t1, ecZ(a, n), ec->f, stack);
534
	// t1 <- t1 xb
535 1
	qrMul(t1, t1, ecX(b), ec->f, stack);
536
	// t2 <- t2 yb
537 1
	qrMul(t2, t2, ecY(b, n), ec->f, stack);
538
	// t1 <- t1 - xa
539 1
	zmSub(t1, t1, ecX(a), ec->f);
540
	// t2 <- t2 - ya
541 1
	zmSub(t2, t2, ecY(a, n), ec->f);
542
	// t1 == 0?
543 1
	if (qrIsZero(t1, ec->f))
544
	{
545
		// t2 == 0 => c <- 2(xb : yb : 1)
546 0
		if (qrIsZero(t2, ec->f))
547 0
			ecpDblAJ(c, b, ec, stack);
548
		// t2 != 0 => c <- O
549
		else
550 0
			qrSetZero(ecZ(c, n), ec->f);
551 0
		return;
552
	}
553
	// zc <- t1 za
554 1
	qrMul(ecZ(c, n), t1, ecZ(a, n), ec->f, stack);
555
	// t3 <- t1^2
556 1
	qrSqr(t3, t1, ec->f, stack);
557
	// t4 <- t1 t3
558 1
	qrMul(t4, t1, t3, ec->f, stack);
559
	// t3 <- t3 xa
560 1
	qrMul(t3, t3, ecX(a), ec->f, stack);
561
	// t1 <- 2 t3
562 1
	gfpDouble(t1, t3, ec->f);
563
	// xc <- t2^2
564 1
	qrSqr(ecX(c), t2, ec->f, stack);
565
	// xc <- xc - t1
566 1
	zmSub(ecX(c), ecX(c), t1, ec->f);
567
	// xc <- xc - t4
568 1
	zmSub(ecX(c), ecX(c), t4, ec->f);
569
	// t3 <- t3 - xc
570 1
	zmSub(t3, t3, ecX(c), ec->f);
571
	// t3 <- t3 t2
572 1
	qrMul(t3, t3, t2, ec->f, stack);
573
	// t4 <- t4 ya
574 1
	qrMul(t4, t4, ecY(a, n), ec->f, stack);
575
	// yc <- t3 - t4
576 1
	zmSub(ecY(c, n), t3, t4, ec->f);
577
}
578

579 1
static size_t ecpAddAJ_deep(size_t n, size_t f_deep)
580
{
581 1
	return O_OF_W(4 * n) +
582 1
		utilMax(2,
583
			f_deep,
584
			ecpDblAJ_deep(n, f_deep));
585
}
586

587
// [3n]c <- [3n]a - [3n]b (P <- P - P)
588 1
void ecpSubJ(word c[], const word a[], const word b[], const ec_o* ec,
589
	void* stack)
590
{
591 1
	const size_t n = ec->f->n;
592
	// переменные в stack
593 1
	word* t = (word*)stack;
594 1
	stack = t + 3 * n;
595
	// pre
596 1
	ASSERT(ecIsOperable(ec) && ec->d == 3);
597 1
	ASSERT(ecpSeemsOn3(a, ec));
598 1
	ASSERT(ecpSeemsOn3(b, ec));
599 1
	ASSERT(wwIsSameOrDisjoint(a, c, 3 * n));
600 1
	ASSERT(wwIsSameOrDisjoint(b, c, 3 * n));
601
	// t <- -b
602 1
	qrCopy(ecX(t), ecX(b), ec->f);
603 1
	zmNeg(ecY(t, n), ecY(b, n), ec->f);
604 1
	qrCopy(ecZ(t, n), ecZ(b, n), ec->f);
605
	// c <- a + t
606 1
	ecpAddJ(c, a, t, ec, stack);
607
}
608

609 1
size_t ecpSubJ_deep(size_t n, size_t f_deep)
610
{
611 1
	return O_OF_W(3 * n) + ecpAddJ_deep(n, f_deep);
612
}
613

614
// [3n]c <- [3n]a - [2n]b (P <- P - A)
615 1
static void ecpSubAJ(word c[], const word a[], const word b[], const ec_o* ec,
616
	void* stack)
617
{
618 1
	const size_t n = ec->f->n;
619
	// переменные в stack
620 1
	word* t = (word*)stack;
621 1
	stack = t + 2 * n;
622
	// pre
623 1
	ASSERT(ecIsOperable(ec) && ec->d == 3);
624 1
	ASSERT(ecpSeemsOn3(a, ec));
625 1
	ASSERT(ecpSeemsOnA(b, ec));
626 1
	ASSERT(wwIsSameOrDisjoint(a,  c, 3 * n));
627 1
	ASSERT(b == c || wwIsDisjoint2(b, 2 * n, c, 3 * n));
628
	// t <- -b
629 1
	qrCopy(ecX(t), ecX(b), ec->f);
630 1
	zmNeg(ecY(t, n), ecY(b, n), ec->f);
631
	// c <- a + t
632 1
	ecpAddAJ(c, a, t, ec, stack);
633
}
634

635 1
static size_t ecpSubAJ_deep(size_t n, size_t f_deep)
636
{
637 1
	return O_OF_W(2 * n) + ecpAddAJ_deep(n, f_deep);
638
}
639

640
// [3n]b <- 3[3n]a (P <- 3P)
641 0
static void ecpTplJ(word b[], const word a[], const ec_o* ec, void* stack)
642
{
643 0
	const size_t n = ec->f->n;
644
	// переменные в stack
645 0
	word* t0 = (word*)stack;
646 0
	word* t1 = t0 + n;
647 0
	word* t2 = t1 + n;
648 0
	word* t3 = t2 + n;
649 0
	word* t4 = t3 + n;
650 0
	word* t5 = t4 + n;
651 0
	word* t6 = t5 + n;
652 0
	word* t7 = t6 + n;
653 0
	stack = t7 + n;
654
	// pre
655 0
	ASSERT(ecIsOperable(ec) && ec->d == 3);
656 0
	ASSERT(ecpSeemsOn3(a, ec));
657 0
	ASSERT(wwIsSameOrDisjoint(a, b, 3 * n));
658
	// t0 <- xa^2 [XX]
659 0
	qrSqr(t0, ecX(a), ec->f, stack);
660
	// t1 <- ya^2 [YY]
661 0
	qrSqr(t1, ecY(a, n), ec->f, stack);
662
	// t2 <- za^2 [ZZ]
663 0
	qrSqr(t2, ecZ(a, n), ec->f, stack);
664
	// t3 <- t1^2 [YYYY]
665 0
	qrSqr(t3, t1, ec->f, stack);
666
	// t4 <- 3 t0 + A t2^2 [M]
667 0
	qrSqr(t4, t2, ec->f, stack);
668 0
	qrMul(t4, t4, ec->A, ec->f, stack);
669 0
	gfpDouble(t5, t0, ec->f);
670 0
	zmAdd(t5, t0, t5, ec->f);
671 0
	zmAdd(t4, t4, t5, ec->f);
672
	// t5 <- t4^2 [MM]
673 0
	qrSqr(t5, t4, ec->f, stack);
674
	// t6 <- 6((xa + t1)^2 - t0 - t3) - t5 [E]
675 0
	zmAdd(t6, ecX(a), t1, ec->f);
676 0
	qrSqr(t6, t6, ec->f, stack);
677 0
	zmSub(t6, t6, t0, ec->f);
678 0
	zmSub(t6, t6, t3, ec->f);
679 0
	gfpDouble(t7, t6, ec->f);
680 0
	zmAdd(t6, t6, t7, ec->f);
681 0
	gfpDouble(t6, t6, ec->f);
682 0
	zmSub(t6, t6, t5, ec->f);
683
	// t7 <- t6^2 [EE]
684 0
	qrSqr(t7, t6, ec->f, stack);
685
	// t3 <- 16 t3 [T]
686 0
	gfpDouble(t3, t3, ec->f);
687 0
	gfpDouble(t3, t3, ec->f);
688 0
	gfpDouble(t3, t3, ec->f);
689 0
	gfpDouble(t3, t3, ec->f);
690
	// zb <- (za + t6)^2 - t2 - t7
691 0
	zmAdd(ecZ(b, n), ecZ(a, n), t6, ec->f);
692 0
	qrSqr(ecZ(b, n), ecZ(b, n), ec->f, stack);
693 0
	zmSub(ecZ(b, n), ecZ(b, n), t2, ec->f);
694 0
	zmSub(ecZ(b, n), ecZ(b, n), t7, ec->f);
695
	// t2 <- (t4 + t6)^2 - t5 - t7 - t3 [U] 
696 0
	zmAdd(t2, t4, t6, ec->f);
697 0
	qrSqr(t2, t2, ec->f, stack);
698 0
	zmSub(t2, t2, t5, ec->f);
699 0
	zmSub(t2, t2, t7, ec->f);
700 0
	zmSub(t2, t2, t3, ec->f);
701
	// yb <- 8 ya (t2(t3 - t2) - t6 t7)
702 0
	zmSub(t3, t3, t2, ec->f);
703 0
	qrMul(t3, t2, t3, ec->f, stack);
704 0
	qrMul(t6, t6, t7, ec->f, stack);
705 0
	zmSub(t3, t3, t6, ec->f);
706 0
	qrMul(ecY(b, n), ecY(a, n), t3, ec->f, stack);
707 0
	gfpDouble(ecY(b, n), ecY(b, n), ec->f);
708 0
	gfpDouble(ecY(b, n), ecY(b, n), ec->f);
709 0
	gfpDouble(ecY(b, n), ecY(b, n), ec->f);
710
	// xb <- 4 (xa t7 - 4 t1 t2)
711 0
	qrMul(t1, t1, t2, ec->f, stack);
712 0
	gfpDouble(t1, t1, ec->f);
713 0
	gfpDouble(t1, t1, ec->f);
714 0
	qrMul(ecX(b), ecX(a), t7, ec->f, stack);
715 0
	zmSub(ecX(b), ecX(b), t1, ec->f);
716 0
	gfpDouble(ecX(b), ecX(b), ec->f);
717 0
	gfpDouble(ecX(b), ecX(b), ec->f);
718
}
719

720 1
static size_t ecpTplJ_deep(size_t n, size_t f_deep)
721
{
722 1
	return O_OF_W(8 * n) + f_deep;
723
}
724

725
// [3n]b <- 3[3n]a (P <- 3P)
726 0
static void ecpTplJA3(word b[], const word a[], const ec_o* ec, void* stack)
727
{
728 0
	const size_t n = ec->f->n;
729
	// переменные в stack
730 0
	word* t1 = (word*)stack;
731 0
	word* t2 = t1 + n;
732 0
	word* t3 = t2 + n;
733 0
	word* t4 = t3 + n;
734 0
	word* t5 = t4 + n;
735 0
	word* t6 = t5 + n;
736 0
	word* t7 = t6 + n;
737 0
	stack = t7 + n;
738
	// pre
739 0
	ASSERT(ecIsOperable(ec) && ec->d == 3);
740 0
	ASSERT(ecpSeemsOn3(a, ec));
741 0
	ASSERT(wwIsSameOrDisjoint(a, b, 3 * n));
742
	// t1 <- ya^2 [YY]
743 0
	qrSqr(t1, ecY(a, n), ec->f, stack);
744
	// t2 <- za^2 [ZZ]
745 0
	qrSqr(t2, ecZ(a, n), ec->f, stack);
746
	// t3 <- t1^2 [YYYY]
747 0
	qrSqr(t3, t1, ec->f, stack);
748
	// t4 <- 3(xa - t2)(xa + t2) [M]
749 0
	zmSub(t4, ecX(a), t2, ec->f);
750 0
	zmAdd(t5, ecX(a), t2, ec->f);
751 0
	qrMul(t4, t4, t5, ec->f, stack);
752 0
	gfpDouble(t5, t4, ec->f);
753 0
	zmAdd(t4, t4, t5, ec->f);
754
	// t5 <- t4^2 [MM]
755 0
	qrSqr(t5, t4, ec->f, stack);
756
	// t6 <- 12 xa t1 - t5 [E]
757 0
	qrMul(t6, ecX(a), t1, ec->f, stack);
758 0
	gfpDouble(t7, t6, ec->f);
759 0
	zmAdd(t6, t6, t7, ec->f);
760 0
	gfpDouble(t6, t6, ec->f);
761 0
	gfpDouble(t6, t6, ec->f);
762 0
	zmSub(t6, t6, t5, ec->f);
763
	// t7 <- t6^2 [EE]
764 0
	qrSqr(t7, t6, ec->f, stack);
765
	// t3 <- 16 t3 [T]
766 0
	gfpDouble(t3, t3, ec->f);
767 0
	gfpDouble(t3, t3, ec->f);
768 0
	gfpDouble(t3, t3, ec->f);
769 0
	gfpDouble(t3, t3, ec->f);
770
	// zb <- (za + t6)^2 - t2 - t7
771 0
	zmAdd(ecZ(b, n), ecZ(a, n), t6, ec->f);
772 0
	qrSqr(ecZ(b, n), ecZ(b, n), ec->f, stack);
773 0
	zmSub(ecZ(b, n), ecZ(b, n), t2, ec->f);
774 0
	zmSub(ecZ(b, n), ecZ(b, n), t7, ec->f);
775
	// t2 <- (t4 + t6)^2 - t5 - t7 - t3 [U] 
776 0
	zmAdd(t2, t4, t6, ec->f);
777 0
	qrSqr(t2, t2, ec->f, stack);
778 0
	zmSub(t2, t2, t5, ec->f);
779 0
	zmSub(t2, t2, t7, ec->f);
780 0
	zmSub(t2, t2, t3, ec->f);
781
	// yb <- 8 ya (t2(t3 - t2) - t6 t7)
782 0
	zmSub(t3, t3, t2, ec->f);
783 0
	qrMul(t3, t2, t3, ec->f, stack);
784 0
	qrMul(t6, t6, t7, ec->f, stack);
785 0
	zmSub(t3, t3, t6, ec->f);
786 0
	qrMul(ecY(b, n), ecY(a, n), t3, ec->f, stack);
787 0
	gfpDouble(ecY(b, n), ecY(b, n), ec->f);
788 0
	gfpDouble(ecY(b, n), ecY(b, n), ec->f);
789 0
	gfpDouble(ecY(b, n), ecY(b, n), ec->f);
790
	// xb <- 4 (xa t7 - 4 t1 t2)
791 0
	qrMul(t1, t1, t2, ec->f, stack);
792 0
	gfpDouble(t1, t1, ec->f);
793 0
	gfpDouble(t1, t1, ec->f);
794 0
	qrMul(ecX(b), ecX(a), t7, ec->f, stack);
795 0
	zmSub(ecX(b), ecX(b), t1, ec->f);
796 0
	gfpDouble(ecX(b), ecX(b), ec->f);
797 0
	gfpDouble(ecX(b), ecX(b), ec->f);
798
}
799

800 1
size_t ecpTplJA3_deep(size_t n, size_t f_deep)
801
{
802 1
	return O_OF_W(7 * n) + f_deep;
803
}
804

805 1
bool_t ecpCreateJ(ec_o* ec, const qr_o* f, const octet A[], const octet B[], 
806
	void* stack)
807
{
808
	register bool_t bA3;
809
	word* t;
810
	// pre
811 1
	ASSERT(memIsValid(ec, sizeof(ec_o)));
812 1
	ASSERT(gfpIsOperable(f));
813 1
	ASSERT(memIsValid(A, f->no)); 
814 1
	ASSERT(memIsValid(B, f->no)); 
815
	// обнулить
816 1
	memSetZero(ec, sizeof(ec_o));
817
	// зафикисровать размерности
818 1
	ec->d = 3;
819
	// запомнить базовое поле
820 1
	ec->f = f;
821
	// сохранить коэффициенты
822 1
	ec->A = (word*)ec->descr;
823 1
	ec->B = ec->A + f->n;
824 1
	if (!qrFrom(ec->A, A, ec->f, stack) || !qrFrom(ec->B, B, ec->f, stack))
825 0
		return FALSE;
826
	// t <- -3
827 1
	t = (word*)stack;
828 1
	gfpDouble(t, f->unity, f);
829 1
	zmAdd(t, t, f->unity, f);
830 1
	zmNeg(t, t, f);
831
	// bA3 <- A == -3?
832 1
	bA3 = qrCmp(t, ec->A, f) == 0;
833
	// подготовить буферы для описания группы точек
834 1
	ec->base = ec->B + f->n;
835 1
	ec->order = ec->base + 2 * f->n;
836
	// настроить интерфейсы
837 1
	ec->froma = ecpFromAJ;
838 1
	ec->toa = ecpToAJ;
839 1
	ec->neg = ecpNegJ;
840 1
	ec->add = ecpAddJ;
841 1
	ec->adda = ecpAddAJ;
842 1
	ec->sub = ecpSubJ;
843 1
	ec->suba = ecpSubAJ;
844 1
	ec->dbl = bA3 ? ecpDblJA3 : ecpDblJ;
845 1
	ec->dbla = ecpDblAJ;
846 1
	ec->tpl = bA3 ? ecpTplJA3 : ecpTplJ;
847 1
	ec->deep = utilMax(9,
848
		ecpToAJ_deep(f->n, f->deep),
849
		ecpAddJ_deep(f->n, f->deep),
850
		ecpAddAJ_deep(f->n, f->deep),
851
		ecpSubJ_deep(f->n, f->deep),
852
		ecpSubAJ_deep(f->n, f->deep),
853 1
		bA3 ? ecpDblJA3_deep(f->n, f->deep) : ecpDblJ_deep(f->n, f->deep),
854
		ecpDblAJ_deep(f->n, f->deep),
855 1
		bA3 ? ecpTplJA3_deep(f->n, f->deep) : ecpTplJ_deep(f->n, f->deep));
856
	// настроить
857 1
	ec->hdr.keep = sizeof(ec_o) + O_OF_W(5 * f->n + 1);
858 1
	ec->hdr.p_count = 6;
859 1
	ec->hdr.o_count = 1;
860
	// все нормально
861 1
	bA3 = 0;
862 1
	return TRUE;
863
}
864

865 1
size_t ecpCreateJ_keep(size_t n)
866
{
867 1
	return sizeof(ec_o) + O_OF_W(5 * n + 1);
868
}
869

870 1
size_t ecpCreateJ_deep(size_t n, size_t f_deep)
871
{
872 1
	return utilMax(11,
873
		O_OF_W(n),
874
		ecpToAJ_deep(n, f_deep),
875
		ecpAddJ_deep(n, f_deep),
876
		ecpAddAJ_deep(n, f_deep),
877
		ecpSubJ_deep(n, f_deep),
878
		ecpSubAJ_deep(n, f_deep),
879
		ecpDblJ_deep(n, f_deep),
880
		ecpDblJA3_deep(n, f_deep),
881
		ecpDblAJ_deep(n, f_deep),
882
		ecpTplJ_deep(n, f_deep),
883
		ecpTplJA3_deep(n, f_deep));
884
}
885

886
/*
887
*******************************************************************************
888
Свойства кривой 
889
*******************************************************************************
890
*/
891

892 1
bool_t ecpIsValid(const ec_o* ec, void* stack)
893
{
894 1
	size_t n = ec->f->n;
895
	// переменные в stack
896 1
	word* t1 = (word*)stack;
897 1
	word* t2 = t1 + n;
898 1
	word* t3 = t2 + n;
899 1
	stack = t3 + n;
900
	// кривая работоспособна? поле ec->f корректно?
901
	// ec->deep >= ec->f->deep?
902
	// A, B \in ec->f?
903 1
	if (!ecIsOperable2(ec) ||
904 1
		!gfpIsValid(ec->f, stack) ||
905 1
		ec->deep < ec->f->deep ||
906 1
		!zmIsIn(ec->A, ec->f) || 
907 1
		!zmIsIn(ec->B, ec->f))
908 0
		return FALSE;
909
	// t1 <- 4 A^3
910 1
	qrSqr(t1, ec->A, ec->f, stack);
911 1
	qrMul(t1, t1, ec->A, ec->f, stack);
912 1
	gfpDouble(t1, t1, ec->f);
913 1
	gfpDouble(t1, t1, ec->f);
914
	// t3 <- 3 B^2
915 1
	qrSqr(t2, ec->B, ec->f, stack);
916 1
	gfpDouble(t3, t2, ec->f);
917 1
	zmAdd(t3, t3, t2, ec->f);
918
	// t2 <- 3 t3 [27 B^2]
919 1
	gfpDouble(t2, t3, ec->f);
920 1
	zmAdd(t2, t3, t2, ec->f);
921
	// t1 <- t1 + t2 [4 A^3 + 27 B^2 -- дискриминант]
922 1
	zmAdd(t1, t1, t2, ec->f);
923
	// t1 == 0 => сингулярная кривая
924 1
	return !qrIsZero(t1, ec->f);
925
}
926

927 1
size_t ecpIsValid_deep(size_t n, size_t f_deep)
928
{
929 1
	return O_OF_W(3 * n) + 
930 1
		utilMax(2,
931
			f_deep,
932
			gfpIsValid_deep(n));
933
}
934

935 1
bool_t ecpSeemsValidGroup(const ec_o* ec, void* stack)
936
{
937 1
	size_t n = ec->f->n;
938
	int cmp;
939
	word w;
940
	// переменные в stack
941 1
	word* t1 = (word*)stack;
942 1
	word* t2 = t1 + n + 1;
943 1
	word* t3 = t2 + n + 2;
944 1
	stack = t3 + 2 * n;
945
	// pre
946 1
	ASSERT(ecIsOperable(ec));
947
	// ecIsOperableGroup(ec) == TRUE? base \in ec?
948 1
	if (!ecIsOperableGroup(ec) ||
949 1
		!ecpIsOnA(ec->base, ec, stack))
950 0
		return FALSE;
951
	// [n + 2]t1 <- order * cofactor
952 1
	t1[n + 1] = zzMulW(t1, ec->order, n + 1, ec->cofactor);
953
	// t1 <- |t1 - (p + 1)|
954 1
	if (zzSubW2(t1, n + 2, 1))
955 0
		return FALSE;
956 1
	if (wwCmp2(t1, n + 2, ec->f->mod, n) >= 0)
957 1
		zzSubW2(t1 + n, 2, zzSub2(t1, ec->f->mod, n));
958
	else
959 1
		zzSub(t1, ec->f->mod, t1, n);
960
	// n <- длина t1
961 1
	n = wwWordSize(t1, n + 2);
962
	// n > ec->f->n => t1^2 > 4 p
963 1
	if (n > ec->f->n)
964 0
		return FALSE;
965
	// [2n]t2 <- ([n]t1)^2
966 1
	zzSqr(t2, t1, n, stack);
967
	// условие Хассе: t2 <= 4 p?
968 1
	w = wwGetBits(t2, 0, 2);
969 1
	wwShLo(t2, 2 * n, 2);
970 1
	cmp = wwCmp2(t2, 2 * n, ec->f->mod, ec->f->n);
971 1
	if (cmp > 0 || cmp == 0 && w != 0)
972 0
		return FALSE;
973
	// все нормально
974 1
	return TRUE;
975
}
976

977 1
size_t ecpSeemsValidGroup_deep(size_t n, size_t f_deep)
978
{
979 1
	return O_OF_W(4 * n + 3) + 
980 1
		utilMax(2, 
981
			ecpIsOnA_deep(n, f_deep),
982
			zzSqr_deep(n));
983
}
984

985 1
bool_t ecpIsSafeGroup(const ec_o* ec, size_t mov_threshold, void* stack)
986
{
987 1
	size_t n1 = ec->f->n + 1;
988
	// переменные в stack
989 1
	word* t1 = (word*)stack;
990 1
	word* t2 = t1 + n1;
991 1
	stack = t2 + n1;
992
	// pre
993 1
	ASSERT(ecIsOperable(ec));
994 1
	ASSERT(ecIsOperableGroup(ec));
995
	// order -- простое?
996 1
	n1 = wwWordSize(ec->order, n1);
997 1
	if (!priIsPrime(ec->order, n1, stack))
998 0
		return FALSE;
999
	// order == p?
1000 1
	if (wwCmp2(ec->f->mod, ec->f->n, ec->order, n1) == 0)
1001 0
		return FALSE;
1002
	// проверка MOV
1003 1
	if (mov_threshold)
1004
	{
1005 1
		zzMod(t1, ec->f->mod, ec->f->n, ec->order, n1, stack);
1006 1
		wwCopy(t2, t1, n1);
1007 1
		if (wwCmpW(t2, n1, 1) == 0)
1008 0
			return FALSE;
1009 1
		while (--mov_threshold)
1010
		{
1011 1
			zzMulMod(t2, t2, t1, ec->order, n1, stack);
1012 1
			if (wwCmpW(t2, n1, 1) == 0)
1013 0
				return FALSE;
1014
		}
1015
	}
1016
	// все нормально
1017 1
	return TRUE;
1018
}
1019

1020 1
size_t ecpIsSafeGroup_deep(size_t n)
1021
{
1022 1
	const size_t n1 = n + 1;
1023 1
	return O_OF_W(2 * n1) +
1024 1
		utilMax(3,
1025
			priIsPrime_deep(n1),
1026
			zzMod_deep(n, n1),
1027
			zzMulMod_deep(n1));
1028
}
1029

1030
/*
1031
*******************************************************************************
1032
Арифметика аффинных точек
1033

1034
Сложение A <- A + A:
1035
	1D + 1S + 1M + 6add \approx 102M
1036

1037
Удвоение A <- 2A:
1038
	1D + 2S + 1M + 5add + 1*3 + 1*2 \approx 103M
1039
*******************************************************************************
1040
*/
1041

1042 1
bool_t ecpIsOnA(const word a[], const ec_o* ec, void* stack)
1043
{
1044 1
	const size_t n = ec->f->n;
1045
	// переменные в stack
1046 1
	word* t1 = (word*)stack;
1047 1
	word* t2 = t1 + n;
1048 1
	stack = t2 + n;
1049
	// pre
1050 1
	ASSERT(ecIsOperable(ec));
1051
	// xa, ya \in ec->f?
1052 1
	if (!zmIsIn(ecX(a), ec->f) || !zmIsIn(ecY(a, n), ec->f))
1053 0
		return FALSE;
1054
	// t1 <- (xa^2 + A)xa + B
1055 1
	qrSqr(t1, ecX(a), ec->f, stack);
1056 1
	zmAdd(t1, t1, ec->A, ec->f);
1057 1
	qrMul(t1, t1, ecX(a), ec->f, stack);
1058 1
	zmAdd(t1, t1, ec->B, ec->f);
1059
	// t2 <- ya^2
1060 1
	qrSqr(t2, ecY(a, n), ec->f, stack);
1061
	// t1 == t2?
1062 1
	return qrCmp(t1, t2, ec->f) == 0;
1063
}
1064

1065 1
size_t ecpIsOnA_deep(size_t n, size_t f_deep)
1066
{
1067 1
	return O_OF_W(2 * n) + f_deep;
1068
}
1069

1070 0
void ecpNegA(word b[], const word a[], const ec_o* ec)
1071
{
1072 0
	const size_t n = ec->f->n;
1073
	// pre
1074 0
	ASSERT(ecIsOperable(ec));
1075 0
	ASSERT(ecpSeemsOnA(a, ec));
1076 0
	ASSERT(wwIsSameOrDisjoint(a, b, 3 * n));
1077
	// (xb, yb) <- (xa, -ya)
1078 0
	qrCopy(ecX(b), ecX(a), ec->f);
1079 0
	zmNeg(ecY(b, n), ecY(a, n), ec->f);
1080
}
1081

1082 0
bool_t ecpAddAA(word c[], const word a[], const word b[], const ec_o* ec,
1083
	void* stack)
1084
{
1085 0
	const size_t n = ec->f->n;
1086
	// переменные в stack
1087 0
	word* t1 = (word*)stack;
1088 0
	word* t2 = t1 + n;
1089 0
	word* t3 = t2 + n;
1090 0
	stack = t3 + n;
1091
	// pre
1092 0
	ASSERT(ecIsOperable(ec));
1093 0
	ASSERT(ecpSeemsOnA(a, ec));
1094 0
	ASSERT(ecpSeemsOnA(b, ec));
1095
	// xa != xb?
1096 0
	if (qrCmp(ecX(a), ecX(b), ec->f) != 0)
1097
	{
1098
		// t1 <- xa - xb
1099 0
		zmSub(t1, ecX(a), ecX(b), ec->f);
1100
		// t2 <- ya - yb
1101 0
		zmSub(t2, ecY(a, n), ecY(b, n), ec->f);
1102
	}
1103
	else
1104
	{
1105
		// ya != yb или yb == 0 => a == -b => a + b == O
1106 0
		if (qrCmp(ecY(a, n), ecY(b, n), ec->f) != 0 || 
1107 0
			qrIsZero(ecY(b, n), ec->f))
1108 0
			return FALSE;
1109
		// t2 <- 3 xa^2 + A
1110 0
		qrSqr(t1, ecX(a), ec->f, stack);
1111 0
		gfpDouble(t2, t1, ec->f);
1112 0
		zmAdd(t2, t2, t1, ec->f);
1113 0
		zmAdd(t2, t2, ec->A, ec->f);
1114
		// t1 <- 2 ya
1115 0
		gfpDouble(t1, ecY(a, n), ec->f);
1116
	}
1117
	// t2 <- t2 / t1 = \lambda
1118 0
	qrDiv(t2, t2, t1, ec->f, stack);
1119
	// t1 <- \lambda^2 - xa - xb = xc
1120 0
	qrSqr(t1, t2, ec->f, stack);
1121 0
	zmSub(t1, t1, ecX(a), ec->f);
1122 0
	zmSub(t1, t1, ecX(b), ec->f);
1123
	// t3 <- xa - xc
1124 0
	zmSub(t3, ecX(a), t1, ec->f);
1125
	// t2 <- \lambda(xa - xc) - ya
1126 0
	qrMul(t2, t2, t3, ec->f, stack);
1127 0
	zmSub(t2, t2, ecY(a, n), ec->f);
1128
	// выгрузить результат
1129 0
	qrCopy(ecX(c), t1, ec->f);
1130 0
	qrCopy(ecY(c, n), t2, ec->f);
1131 0
	return TRUE;
1132
}
1133

1134 0
size_t ecpAddAA_deep(size_t n, size_t f_deep)
1135
{
1136 0
	return O_OF_W(3 * n) + f_deep;
1137
}
1138

1139 1
bool_t ecpSubAA(word c[], const word a[], const word b[], const ec_o* ec,
1140
	void* stack)
1141
{
1142 1
	const size_t n = ec->f->n;
1143
	// переменные в stack
1144 1
	word* t1 = (word*)stack;
1145 1
	word* t2 = t1 + n;
1146 1
	word* t3 = t2 + n;
1147 1
	stack = t3 + n;
1148
	// pre
1149 1
	ASSERT(ecIsOperable(ec));
1150 1
	ASSERT(ecpSeemsOnA(a, ec));
1151 1
	ASSERT(ecpSeemsOnA(b, ec));
1152
	// xa != xb?
1153 1
	if (qrCmp(ecX(a), ecX(b), ec->f) != 0)
1154
	{
1155
		// t1 <- xa - xb
1156 1
		zmSub(t1, ecX(a), ecX(b), ec->f);
1157
		// t2 <- ya + yb
1158 1
		zmAdd(t2, ecY(a, n), ecY(b, n), ec->f);
1159
	}
1160
	else
1161
	{
1162
		// ya == yb => a == b => a - b == O
1163 0
		if (qrCmp(ecY(a, n), ecY(b, n), ec->f) == 0)
1164 0
			return FALSE;
1165
		// здесь должно быть a == -b => a - b = 2a
1166
		// t2 <- 3 xa^2 + A
1167 0
		qrSqr(t1, ecX(a), ec->f, stack);
1168 0
		gfpDouble(t2, t1, ec->f);
1169 0
		zmAdd(t2, t2, t1, ec->f);
1170 0
		zmAdd(t2, t2, ec->A, ec->f);
1171
		// t1 <- 2 ya
1172 0
		gfpDouble(t1, ecY(a, n), ec->f);
1173
	}
1174
	// t2 <- t2 / t1 = \lambda
1175 1
	qrDiv(t2, t2, t1, ec->f, stack);
1176
	// t1 <- \lambda^2 - xa - xb = xc
1177 1
	qrSqr(t1, t2, ec->f, stack);
1178 1
	zmSub(t1, t1, ecX(a), ec->f);
1179 1
	zmSub(t1, t1, ecX(b), ec->f);
1180
	// t3 <- xa - xc
1181 1
	zmSub(t3, ecX(a), t1, ec->f);
1182
	// t2 <- \lambda(xa - xc) - ya
1183 1
	qrMul(t2, t2, t3, ec->f, stack);
1184 1
	zmSub(t2, t2, ecY(a, n), ec->f);
1185
	// выгрузить результат
1186 1
	qrCopy(ecX(c), t1, ec->f);
1187 1
	qrCopy(ecY(c, n), t2, ec->f);
1188 1
	return TRUE;
1189
}
1190

1191 1
size_t ecpSubAA_deep(size_t n, size_t f_deep)
1192
{
1193 1
	return O_OF_W(3 * n) + f_deep;
1194
}
1195

1196
/*
1197
*******************************************************************************
1198
Алгоритм SWU
1199

1200
\todo Регуляризировать (qrIsUnitySafe).
1201
*******************************************************************************
1202
*/
1203

1204 1
void ecpSWU(word b[], const word a[], const ec_o* ec, void* stack)
1205
{
1206 1
	const size_t n = ec->f->n;
1207
	register size_t mask;
1208
	// переменные в stack [x2 после x1, s после y!]
1209 1
	word* t = (word*)stack;
1210 1
	word* x1 = t + n;
1211 1
	word* x2 = x1 + n;
1212 1
	word* y = x2 + n;
1213 1
	word* s = y + n; 
1214 1
	stack = s + n;
1215
	// pre
1216 1
	ASSERT(ecIsOperable(ec));
1217 1
	ASSERT(zmIsIn(a, ec->f));
1218 1
	ASSERT(wwGetBits(ec->f->mod, 0, 2) == 3);
1219 1
	ASSERT(!qrIsZero(ec->A, ec->f) && !qrIsZero(ec->B, ec->f));
1220
	// t <- -a^2
1221 1
	qrSqr(t, a, ec->f, stack);
1222 1
	zmNeg(t, t, ec->f);
1223
	// s <- p - 2
1224 1
	wwCopy(s, ec->f->mod, n);
1225 1
	zzSubW2(s, n, 2);
1226
	// x1 <- -B(1 + t + t^2)(A(t + t^2))^{p - 2} 
1227 1
	qrSqr(x2, t, ec->f, stack);
1228 1
	qrAdd(x2, x2, t, ec->f);
1229 1
	qrMul(x1, x2, ec->A, ec->f, stack);
1230 1
	qrPower(x1, x1, s, n, ec->f, stack);
1231 1
	qrAddUnity(x2, x2, ec->f);
1232 1
	qrMul(x1, x1, x2, ec->f, stack);
1233 1
	qrMul(x1, x1, ec->B, ec->f, stack);
1234 1
	zmNeg(x1, x1, ec->f);
1235
	// y <- (x1)^3 + A x1 + B
1236 1
	qrSqr(x2, x1, ec->f, stack);
1237 1
	qrMul(x2, x2, x1, ec->f, stack);
1238 1
	qrMul(y, x1, ec->A, ec->f, stack);
1239 1
	qrAdd(y, y, x2, ec->f);
1240 1
	qrAdd(y, y, ec->B, ec->f);
1241
	// x2 <- x1 t
1242 1
	qrMul(x2, x1, t, ec->f, stack);
1243
	// t <- y^{(p - 1) - (p + 1) / 4} = y^{s - (p - 3) / 4}
1244 1
	wwCopy(t, ec->f->mod, n);
1245 1
	wwShLo(t, n, 2);
1246 1
	zzSub(s, s, t, n);
1247 1
	qrPower(t, y, s, n, ec->f, stack);
1248
	// s <- a^3 y
1249 1
	qrSqr(s, a, ec->f, stack);
1250 1
	qrMul(s, s, a, ec->f, stack);
1251 1
	qrMul(s, s, y, ec->f, stack);
1252
	// mask <- t^2 y == 1 ? 0 : -1;
1253 1
	qrSqr(b, t, ec->f, stack);
1254 1
	qrMul(b, b, y, ec->f, stack);
1255 1
	mask = qrIsUnity(b, ec->f) - SIZE_1;
1256
	// b <- mask == 0 ? (x1, t y) : (x2, t s)
1257 1
	qrCopy(ecX(b), x1 + (mask & n), ec->f);
1258 1
	qrMul(ecY(b, n), t, y + (mask & n), ec->f, stack);
1259
	// очистка
1260 1
	mask = 0;
1261
}
1262

1263 1
size_t ecpSWU_deep(size_t n, size_t f_deep)
1264
{
1265 1
	return O_OF_W(5 * n) + 
1266 1
		utilMax(2,
1267
			f_deep,
1268
			qrPower_deep(n, n, f_deep));
1269
}

Read our documentation on viewing source code .

Loading