1
/*
2
*******************************************************************************
3
\file zz_gcd.c
4
\brief Multiple-precision unsigned integers: Euclidian gcd algorithms
5
\project bee2 [cryptographic library]
6
\author (C) Sergey Agievich [agievich@{bsu.by|gmail.com}]
7
\created 2012.04.22
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/util.h"
16
#include "bee2/math/ww.h"
17
#include "bee2/math/zz.h"
18

19
/*
20
*******************************************************************************
21
Алгоритм Евклида
22

23
В функциях zzGCD(), zzExGCD() реализованы бинарные алгоритмы,
24
не требующие прямых делений.
25

26
В функции zzExGCD() пересчитываются числа da, db, da0, db0 такие, что
27
	da0 * aa - db0 * bb = (-1)^sign0 u,
28
	da * aa - db * bb = (-1)^sign v,
29
где aa = a / 2^s, bb = b / 2^s, s -- max целое т.ч. 2^s | a и 2^s | b.
30

31
Числа u и v поддерживают вычисление НОД(aa, bb). Если u >= v, то u
32
заменяется на u - v, а если u < v, то v заменяется на v - u.
33
Как только u == 0 вычисления останавливаются и возвращается тройка
34
(2^s * v, da, db).
35

36
В функции zzExGCD() реализован алгоритм:
37
	u <- aa
38
	da0 <- 1, db0 <- 0, sign0 <- 0
39
	v <- bb
40
	da <- 0, db <- 1, sign <- 0
41
	пока (u != 0)
42
	{
43
		пока (u -- четное)
44
			u <- u / 2
45
			если (da0 и db0 -- четные)
46
				da0 <- da0 / 2, db0 <- db0 / 2
47
			иначе 
48
[
49
	Пусть da0 -- четное, db0 -- нечетное. Поскольку da0 aa + db0 bb -- четное,
50
	bb -- четное. Но aa и bb не могут быть одновременно четными. 
51
	Поэтому aa -- нечетное. В конце концов, da0 + bb и db0 + aa -- четные.
52
	Аналогично рассматриваются другие варианты четности da0 и db0.
53
]
54
				da0 <- (da0 + bb) / 2, db0 <- (db0 + aa) / 2
55
		пока (v -- четное)
56
			v <- v / 2
57
			если (da и db -- четные)
58
				da <- da / 2, db <- db / 2
59
			иначе
60
				da <- (da + bb) / 2, db <- (db + aa) / 2
61
		если (u >= v)
62
			u <- u - v
63
			если (sign0 == sign)
64
				da0 <- da0 - da, db0 <- db0 - db
65
				если (da0 < 0)
66
					da0 <- -da0, db0 <- -db0, sign0 <- 1 - sign0
67
			иначе // sign0 != sign
68
				da0 <- da0 + da, db0 <- db0 + db
69
				если (da0 > bb)
70
					da0 <- da0 - bb, db0 <- db0 - aa			(*)
71
		иначе // u < v
72
			v <- v - u
73
			если (sign0 == sign)
74
				da <- da - da0, db <- db - db0
75
				если (da < 0)
76
					da <- -da, db <- -db, sign <- 1 - sign
77
			иначе // sign0 != sign
78
				da <- da + da0, db <- db + db0
79
				если (da > bb)
80
					da <- da - bb, db <- db - aa				(**)
81
	}
82
	Корректировки (*), (**) гарантируют, что da0, da < bb и db0, db < aa.
83

84
\remark Эксперименты с различными реализациями алгоритма Евклида:
85
1) бинарный алгоритм опережает обычный (полные деления) примерно в 2 раза и 
86
2) опережает смешанный (полные деления и деления на 2) примерно в 1.5 раза.
87

88
\todo Регуляризация?
89
*******************************************************************************
90
*/
91

92 1
void zzGCD(word d[], const word a[], size_t n, const word b[], size_t m,
93
	void* stack)
