1
/*
2
*******************************************************************************
3
\file zz_mul.c
4
\brief Multiple-precision unsigned integers: multiplicative operations
5
\project bee2 [cryptographic library]
6
\author (C) Sergey Agievich [agievich@{bsu.by|gmail.com}]
7
\created 2012.04.22
8
\version 2019.06.26
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/core/word.h"
17
#include "bee2/math/ww.h"
18
#include "bee2/math/zz.h"
19
#include "zz_lcl.h"
20

21
/*
22
*******************************************************************************
23
Умножение / возведение в квадрат
24

25
\todo Возведение в квадрат за один проход (?), сначала с квадратов (?).
26
\todo Умножение Карацубы.
27
*******************************************************************************
28
*/
29

30 1
word zzMulW(word b[], const word a[], size_t n, register word w)
31
{
32 1
	register word carry = 0;
33
	register dword prod;
34
	size_t i;
35 1
	ASSERT(wwIsSameOrDisjoint(a, b, n));
36 1
	for (i = 0; i < n; ++i)
37
	{
38 1
		_MUL(prod, w, a[i]);
39 1
		prod += carry;
40 1
		b[i] = (word)prod;
41 1
		carry = (word)(prod >> B_PER_W);
42
	}
43 1
	prod = 0, w = 0;
44 1
	return carry;
45
}
46

47 1
word zzAddMulW(word b[], const word a[], size_t n, register word w)
48
{
49 1
	register word carry = 0;
50
	register dword prod;
51
	size_t i;
52 1
	ASSERT(wwIsSameOrDisjoint(a, b, n));
53 1
	for (i = 0; i < n; ++i)
54
	{
55 1
		_MUL(prod, w, a[i]);
56 1
		prod += carry;
57 1
		prod += b[i];
58 1
		b[i] = (word)prod;
59 1
		carry = (word)(prod >> B_PER_W);
60
	}
61 1
	prod = 0, w = 0;
62 1
	return carry;
63
}
64

65 1
word zzSubMulW(word b[], const word a[], size_t n, register word w)
66
{
67 1
	register word borrow = 0;
68
	register dword prod;
69
	size_t i;
70 1
	ASSERT(wwIsSameOrDisjoint(a, b, n));
71 1
	for (i = 0; i < n; ++i)
72
	{
73 1
		_MUL(prod, w, a[i]);
74 1
		prod = (dword)0 - prod;
75 1
		prod += b[i];
76 1
		prod -= borrow;
77 1
		b[i] = (word)prod;
78 1
		borrow = WORD_0 - (word)(prod >> B_PER_W);
79
	}
80 1
	prod = 0, w = 0;
81 1
	return borrow;
82
}
83

84 1
void zzMul(word c[], const word a[], size_t n, const word b[], size_t m, 
85
	void* stack)
86
{
87 1
	register word carry = 0;
88
	register dword prod;
89
	size_t i, j;
90 1
	ASSERT(wwIsDisjoint2(a, n, c, n + m));
91 1
	ASSERT(wwIsDisjoint2(b, m, c, n + m));
92 1
	wwSetZero(c, n + m);
93 1
	for (i = 0; i < n; ++i)
94
	{
95 1
		for (j = 0; j < m; ++j)
96
		{
97 1
			_MUL(prod, a[i], b[j]);
98 1
			prod += carry;
99 1
			prod += c[i + j];
100 1
			c[i + j] = (word)prod;
101 1
			carry = (word)(prod >> B_PER_W);
102
		}
103 1
		c[i + j] = carry;
104 1
		carry = 0;
105
	}
106 1
	prod = 0;
107
}
108

109 1
size_t zzMul_deep(size_t n, size_t m)
110
{
111 1
	return 0;
112
}
113

114 1
void zzSqr(word b[], const word a[], size_t n, void* stack)
115
{
116 1
	register word carry = 0;
117
	register word carry1;
118
	register dword prod;
119
	size_t i, j;
120 1
	ASSERT(wwIsDisjoint2(a, n, b, n + n));
121
	// b <- \sum_{i < j} a_i a_j B^{i + j}
122 1
	wwSetZero(b, n + n);
123 1
	for (i = 0; i < n; ++i)
124
	{
125 1
		for (j = i + 1; j < n; ++j)
126
		{
127 1
			_MUL(prod, a[i], a[j]);
128 1
			prod += carry;
129 1
			prod += b[i + j];
130 1
			b[i + j] = (word)prod;
131 1
			carry = (word)(prod >> B_PER_W);
132
		}
133 1
		b[i + j] = carry;
134 1
		carry = 0;
135
	}
136
	// b <- 2 b
137 1
	for (i = 0; i < n + n; ++i)
138
	{
139 1
		carry1 = b[i] >> (B_PER_W - 1);
140 1
		b[i] = (b[i] << 1) | carry;
141 1
		carry = carry1;
142
	}
143
	// b <- b + \sum_i a_i^2 B^{i + i}
144 1
	for (i = 0; i < n; ++i)
145
	{
146 1
		_MUL(prod, a[i], a[i]);
147 1
		prod += carry;
148 1
		prod += b[i + i];
149 1
		b[i + i] = (word)prod;
150 1
		prod >>= B_PER_W;
151 1
		prod += b[i + i + 1];
152 1
		b[i + i + 1] = (word)prod;
153 1
		carry = (word)(prod >> B_PER_W);
154
	}
155 1
	prod = 0;
156 1
	carry = carry1 = 0;
157
}
158

