aboutsummaryrefslogtreecommitdiffstats
path: root/test/monniaux/BearSSL/src/int/i15_moddiv.c
blob: 45af756d346039fafe4cbe57164b6c1bb0a11e9e (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
/*
 * Copyright (c) 2018 Thomas Pornin <pornin@bolet.org>
 *
 * Permission is hereby granted, free of charge, to any person obtaining 
 * a copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject to
 * the following conditions:
 *
 * The above copyright notice and this permission notice shall be 
 * included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 
 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
 * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
 * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#include "inner.h"

/*
 * In this file, we handle big integers with a custom format, i.e.
 * without the usual one-word header. Value is split into 15-bit words,
 * each stored in a 16-bit slot (top bit is zero) in little-endian
 * order. The length (in words) is provided explicitly. In some cases,
 * the value can be negative (using two's complement representation). In
 * some cases, the top word is allowed to have a 16th bit.
 */

/*
 * Negate big integer conditionally. The value consists of 'len' words,
 * with 15 bits in each word (the top bit of each word should be 0,
 * except possibly for the last word). If 'ctl' is 1, the negation is
 * computed; otherwise, if 'ctl' is 0, then the value is unchanged.
 */
static void
cond_negate(uint16_t *a, size_t len, uint32_t ctl)
{
	size_t k;
	uint32_t cc, xm;

	cc = ctl;
	xm = 0x7FFF & -ctl;
	for (k = 0; k < len; k ++) {
		uint32_t aw;

		aw = a[k];
		aw = (aw ^ xm) + cc;
		a[k] = aw & 0x7FFF;
		cc = (aw >> 15) & 1;
	}
}

/*
 * Finish modular reduction. Rules on input parameters:
 *
 *   if neg = 1, then -m <= a < 0
 *   if neg = 0, then 0 <= a < 2*m
 *
 * If neg = 0, then the top word of a[] may use 16 bits.
 *
 * Also, modulus m must be odd.
 */
static void
finish_mod(uint16_t *a, size_t len, const uint16_t *m, uint32_t neg)
{
	size_t k;
	uint32_t cc, xm, ym;

	/*
	 * First pass: compare a (assumed nonnegative) with m.
	 */
	cc = 0;
	for (k = 0; k < len; k ++) {
		uint32_t aw, mw;

		aw = a[k];
		mw = m[k];
		cc = (aw - mw - cc) >> 31;
	}

	/*
	 * At this point:
	 *   if neg = 1, then we must add m (regardless of cc)
	 *   if neg = 0 and cc = 0, then we must subtract m
	 *   if neg = 0 and cc = 1, then we must do nothing
	 */
	xm = 0x7FFF & -neg;
	ym = -(neg | (1 - cc));
	cc = neg;
	for (k = 0; k < len; k ++) {
		uint32_t aw, mw;

		aw = a[k];
		mw = (m[k] ^ xm) & ym;
		aw = aw - mw - cc;
		a[k] = aw & 0x7FFF;
		cc = aw >> 31;
	}
}

/*
 * Compute:
 *   a <- (a*pa+b*pb)/(2^15)
 *   b <- (a*qa+b*qb)/(2^15)
 * The division is assumed to be exact (i.e. the low word is dropped).
 * If the final a is negative, then it is negated. Similarly for b.
 * Returned value is the combination of two bits:
 *   bit 0: 1 if a had to be negated, 0 otherwise
 *   bit 1: 1 if b had to be negated, 0 otherwise
 *
 * Factors pa, pb, qa and qb must be at most 2^15 in absolute value.
 * Source integers a and b must be nonnegative; top word is not allowed
 * to contain an extra 16th bit.
 */
static uint32_t
co_reduce(uint16_t *a, uint16_t *b, size_t len,
	int32_t pa, int32_t pb, int32_t qa, int32_t qb)
{
	size_t k;
	int32_t cca, ccb;
	uint32_t nega, negb;

	cca = 0;
	ccb = 0;
	for (k = 0; k < len; k ++) {
		uint32_t wa, wb, za, zb;
		uint16_t tta, ttb;

		/*
		 * Since:
		 *   |pa| <= 2^15
		 *   |pb| <= 2^15
		 *   0 <= wa <= 2^15 - 1
		 *   0 <= wb <= 2^15 - 1
		 *   |cca| <= 2^16 - 1
		 * Then:
		 *   |za| <= (2^15-1)*(2^16) + (2^16-1) = 2^31 - 1
		 *
		 * Thus, the new value of cca is such that |cca| <= 2^16 - 1.
		 * The same applies to ccb.
		 */
		wa = a[k];
		wb = b[k];
		za = wa * (uint32_t)pa + wb * (uint32_t)pb + (uint32_t)cca;
		zb = wa * (uint32_t)qa + wb * (uint32_t)qb + (uint32_t)ccb;
		if (k > 0) {
			a[k - 1] = za & 0x7FFF;
			b[k - 1] = zb & 0x7FFF;
		}
		tta = za >> 15;
		ttb = zb >> 15;
		cca = *(int16_t *)&tta;
		ccb = *(int16_t *)&ttb;
	}
	a[len - 1] = (uint16_t)cca;
	b[len - 1] = (uint16_t)ccb;
	nega = (uint32_t)cca >> 31;
	negb = (uint32_t)ccb >> 31;
	cond_negate(a, len, nega);
	cond_negate(b, len, negb);
	return nega | (negb << 1);
}

/*
 * Compute:
 *   a <- (a*pa+b*pb)/(2^15) mod m
 *   b <- (a*qa+b*qb)/(2^15) mod m
 *
 * m0i is equal to -1/m[0] mod 2^15.
 *
 * Factors pa, pb, qa and qb must be at most 2^15 in absolute value.
 * Source integers a and b must be nonnegative; top word is not allowed
 * to contain an extra 16th bit.
 */
static void
co_reduce_mod(uint16_t *a, uint16_t *b, size_t len,
	int32_t pa, int32_t pb, int32_t qa, int32_t qb,
	const uint16_t *m, uint16_t m0i)
{
	size_t k;
	int32_t cca, ccb, fa, fb;

	cca = 0;
	ccb = 0;
	fa = ((a[0] * (uint32_t)pa + b[0] * (uint32_t)pb) * m0i) & 0x7FFF;
	fb = ((a[0] * (uint32_t)qa + b[0] * (uint32_t)qb) * m0i) & 0x7FFF;
	for (k = 0; k < len; k ++) {
		uint32_t wa, wb, za, zb;
		uint32_t tta, ttb;

		/*
		 * In this loop, carries 'cca' and 'ccb' always fit on
		 * 17 bits (in absolute value).
		 */
		wa = a[k];
		wb = b[k];
		za = wa * (uint32_t)pa + wb * (uint32_t)pb
			+ m[k] * (uint32_t)fa + (uint32_t)cca;
		zb = wa * (uint32_t)qa + wb * (uint32_t)qb
			+ m[k] * (uint32_t)fb + (uint32_t)ccb;
		if (k > 0) {
			a[k - 1] = za & 0x7FFF;
			b[k - 1] = zb & 0x7FFF;
		}

		/*
		 * The XOR-and-sub construction below does an arithmetic
		 * right shift in a portable way (technically, right-shifting
		 * a negative signed value is implementation-defined in C).
		 */
#define M   ((uint32_t)1 << 16)
		tta = za >> 15;
		ttb = zb >> 15;
		tta = (tta ^ M) - M;
		ttb = (ttb ^ M) - M;
		cca = *(int32_t *)&tta;
		ccb = *(int32_t *)&ttb;
#undef M
	}
	a[len - 1] = (uint32_t)cca;
	b[len - 1] = (uint32_t)ccb;

	/*
	 * At this point:
	 *   -m <= a < 2*m
	 *   -m <= b < 2*m
	 * (this is a case of Montgomery reduction)
	 * The top word of 'a' and 'b' may have a 16-th bit set.
	 * We may have to add or subtract the modulus.
	 */
	finish_mod(a, len, m, (uint32_t)cca >> 31);
	finish_mod(b, len, m, (uint32_t)ccb >> 31);
}

/* see inner.h */
uint32_t
br_i15_moddiv(uint16_t *x, const uint16_t *y, const uint16_t *m, uint16_t m0i,
	uint16_t *t)
{
	/*
	 * Algorithm is an extended binary GCD. We maintain four values
	 * a, b, u and v, with the following invariants:
	 *
	 *   a * x = y * u mod m
	 *   b * x = y * v mod m
	 *
	 * Starting values are:
	 *
	 *   a = y
	 *   b = m
	 *   u = x
	 *   v = 0
	 *
	 * The formal definition of the algorithm is a sequence of steps:
	 *
	 *   - If a is even, then a <- a/2 and u <- u/2 mod m.
	 *   - Otherwise, if b is even, then b <- b/2 and v <- v/2 mod m.
	 *   - Otherwise, if a > b, then a <- (a-b)/2 and u <- (u-v)/2 mod m.
	 *   - Otherwise, b <- (b-a)/2 and v <- (v-u)/2 mod m.
	 *
	 * Algorithm stops when a = b. At that point, they both are equal
	 * to GCD(y,m); the modular division succeeds if that value is 1.
	 * The result of the modular division is then u (or v: both are
	 * equal at that point).
	 *
	 * Each step makes either a or b shrink by at least one bit; hence,
	 * if m has bit length k bits, then 2k-2 steps are sufficient.
	 *
	 *
	 * Though complexity is quadratic in the size of m, the bit-by-bit
	 * processing is not very efficient. We can speed up processing by
	 * remarking that the decisions are taken based only on observation
	 * of the top and low bits of a and b.
	 *
	 * In the loop below, at each iteration, we use the two top words
	 * of a and b, and the low words of a and b, to compute reduction
	 * parameters pa, pb, qa and qb such that the new values for a
	 * and b are:
	 *
	 *   a' = (a*pa + b*pb) / (2^15)
	 *   b' = (a*qa + b*qb) / (2^15)
	 *
	 * the division being exact.
	 *
	 * Since the choices are based on the top words, they may be slightly
	 * off, requiring an optional correction: if a' < 0, then we replace
	 * pa with -pa, and pb with -pb. The total length of a and b is
	 * thus reduced by at least 14 bits at each iteration.
	 *
	 * The stopping conditions are still the same, though: when a
	 * and b become equal, they must be both odd (since m is odd,
	 * the GCD cannot be even), therefore the next operation is a
	 * subtraction, and one of the values becomes 0. At that point,
	 * nothing else happens, i.e. one value is stuck at 0, and the
	 * other one is the GCD.
	 */
	size_t len, k;
	uint16_t *a, *b, *u, *v;
	uint32_t num, r;

	len = (m[0] + 15) >> 4;
	a = t;
	b = a + len;
	u = x + 1;
	v = b + len;
	memcpy(a, y + 1, len * sizeof *y);
	memcpy(b, m + 1, len * sizeof *m);
	memset(v, 0, len * sizeof *v);

	/*
	 * Loop below ensures that a and b are reduced by some bits each,
	 * for a total of at least 14 bits.
	 */
	for (num = ((m[0] - (m[0] >> 4)) << 1) + 14; num >= 14; num -= 14) {
		size_t j;
		uint32_t c0, c1;
		uint32_t a0, a1, b0, b1;
		uint32_t a_hi, b_hi, a_lo, b_lo;
		int32_t pa, pb, qa, qb;
		int i;

		/*
		 * Extract top words of a and b. If j is the highest
		 * index >= 1 such that a[j] != 0 or b[j] != 0, then we want
		 * (a[j] << 15) + a[j - 1], and (b[j] << 15) + b[j - 1].
		 * If a and b are down to one word each, then we use a[0]
		 * and b[0].
		 */
		c0 = (uint32_t)-1;
		c1 = (uint32_t)-1;
		a0 = 0;
		a1 = 0;
		b0 = 0;
		b1 = 0;
		j = len;
		while (j -- > 0) {
			uint32_t aw, bw;

			aw = a[j];
			bw = b[j];
			a0 ^= (a0 ^ aw) & c0;
			a1 ^= (a1 ^ aw) & c1;
			b0 ^= (b0 ^ bw) & c0;
			b1 ^= (b1 ^ bw) & c1;
			c1 = c0;
			c0 &= (((aw | bw) + 0xFFFF) >> 16) - (uint32_t)1;
		}

		/*
		 * If c1 = 0, then we grabbed two words for a and b.
		 * If c1 != 0 but c0 = 0, then we grabbed one word. It
		 * is not possible that c1 != 0 and c0 != 0, because that
		 * would mean that both integers are zero.
		 */
		a1 |= a0 & c1;
		a0 &= ~c1;
		b1 |= b0 & c1;
		b0 &= ~c1;
		a_hi = (a0 << 15) + a1;
		b_hi = (b0 << 15) + b1;
		a_lo = a[0];
		b_lo = b[0];

		/*
		 * Compute reduction factors:
		 *
		 *   a' = a*pa + b*pb
		 *   b' = a*qa + b*qb
		 *
		 * such that a' and b' are both multiple of 2^15, but are
		 * only marginally larger than a and b.
		 */
		pa = 1;
		pb = 0;
		qa = 0;
		qb = 1;
		for (i = 0; i < 15; i ++) {
			/*
			 * At each iteration:
			 *
			 *   a <- (a-b)/2 if: a is odd, b is odd, a_hi > b_hi
			 *   b <- (b-a)/2 if: a is odd, b is odd, a_hi <= b_hi
			 *   a <- a/2 if: a is even
			 *   b <- b/2 if: a is odd, b is even
			 *
			 * We multiply a_lo and b_lo by 2 at each
			 * iteration, thus a division by 2 really is a
			 * non-multiplication by 2.
			 */
			uint32_t r, oa, ob, cAB, cBA, cA;

			/*
			 * cAB = 1 if b must be subtracted from a
			 * cBA = 1 if a must be subtracted from b
			 * cA = 1 if a is divided by 2, 0 otherwise
			 *
			 * Rules:
			 *
			 *   cAB and cBA cannot be both 1.
			 *   if a is not divided by 2, b is.
			 */
			r = GT(a_hi, b_hi);
			oa = (a_lo >> i) & 1;
			ob = (b_lo >> i) & 1;
			cAB = oa & ob & r;
			cBA = oa & ob & NOT(r);
			cA = cAB | NOT(oa);

			/*
			 * Conditional subtractions.
			 */
			a_lo -= b_lo & -cAB;
			a_hi -= b_hi & -cAB;
			pa -= qa & -(int32_t)cAB;
			pb -= qb & -(int32_t)cAB;
			b_lo -= a_lo & -cBA;
			b_hi -= a_hi & -cBA;
			qa -= pa & -(int32_t)cBA;
			qb -= pb & -(int32_t)cBA;

			/*
			 * Shifting.
			 */
			a_lo += a_lo & (cA - 1);
			pa += pa & ((int32_t)cA - 1);
			pb += pb & ((int32_t)cA - 1);
			a_hi ^= (a_hi ^ (a_hi >> 1)) & -cA;
			b_lo += b_lo & -cA;
			qa += qa & -(int32_t)cA;
			qb += qb & -(int32_t)cA;
			b_hi ^= (b_hi ^ (b_hi >> 1)) & (cA - 1);
		}

		/*
		 * Replace a and b with new values a' and b'.
		 */
		r = co_reduce(a, b, len, pa, pb, qa, qb);
		pa -= pa * ((r & 1) << 1);
		pb -= pb * ((r & 1) << 1);
		qa -= qa * (r & 2);
		qb -= qb * (r & 2);
		co_reduce_mod(u, v, len, pa, pb, qa, qb, m + 1, m0i);
	}

	/*
	 * Now one of the arrays should be 0, and the other contains
	 * the GCD. If a is 0, then u is 0 as well, and v contains
	 * the division result.
	 * Result is correct if and only if GCD is 1.
	 */
	r = (a[0] | b[0]) ^ 1;
	u[0] |= v[0];
	for (k = 1; k < len; k ++) {
		r |= a[k] | b[k];
		u[k] |= v[k];
	}
	return EQ0(r);
}