94
{
95
	register size_t s;
96
	// переменные в stack
97 1
	word* u = (word*)stack;
98 1
	word* v = u + n;
99 1
	stack = v + m;
100
	// pre
101 1
	ASSERT(wwIsDisjoint2(a, n, d, MIN2(n, m)));
102 1
	ASSERT(wwIsDisjoint2(b, m, d, MIN2(n, m)));
103 1
	ASSERT(!wwIsZero(a, n) && !wwIsZero(b, m));
104
	// d <- 0
105 1
	wwSetZero(d, MIN2(n, m));
106
	// u <- a, v <- b
107 1
	wwCopy(u, a, n);
108 1
	wwCopy(v, b, m);
109
	// найти максимальное s т.ч. 2^s | u и 2^s | v
110 1
	s = utilMin(2, wwLoZeroBits(u, n), wwLoZeroBits(v, m));
111
	// u <- u / 2^s, v <- v / 2^s
112 1
	wwShLo(u, n, s);
113 1
	n = wwWordSize(u, n);
114 1
	wwShLo(v, m, s);
115 1
	m = wwWordSize(v, m);
116
	// итерации
117
	do
118
	{
119 1
		wwShLo(u, n, wwLoZeroBits(u, n));
120 1
		n = wwWordSize(u, n);
121 1
		wwShLo(v, m, wwLoZeroBits(v, m));
122 1
		m = wwWordSize(v, m);
123
		// u >= v?
124 1
		if (wwCmp2(u, n, v, m) >= 0)
125
			// u <- u - v
126 1
			zzSubW2(u + m, n - m, zzSub2(u, v, m));
127
		else
128
			// v <- v - u
129 1
			zzSubW2(v + n, m - n, zzSub2(v, u, n));
130
	}
131 1
	while (!wwIsZero(u, n));
132
	// d <- v
133 1
	wwCopy(d, v, m);
134
	// d <- d * 2^s
135 1
	wwShHi(d, W_OF_B(wwBitSize(d, m) + s), s);
136
	// очистка
137 1
	s = 0;
138
}
139

140 1
size_t zzGCD_deep(size_t n, size_t m)
141
{
142 1
	return O_OF_W(n + m);
143
}
144

145 1
bool_t zzIsCoprime(const word a[], size_t n, const word b[], size_t m, void* stack)
146
{
147 1
	word* d = (word*)stack;
148 1
	stack = d + MIN2(n, m);
149
	// a == 0 => (a, b) = b
150 1
	if (wwIsZero(a, n))
151 0
		return wwIsW(b, m, 1);
152
	// b == 0 => (a, b) = a
153 1
	if (wwIsZero(b, m))
154 0
		return wwIsW(a, n, 1);
155
	// d <- (a, b), d == 1?
156 1
	zzGCD(d, a, n, b, m, stack);
157 1
	return wwIsW(d, MIN2(n, m), 1);
158
}
159

160 1
size_t zzIsCoprime_deep(size_t n, size_t m)
161
{
162 1
	return O_OF_W(MIN2(n, m)) + zzGCD_deep(n, m);
163
}
164

165 1
void zzLCM(word d[], const word a[], size_t n, const word b[], size_t m,
166
	void* stack)
167
{
168
	// переменные в stack
169 1
	word* prod = (word*)stack;
170 1
	word* gcd = prod + n + m;
171 1
	word* r = gcd + MIN2(n, m);
172 1
	stack = r + MIN2(n, m);
173
	// pre
174 1
	ASSERT(wwIsDisjoint2(a, n, d, MAX2(n, m)));
175 1
	ASSERT(wwIsDisjoint2(b, m, d, MAX2(n, m)));
176 1
	ASSERT(!wwIsZero(a, n) && !wwIsZero(b, m));
177
	// d <- 0
178 1
	wwSetZero(d, n + m);
179
	// нормализация
180 1
	n = wwWordSize(a, n);
181 1
	m = wwWordSize(b, m);
182
	// prod <- a * b
183 1
	zzMul(prod, a, n, b, m, stack);
184
	// gcd <- (a, b)
185 1
	zzGCD(gcd, a, n, b, m, stack);
186
	// (n, m) <- (|prod|, |gcd|)
187 1
	if (n < m)
188 1
		SWAP(n, m);
189 1
	n += m;
190 1
	m = wwWordSize(gcd, m);
191
	// d <- prod \div gcd
192 1
	zzDiv(d, r, prod, n, gcd, m, stack);
193
}
194

195 1
size_t zzLCM_deep(size_t n, size_t m)
196
{
197 1
	return O_OF_W(n + m + MIN2(n, m) + MIN2(n, m)) +
198 1
		utilMax(3, 
199
			zzMul_deep(n, m), 
200
			zzGCD_deep(n, m), 
201
			zzMod_deep(n + m, MIN2(n, m)));
202
}
203