159 1
size_t zzSqr_deep(size_t n)
160
{
161 1
	return 0;
162
}
163

164
/*
165
*******************************************************************************
166
Деление на машинное слово
167

168
В функции zzModW2() сначала определяется значение (b = B \mod mod):
169
	r = \sum_i a[i] b^i \equiv \sum_i a[i] B^i = a \mod mod,
170
которое затем приводится \mod mod.
171
Используется следующий алгоритм:
172
	r = (r1 r0) <- 0
173
	for i = n - 1,..., 0:
174
		r <- (r1 b + r0)b + a[i]	(*)
175
	while (r1 != 0)
176
		r <- r1 b + (r0 % mod)		(**)
177
	return r0 % mod
178
После каждой итерации (*):
179
	r <= (B - 1)(1 + b + b^2) <= (B - 1)(mod^2 - mod + 1)
180
	  <= (B - 1)(B + 1) < B^2.
181
По окончании первой итерации (**):
182
	r <= (B - 1)(mod - 1) + (mod - 1) = B(mod - 1).
183
По окончании второй итерации (**):
184
	r <= (mod - 1)(mod - 1) + (mod - 1) = mod(mod - 1) < B.
185
Таким образом, r \mod mod = r0 \mod mod.
186
*******************************************************************************
187
*/
188

189 1
word zzDivW(word q[], const word a[], size_t n, register word w)
190
{
191 1
	register word r = 0;
192
	register dword divisor;
193 1
	ASSERT(w > 0);
194 1
	ASSERT(wwIsSameOrDisjoint(a, q, n));
195 1
	while (n--)
196
	{
197 1
		divisor = r;
198 1
		divisor <<= B_PER_W;
199 1
		divisor |= a[n];
200 1
		q[n] = (word)(divisor / w);
201 1
		r = (word)(divisor % w);
202
	}
203 1
	divisor = 0, w = 0;
204 1
	return r;
205
}
206

207 1
word zzModW(const word a[], size_t n, register word w)
208
{
209 1
	register word r = 0;
210
	register dword divisor;
211 1
	ASSERT(w > 0);
212 1
	ASSERT(wwIsValid(a, n));
213 1
	while (n--)
214
	{
215 1
		divisor = r;
216 1
		divisor <<= B_PER_W;
217 1
		divisor |= a[n];
218 1
		r = (word)(divisor % w);
219
	}
220 1
	divisor = 0, w = 0;
221 1
	return r;
222
}
223

224 1
word zzModW2(const word a[], size_t n, register word w)
225
{
226 1
	register word r0 = 0;
227 1
	register dword r1 = 0;
228
	register word b;
229
	// pre
230 1
	ASSERT(w > 0);
231 1
	ASSERT(w <= WORD_BIT_HALF);
232 1
	ASSERT(wwIsValid(a, n));
233
	// b <- B \mod mod
234 1
	b = (WORD_MAX - w + 1) % w;
235
	// (r1 r0) <- \sum_i a[i] b^i
236 1
	while (n--)
237
	{
238 1
		r1 *= b;
239 1
		r1 += r0;
240 1
		r1 *= b;
241 1
		r1 += a[n];
242 1
		r0 = (word)r1;
243 1
		r1 >>= B_PER_W;
244
	}
245
	// нормализация
246
#ifdef SAFE_FAST
247
	while (r1 != 0)
248
	{
249
		r1 *= b;
250
		r1 += r0 % w;
251
		r0 = (word)r1;
252
		r1 >>= B_PER_W;
253
	}
254
	r0 %= w;
255
#else
256 1
	r1 *= b;
257 1
	r1 += r0 % w;
258 1
	r0 = (word)r1;
259 1
	r1 >>= B_PER_W;
260 1
	r1 *= b;
261 1
	r1 += r0 % w;
262 1
	r0 = (word)r1 % w;
263
#endif
264
	// очистка и возврат
265 1
	r1 = 0, b = w = 0;
266 1
	return r0;
267
}
268

