1/*
2 * Copyright (c) 2007, 2011, Oracle and/or its affiliates. All rights reserved.
3 * Use is subject to license terms.
4 *
5 * This library is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU Lesser General Public
7 * License as published by the Free Software Foundation; either
8 * version 2.1 of the License, or (at your option) any later version.
9 *
10 * This library is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * Lesser General Public License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public License
16 * along with this library; if not, write to the Free Software Foundation,
17 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 *
19 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
20 * or visit www.oracle.com if you need additional information or have any
21 * questions.
22 */
23
24/* *********************************************************************
25 *
26 * The Original Code is the Multi-precision Binary Polynomial Arithmetic Library.
27 *
28 * The Initial Developer of the Original Code is
29 * Sun Microsystems, Inc.
30 * Portions created by the Initial Developer are Copyright (C) 2003
31 * the Initial Developer. All Rights Reserved.
32 *
33 * Contributor(s):
34 * Sheueling Chang Shantz <sheueling.chang@sun.com> and
35 * Douglas Stebila <douglas@stebila.ca> of Sun Laboratories.
36 *
37 *********************************************************************** */
38
39#include "mp_gf2m.h"
40#include "mp_gf2m-priv.h"
41#include "mplogic.h"
42#include "mpi-priv.h"
43
44const mp_digit mp_gf2m_sqr_tb[16] =
45{
46 0, 1, 4, 5, 16, 17, 20, 21,
47 64, 65, 68, 69, 80, 81, 84, 85
48};
49
50/* Multiply two binary polynomials mp_digits a, b.
51 * Result is a polynomial with degree < 2 * MP_DIGIT_BITS - 1.
52 * Output in two mp_digits rh, rl.
53 */
54#if MP_DIGIT_BITS == 32
55void
56s_bmul_1x1(mp_digit *rh, mp_digit *rl, const mp_digit a, const mp_digit b)
57{
58 register mp_digit h, l, s;
59 mp_digit tab[8], top2b = a >> 30;
60 register mp_digit a1, a2, a4;
61
62 a1 = a & (0x3FFFFFFF); a2 = a1 << 1; a4 = a2 << 1;
63
64 tab[0] = 0; tab[1] = a1; tab[2] = a2; tab[3] = a1^a2;
65 tab[4] = a4; tab[5] = a1^a4; tab[6] = a2^a4; tab[7] = a1^a2^a4;
66
67 s = tab[b & 0x7]; l = s;
68 s = tab[b >> 3 & 0x7]; l ^= s << 3; h = s >> 29;
69 s = tab[b >> 6 & 0x7]; l ^= s << 6; h ^= s >> 26;
70 s = tab[b >> 9 & 0x7]; l ^= s << 9; h ^= s >> 23;
71 s = tab[b >> 12 & 0x7]; l ^= s << 12; h ^= s >> 20;
72 s = tab[b >> 15 & 0x7]; l ^= s << 15; h ^= s >> 17;
73 s = tab[b >> 18 & 0x7]; l ^= s << 18; h ^= s >> 14;
74 s = tab[b >> 21 & 0x7]; l ^= s << 21; h ^= s >> 11;
75 s = tab[b >> 24 & 0x7]; l ^= s << 24; h ^= s >> 8;
76 s = tab[b >> 27 & 0x7]; l ^= s << 27; h ^= s >> 5;
77 s = tab[b >> 30 ]; l ^= s << 30; h ^= s >> 2;
78
79 /* compensate for the top two bits of a */
80
81 if (top2b & 01) { l ^= b << 30; h ^= b >> 2; }
82 if (top2b & 02) { l ^= b << 31; h ^= b >> 1; }
83
84 *rh = h; *rl = l;
85}
86#else
87void
88s_bmul_1x1(mp_digit *rh, mp_digit *rl, const mp_digit a, const mp_digit b)
89{
90 register mp_digit h, l, s;
91 mp_digit tab[16], top3b = a >> 61;
92 register mp_digit a1, a2, a4, a8;
93
94 a1 = a & (0x1FFFFFFFFFFFFFFFULL); a2 = a1 << 1;
95 a4 = a2 << 1; a8 = a4 << 1;
96 tab[ 0] = 0; tab[ 1] = a1; tab[ 2] = a2; tab[ 3] = a1^a2;
97 tab[ 4] = a4; tab[ 5] = a1^a4; tab[ 6] = a2^a4; tab[ 7] = a1^a2^a4;
98 tab[ 8] = a8; tab[ 9] = a1^a8; tab[10] = a2^a8; tab[11] = a1^a2^a8;
99 tab[12] = a4^a8; tab[13] = a1^a4^a8; tab[14] = a2^a4^a8; tab[15] = a1^a2^a4^a8;
100
101 s = tab[b & 0xF]; l = s;
102 s = tab[b >> 4 & 0xF]; l ^= s << 4; h = s >> 60;
103 s = tab[b >> 8 & 0xF]; l ^= s << 8; h ^= s >> 56;
104 s = tab[b >> 12 & 0xF]; l ^= s << 12; h ^= s >> 52;
105 s = tab[b >> 16 & 0xF]; l ^= s << 16; h ^= s >> 48;
106 s = tab[b >> 20 & 0xF]; l ^= s << 20; h ^= s >> 44;
107 s = tab[b >> 24 & 0xF]; l ^= s << 24; h ^= s >> 40;
108 s = tab[b >> 28 & 0xF]; l ^= s << 28; h ^= s >> 36;
109 s = tab[b >> 32 & 0xF]; l ^= s << 32; h ^= s >> 32;
110 s = tab[b >> 36 & 0xF]; l ^= s << 36; h ^= s >> 28;
111 s = tab[b >> 40 & 0xF]; l ^= s << 40; h ^= s >> 24;
112 s = tab[b >> 44 & 0xF]; l ^= s << 44; h ^= s >> 20;
113 s = tab[b >> 48 & 0xF]; l ^= s << 48; h ^= s >> 16;
114 s = tab[b >> 52 & 0xF]; l ^= s << 52; h ^= s >> 12;
115 s = tab[b >> 56 & 0xF]; l ^= s << 56; h ^= s >> 8;
116 s = tab[b >> 60 ]; l ^= s << 60; h ^= s >> 4;
117
118 /* compensate for the top three bits of a */
119
120 if (top3b & 01) { l ^= b << 61; h ^= b >> 3; }
121 if (top3b & 02) { l ^= b << 62; h ^= b >> 2; }
122 if (top3b & 04) { l ^= b << 63; h ^= b >> 1; }
123
124 *rh = h; *rl = l;
125}
126#endif
127
128/* Compute xor-multiply of two binary polynomials (a1, a0) x (b1, b0)
129 * result is a binary polynomial in 4 mp_digits r[4].
130 * The caller MUST ensure that r has the right amount of space allocated.
131 */
132void
133s_bmul_2x2(mp_digit *r, const mp_digit a1, const mp_digit a0, const mp_digit b1,
134 const mp_digit b0)
135{
136 mp_digit m1, m0;
137 /* r[3] = h1, r[2] = h0; r[1] = l1; r[0] = l0 */
138 s_bmul_1x1(r+3, r+2, a1, b1);
139 s_bmul_1x1(r+1, r, a0, b0);
140 s_bmul_1x1(&m1, &m0, a0 ^ a1, b0 ^ b1);
141 /* Correction on m1 ^= l1 ^ h1; m0 ^= l0 ^ h0; */
142 r[2] ^= m1 ^ r[1] ^ r[3]; /* h0 ^= m1 ^ l1 ^ h1; */
143 r[1] = r[3] ^ r[2] ^ r[0] ^ m1 ^ m0; /* l1 ^= l0 ^ h0 ^ m0; */
144}
145
146/* Compute xor-multiply of two binary polynomials (a2, a1, a0) x (b2, b1, b0)
147 * result is a binary polynomial in 6 mp_digits r[6].
148 * The caller MUST ensure that r has the right amount of space allocated.
149 */
150void
151s_bmul_3x3(mp_digit *r, const mp_digit a2, const mp_digit a1, const mp_digit a0,
152 const mp_digit b2, const mp_digit b1, const mp_digit b0)
153{
154 mp_digit zm[4];
155
156 s_bmul_1x1(r+5, r+4, a2, b2); /* fill top 2 words */
157 s_bmul_2x2(zm, a1, a2^a0, b1, b2^b0); /* fill middle 4 words */
158 s_bmul_2x2(r, a1, a0, b1, b0); /* fill bottom 4 words */
159
160 zm[3] ^= r[3];
161 zm[2] ^= r[2];
162 zm[1] ^= r[1] ^ r[5];
163 zm[0] ^= r[0] ^ r[4];
164
165 r[5] ^= zm[3];
166 r[4] ^= zm[2];
167 r[3] ^= zm[1];
168 r[2] ^= zm[0];
169}
170
171/* Compute xor-multiply of two binary polynomials (a3, a2, a1, a0) x (b3, b2, b1, b0)
172 * result is a binary polynomial in 8 mp_digits r[8].
173 * The caller MUST ensure that r has the right amount of space allocated.
174 */
175void s_bmul_4x4(mp_digit *r, const mp_digit a3, const mp_digit a2, const mp_digit a1,
176 const mp_digit a0, const mp_digit b3, const mp_digit b2, const mp_digit b1,
177 const mp_digit b0)
178{
179 mp_digit zm[4];
180
181 s_bmul_2x2(r+4, a3, a2, b3, b2); /* fill top 4 words */
182 s_bmul_2x2(zm, a3^a1, a2^a0, b3^b1, b2^b0); /* fill middle 4 words */
183 s_bmul_2x2(r, a1, a0, b1, b0); /* fill bottom 4 words */
184
185 zm[3] ^= r[3] ^ r[7];
186 zm[2] ^= r[2] ^ r[6];
187 zm[1] ^= r[1] ^ r[5];
188 zm[0] ^= r[0] ^ r[4];
189
190 r[5] ^= zm[3];
191 r[4] ^= zm[2];
192 r[3] ^= zm[1];
193 r[2] ^= zm[0];
194}
195
196/* Compute addition of two binary polynomials a and b,
197 * store result in c; c could be a or b, a and b could be equal;
198 * c is the bitwise XOR of a and b.
199 */
200mp_err
201mp_badd(const mp_int *a, const mp_int *b, mp_int *c)
202{
203 mp_digit *pa, *pb, *pc;
204 mp_size ix;
205 mp_size used_pa, used_pb;
206 mp_err res = MP_OKAY;
207
208 /* Add all digits up to the precision of b. If b had more
209 * precision than a initially, swap a, b first
210 */
211 if (MP_USED(a) >= MP_USED(b)) {
212 pa = MP_DIGITS(a);
213 pb = MP_DIGITS(b);
214 used_pa = MP_USED(a);
215 used_pb = MP_USED(b);
216 } else {
217 pa = MP_DIGITS(b);
218 pb = MP_DIGITS(a);
219 used_pa = MP_USED(b);
220 used_pb = MP_USED(a);
221 }
222
223 /* Make sure c has enough precision for the output value */
224 MP_CHECKOK( s_mp_pad(c, used_pa) );
225
226 /* Do word-by-word xor */
227 pc = MP_DIGITS(c);
228 for (ix = 0; ix < used_pb; ix++) {
229 (*pc++) = (*pa++) ^ (*pb++);
230 }
231
232 /* Finish the rest of digits until we're actually done */
233 for (; ix < used_pa; ++ix) {
234 *pc++ = *pa++;
235 }
236
237 MP_USED(c) = used_pa;
238 MP_SIGN(c) = ZPOS;
239 s_mp_clamp(c);
240
241CLEANUP:
242 return res;
243}
244
245#define s_mp_div2(a) MP_CHECKOK( mpl_rsh((a), (a), 1) );
246
247/* Compute binary polynomial multiply d = a * b */
248static void
249s_bmul_d(const mp_digit *a, mp_size a_len, mp_digit b, mp_digit *d)
250{
251 mp_digit a_i, a0b0, a1b1, carry = 0;
252 while (a_len--) {
253 a_i = *a++;
254 s_bmul_1x1(&a1b1, &a0b0, a_i, b);
255 *d++ = a0b0 ^ carry;
256 carry = a1b1;
257 }
258 *d = carry;
259}
260
261/* Compute binary polynomial xor multiply accumulate d ^= a * b */
262static void
263s_bmul_d_add(const mp_digit *a, mp_size a_len, mp_digit b, mp_digit *d)
264{
265 mp_digit a_i, a0b0, a1b1, carry = 0;
266 while (a_len--) {
267 a_i = *a++;
268 s_bmul_1x1(&a1b1, &a0b0, a_i, b);
269 *d++ ^= a0b0 ^ carry;
270 carry = a1b1;
271 }
272 *d ^= carry;
273}
274
275/* Compute binary polynomial xor multiply c = a * b.
276 * All parameters may be identical.
277 */
278mp_err
279mp_bmul(const mp_int *a, const mp_int *b, mp_int *c)
280{
281 mp_digit *pb, b_i;
282 mp_int tmp;
283 mp_size ib, a_used, b_used;
284 mp_err res = MP_OKAY;
285
286 MP_DIGITS(&tmp) = 0;
287
288 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
289
290 if (a == c) {
291 MP_CHECKOK( mp_init_copy(&tmp, a) );
292 if (a == b)
293 b = &tmp;
294 a = &tmp;
295 } else if (b == c) {
296 MP_CHECKOK( mp_init_copy(&tmp, b) );
297 b = &tmp;
298 }
299
300 if (MP_USED(a) < MP_USED(b)) {
301 const mp_int *xch = b; /* switch a and b if b longer */
302 b = a;
303 a = xch;
304 }
305
306 MP_USED(c) = 1; MP_DIGIT(c, 0) = 0;
307 MP_CHECKOK( s_mp_pad(c, USED(a) + USED(b)) );
308
309 pb = MP_DIGITS(b);
310 s_bmul_d(MP_DIGITS(a), MP_USED(a), *pb++, MP_DIGITS(c));
311
312 /* Outer loop: Digits of b */
313 a_used = MP_USED(a);
314 b_used = MP_USED(b);
315 MP_USED(c) = a_used + b_used;
316 for (ib = 1; ib < b_used; ib++) {
317 b_i = *pb++;
318
319 /* Inner product: Digits of a */
320 if (b_i)
321 s_bmul_d_add(MP_DIGITS(a), a_used, b_i, MP_DIGITS(c) + ib);
322 else
323 MP_DIGIT(c, ib + a_used) = b_i;
324 }
325
326 s_mp_clamp(c);
327
328 SIGN(c) = ZPOS;
329
330CLEANUP:
331 mp_clear(&tmp);
332 return res;
333}
334
335
336/* Compute modular reduction of a and store result in r.
337 * r could be a.
338 * For modular arithmetic, the irreducible polynomial f(t) is represented
339 * as an array of int[], where f(t) is of the form:
340 * f(t) = t^p[0] + t^p[1] + ... + t^p[k]
341 * where m = p[0] > p[1] > ... > p[k] = 0.
342 */
343mp_err
344mp_bmod(const mp_int *a, const unsigned int p[], mp_int *r)
345{
346 int j, k;
347 int n, dN, d0, d1;
348 mp_digit zz, *z, tmp;
349 mp_size used;
350 mp_err res = MP_OKAY;
351
352 /* The algorithm does the reduction in place in r,
353 * if a != r, copy a into r first so reduction can be done in r
354 */
355 if (a != r) {
356 MP_CHECKOK( mp_copy(a, r) );
357 }
358 z = MP_DIGITS(r);
359
360 /* start reduction */
361 dN = p[0] / MP_DIGIT_BITS;
362 used = MP_USED(r);
363
364 for (j = used - 1; j > dN;) {
365
366 zz = z[j];
367 if (zz == 0) {
368 j--; continue;
369 }
370 z[j] = 0;
371
372 for (k = 1; p[k] > 0; k++) {
373 /* reducing component t^p[k] */
374 n = p[0] - p[k];
375 d0 = n % MP_DIGIT_BITS;
376 d1 = MP_DIGIT_BITS - d0;
377 n /= MP_DIGIT_BITS;
378 z[j-n] ^= (zz>>d0);
379 if (d0)
380 z[j-n-1] ^= (zz<<d1);
381 }
382
383 /* reducing component t^0 */
384 n = dN;
385 d0 = p[0] % MP_DIGIT_BITS;
386 d1 = MP_DIGIT_BITS - d0;
387 z[j-n] ^= (zz >> d0);
388 if (d0)
389 z[j-n-1] ^= (zz << d1);
390
391 }
392
393 /* final round of reduction */
394 while (j == dN) {
395
396 d0 = p[0] % MP_DIGIT_BITS;
397 zz = z[dN] >> d0;
398 if (zz == 0) break;
399 d1 = MP_DIGIT_BITS - d0;
400
401 /* clear up the top d1 bits */
402 if (d0) z[dN] = (z[dN] << d1) >> d1;
403 *z ^= zz; /* reduction t^0 component */
404
405 for (k = 1; p[k] > 0; k++) {
406 /* reducing component t^p[k]*/
407 n = p[k] / MP_DIGIT_BITS;
408 d0 = p[k] % MP_DIGIT_BITS;
409 d1 = MP_DIGIT_BITS - d0;
410 z[n] ^= (zz << d0);
411 tmp = zz >> d1;
412 if (d0 && tmp)
413 z[n+1] ^= tmp;
414 }
415 }
416
417 s_mp_clamp(r);
418CLEANUP:
419 return res;
420}
421
422/* Compute the product of two polynomials a and b, reduce modulo p,
423 * Store the result in r. r could be a or b; a could be b.
424 */
425mp_err
426mp_bmulmod(const mp_int *a, const mp_int *b, const unsigned int p[], mp_int *r)
427{
428 mp_err res;
429
430 if (a == b) return mp_bsqrmod(a, p, r);
431 if ((res = mp_bmul(a, b, r) ) != MP_OKAY)
432 return res;
433 return mp_bmod(r, p, r);
434}
435
436/* Compute binary polynomial squaring c = a*a mod p .
437 * Parameter r and a can be identical.
438 */
439
440mp_err
441mp_bsqrmod(const mp_int *a, const unsigned int p[], mp_int *r)
442{
443 mp_digit *pa, *pr, a_i;
444 mp_int tmp;
445 mp_size ia, a_used;
446 mp_err res;
447
448 ARGCHK(a != NULL && r != NULL, MP_BADARG);
449 MP_DIGITS(&tmp) = 0;
450
451 if (a == r) {
452 MP_CHECKOK( mp_init_copy(&tmp, a) );
453 a = &tmp;
454 }
455
456 MP_USED(r) = 1; MP_DIGIT(r, 0) = 0;
457 MP_CHECKOK( s_mp_pad(r, 2*USED(a)) );
458
459 pa = MP_DIGITS(a);
460 pr = MP_DIGITS(r);
461 a_used = MP_USED(a);
462 MP_USED(r) = 2 * a_used;
463
464 for (ia = 0; ia < a_used; ia++) {
465 a_i = *pa++;
466 *pr++ = gf2m_SQR0(a_i);
467 *pr++ = gf2m_SQR1(a_i);
468 }
469
470 MP_CHECKOK( mp_bmod(r, p, r) );
471 s_mp_clamp(r);
472 SIGN(r) = ZPOS;
473
474CLEANUP:
475 mp_clear(&tmp);
476 return res;
477}
478
479/* Compute binary polynomial y/x mod p, y divided by x, reduce modulo p.
480 * Store the result in r. r could be x or y, and x could equal y.
481 * Uses algorithm Modular_Division_GF(2^m) from
482 * Chang-Shantz, S. "From Euclid's GCD to Montgomery Multiplication to
483 * the Great Divide".
484 */
485int
486mp_bdivmod(const mp_int *y, const mp_int *x, const mp_int *pp,
487 const unsigned int p[], mp_int *r)
488{
489 mp_int aa, bb, uu;
490 mp_int *a, *b, *u, *v;
491 mp_err res = MP_OKAY;
492
493 MP_DIGITS(&aa) = 0;
494 MP_DIGITS(&bb) = 0;
495 MP_DIGITS(&uu) = 0;
496
497 MP_CHECKOK( mp_init_copy(&aa, x) );
498 MP_CHECKOK( mp_init_copy(&uu, y) );
499 MP_CHECKOK( mp_init_copy(&bb, pp) );
500 MP_CHECKOK( s_mp_pad(r, USED(pp)) );
501 MP_USED(r) = 1; MP_DIGIT(r, 0) = 0;
502
503 a = &aa; b= &bb; u=&uu; v=r;
504 /* reduce x and y mod p */
505 MP_CHECKOK( mp_bmod(a, p, a) );
506 MP_CHECKOK( mp_bmod(u, p, u) );
507
508 while (!mp_isodd(a)) {
509 s_mp_div2(a);
510 if (mp_isodd(u)) {
511 MP_CHECKOK( mp_badd(u, pp, u) );
512 }
513 s_mp_div2(u);
514 }
515
516 do {
517 if (mp_cmp_mag(b, a) > 0) {
518 MP_CHECKOK( mp_badd(b, a, b) );
519 MP_CHECKOK( mp_badd(v, u, v) );
520 do {
521 s_mp_div2(b);
522 if (mp_isodd(v)) {
523 MP_CHECKOK( mp_badd(v, pp, v) );
524 }
525 s_mp_div2(v);
526 } while (!mp_isodd(b));
527 }
528 else if ((MP_DIGIT(a,0) == 1) && (MP_USED(a) == 1))
529 break;
530 else {
531 MP_CHECKOK( mp_badd(a, b, a) );
532 MP_CHECKOK( mp_badd(u, v, u) );
533 do {
534 s_mp_div2(a);
535 if (mp_isodd(u)) {
536 MP_CHECKOK( mp_badd(u, pp, u) );
537 }
538 s_mp_div2(u);
539 } while (!mp_isodd(a));
540 }
541 } while (1);
542
543 MP_CHECKOK( mp_copy(u, r) );
544
545CLEANUP:
546 /* XXX this appears to be a memory leak in the NSS code */
547 mp_clear(&aa);
548 mp_clear(&bb);
549 mp_clear(&uu);
550 return res;
551
552}
553
554/* Convert the bit-string representation of a polynomial a into an array
555 * of integers corresponding to the bits with non-zero coefficient.
556 * Up to max elements of the array will be filled. Return value is total
557 * number of coefficients that would be extracted if array was large enough.
558 */
559int
560mp_bpoly2arr(const mp_int *a, unsigned int p[], int max)
561{
562 int i, j, k;
563 mp_digit top_bit, mask;
564
565 top_bit = 1;
566 top_bit <<= MP_DIGIT_BIT - 1;
567
568 for (k = 0; k < max; k++) p[k] = 0;
569 k = 0;
570
571 for (i = MP_USED(a) - 1; i >= 0; i--) {
572 mask = top_bit;
573 for (j = MP_DIGIT_BIT - 1; j >= 0; j--) {
574 if (MP_DIGITS(a)[i] & mask) {
575 if (k < max) p[k] = MP_DIGIT_BIT * i + j;
576 k++;
577 }
578 mask >>= 1;
579 }
580 }
581
582 return k;
583}
584
585/* Convert the coefficient array representation of a polynomial to a
586 * bit-string. The array must be terminated by 0.
587 */
588mp_err
589mp_barr2poly(const unsigned int p[], mp_int *a)
590{
591
592 mp_err res = MP_OKAY;
593 int i;
594
595 mp_zero(a);
596 for (i = 0; p[i] > 0; i++) {
597 MP_CHECKOK( mpl_set_bit(a, p[i], 1) );
598 }
599 MP_CHECKOK( mpl_set_bit(a, 0, 1) );
600
601CLEANUP:
602 return res;
603}
604