204 1
int zzExGCD(word d[], word da[], word db[], const word a[], size_t n,
205
	const word b[], size_t m, void* stack)
206
{
207
	register size_t s;
208
	register size_t nu, mv;
209 1
	register int sign0 = 0, sign = 1;
210
	// переменные в stack
211 1
	word* aa = (word*)stack;
212 1
	word* bb = aa + n;
213 1
	word* u = bb + m;
214 1
	word* v = u + n;
215 1
	word* da0 = v + m;
216 1
	word* db0 = da0 + m;
217 1
	stack = db0 + n;
218
	// pre
219 1
	ASSERT(wwIsDisjoint3(da, m, db, n, d, MIN2(n, m)));
220 1
	ASSERT(wwIsDisjoint2(a, n, d, MIN2(n, m)));
221 1
	ASSERT(wwIsDisjoint2(b, m, d, MIN2(n, m)));
222 1
	ASSERT(wwIsDisjoint2(a, n, da, m));
223 1
	ASSERT(wwIsDisjoint2(b, m, da, m));
224 1
	ASSERT(wwIsDisjoint2(a, n, db, n));
225 1
	ASSERT(wwIsDisjoint2(b, m, db, n));
226 1
	ASSERT(!wwIsZero(a, n) && !wwIsZero(b, m));
227
	// d <- 0, da0 <- 1, db0 <- 0, da <- 0, db <- 1
228 1
	wwSetZero(d, MIN2(n, m));
229 1
	wwSetW(da0, m, 1);
230 1
	wwSetZero(db0, n);
231 1
	wwSetZero(da, m);
232 1
	wwSetW(db, n, 1);
233
	// найти максимальное s т.ч. 2^s | a и 2^s | b
234 1
	s = utilMin(2, wwLoZeroBits(a, n), wwLoZeroBits(b, m));
235
	// aa <- a / 2^s, bb <- b / 2^s
236 1
	wwCopy(aa, a, n), wwShLo(aa, n, s), n = wwWordSize(aa, n);
237 1
	wwCopy(bb, b, m), wwShLo(bb, m, s), m = wwWordSize(bb, m);
238
	// u <- aa, v <- bb
239 1
	wwCopy(u, aa, n);
240 1
	wwCopy(v, bb, m);
241 1
	nu = n, mv = m;
242
	// итерации
243
	do
244
	{
245
		// пока u четное
246 1
		for (; u[0] % 2 == 0; wwShLo(u, nu, 1))
247 1
			if (da0[0] % 2 == 0 && db0[0] % 2 == 0)
248
			{
249
				// da0 <- da0 / 2, db0 <- db0 / 2
250 1
				wwShLo(da0, m, 1);
251 1
				wwShLo(db0, n, 1);
252
			}
253
			else
254
			{
255 1
				ASSERT((da0[0] + bb[0]) % 2 == 0);
256 1
				ASSERT((db0[0] + aa[0]) % 2 == 0);
257
				// da0 <- (da0 + bb) / 2, db0 <- (db0 + aa) / 2
258 1
				wwShLoCarry(da0, m, 1, zzAdd2(da0, bb, m));
259 1
				wwShLoCarry(db0, n, 1, zzAdd2(db0, aa, n));
260
			}
261
		// пока v четное
262 1
		for (; v[0] % 2 == 0; wwShLo(v, mv, 1))
263 1
			if (da[0] % 2 == 0 && db[0] % 2 == 0)
264
			{
265
				// da <- da / 2, db <- db / 2
266 1
				wwShLo(da, m, 1);
267 1
				wwShLo(db, n, 1);
268
			}
269
			else
270
			{
271 1
				ASSERT((da[0] + bb[0]) % 2 == 0);
272 1
				ASSERT((db[0] + aa[0]) % 2 == 0);
273
				// da <- (da + bb) / 2, db <- (db + aa) / 2
274 1
				wwShLoCarry(da, m, 1, zzAdd2(da, bb, m));
275 1
				wwShLoCarry(db, n, 1, zzAdd2(db, aa, n));
276
			}
277
		// нормализация
278 1
		nu = wwWordSize(u, nu);
279 1
		mv = wwWordSize(v, mv);
280
		// u >= v?
281 1
		if (wwCmp2(u, nu, v, mv) >= 0)
282
		{
283
			// u <- u - v
284 1
			zzSubW2(u + mv, nu - mv, zzSub2(u, v, mv));
285 1
			if (sign0 != sign)
286
			{
287 1
				if (zzAdd2(da0, da, m) || wwCmp(da0, bb, m) >= 0)
288 1
					zzSub2(da0, bb, m);
289 1
				if (zzAdd2(db0, db, n) || wwCmp(db0, aa, n) >= 0)
290 1
					zzSub2(db0, aa, n);
291
			}
292 0
			else if (wwCmp(da0, da, m) >= 0)
293
			{
294 0
				ASSERT(wwCmp(db0, db, n) >= 0);
295 0
				zzSub2(da0, da, m);
296 0
				zzSub2(db0, db, n);
297
			}
298
			else
299
			{
300 0
				ASSERT(wwCmp(db0, db, n) < 0);
301 0
				zzSub(da0, da, da0, m);
302 0
				zzSub(db0, db, db0, n);
303 0
				sign0 = 1 - sign0;
304
			}
305
		}
306
		else
307
		{
308
			// v <- v - u
309 1
			zzSubW2(v + nu, mv - nu, zzSub2(v, u, nu));
310 1
			if (sign0 != sign)
311
			{
312 1
				if (zzAdd2(da, da0, m) || wwCmp(da, bb, m) >= 0)
313 1
					zzSub2(da, bb, m);
314 1
				if (zzAdd2(db, db0, n) || wwCmp(db, aa, n) >= 0)
315 1
					zzSub2(db, aa, n);
316
			}
317 0
			else if (wwCmp(da, da0, m) >= 0)
318
			{
319 0
				ASSERT(wwCmp(db, db0, n) >= 0);
320 0
				zzSub2(da, da0, m);
321 0
				zzSub2(db, db0, n);
322
			}
323
			else
324
			{
325 0
				ASSERT(wwCmp(db, db0, n) < 0);
326 0
				zzSub(da, da0, da, m);
327 0
				zzSub(db, db0, db, n);
328 0
				sign = 1 - sign;
329
			}
330
		}
331
	}
332 1
	while (!wwIsZero(u, nu));
333
	// d <- v
334 1
	wwCopy(d, v, m);
335
	// d <- d * 2^s
336 1
	wwShHi(d, W_OF_B(wwBitSize(d, m) + s), s);
337
	// очистка
338 1
	s = 0;
339 1
	sign0 = sign = 0;
340 1
	nu = mv = 0;
341
	// возврат
342 1
	return sign;
343
}
344

