1 /* $OpenBSD: math_2n.c,v 1.23 2005/05/03 13:50:44 moritz Exp $ */
2 /* $EOM: math_2n.c,v 1.15 1999/04/20 09:23:30 niklas Exp $ */
3
4 /*
5 * Copyright (c) 1998 Niels Provos. All rights reserved.
6 * Copyright (c) 1999 Niklas Hallqvist. All rights reserved.
7 *
8 * Redistribution and use in source and binary forms, with or without
9 * modification, are permitted provided that the following conditions
10 * are met:
11 * 1. Redistributions of source code must retain the above copyright
12 * notice, this list of conditions and the following disclaimer.
13 * 2. Redistributions in binary form must reproduce the above copyright
14 * notice, this list of conditions and the following disclaimer in the
15 * documentation and/or other materials provided with the distribution.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
18 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
19 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
20 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
21 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
22 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
23 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
24 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
26 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 */
28
29 /*
30 * This code was written under funding by Ericsson Radio Systems.
31 */
32
33 /*
34 * B2N is a module for doing arithmetic on the Field GF(2**n) which is
35 * isomorph to ring of polynomials GF(2)[x]/p(x) where p(x) is an
36 * irreduciable polynomial over GF(2)[x] with grade n.
37 *
38 * First we need functions which operate on GF(2)[x], operation
39 * on GF(2)[x]/p(x) can be done as for Z_p then.
40 */
41
42 #include <stdlib.h>
43 #include <string.h>
44 #include <stdio.h>
45
46 #include "math_2n.h"
47 #include "util.h"
48
49 static u_int8_t hex2int(char);
50
51 CHUNK_TYPE b2n_mask[CHUNK_BITS] = {
52 0x01, 0x02, 0x04, 0x08, 0x10, 0x20, 0x40, 0x80,
53 #if CHUNK_BITS > 8
54 0x0100, 0x0200, 0x0400, 0x0800, 0x1000, 0x2000, 0x4000, 0x8000,
55 #if CHUNK_BITS > 16
56 0x00010000, 0x00020000, 0x00040000, 0x00080000,
57 0x00100000, 0x00200000, 0x00400000, 0x00800000,
58 0x01000000, 0x02000000, 0x04000000, 0x08000000,
59 0x10000000, 0x20000000, 0x40000000, 0x80000000,
60 #endif
61 #endif
62 };
63
64 /* Convert a hex character to its integer value. */
65 static u_int8_t
hex2int(char c)66 hex2int(char c)
67 {
68 if (c <= '9')
69 return c - '0';
70 if (c <= 'f')
71 return 10 + c - 'a';
72
73 return 0;
74 }
75
76 int
b2n_random(b2n_ptr n,u_int32_t bits)77 b2n_random(b2n_ptr n, u_int32_t bits)
78 {
79 if (b2n_resize(n, (CHUNK_MASK + bits) >> CHUNK_SHIFTS))
80 return -1;
81
82 getrandom((u_int8_t *) n->limp, CHUNK_BYTES * n->chunks);
83
84 /* Get the number of significant bits right */
85 if (bits & CHUNK_MASK) {
86 CHUNK_TYPE m =
87 (((1 << ((bits & CHUNK_MASK) - 1)) - 1) << 1) | 1;
88 n->limp[n->chunks - 1] &= m;
89 }
90 n->dirty = 1;
91 return 0;
92 }
93
94 /* b2n management functions */
95
96 void
b2n_init(b2n_ptr n)97 b2n_init(b2n_ptr n)
98 {
99 n->chunks = 0;
100 n->limp = 0;
101 }
102
103 void
b2n_clear(b2n_ptr n)104 b2n_clear(b2n_ptr n)
105 {
106 if (n->limp)
107 free(n->limp);
108 }
109
110 int
b2n_resize(b2n_ptr n,unsigned int chunks)111 b2n_resize(b2n_ptr n, unsigned int chunks)
112 {
113 size_t old = n->chunks;
114 size_t size;
115 CHUNK_TYPE *new;
116
117 if (chunks == 0)
118 chunks = 1;
119
120 if (chunks == old)
121 return 0;
122
123 size = CHUNK_BYTES * chunks;
124
125 new = realloc(n->limp, size);
126 if (!new)
127 return -1;
128
129 n->limp = new;
130 n->chunks = chunks;
131 n->bits = chunks << CHUNK_SHIFTS;
132 n->dirty = 1;
133
134 if (chunks > old)
135 bzero(n->limp + old, size - CHUNK_BYTES * old);
136
137 return 0;
138 }
139
140 /* Simple assignment functions. */
141
142 int
b2n_set(b2n_ptr d,b2n_ptr s)143 b2n_set(b2n_ptr d, b2n_ptr s)
144 {
145 if (d == s)
146 return 0;
147
148 b2n_sigbit(s);
149 if (b2n_resize(d, (CHUNK_MASK + s->bits) >> CHUNK_SHIFTS))
150 return -1;
151 memcpy(d->limp, s->limp, CHUNK_BYTES * d->chunks);
152 d->bits = s->bits;
153 d->dirty = s->dirty;
154 return 0;
155 }
156
157 int
b2n_set_null(b2n_ptr n)158 b2n_set_null(b2n_ptr n)
159 {
160 if (b2n_resize(n, 1))
161 return -1;
162 n->limp[0] = n->bits = n->dirty = 0;
163 return 0;
164 }
165
166 int
b2n_set_ui(b2n_ptr n,unsigned int val)167 b2n_set_ui(b2n_ptr n, unsigned int val)
168 {
169 #if CHUNK_BITS < 32
170 int i, chunks;
171
172 chunks = (CHUNK_BYTES - 1 + sizeof(val)) / CHUNK_BYTES;
173
174 if (b2n_resize(n, chunks))
175 return -1;
176
177 for (i = 0; i < chunks; i++) {
178 n->limp[i] = val & CHUNK_BMASK;
179 val >>= CHUNK_BITS;
180 }
181 #else
182 if (b2n_resize(n, 1))
183 return -1;
184 n->limp[0] = val;
185 #endif
186 n->dirty = 1;
187 return 0;
188 }
189
190 /* XXX This one only takes hex at the moment. */
191 int
b2n_set_str(b2n_ptr n,char * str)192 b2n_set_str(b2n_ptr n, char *str)
193 {
194 int i, j, w, len, chunks;
195 CHUNK_TYPE tmp;
196
197 if (strncasecmp(str, "0x", 2))
198 return -1;
199
200 /* Make the hex string even lengthed */
201 len = strlen(str) - 2;
202 if (len & 1) {
203 len++;
204 str++;
205 } else
206 str += 2;
207
208 len /= 2;
209
210 chunks = (CHUNK_BYTES - 1 + len) / CHUNK_BYTES;
211 if (b2n_resize(n, chunks))
212 return -1;
213 bzero(n->limp, CHUNK_BYTES * n->chunks);
214
215 for (w = 0, i = 0; i < chunks; i++) {
216 tmp = 0;
217 for (j = (i == 0 ?
218 ((len - 1) % CHUNK_BYTES) + 1 : CHUNK_BYTES);
219 j > 0; j--) {
220 tmp <<= 8;
221 tmp |= (hex2int(str[w]) << 4) | hex2int(str[w + 1]);
222 w += 2;
223 }
224 n->limp[chunks - 1 - i] = tmp;
225 }
226
227 n->dirty = 1;
228 return 0;
229 }
230
231 /* Arithmetic functions. */
232
233 u_int32_t
b2n_sigbit(b2n_ptr n)234 b2n_sigbit(b2n_ptr n)
235 {
236 int i, j;
237
238 if (!n->dirty)
239 return n->bits;
240
241 for (i = n->chunks - 1; i > 0; i--)
242 if (n->limp[i])
243 break;
244
245 if (!n->limp[i])
246 return 0;
247
248 for (j = CHUNK_MASK; j > 0; j--)
249 if (n->limp[i] & b2n_mask[j])
250 break;
251
252 n->bits = (i << CHUNK_SHIFTS) + j + 1;
253 n->dirty = 0;
254 return n->bits;
255 }
256
257 /* Addition on GF(2)[x] is nice, its just an XOR. */
258 int
b2n_add(b2n_ptr d,b2n_ptr a,b2n_ptr b)259 b2n_add(b2n_ptr d, b2n_ptr a, b2n_ptr b)
260 {
261 int i;
262 b2n_ptr bmin, bmax;
263
264 if (!b2n_cmp_null(a))
265 return b2n_set(d, b);
266
267 if (!b2n_cmp_null(b))
268 return b2n_set(d, a);
269
270 bmin = B2N_MIN(a, b);
271 bmax = B2N_MAX(a, b);
272
273 if (b2n_resize(d, bmax->chunks))
274 return -1;
275
276 for (i = 0; i < bmin->chunks; i++)
277 d->limp[i] = bmax->limp[i] ^ bmin->limp[i];
278
279 /*
280 * If d is not bmax, we have to copy the rest of the bytes, and also
281 * need to adjust to number of relevant bits.
282 */
283 if (d != bmax) {
284 for (; i < bmax->chunks; i++)
285 d->limp[i] = bmax->limp[i];
286
287 d->bits = bmax->bits;
288 }
289 /*
290 * Help to converse memory. When the result of the addition is zero
291 * truncate the used amount of memory.
292 */
293 if (d != bmax && !b2n_cmp_null(d))
294 return b2n_set_null(d);
295 else
296 d->dirty = 1;
297 return 0;
298 }
299
300 /* Compare two polynomials. */
301 int
b2n_cmp(b2n_ptr n,b2n_ptr m)302 b2n_cmp(b2n_ptr n, b2n_ptr m)
303 {
304 int sn, sm;
305 int i;
306
307 sn = b2n_sigbit(n);
308 sm = b2n_sigbit(m);
309
310 if (sn > sm)
311 return 1;
312 if (sn < sm)
313 return -1;
314
315 for (i = n->chunks - 1; i >= 0; i--)
316 if (n->limp[i] > m->limp[i])
317 return 1;
318 else if (n->limp[i] < m->limp[i])
319 return -1;
320
321 return 0;
322 }
323
324 int
b2n_cmp_null(b2n_ptr a)325 b2n_cmp_null(b2n_ptr a)
326 {
327 int i = 0;
328
329 do {
330 if (a->limp[i])
331 return 1;
332 } while (++i < a->chunks);
333
334 return 0;
335 }
336
337 /* Left shift, needed for polynomial multiplication. */
338 int
b2n_lshift(b2n_ptr d,b2n_ptr n,unsigned int s)339 b2n_lshift(b2n_ptr d, b2n_ptr n, unsigned int s)
340 {
341 int i, maj, min, chunks;
342 u_int16_t bits = b2n_sigbit(n), add;
343 CHUNK_TYPE *p, *op;
344
345 if (!s)
346 return b2n_set(d, n);
347
348 maj = s >> CHUNK_SHIFTS;
349 min = s & CHUNK_MASK;
350
351 add = (!(bits & CHUNK_MASK) ||
352 ((bits & CHUNK_MASK) + min) > CHUNK_MASK) ? 1 : 0;
353 chunks = n->chunks;
354 if (b2n_resize(d, chunks + maj + add))
355 return -1;
356 memmove(d->limp + maj, n->limp, CHUNK_BYTES * chunks);
357
358 if (maj)
359 bzero(d->limp, CHUNK_BYTES * maj);
360 if (add)
361 d->limp[d->chunks - 1] = 0;
362
363 /* If !min there are no bit shifts, we are done */
364 if (!min)
365 return 0;
366
367 op = p = &d->limp[d->chunks - 1];
368 for (i = d->chunks - 2; i >= maj; i--) {
369 op--;
370 *p = (*p << min) | (*op >> (CHUNK_BITS - min));
371 p--;
372 }
373 *p <<= min;
374
375 d->dirty = 0;
376 d->bits = bits + (maj << CHUNK_SHIFTS) + min;
377 return 0;
378 }
379
380 /* Right shift, needed for polynomial division. */
381 int
b2n_rshift(b2n_ptr d,b2n_ptr n,unsigned int s)382 b2n_rshift(b2n_ptr d, b2n_ptr n, unsigned int s)
383 {
384 int maj, min, size = n->chunks, newsize;
385 b2n_ptr tmp;
386
387 if (!s)
388 return b2n_set(d, n);
389
390 maj = s >> CHUNK_SHIFTS;
391
392 newsize = size - maj;
393
394 if (size < maj)
395 return b2n_set_null(d);
396
397 min = (CHUNK_BITS - (s & CHUNK_MASK)) & CHUNK_MASK;
398 if (min) {
399 if ((b2n_sigbit(n) & CHUNK_MASK) > (u_int32_t) min)
400 newsize++;
401
402 if (b2n_lshift(d, n, min))
403 return -1;
404 tmp = d;
405 } else
406 tmp = n;
407
408 memmove(d->limp, tmp->limp + maj + (min ? 1 : 0),
409 CHUNK_BYTES * newsize);
410 if (b2n_resize(d, newsize))
411 return -1;
412
413 d->bits = tmp->bits - ((maj + (min ? 1 : 0)) << CHUNK_SHIFTS);
414 return 0;
415 }
416
417 /* Normal polynomial multiplication. */
418 int
b2n_mul(b2n_ptr d,b2n_ptr n,b2n_ptr m)419 b2n_mul(b2n_ptr d, b2n_ptr n, b2n_ptr m)
420 {
421 int i, j;
422 b2n_t tmp, tmp2;
423
424 if (!b2n_cmp_null(m) || !b2n_cmp_null(n))
425 return b2n_set_null(d);
426
427 if (b2n_sigbit(m) == 1)
428 return b2n_set(d, n);
429
430 if (b2n_sigbit(n) == 1)
431 return b2n_set(d, m);
432
433 b2n_init(tmp);
434 b2n_init(tmp2);
435
436 if (b2n_set(tmp, B2N_MAX(n, m)))
437 goto fail;
438 if (b2n_set(tmp2, B2N_MIN(n, m)))
439 goto fail;
440
441 if (b2n_set_null(d))
442 goto fail;
443
444 for (i = 0; i < tmp2->chunks; i++)
445 if (tmp2->limp[i])
446 for (j = 0; j < CHUNK_BITS; j++) {
447 if (tmp2->limp[i] & b2n_mask[j])
448 if (b2n_add(d, d, tmp))
449 goto fail;
450
451 if (b2n_lshift(tmp, tmp, 1))
452 goto fail;
453 }
454 else if (b2n_lshift(tmp, tmp, CHUNK_BITS))
455 goto fail;
456
457 b2n_clear(tmp);
458 b2n_clear(tmp2);
459 return 0;
460
461 fail:
462 b2n_clear(tmp);
463 b2n_clear(tmp2);
464 return -1;
465 }
466
467 /*
468 * Squaring in this polynomial ring is more efficient than normal
469 * multiplication.
470 */
471 int
b2n_square(b2n_ptr d,b2n_ptr n)472 b2n_square(b2n_ptr d, b2n_ptr n)
473 {
474 int i, j, maj, min, bits, chunk;
475 b2n_t t;
476
477 maj = b2n_sigbit(n);
478 min = maj & CHUNK_MASK;
479 maj = (maj + CHUNK_MASK) >> CHUNK_SHIFTS;
480
481 b2n_init(t);
482 if (b2n_resize(t,
483 2 * maj + ((CHUNK_MASK + 2 * min) >> CHUNK_SHIFTS))) {
484 b2n_clear(t);
485 return -1;
486 }
487 chunk = 0;
488 bits = 0;
489
490 for (i = 0; i < maj; i++)
491 if (n->limp[i])
492 for (j = 0; j < CHUNK_BITS; j++) {
493 if (n->limp[i] & b2n_mask[j])
494 t->limp[chunk] ^= b2n_mask[bits];
495
496 bits += 2;
497 if (bits >= CHUNK_BITS) {
498 chunk++;
499 bits &= CHUNK_MASK;
500 }
501 }
502 else
503 chunk += 2;
504
505 t->dirty = 1;
506 B2N_SWAP(d, t);
507 b2n_clear(t);
508 return 0;
509 }
510
511 /*
512 * Normal polynomial division.
513 * These functions are far from optimal in speed.
514 */
515 int
b2n_div_q(b2n_ptr d,b2n_ptr n,b2n_ptr m)516 b2n_div_q(b2n_ptr d, b2n_ptr n, b2n_ptr m)
517 {
518 b2n_t r;
519 int rv;
520
521 b2n_init(r);
522 rv = b2n_div(d, r, n, m);
523 b2n_clear(r);
524 return rv;
525 }
526
527 int
b2n_div_r(b2n_ptr r,b2n_ptr n,b2n_ptr m)528 b2n_div_r(b2n_ptr r, b2n_ptr n, b2n_ptr m)
529 {
530 b2n_t q;
531 int rv;
532
533 b2n_init(q);
534 rv = b2n_div(q, r, n, m);
535 b2n_clear(q);
536 return rv;
537 }
538
539 int
b2n_div(b2n_ptr q,b2n_ptr r,b2n_ptr n,b2n_ptr m)540 b2n_div(b2n_ptr q, b2n_ptr r, b2n_ptr n, b2n_ptr m)
541 {
542 int i, j, len, bits;
543 u_int32_t sm, sn;
544 b2n_t nenn, div, shift, mask;
545
546 /* If Teiler > Zaehler, the result is 0 */
547 if ((sm = b2n_sigbit(m)) > (sn = b2n_sigbit(n))) {
548 if (b2n_set_null(q))
549 return -1;
550 return b2n_set(r, n);
551 }
552 if (sm == 0)
553 /* Division by Zero */
554 return -1;
555 else if (sm == 1) {
556 /* Division by the One-Element */
557 if (b2n_set(q, n))
558 return -1;
559 return b2n_set_null(r);
560 }
561 b2n_init(nenn);
562 b2n_init(div);
563 b2n_init(shift);
564 b2n_init(mask);
565
566 if (b2n_set(nenn, n))
567 goto fail;
568 if (b2n_set(div, m))
569 goto fail;
570 if (b2n_set(shift, m))
571 goto fail;
572 if (b2n_set_ui(mask, 1))
573 goto fail;
574
575 if (b2n_resize(q, (sn - sm + CHUNK_MASK) >> CHUNK_SHIFTS))
576 goto fail;
577 bzero(q->limp, CHUNK_BYTES * q->chunks);
578
579 if (b2n_lshift(shift, shift, sn - sm))
580 goto fail;
581 if (b2n_lshift(mask, mask, sn - sm))
582 goto fail;
583
584 /* Number of significant octets */
585 len = (sn - 1) >> CHUNK_SHIFTS;
586 /* The first iteration is done over the relevant bits */
587 bits = (CHUNK_MASK + sn) & CHUNK_MASK;
588 for (i = len; i >= 0 && b2n_sigbit(nenn) >= sm; i--)
589 for (j = (i == len ? bits : CHUNK_MASK); j >= 0 &&
590 b2n_sigbit(nenn) >= sm; j--) {
591 if (nenn->limp[i] & b2n_mask[j]) {
592 if (b2n_sub(nenn, nenn, shift))
593 goto fail;
594 if (b2n_add(q, q, mask))
595 goto fail;
596 }
597 if (b2n_rshift(shift, shift, 1))
598 goto fail;
599 if (b2n_rshift(mask, mask, 1))
600 goto fail;
601 }
602
603 B2N_SWAP(r, nenn);
604
605 b2n_clear(nenn);
606 b2n_clear(div);
607 b2n_clear(shift);
608 b2n_clear(mask);
609 return 0;
610
611 fail:
612 b2n_clear(nenn);
613 b2n_clear(div);
614 b2n_clear(shift);
615 b2n_clear(mask);
616 return -1;
617 }
618
619 /* Functions for Operation on GF(2**n) ~= GF(2)[x]/p(x). */
620 int
b2n_mod(b2n_ptr m,b2n_ptr n,b2n_ptr p)621 b2n_mod(b2n_ptr m, b2n_ptr n, b2n_ptr p)
622 {
623 int bits, size;
624
625 if (b2n_div_r(m, n, p))
626 return -1;
627
628 bits = b2n_sigbit(m);
629 size = ((CHUNK_MASK + bits) >> CHUNK_SHIFTS);
630 if (size == 0)
631 size = 1;
632 if (m->chunks > size)
633 if (b2n_resize(m, size))
634 return -1;
635
636 m->bits = bits;
637 m->dirty = 0;
638 return 0;
639 }
640
641 int
b2n_gcd(b2n_ptr e,b2n_ptr go,b2n_ptr ho)642 b2n_gcd(b2n_ptr e, b2n_ptr go, b2n_ptr ho)
643 {
644 b2n_t g, h;
645
646 b2n_init(g);
647 b2n_init(h);
648 if (b2n_set(g, go))
649 goto fail;
650 if (b2n_set(h, ho))
651 goto fail;
652
653 while (b2n_cmp_null(h)) {
654 if (b2n_mod(g, g, h))
655 goto fail;
656 B2N_SWAP(g, h);
657 }
658
659 B2N_SWAP(e, g);
660
661 b2n_clear(g);
662 b2n_clear(h);
663 return 0;
664
665 fail:
666 b2n_clear(g);
667 b2n_clear(h);
668 return -1;
669 }
670
671 int
b2n_mul_inv(b2n_ptr ga,b2n_ptr be,b2n_ptr p)672 b2n_mul_inv(b2n_ptr ga, b2n_ptr be, b2n_ptr p)
673 {
674 b2n_t a;
675
676 b2n_init(a);
677 if (b2n_set_ui(a, 1))
678 goto fail;
679
680 if (b2n_div_mod(ga, a, be, p))
681 goto fail;
682
683 b2n_clear(a);
684 return 0;
685
686 fail:
687 b2n_clear(a);
688 return -1;
689 }
690
691 int
b2n_div_mod(b2n_ptr ga,b2n_ptr a,b2n_ptr be,b2n_ptr p)692 b2n_div_mod(b2n_ptr ga, b2n_ptr a, b2n_ptr be, b2n_ptr p)
693 {
694 b2n_t s0, s1, s2, q, r0, r1;
695
696 /* There is no multiplicative inverse to Null. */
697 if (!b2n_cmp_null(be))
698 return b2n_set_null(ga);
699
700 b2n_init(s0);
701 b2n_init(s1);
702 b2n_init(s2);
703 b2n_init(r0);
704 b2n_init(r1);
705 b2n_init(q);
706
707 if (b2n_set(r0, p))
708 goto fail;
709 if (b2n_set(r1, be))
710 goto fail;
711
712 if (b2n_set_null(s0))
713 goto fail;
714 if (b2n_set(s1, a))
715 goto fail;
716
717 while (b2n_cmp_null(r1)) {
718 if (b2n_div(q, r0, r0, r1))
719 goto fail;
720 B2N_SWAP(r0, r1);
721
722 if (b2n_mul(s2, q, s1))
723 goto fail;
724 if (b2n_mod(s2, s2, p))
725 goto fail;
726 if (b2n_sub(s2, s0, s2))
727 goto fail;
728
729 B2N_SWAP(s0, s1);
730 B2N_SWAP(s1, s2);
731 }
732 B2N_SWAP(ga, s0);
733
734 b2n_clear(s0);
735 b2n_clear(s1);
736 b2n_clear(s2);
737 b2n_clear(r0);
738 b2n_clear(r1);
739 b2n_clear(q);
740 return 0;
741
742 fail:
743 b2n_clear(s0);
744 b2n_clear(s1);
745 b2n_clear(s2);
746 b2n_clear(r0);
747 b2n_clear(r1);
748 b2n_clear(q);
749 return -1;
750 }
751
752 /*
753 * The trace tells us if there do exist any square roots
754 * for 'a' in GF(2)[x]/p(x). The number of square roots is
755 * 2 - 2*Trace.
756 * If z is a square root, z + 1 is the other.
757 */
758 int
b2n_trace(b2n_ptr ho,b2n_ptr a,b2n_ptr p)759 b2n_trace(b2n_ptr ho, b2n_ptr a, b2n_ptr p)
760 {
761 int i, m = b2n_sigbit(p) - 1;
762 b2n_t h;
763
764 b2n_init(h);
765 if (b2n_set(h, a))
766 goto fail;
767
768 for (i = 0; i < m - 1; i++) {
769 if (b2n_square(h, h))
770 goto fail;
771 if (b2n_mod(h, h, p))
772 goto fail;
773
774 if (b2n_add(h, h, a))
775 goto fail;
776 }
777 B2N_SWAP(ho, h);
778
779 b2n_clear(h);
780 return 0;
781
782 fail:
783 b2n_clear(h);
784 return -1;
785 }
786
787 /*
788 * The halftrace yields the square root if the degree of the
789 * irreduceable polynomial is odd.
790 */
791 int
b2n_halftrace(b2n_ptr ho,b2n_ptr a,b2n_ptr p)792 b2n_halftrace(b2n_ptr ho, b2n_ptr a, b2n_ptr p)
793 {
794 int i, m = b2n_sigbit(p) - 1;
795 b2n_t h;
796
797 b2n_init(h);
798 if (b2n_set(h, a))
799 goto fail;
800
801 for (i = 0; i < (m - 1) / 2; i++) {
802 if (b2n_square(h, h))
803 goto fail;
804 if (b2n_mod(h, h, p))
805 goto fail;
806 if (b2n_square(h, h))
807 goto fail;
808 if (b2n_mod(h, h, p))
809 goto fail;
810
811 if (b2n_add(h, h, a))
812 goto fail;
813 }
814
815 B2N_SWAP(ho, h);
816
817 b2n_clear(h);
818 return 0;
819
820 fail:
821 b2n_clear(h);
822 return -1;
823 }
824
825 /*
826 * Solving the equation: y**2 + y = b in GF(2**m) where ip is the
827 * irreduceable polynomial. If m is odd, use the half trace.
828 */
829 int
b2n_sqrt(b2n_ptr zo,b2n_ptr b,b2n_ptr ip)830 b2n_sqrt(b2n_ptr zo, b2n_ptr b, b2n_ptr ip)
831 {
832 int i, m = b2n_sigbit(ip) - 1;
833 b2n_t w, p, temp, z;
834
835 if (!b2n_cmp_null(b))
836 return b2n_set_null(z);
837
838 if (m & 1)
839 return b2n_halftrace(zo, b, ip);
840
841 b2n_init(z);
842 b2n_init(w);
843 b2n_init(p);
844 b2n_init(temp);
845
846 do {
847 if (b2n_random(p, m))
848 goto fail;
849 if (b2n_set_null(z))
850 goto fail;
851 if (b2n_set(w, p))
852 goto fail;
853
854 for (i = 1; i < m; i++) {
855 if (b2n_square(z, z)) /* z**2 */
856 goto fail;
857 if (b2n_mod(z, z, ip))
858 goto fail;
859
860 if (b2n_square(w, w)) /* w**2 */
861 goto fail;
862 if (b2n_mod(w, w, ip))
863 goto fail;
864
865 if (b2n_mul(temp, w, b)) /* w**2 * b */
866 goto fail;
867 if (b2n_mod(temp, temp, ip))
868 goto fail;
869 if (b2n_add(z, z, temp)) /* z**2 + w**2 + b */
870 goto fail;
871
872 if (b2n_add(w, w, p)) /* w**2 + p */
873 goto fail;
874 }
875 } while (!b2n_cmp_null(w));
876
877 B2N_SWAP(zo, z);
878
879 b2n_clear(w);
880 b2n_clear(p);
881 b2n_clear(temp);
882 b2n_clear(z);
883 return 0;
884
885 fail:
886 b2n_clear(w);
887 b2n_clear(p);
888 b2n_clear(temp);
889 b2n_clear(z);
890 return -1;
891 }
892
893 /* Exponentiation modulo a polynomial. */
894 int
b2n_exp_mod(b2n_ptr d,b2n_ptr b0,u_int32_t e,b2n_ptr p)895 b2n_exp_mod(b2n_ptr d, b2n_ptr b0, u_int32_t e, b2n_ptr p)
896 {
897 b2n_t u, b;
898
899 b2n_init(u);
900 b2n_init(b);
901 if (b2n_set_ui(u, 1))
902 goto fail;
903 if (b2n_mod(b, b0, p))
904 goto fail;
905
906 while (e) {
907 if (e & 1) {
908 if (b2n_mul(u, u, b))
909 goto fail;
910 if (b2n_mod(u, u, p))
911 goto fail;
912 }
913 if (b2n_square(b, b))
914 goto fail;
915 if (b2n_mod(b, b, p))
916 goto fail;
917 e >>= 1;
918 }
919
920 B2N_SWAP(d, u);
921
922 b2n_clear(u);
923 b2n_clear(b);
924 return 0;
925
926 fail:
927 b2n_clear(u);
928 b2n_clear(b);
929 return -1;
930 }
931
932 /*
933 * Low-level function to speed up scalar multiplication with
934 * elliptic curves.
935 * Multiplies a normal number by 3.
936 */
937
938 /* Normal addition behaves as Z_{2**n} and not F_{2**n}. */
939 int
b2n_nadd(b2n_ptr d0,b2n_ptr a0,b2n_ptr b0)940 b2n_nadd(b2n_ptr d0, b2n_ptr a0, b2n_ptr b0)
941 {
942 int i, carry;
943 b2n_ptr a, b;
944 b2n_t d;
945
946 if (!b2n_cmp_null(a0))
947 return b2n_set(d0, b0);
948
949 if (!b2n_cmp_null(b0))
950 return b2n_set(d0, a0);
951
952 b2n_init(d);
953 a = B2N_MAX(a0, b0);
954 b = B2N_MIN(a0, b0);
955
956 if (b2n_resize(d, a->chunks + 1)) {
957 b2n_clear(d);
958 return -1;
959 }
960 for (carry = i = 0; i < b->chunks; i++) {
961 d->limp[i] = a->limp[i] + b->limp[i] + carry;
962 carry = (d->limp[i] < a->limp[i] ? 1 : 0);
963 }
964
965 for (; i < a->chunks && carry; i++) {
966 d->limp[i] = a->limp[i] + carry;
967 carry = (d->limp[i] < a->limp[i] ? 1 : 0);
968 }
969
970 if (i < a->chunks)
971 memcpy(d->limp + i, a->limp + i,
972 CHUNK_BYTES * (a->chunks - i));
973
974 d->dirty = 1;
975 B2N_SWAP(d0, d);
976
977 b2n_clear(d);
978 return 0;
979 }
980
981 /* Very special sub, a > b. */
982 int
b2n_nsub(b2n_ptr d0,b2n_ptr a,b2n_ptr b)983 b2n_nsub(b2n_ptr d0, b2n_ptr a, b2n_ptr b)
984 {
985 int i, carry;
986 b2n_t d;
987
988 if (b2n_cmp(a, b) <= 0)
989 return b2n_set_null(d0);
990
991 b2n_init(d);
992 if (b2n_resize(d, a->chunks)) {
993 b2n_clear(d);
994 return -1;
995 }
996 for (carry = i = 0; i < b->chunks; i++) {
997 d->limp[i] = a->limp[i] - b->limp[i] - carry;
998 carry = (d->limp[i] > a->limp[i] ? 1 : 0);
999 }
1000
1001 for (; i < a->chunks && carry; i++) {
1002 d->limp[i] = a->limp[i] - carry;
1003 carry = (d->limp[i] > a->limp[i] ? 1 : 0);
1004 }
1005
1006 if (i < a->chunks)
1007 memcpy(d->limp + i, a->limp + i,
1008 CHUNK_BYTES * (a->chunks - i));
1009
1010 d->dirty = 1;
1011
1012 B2N_SWAP(d0, d);
1013
1014 b2n_clear(d);
1015 return 0;
1016 }
1017
1018 int
b2n_3mul(b2n_ptr d0,b2n_ptr e)1019 b2n_3mul(b2n_ptr d0, b2n_ptr e)
1020 {
1021 b2n_t d;
1022
1023 b2n_init(d);
1024 if (b2n_lshift(d, e, 1))
1025 goto fail;
1026
1027 if (b2n_nadd(d0, d, e))
1028 goto fail;
1029
1030 b2n_clear(d);
1031 return 0;
1032
1033 fail:
1034 b2n_clear(d);
1035 return -1;
1036 }
1037