269
/*
270
*******************************************************************************
271
Общее деление
272

273
\opt При делении слов определять остаток по частному: (/,*) вместо (/, %).
274

275
\todo Убрать ограничение n >= m в zzDiv().
276

277
\todo T. Jabelean. An Algorithm for exact division. J. of Symb. Computations, 
278
15 (2): 169-180, 1993.
279

280
\todo: В zzMod() отказаться от divident.
281

282
В функциях zzDiv(), zzMod() делимое a = a[n - 1]...a[0] и делитель
283
b = b[m - 1]...b[0] предварительно нормализуются:
284
	a = a[n]...a[0] <- a * 2^shift;
285
	b = b[m - 1]...b[0] <- b * 2^shift.
286
Здесь shift --- минимальное натуральное т.ч. старший бит b[m - 1] * 2^shift
287
равняется 1.
288

289
Деление выполняется по алгоритму 14.20 из [Menezes A., van Oorschot P.,
290
Vanstone S. Handbook of Applied Cryptography]:
291
	for i = n, n - 1, ...., m:
292
		if a[i] == b[m - 1]											(#)
293
			q[i - m] <- B - 1
294
		else
295
			q[i - m] <- a[i]a[i - 1] div b[m - 1]
296
		while (q[i - m] * b[m - 1]b[m - 2] > a[i]a[i - 1]a[i - 2])	(##)
297
			q[i - m]--
298
		a <- a - q[i - m] * b * B^{i - m}
299
		if (a < 0)
300
			a += b * B^{i - m}, q[i - m]--
301
	return q = q[n - m]...q[0] --- частное и a --- остаток
302

303
В реализации вместо (#):
304
	d <- a[i]a[i - 1] div b[m - 1]
305
	if d >= B
306
		d <- B - 1
307
	q[i - m] <- d
308

309
\opt Если a[i] == b[m - 1] в (#), то цикл (##) можно не выполнять:
310
	q[i - m] * b[m - 1]b[m - 2] <=
311
		(B - 1) * (a[i] * B + (B - 1)) =
312
		B^2 * a[i] + B^2 - 1 - a[i] * B < a[i]a[i - 1]a[i - 2]
313

314
\opt Если известен остаток d = a[i]a[i - 1] mod b[m - 1], то (##) можно
315
	заменить на
316
		while (q[i - m] * b[m - 2] > d * B + a[i - 2])
317
			q[i - m]--, d += b[m - 1]
318
*******************************************************************************
319
*/
320

321 1
void zzDiv(word q[], word r[], const word a[], size_t n, const word b[],
322
	size_t m, void* stack)
323
{
324
	register dword dividentHi;
325
	register word borrow;
326
	register size_t shift;
327
	size_t i;
328
	// переменные в stack
329
	word* divident;		/*< нормализованное делимое (n + 1 слово) */
330
	word* divisor;		/*< нормализованный делитель (m слов) */
331
	word* mul;			/*< вспомогательное произведение (3 слова) */
332
	// pre
333 1
	ASSERT(n >= m);
334 1
	ASSERT(wwIsValid(a, n) && wwIsValid(b, m));
335 1
	ASSERT(m > 0 && b[m - 1] > 0);
336 1
	ASSERT(wwIsDisjoint2(q, n + 1 - m, r, m));
337 1
	ASSERT(a == r || wwIsDisjoint2(a, n, r, m));
338
	// a < b?
339 1
	if (wwCmp2(a, n, b, m) < 0)
340
	{
341
		// q <- 0, r <- a
342 1
		wwSetZero(q, n - m + 1);
343 1
		wwCopy(r, a, m);
344 1
		return;
345
	}
346
	// делим на одноразрядное число?
347 1
	if (m == 1)
348
	{
349 1
		r[0] = zzDivW(q, a, n, b[0]);
350 1
		return;
351
	}
352
	// резервируем переменные в stack
353 1
	divident = (word*)stack;
354 1
	divisor = divident + n + 1;
355 1
	mul = divisor + m;
356 1
	stack = mul + 3;
357
	// divident <- a
358 1
	wwCopy(divident, a, n);
359 1
	divident[n] = 0;
360
	// divisor <- b
361 1
	wwCopy(divisor, b, m);
362
	// нормализация
363 1
	shift = wordCLZ(b[m - 1]);
364 1
	wwShHi(divident, n + 1, shift);
365 1
	wwShHi(divisor, m, shift);
366
	// цикл по разрядам делимого
367 1
	for (i = n; i >= m; --i)
368
	{
369
		// вычислить пробное частное
370 1
		dividentHi = divident[i];
371 1
		dividentHi <<= B_PER_W;
372 1
		dividentHi |= divident[i - 1];
373 1
		dividentHi /= divisor[m - 1];
374 1
		if (dividentHi > WORD_MAX)
375 0
			q[i - m] = WORD_MAX;
376
		else
377 1
			q[i - m] = (word)dividentHi;
378
		// уточнить пробное частное
379 1
		wwCopy(mul, divisor + m - 2, 2);
380 1
		mul[2] = zzMulW(mul, mul, 2, q[i - m]);
381 1
		while (wwCmp2(mul, 3, divident + i - 2, 3) > 0)
382
		{
383 1
			q[i - m]--;
384 1
			mul[2] -= zzSub2(mul, divisor + m - 2, 2);
385
		}
386
		// учесть пробное частное
387 1
		borrow = zzSubMulW(divident + i - m, divisor, m, q[i - m]);
388 1
		divident[i] -= borrow;
389 1
		if (divident[i] > (word)~borrow)
390
		{
391
			// окончательно подправить пробное частное
392 0
			q[i - m]--;
393
			// корректирующее сложение
394 0
			divident[i] += zzAdd2(divident + i - m, divisor, m);
395
		}
396
	}
397
	// денормализация
398 1
	wwShLo(divident, n + 1, shift);
399
	// сохранить остаток
400 1
	wwCopy(r, divident, m);
401
	// очистить регистровые переменные
402 1
	shift = 0;
403 1
	borrow = 0;
404 1
	dividentHi = 0;
405
}
406