345 1
size_t zzExGCD_deep(size_t n, size_t m)
346
{
347 1
	return O_OF_W(3 * n + 3 * m);
348
}
349

350

351
/*
352
*******************************************************************************
353
Деление по модулю
354

355
В zzDivMod() реализован упрощенный вариант zzExGCD(): рассчитываются
356
только da0, da, причем da0 = divident (а не 1).
357

358
\todo Реализовать в zzDivMod() случай произвольного (а не только 
359
нечетного) mod.
360
*******************************************************************************
361
*/
362

363 1
void zzDivMod(word b[], const word divident[], const word a[],
364
	const word mod[], size_t n, void* stack)
365
{
366
	register size_t nu, nv;
367 1
	register int sign0 = 0, sign = 1;
368
	// переменные в stack
369 1
	word* u = (word*)stack;
370 1
	word* v = u + n;
371 1
	word* da0 = v + n;
372 1
	word* da = da0 + n;
373 1
	stack = da + n;
374
	// pre
375 1
	ASSERT(wwCmp(a, mod, n) < 0);
376 1
	ASSERT(wwCmp(divident, mod, n) < 0);
377 1
	ASSERT(wwIsDisjoint(b, mod, n));
378 1
	ASSERT(zzIsOdd(mod, n) && mod[n - 1] != 0);
379
	// da0 <- divident, da <- 0
380 1
	wwCopy(da0, divident, n);
381 1
	wwSetZero(da, n);
382
	// u <- a, v <- mod
383 1
	wwCopy(u, a, n);
384 1
	wwCopy(v, mod, n);
385 1
	nu = wwWordSize(u, n);
386 1
	nv = n;
387
	// итерации со следующими инвариантами:
388
	//	da0 * a \equiv (-1)^sign0 * divident * u \mod mod
389
	//	da * a \equiv (-1)^sign * divident * v \mod mod
390 1
	while (!wwIsZero(u, nu))
391
	{
392
		// пока u -- четное
393 1
		for (; u[0] % 2 == 0; wwShLo(u, nu, 1))
394 1
			if (da0[0] % 2 == 0)
395
				// da0 <- da0 / 2
396 1
				wwShLo(da0, n, 1);
397
			else
398
				// da0 <- (da0 + mod) / 2
399 1
				wwShLoCarry(da0, n, 1, zzAdd2(da0, mod, n));
400
		// пока v -- четное
401 1
		for (; v[0] % 2 == 0; wwShLo(v, nv, 1))
402 1
			if (da[0] % 2 == 0)
403
				// da <- da / 2
404 1
				wwShLo(da, n, 1);
405
			else
406
				// da <- (da + mod) / 2
407 1
				wwShLoCarry(da, n, 1, zzAdd2(da, mod, n));
408
		// нормализация
409 1
		nu = wwWordSize(u, nu);
410 1
		nv = wwWordSize(v, nv);
411
		// u >= v?
412 1
		if (wwCmp2(u, nu, v, nv) >= 0)
413
		{
414
			// u <- u - v
415 1
			zzSubW2(u + nv, nu - nv, zzSub2(u, v, nv));
416 1
			if (sign0 != sign)
417
			{
418 1
				if (zzAdd2(da0, da, n) || wwCmp(da0, mod, n) >= 0)
419 1
					zzSub2(da0, mod, n);
420
			}
421 0
			else if (wwCmp(da0, da, n) >= 0)
422 0
				zzSub2(da0, da, n);
423
			else
424 0
				zzSub(da0, da, da0, n),
425 0
				sign0 = 1 - sign0;
426
		}
427
		else
428
		{
429
			// v <- v - u
430 1
			zzSubW2(v + nu, nv - nu, zzSub2(v, u, nu));
431 1
			if (sign0 != sign)
432
			{
433 1
				if (zzAdd2(da, da0, n) || wwCmp(da, mod, n) >= 0)
434 1
					zzSub2(da, mod, n);
435
			}
436 0
			else if (wwCmp(da, da0, n) >= 0)
437 0
				zzSub2(da, da0, n);
438
			else
439 0
				zzSub(da, da0, da, n),
440 0
				sign = 1 - sign;
441
		}
442
	}
443
	// здесь v == (a, mod)
444
	EXPECT(wwIsW(v, nv, 1));
445
	// \gcd(a, mod) != 1? b <- 0
446 1
	if (!wwIsW(v, nv, 1))
447 0
		wwSetZero(b, n);
448
	// здесь da * a \equiv (-1)^sign * divident \mod mod
449
	// если sign == 1, то b <- mod - da, иначе b <- da
450 1
	else if (sign == 1)
451 1
		zzSub(b, mod, da, n);
452
	else
453 0
		wwCopy(b, da, n);
454
	// очистка
455 1
	sign0 = sign = 0;
456 1
	nu = nv = 0;
457
}
458

