|
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
|
|
}
|
agievich