407 1
size_t zzDiv_deep(size_t n, size_t m)
408
{
409 1
	return O_OF_W(n + m + 4);
410
}
411

412 1
void zzMod(word r[], const word a[], size_t n, const word b[], size_t m, void* stack)
413
{
414
	register dword dividentHi;
415
	register word temp;
416
	register size_t shift;
417
	size_t i;
418
	// переменные в stack
419
	word* divident;		/*< нормализованное делимое (n + 1 слово) */
420
	word* divisor;		/*< нормализованный делитель (m слов) */
421
	word* mul;			/*< вспомогательное произведение (3 слова) */
422
	// pre
423 1
	ASSERT(wwIsValid(a, n) && wwIsValid(b, m));
424 1
	ASSERT(m > 0 && b[m - 1] > 0);
425 1
	ASSERT(a == r || wwIsDisjoint2(a, n, r, m));
426
	// a < b?
427 1
	if (wwCmp2(a, n, b, m) < 0)
428
	{
429
		// r <- a
430 1
		if (n < m)
431 0
			wwSetZero(r + n, m - n), m = n;
432 1
		wwCopy(r, a, m);
433 1
		return;
434
	}
435
	// делим на одноразрядное число?
436 1
	if (m == 1)
437
	{
438 1
		r[0] = zzModW(a, n, b[0]);
439 1
		return;
440
	}
441
	// резервируем переменные в stack
442 1
	divident = (word*)stack;
443 1
	divisor = divident + n + 1;
444 1
	mul = divisor + m;
445 1
	stack = mul + 3;
446
	// divident <- a
447 1
	wwCopy(divident, a, n);
448 1
	divident[n] = 0;
449
	// divisor <- b
450 1
	wwCopy(divisor, b, m);
451
	// нормализация
452 1
	shift = wordCLZ(b[m - 1]);
453 1
	wwShHi(divident, n + 1, shift);
454 1
	wwShHi(divisor, m, shift);
455
	// цикл по разрядам делимого
456 1
	for (i = n; i >= m; --i)
457
	{
458
		// вычислить пробное частное
459 1
		dividentHi = divident[i];
460 1
		dividentHi <<= B_PER_W;
461 1
		dividentHi |= divident[i - 1];
462 1
		dividentHi /= divisor[m - 1];
463 1
		if (dividentHi > WORD_MAX)
464 1
			temp = WORD_MAX;
465
		else
466 1
			temp = (word)dividentHi;
467
		// уточнить пробное частное
468 1
		wwCopy(mul, divisor + m - 2, 2);
469 1
		mul[2] = zzMulW(mul, mul, 2, temp);
470 1
		while (wwCmp2(mul, 3, divident + i - 2, 3) > 0)
471
		{
472 1
			temp--;
473 1
			mul[2] -= zzSub2(mul, divisor + m - 2, 2);
474
		}
475
		// учесть пробное частное
476 1
		temp = zzSubMulW(divident + i - m, divisor, m, temp);
477 1
		divident[i] -= temp;
478 1
		if (divident[i] > (word)~temp)
479
			// корректирующее сложение
480 1
			divident[i] += zzAdd2(divident + i - m, divisor, m);
481
	}
482
	// денормализация
483 1
	wwShLo(divident, n + 1, shift);
484
	// сохранить остаток
485 1
	wwCopy(r, divident, m);
486
	// очистить регистровые переменные
487 1
	shift = 0;
488 1
	temp = 0;
489 1
	dividentHi = 0;
490
}
491

492 1
size_t zzMod_deep(size_t n, size_t m)
493
{
494 1
	return O_OF_W(n + m + 4);
495
}

Read our documentation on viewing source code .

Loading