459 1
size_t zzDivMod_deep(size_t n)
460
{
461 1
	return O_OF_W(4 * n);
462
}
463

464

465
/*
466
*******************************************************************************
467
Почти обращение по модулю
468

469
В zzAlmostDivMod() реализован алгоритм Калиски [B.S.Kaliski Jr. The Montgomery 
470
inverse and its applications. IEEE Transactions on Computers, 44(8):1064–1065, 
471
1995]:
472
	u <- a
473
	da0 <- 1
474
	v <- mod
475
	da <- 0
476
	k <- 0
477
	пока (u != 0)
478
	{
479
		если (v -- четное)
480
			v <- v / 2, da0 <- da0 * 2
481
		иначе если (u -- четное)
482
			u <- u / 2, da <- da * 2
483
		иначе если (v > u)
484
			v <- (v - u) / 2, da <- da + da0, da0 <- da0 * 2
485
		иначе // если (u >= v)
486
			u <- (u - v) / 2, da0 <- da0 + da, da <- da * 2
487
		k <- k + 1
488
	}
489
	если (da >= mod)
490
		da <- da - mod
491
	da <- mod - da
492
	return (da, k)
493

494
Инварианты на итерациях zzAlmostDivMod():
495
	mod = v * da0 + u * da
496
	a * da = -v (\mod mod)
497

498
В оригинальной статье Калиски доказано, что при 0 < a < mod: 
499
-	числа da, da0 лежат в интервале [0, 2 * mod - 1];
500
-	wwBitSize(mod) <= k <= 2 * wwBitSize(mod).
501

502
\remark В [E. Savas, K. Koc. The Montgomery Modular Inverse -- Revisited. 
503
IEEE Transactions on Computers, 49(7):763–766, 2000] рассмотрен случай, когда 
504
a и mod < 2^m, причем условие a < mod может нарушаться. Доказано, что в этом 
505
случае
506
-	числа da, da0 лежат в интервале [0, 2 * mod - 1];
507
-	wwBitSize(mod) <= k <= m + wwBitSize(mod).
508

509
\todo В перечисленных статьях предполагается, что mod -- простое число.
510
Проверить, что результаты можно распространить на случай произвольного 
511
нечетного mod. Можно ли сузить интервал [0, 2 * mod - 1]?
512
*******************************************************************************
513
*/
514

515 1
size_t zzAlmostInvMod(word b[], const word a[], const word mod[], size_t n,
516
	void* stack)
517
{
518 1
	register size_t k = 0;
519
	size_t nu, nv;
520
	// переменные в stack
521 1
	word* u = (word*)stack;
522 1
	word* v = u + n;
523 1
	word* da0 = v + n;
524 1
	word* da = da0 + n + 1;
525 1
	stack = da + n + 1;
526
	// pre
527 1
	ASSERT(!wwIsZero(a, n));
528 1
	ASSERT(wwCmp(a, mod, n) < 0);
529 1
	ASSERT(wwIsDisjoint(b, mod, n));
530 1
	ASSERT(zzIsOdd(mod, n) && mod[n - 1] != 0);
531
	// da0 <- 1, da <- 0
532 1
	wwSetW(da0, n + 1, 1);
533 1
	wwSetZero(da, n + 1);
534
	// u <- a, v <- mod
535 1
	wwCopy(u, a, n);
536 1
	wwCopy(v, mod, n);
537 1
	nu = wwWordSize(u, n);
538 1
	nv = n;
539
	// пока (u != 0)
540
	do
541
	{
542
		// v -- четное?
543 1
		if (zzIsEven(v, nv))
544
		{
545 1
			wwShLo(v, nv, 1);
546 1
			nv = wwWordSize(v, nv);
547 1
			wwShHi(da0, n + 1, 1);
548
		}
549
		// u -- четное?
550 1
		else if (zzIsEven(u, nu))
551
		{
552 1
			wwShLo(u, nu, 1);
553 1
			nu = wwWordSize(u, nu);
554 1
			wwShHi(da, n + 1, 1);
555
		}
556
		// v > u?
557 1
		else if (wwCmp2(v, nv, u, nu) > 0)
558
		{
559 1
			ASSERT(nv >= nu);
560 1
			zzSubW2(v + nu, nv - nu, zzSub2(v, u, nu));
561 1
			wwShLo(v, nv, 1); 
562 1
			nv = wwWordSize(v, nv);
563 1
			zzAdd2(da, da0, n + 1);
564 1
			wwShHi(da0, n + 1, 1);
565
		}
566
		// u >= v?
567
		else
568
		{
569 1
			ASSERT(nu >= nv);
570 1
			zzSubW2(u + nv, nu - nv, zzSub2(u, v, nv));
571 1
			wwShLo(u, nu, 1); 
572 1
			nu = wwWordSize(u, nu);
573 1
			zzAdd2(da0, da, n + 1);
574 1
			wwShHi(da, n + 1, 1);
575
		}
576
		// k <- k + 1
577 1
		k = k + 1;
578
	}
579 1
	while (!wwIsZero(u, nu));
580
	// здесь v == (a, mod)
581
	EXPECT(wwIsW(v, nv, 1));
582
	// \gcd(a, mod) != 1? b <- 0
583 1
	if (!wwIsW(v, nv, 1))
584 0
		wwSetZero(b, n);
585
	// da >= mod => da -= mod
586 1
	if (wwCmp2(da, n + 1, mod, n) >= 0)
587 1
		da[n] -= zzSub2(da, mod, n);
588 1
	ASSERT(wwCmp2(da, n + 1, mod, n) < 0);
589
	// b <- mod - da
590 1
	zzNegMod(b, da, mod, n);
591
	// возврат
592 1
	return k;
593
}
594

595 1
size_t zzAlmostInvMod_deep(size_t n)
596
{
597 1
	return O_OF_W(4 * n + 2);
598
}

Read our documentation on viewing source code .

Loading