projects / u / mdw / catacomb / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
Initial import.
[u/mdw/catacomb] / mp.c
1 /* -*-c-*-
2 *
3 * $Id: mp.c,v 1.1 1999/09/03 08:41:12 mdw Exp $
4 *
5 * Multiprecision arithmetic
6 *
7 * (c) 1998 Mark Wooding
8 */
9
10 /*----- Licensing notice --------------------------------------------------*
11 *
12 * This file is part of Catacomb.
13 *
14 * Catacomb is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU Library General Public License as
16 * published by the Free Software Foundation; either version 2 of the
17 * License, or (at your option) any later version.
18 *
19 * Catacomb is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU Library General Public License for more details.
23 *
24 * You should have received a copy of the GNU Library General Public
25 * License along with Catacomb; if not, write to the Free
26 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
27 * MA 02111-1307, USA.
28 */
29
30 /*----- Revision history --------------------------------------------------*
31 *
32 * $Log: mp.c,v $
33 * Revision 1.1 1999/09/03 08:41:12 mdw
34 * Initial import.
35 *
36 */
37
38 /*----- Header files ------------------------------------------------------*/
39
40 #include <limits.h>
41 #include <stddef.h>
42 #include <stdlib.h>
43 #include <string.h>
44
45 #include "alloc.c"
46 #include "bits.h"
47 #include "exc.h"
48
49 /*----- Important types ---------------------------------------------------*/
50
51 /* For now... */
52
53 #ifndef TEST_RIG
54
55 typedef unsigned short mpw;
56 #define MPW_BITS 16
57 #define MPW_MAX 0xffffu
58 typedef unsigned int mpd;
59 #define MPD_BITS 32
60 #define MPD_MAX 0xffffffffu
61
62 #else
63
64 typedef unsigned long mpw;
65 static unsigned mpw_bits;
66 static mpw mpw_max;
67 typedef unsigned long mpd;
68
69 #define MPW_BITS mpw_bits
70 #define MPW_MAX mpw_max
71
72 #endif
73
74 /*----- Data structures ---------------------------------------------------*/
75
76 /*----- Some basic values which are important -----------------------------*/
77
78 static mpw mp__v[] = { 1, 2, 3, 4, 5, 10 };
79 const mp mp_std[] = {
80 { 0, 0, 0, 0 }, /* 0 */
81 { mp__v + 0, 0, 0, 1 }, /* 1 */
82 { mp__v + 1, 0, 0, 1 }, /* 2 */
83 { mp__v + 2, 0, 0, 1 }, /* 3 */
84 { mp__v + 3, 0, 0, 1 }, /* 4 */
85 { mp__v + 4, 0, 0, 1 }, /* 5 */
86 { mp__v + 5, 0, 0, 1 }, /* 10 */
87 { mp__v + 0, MPF_SIGN, 0, 1 } /* -1 */
88 };
89
90 /*----- Memory management -------------------------------------------------*/
91
92 /* --- @mp_create@ --- *
93 *
94 * Arguments @mp *x@ = pointer to MP head
95 *
96 * Returns: ---
97 *
98 * Use: Initializes an MP ready for use. The initial value is zero.
99 */
100
101 void mp_create(mp *x)
102 {
103 x->v = 0;
104 x->f = 0;
105 x->sz = 0;
106 x->len = 0;
107 }
108
109 /* --- @mp_destroy@ --- *
110 *
111 * Arguments: @mp *x@ = pointer to MP head
112 *
113 * Returns: ---
114 *
115 * Use: Releases the memory used by an MP.
116 */
117
118 void mp_destroy(mp *x) { MP_DESTROY(x); }
119
120 /* --- @mp_resize@ --- *
121 *
122 * Arguments: @mp *x@ = pointer to MP head
123 * @size_t sz@ = size required (in words)
124 *
125 * Returns: ---
126 *
127 * Use: Resizes the MP so that its word vector has space for
128 * exactly @sz@ words.
129 */
130
131 void mp_resize(mp *x, size_t sz)
132 {
133 if (sz == 0)
134 mp_destroy(x);
135 else {
136 size_t min = sz > x->sz ? x->sz : sz;
137 mpw *v = xmalloc(sz * sizeof(mpw));
138 if (min)
139 memcpy(v, x->v, min * sizeof(mpw));
140 MP_BURN(x);
141 if (x->v)
142 free(x->v);
143 if (sz > min)
144 memset(v + min, 0, (sz - min) * sizeof(mpw));
145 x->v = v;
146 x->sz = sz;
147 }
148 }
149
150 /* --- @mp_norm@ --- *
151 *
152 * Arguments: @mp *x@ = pointer to MP head
153 *
154 * Returns: ---
155 *
156 * Use: `Normalizes' an MP. Fixes the @len@ field so that it's
157 * correct. Assumes that @len@ is either already correct or
158 * too large.
159 */
160
161 void mp_norm(mp *x) { MP_NORM(x); }
162
163 /* --- @mp_dump@ --- *
164 *
165 * Arguments: @mp *x@ = pointer to MP head
166 * @FILE *fp@ = pointer to stream to write on
167 *
168 * Returns: ---
169 *
170 * Use: Dumps an MP to a stream.
171 */
172
173 void mp_dump(mp *x, FILE *fp)
174 {
175 if (x->v) {
176 mpw *v = x->v, *vl = v + x->len;
177 while (v < vl) {
178 fprintf(fp, "%lx", (unsigned long)*v++);
179 if (v < vl) fputc('-', fp);
180 }
181 } else
182 fputs("<0>", fp);
183 }
184
185 /* --- @mp_copy@ --- *
186 *
187 * Arguments: @mp *d@ = pointer to MP head for destination
188 * @const mp *s@ = pointer to MP head for source
189 *
190 * Returns: ---
191 *
192 * Use: Copies an MP.
193 */
194
195 void mp_copy(mp *d, const mp *s)
196 {
197 if (d == s) {
198 MP_NORM(d);
199 return;
200 }
201 MP_BURN(d);
202 MP_ENSURE(d, s->len);
203 memcpy(d->v, s->v, s->len * sizeof(mpw));
204 if (d->f & MPF_BURN && d->sz > s->len)
205 memset(d->v + s->len, 0, (d->sz - s->len) * sizeof(mpw));
206 d->len = s->len;
207 d->f = s->f;
208 MP_NORM(d);
209 }
210
211 /* --- @mp_bits@ --- *
212 *
213 * Arguments: @mp *x@ = pointer to MP head
214 *
215 * Returns: Length of the number in bits.
216 *
217 * Use: Calculates the number of bits required to represent a number.
218 * The number must be normalized.
219 */
220
221 unsigned long mp_bits(mp *x)
222 {
223 if (!x->v)
224 return (0);
225 else {
226 unsigned long bits = MPW_BITS * (x->len - 1);
227 mpw w = x->v[x->len - 1];
228 while (w) {
229 bits++;
230 w >>= 1;
231 }
232 return (bits);
233 }
234 }
235
236 /* --- @mp_octets@ --- *
237 *
238 * Arguments: @mp *x@ = pointer to MP head
239 *
240 * Returns: Length of number in octets.
241 *
242 * Use: Calculates the number of octets required to represent a
243 * number. The number must be normalized.
244 */
245
246 size_t mp_octets(mp *x)
247 {
248 return ((mp_bits(x) + 7) & 7);
249 }
250
251 /*----- Loading and storing as binary data --------------------------------*/
252
253 /* --- @mp_storel@ --- *
254 *
255 * Arguments: @mp *x@ = pointer to MP head
256 * @octet *p@ = pointer to octet array
257 * @size_t sz@ = size of octet array
258 *
259 * Returns: ---
260 *
261 * Use: Stores an MP in an octet array, least significant octet
262 * first. High-end octets are silently discarded if there
263 * isn't enough space for them.
264 */
265
266 void mp_storel(mp *x, octet *p, size_t sz)
267 {
268 mpw *v = x->v, *vl = x->v + x->len;
269 mpw n, w = 0;
270 octet *q = p + sz;
271 unsigned bits = 0;
272
273 while (p < q) {
274 if (bits < 8) {
275 n = (v >= vl) ? 0 : *v++;
276 *p++ = (w | n << bits) & MASK8;
277 w = n >> (8 - bits);
278 bits += MPW_BITS - 8;
279 } else {
280 *p++ = w & MASK8;
281 w >>= 8;
282 bits -= 8;
283 }
284 }
285 }
286
287 /* --- @mp_loadl@ --- *
288 *
289 * Arguments: @mp *x@ = pointer to MP head
290 * @const octet *p@ = pointer to octet array
291 * @size_t sz@ = size of octet array
292 *
293 * Returns: ---
294 *
295 * Use: Loads an MP in an octet array, least significant octet
296 * first.
297 */
298
299 void mp_loadl(mp *x, const octet *p, size_t sz)
300 {
301 mpw *v;
302 mpw w = 0;
303 unsigned n;
304 const octet *q = p + sz;
305 unsigned bits = 0;
306
307 MP_BURN(x);
308 MP_ENSURE(x, ((sz * 8 + MPW_BITS - 1) * MPW_BITS) / MPW_BITS);
309 v = x->v;
310
311 while (p < q) {
312 n = *p++ & MASK8;
313 w |= n << bits;
314 bits += 8;
315 if (bits >= MPW_BITS) {
316 *v++ = w & MPW_MAX;
317 w = n >> (MPW_BITS - bits + 8);
318 bits -= MPW_BITS;
319 }
320 }
321 *v++ = w;
322 x->len = v - x->v;
323 x->f = 0;
324 MP_NORM(x);
325 }
326
327 /* --- @mp_storeb@ --- *
328 *
329 * Arguments: @mp *x@ = pointer to MP head
330 * @octet *p@ = pointer to octet array
331 * @size_t sz@ = size of octet array
332 *
333 * Returns: ---
334 *
335 * Use: Stores an MP in an octet array, most significant octet
336 * first. High-end octets are silently discarded if there
337 * isn't enough space for them.
338 */
339
340 void mp_storeb(mp *x, octet *p, size_t sz)
341 {
342 mpw *v = x->v, *vl = x->v + x->len;
343 mpw n, w = 0;
344 octet *q = p + sz;
345 unsigned bits = 0;
346
347 while (q > p) {
348 if (bits < 8) {
349 n = (v >= vl) ? 0 : *v++;
350 *--q = (w | n << bits) & MASK8;
351 w = n >> (8 - bits);
352 bits += MPW_BITS - 8;
353 } else {
354 *--q= w & MASK8;
355 w >>= 8;
356 bits -= 8;
357 }
358 }
359 }
360
361 /* --- @mp_loadb@ --- *
362 *
363 * Arguments: @mp *x@ = pointer to MP head
364 * @const octet *p@ = pointer to octet array
365 * @size_t sz@ = size of octet array
366 *
367 * Returns: ---
368 *
369 * Use: Loads an MP in an octet array, most significant octet
370 * first.
371 */
372
373 void mp_loadb(mp *x, const octet *p, size_t sz)
374 {
375 mpw *v;
376 mpw w = 0;
377 unsigned n;
378 const octet *q = p + sz;
379 unsigned bits = 0;
380
381 MP_BURN(x);
382 MP_ENSURE(x, ((sz * 8 + MPW_BITS - 1) * MPW_BITS) / MPW_BITS);
383 v = x->v;
384
385 while (q > p) {
386 n = *--q & MASK8;
387 w |= n << bits;
388 bits += 8;
389 if (bits >= MPW_BITS) {
390 *v++ = w & MPW_MAX;
391 w = n >> (MPW_BITS - bits + 8);
392 bits -= MPW_BITS;
393 }
394 }
395 *v++ = w;
396 x->len = v - x->v;
397 x->f = 0;
398 MP_NORM(x);
399 }
400
401 /*----- Iterating through bits --------------------------------------------*/
402
403 /* --- @mp_mkbitscan@ --- *
404 *
405 * Arguments: @mp_bitscan *sc@ = pointer to bitscan object
406 * @const mp *x@ = pointer to MP head
407 *
408 * Returns: ---
409 *
410 * Use: Initializes a bitscan object.
411 */
412
413 void mp_mkbitscan(mp_bitscan *sc, const mp *x)
414 {
415 sc->x = x;
416 sc->bits = 0;
417 sc->i = 0;
418 }
419
420 /* --- @mp_bstep@ --- *
421 *
422 * Arguments: @mp_bitscan *sc@ = pointer to bitscanner object
423 *
424 * Returns: Nonzero if there is another bit to read.
425 *
426 * Use: Steps on to the next bit, and tells the caller whether one
427 * exists.
428 */
429
430 int mp_bstep(mp_bitscan *sc)
431 {
432 if (sc->bits) {
433 sc->w >>= 1;
434 sc->bits--;
435 return (1);
436 }
437 if (sc->i >= sc->x->len)
438 return (0);
439 sc->w = sc->x->v[sc->i++];
440 sc->bits = MPW_BITS - 1;
441 return (1);
442 }
443
444 /* --- @mp_bit@ --- *
445 *
446 * Arguments: @const mp_bitscan *sc@ = pointer to bitscanner
447 *
448 * Returns: Current bit value.
449 *
450 * Use: Returns the value of the current bit.
451 */
452
453 int mp_bit(const mp_bitscan *sc) { return (MP_BIT(sc)); }
454
455 /*----- Shifting ----------------------------------------------------------*/
456
457 /* --- @mp_lsl@ --- *
458 *
459 * Arguments: @mp *d@ = pointer to MP head of destination
460 * @const mp *x@ = pointer to MP head of source
461 * @size_t n@ = number of bits to shift
462 *
463 * Returns: ---
464 *
465 * Use: Shifts a number left by a given number of bit positions.
466 */
467
468 void mp_lsl(mp *d, const mp *x, size_t n)
469 {
470 size_t nw = n / MPW_BITS;
471 unsigned nb = n % MPW_BITS;
472 unsigned nr = MPW_BITS - nb;
473 size_t req;
474
475 /* --- Trivial special case --- */
476
477 if (n == 0 || mp_ucmp(x, MP_ZERO) == 0) {
478 mp_copy(d, x);
479 goto copied;
480 }
481
482 /* --- Decide on the about of memory needed in the destination --- */
483
484 req = x->len + nw;
485 if (nb && (x->v[x->len - 1] >> nr) != 0)
486 req++;
487 if (d != x)
488 MP_BURN(d);
489 MP_ENSURE(d, req);
490
491 /* --- Handle single bit shifting --- */
492
493 if (n == 1) {
494 mpw *v = d->v;
495 const mpw *vx = x->v;
496 const mpw *vl;
497 mpw w = 0;
498
499 vl = vx + x->len;
500 while (vx < vl) {
501 mpw t = *vx++;
502 *v++ = ((t << 1) | w) & MPW_MAX;
503 w = t >> (MPW_BITS - 1);
504 }
505 if (v < d->v + req)
506 *v++ = w & MPW_MAX;
507 goto done;
508 }
509
510 /* --- Handle shifts by a multiple of the word size --- *
511 *
512 * This would be easy, except that C is irritating. Shifting an integer
513 * by an amount equal to the type's width yields undefined behaviour;
514 * in particular, under Intel it's a no-op.
515 */
516
517 if (!nb) {
518 mpw *v = d->v;
519 const mpw *vx = x->v;
520 memmove(v + nw, vx, (req - nw) * sizeof(mpw));
521 memset(v, 0, nw * sizeof(mpw));
522 goto done;
523 }
524
525 /* --- Now do the difficult version --- */
526
527 {
528 mpw *v = d->v + req;
529 const mpw *vx = x->v + x->len;
530 const mpw *vl;
531 mpw w = *--vx;
532
533 /* --- First shift the data over --- */
534
535 if (req > x->len + nw)
536 *--v = w >> nr;
537 vl = x->v;
538 while (vx > vl) {
539 mpw t = *--vx;
540 *--v = ((w << nb) | (t >> nr)) & MPW_MAX;
541 w = t;
542 }
543
544 /* --- Deal with tail-end data --- */
545
546 vl = d->v;
547 if (v > vl) {
548 *--v = w << nb;
549 while (v > vl)
550 *--v = 0;
551 }
552 }
553
554 /* --- Common end code --- */
555
556 done:
557 d->len = req;
558 d->f = x->f;
559 copied:
560 MP_NORM(d);
561 }
562
563 /* --- @mp_lsr@ --- *
564 *
565 * Arguments: @mp *d@ = pointer to MP head of destination
566 * @const mp *x@ = pointer to MP head of source
567 * @size_t n@ = number of bits to shift
568 *
569 * Returns: ---
570 *
571 * Use: Shifts a number right by a given number of bit positions.
572 */
573
574 void mp_lsr(mp *d, const mp *x, size_t n)
575 {
576 size_t nw = n / MPW_BITS;
577 unsigned nb = n % MPW_BITS;
578 unsigned nr = MPW_BITS - nb;
579 size_t req;
580
581 /* --- Trivial special case --- */
582
583 if (n == 0 || mp_ucmp(x, MP_ZERO) == 0) {
584 mp_copy(d, x);
585 goto copied;
586 }
587
588 /* --- Decide on the about of memory needed in the destination --- */
589
590 req = x->len - nw;
591 if ((x->v[x->len - 1] >> nb) == 0)
592 req--;
593 if (d != x)
594 MP_BURN(d);
595 MP_ENSURE(d, req);
596
597 /* --- Handle single bit shifting --- */
598
599 if (n == 1) {
600 mpw *v = d->v + req;
601 const mpw *vx = x->v + x->len;
602 const mpw *vl;
603 mpw w = *--vx;
604
605 if (req == x->len)
606 *--v = (w >> 1) & MPW_MAX;
607 vl = x->v;
608 while (vx > vl) {
609 mpw t = *--vx;
610 *--v = (w << (MPW_BITS - 1) | (t >> 1)) & MPW_MAX;
611 w = t;
612 }
613 goto done;
614 }
615
616 /* --- Handle shifts by a multiple of the word size --- *
617 *
618 * This would be easy, except that C is irritating. Shifting an integer
619 * by an amount equal to the type's width yields undefined behaviour;
620 * in particular, under Intel it's a no-op.
621 */
622
623 if (!nb) {
624 mpw *v = d->v;
625 const mpw *vx = x->v;
626 memmove(v, vx + nw, (x->len - nw) * sizeof(mpw));
627 goto done;
628 }
629
630 /* --- Now do the difficult version --- */
631
632 {
633 mpw *v = d->v;
634 const mpw *vx = x->v + nw;
635 const mpw *vl;
636 mpw w = *vx++;
637
638 /* --- First shift the data over --- */
639
640 vl = x->v + x->len;
641 while (vx < vl) {
642 mpw t = *vx++;
643 *v++ = ((w >> nb) | (t << nr)) & MPW_MAX;
644 w = t;
645 }
646 if (req == x->len - nw) {
647 *v++ = (w >> nb) & MPW_MAX;
648 }
649 }
650
651 /* --- Common end code --- */
652
653 done:
654 d->len = req;
655 d->f = x->f;
656 copied:
657 MP_NORM(d);
658 }
659
660 /*----- Adding and subtracting --------------------------------------------*/
661
662 /* --- @mp_uadd@ --- *
663 *
664 * Arguments: @const mp *d@ = pointers to MP head of destination
665 * @const mp *x, *y@ = pointers to MP heads of operands
666 *
667 * Returns: ---
668 *
669 * Use: Performs unsigned MP addition.
670 */
671
672 void mp_uadd(mp *d, const mp *x, const mp *y)
673 {
674 mpd c;
675 mpw *v;
676 const mpw *vx, *vy;
677 const mpw *vxl, *vyl;
678
679 /* --- Some trivial initialization --- */
680
681 if (d != x && d != y)
682 MP_BURN(d);
683 MP_ENSURE(d, (x->len > y->len ? x->len : y->len) + 1);
684 vx = x->v; vxl = vx + x->len;
685 vy = y->v; vyl = vy + y->len;
686 v = d->v;
687 c = 0;
688
689 /* --- Start on the work --- */
690
691 while (vx < vxl || vy < vyl) {
692 if (vx < vxl) c += *vx++;
693 if (vy < vyl) c += *vy++;
694 *v++ = c & MPW_MAX;
695 c >>= MPW_BITS;
696 }
697 if (c)
698 *v++ = c & MPW_MAX;
699
700 /* --- Tidy up --- */
701
702 d->len = v - d->v;
703 d->f = 0;
704 if (x->f & MPF_BURN || y->f & MPF_BURN)
705 MP_PARANOID(d);
706 MP_NORM(d);
707 }
708
709 /* --- @mp_usub@ --- *
710 *
711 * Arguments: @const mp *d@ = pointers to MP head of destination
712 * @const mp *x, *y@ = pointers to MP heads of operands
713 *
714 * Returns: ---
715 *
716 * Use: Performs unsigned MP subtraction.
717 */
718
719 void mp_usub(mp *d, const mp *x, const mp *y)
720 {
721 mpd c;
722 mpw *v;
723 const mpw *vx, *vy;
724 const mpw *vxl, *vyl;
725
726 /* --- Some trivial initialization --- */
727
728 if (d != x && d != y)
729 MP_BURN(d);
730 MP_ENSURE(d, x->len);
731 vx = x->v; vxl = vx + x->len;
732 vy = y->v; vyl = vy + y->len;
733 v = d->v;
734 c = 0;
735
736 /* --- Start on the work --- */
737
738 while (vx < vxl) {
739 c += *vx++;
740 if (vy < vyl) c -= *vy++;
741 *v++ = c & MPW_MAX;
742 if (c & ~MPW_MAX)
743 c = ~0ul;
744 else
745 c = 0;
746 }
747
748 /* --- Tidy up --- */
749
750 d->len = v - d->v;
751 d->f = 0;
752 if (x->f & MPF_BURN || y->f & MPF_BURN)
753 MP_PARANOID(d);
754 MP_NORM(d);
755 }
756
757 /* --- @mp_ucmp@ --- *
758 *
759 * Arguments: @const mp *x, *y@ = pointers to MP heads of operands
760 *
761 * Returns: Less than, equal to, or greater than zero.
762 *
763 * Use: Performs unsigned MP comparison.
764 */
765
766 int mp_ucmp(const mp *x, const mp *y)
767 {
768 int i;
769 mpw a, b;
770
771 /* --- Decide which to examine --- */
772
773 if (x->len > y->len)
774 i = x->len;
775 else
776 i = y->len;
777
778 /* --- Loop through the data --- */
779
780 while (i > 0) {
781 i--;
782 a = (i < x->len ? x->v[i] : 0);
783 b = (i < y->len ? y->v[i] : 0);
784 if (a < b)
785 return (-1);
786 else if (a > b)
787 return (1);
788 }
789
790 /* --- Finished --- */
791
792 return (0);
793 }
794
795 /*----- Multiplying and dividing ------------------------------------------*/
796
797 /* --- @mp_umul@ --- *
798 *
799 * Arguments: @mp *d@ = pointer to MP head of destination
800 * @const mp *x, *y@ = pointes to MP heads of operands
801 *
802 * Returns: ---
803 *
804 * Use: Performs unsigned MP multiplication.
805 */
806
807 void mp_umul(mp *d, const mp *x, const mp *y)
808 {
809 mpd c, f;
810 mpw *v = 0, *vbase;
811 const mpw *vx, *vy;
812 const mpw *vxl, *vyl;
813
814 /* --- Check for special cases --- */
815
816 if (mp_ucmp(x, MP_ZERO) == 0 || mp_ucmp(y, MP_ZERO) == 0) {
817 mp_copy(d, MP_ZERO);
818 return;
819 }
820
821 /* --- Some trivial initialization --- */
822
823 MP_BURN(d);
824 MP_ENSURE(d, x->len + y->len);
825 vxl = x->v + x->len;
826 vyl = y->v + y->len;
827 vbase = d->v;
828 memset(v, 0, (x->len + y->len) * sizeof(mpw));
829
830 /* --- Main loop --- */
831
832 for (vx = x->v; vx < vxl; vx++) {
833 c = 0;
834 f = *vx;
835 v = vbase++;
836 for (vy = y->v; vy < vyl; vy++) {
837 c += *v + f * *vy;
838 *v++ = c & MPW_MAX;
839 c >>= MPW_BITS;
840 }
841 *v++ = c & MPW_MAX;
842 }
843
844 /* --- Tidying up --- */
845
846 d->len = v - d->v;
847 d->f = 0;
848 if (x->f & MPF_BURN || y->f & MPF_BURN)
849 MP_PARANOID(d);
850 MP_NORM(d);
851 }
852
853 /* --- @mp_udiv@ --- *
854 *
855 * Arguments: @mp *q, *r@ = pointers to MP heads for quotient, remainder
856 * @const mp *x, *y@ = pointers to MP heads for operands
857 *
858 * Returns: ---
859 *
860 * Use: Performs unsigned MP division.
861 */
862
863 void mp_udiv(mp *q, mp *r, const mp *x, const mp *y)
864 {
865 size_t n = x->len, t = y->len;
866 const mpw *vy, *ovy;
867 mpw *vx, *vq, *vxl, *va, *vb;
868 mpw yt, ytt;
869 mpd yd;
870 mp nx = MP_INIT, ny = MP_INIT, nq = MP_INIT;
871 mp tt = MP_INIT;
872 unsigned int shift;
873
874 /* --- Fiddle with some pointers --- */
875
876 if (!r)
877 r = &nx;
878 if (!q)
879 q = &nq;
880
881 /* --- Find the top word in @y@ --- */
882
883 vy = y->v;
884 for (;;) {
885 if (!t)
886 THROW(EXC_FAIL, "mp_udiv detected divide-by-zero");
887 if ((yt = vy[t - 1]) != 0)
888 break;
889 t--;
890 }
891
892 /* --- Check for a zero dividend --- */
893
894 if (mp_ucmp(x, MP_ZERO)) {
895 if (q) mp_copy(q, MP_ZERO);
896 if (r) mp_copy(r, MP_ZERO);
897 return;
898 }
899
900 /* --- Test for some other trivial cases --- */
901
902 {
903 int i = mp_ucmp(x, y);
904 if (i < 0) {
905 if (r) mp_copy(r, x);
906 if (q) mp_copy(q, MP_ZERO);
907 return;
908 } else if (i == 0) {
909 if (r) mp_copy(r, MP_ZERO);
910 if (q) mp_copy(q, MP_ONE);
911 return;
912 }
913 }
914
915 /* --- Normalize the divisor --- *
916 *
917 * Cheat. The original algorithm wants two-word values at least, so I just
918 * shift everything up by a word if necessary.
919 */
920
921 shift = 0;
922 while (yt < MPW_MAX / 2) {
923 yt <<= 1;
924 shift++;
925 }
926 if (t <= 1)
927 shift += MPW_BITS;
928 mp_lsl(r, x, shift);
929 mp_lsl(&ny, y, shift);
930
931 n = r->len;
932 t = ny.len;
933 ytt = ny.v[t - 2];
934 yd = ((mpd)yt << MPW_BITS) | (mpd)ytt;
935
936 /* --- Initialize the quotient --- */
937
938 MP_ENSURE(q, n - t + 1);
939 vq = q->v + n - t;
940 memset(vq, 0, (n - t) * sizeof(mpw));
941
942 /* --- Shift the divisor up to match the dividend --- */
943
944 mp_lsl(&ny, (n - t) * MPW_BITS);
945 ovy = ny.v;
946
947 /* --- Get the most significant quotient digit --- *
948 *
949 * Because of the normalization, this should only happen once.
950 */
951
952 while (mp_ucmp(x, y) >= 0) {
953 mp_usub(x, x, y);
954 (*vq)++;
955 }
956
957 /* --- Now do the main loop --- */
958
959 vx = x->v + n;
960 vxl = x->v + t;
961
962 while (vx > vxl) {
963 mpw xi, xii, xiii;
964 mpd qi;
965 mpd xd, yhi, ylo;
966
967 /* --- Fetch the top words from @x@ --- */
968
969 vx--;
970 xi = vx[0];
971 xii = vx[-1];
972 xiii = vx[-2];
973 xd = ((mpd)xi << MPW_BITS) | (mpd)xii;
974
975 /* --- Get an approximation for @qi@ --- */
976
977 if (xi == yt)
978 qi = MPW_MAX;
979 else
980 qi = xd / yt;
981
982 /* --- Work out how close the approximation is --- *
983 *
984 * This is more than a little ugly.
985 */
986
987 yhi = (mpd)yt * qi;
988 ylo = (mpd)ytt * qi;
989 yhi += ylo >> MPW_BITS;
990 ylo &= MPW_MAX;
991
992 /* --- Now fix the approximation --- *
993 *
994 * Again, the normalization helps; this is never done more than twice.
995 */
996
997 while (!(yhi <= xd && ylo <= xiii)) {
998 qi--;
999 yhi -= yt;
1000 if (ylo < ytt)
1001 yhi--;
1002 ylo = (ylo - ytt) & MPW_MAX;
1003 }
1004
1005 /* --- Subtract off a goodly big chunk --- */
1006
1007
1008 }
1009
1010 /*----- Test rig ----------------------------------------------------------*/
1011
1012 #ifdef TEST_RIG
1013
1014 #include <stdio.h>
1015 #include "dstr.h"
1016 #include "testrig.h"
1017
1018 /* --- Loading and storing numbers --- *
1019 *
1020 * The most reliable way I can think of for doing this is to convert both
1021 * numbers (the MP and the original hexgorp) into binary text, and compare.
1022 * This ought to be reliable, and it's sufficiently different from what the
1023 * actual load and store routines are doing not to have any common bugs.
1024 */
1025
1026 static void mpbindump(const mp *x, char *p)
1027 {
1028 mp_bitscan sc;
1029 for (mp_mkbitscan(&sc, x); mp_bstep(&sc); )
1030 *p++ = '0' + MP_BIT(&sc);
1031 *p++ = 0;
1032 }
1033
1034 static void bindumpl(const octet *q, size_t sz, char *p)
1035 {
1036 unsigned w;
1037 int bits = 0;
1038
1039 while (sz) {
1040 w = *q++;
1041 for (bits = 8; bits > 0; bits--) {
1042 *p++ = '0' + (w & 1);
1043 w >>= 1;
1044 }
1045 sz--;
1046 }
1047 *p++ = 0;
1048 }
1049
1050 static void bindumpb(const octet *q, size_t sz, char *p)
1051 {
1052 unsigned w;
1053 int bits = 0;
1054
1055 q += sz;
1056 while (sz) {
1057 w = *--q;
1058 for (bits = 8; bits > 0; bits--) {
1059 *p++ = '0' + (w & 1);
1060 w >>= 1;
1061 }
1062 sz--;
1063 }
1064 *p++ = 0;
1065 }
1066
1067 static int bincmp(const char *p, const char *q)
1068 {
1069 for (;;) {
1070 if (!*p && !*q)
1071 return (0);
1072 else if (!*p && *q == '0')
1073 q++;
1074 else if (*p == '0' && !*q)
1075 p++;
1076 else if (*p == *q)
1077 p++, q++;
1078 else
1079 return (-1);
1080 }
1081 }
1082
1083 static int lscheck(const char *reason, int w, dstr *d, mp *x,
1084 const char *bufa, const char *bufb)
1085 {
1086 if (bincmp(bufa, bufb)) {
1087 printf("\nmismatch in %s (width = %i):\n"
1088 "\tvalue = ", reason, w);
1089 type_hex.dump(d, stdout);
1090 printf("\n\tmp = ");
1091 mp_dump(x, stdout);
1092 printf("\n\texpected = %s\n"
1093 "\tcalculated = %s\n",
1094 bufb, bufa);
1095 return (1);
1096 }
1097 return (0);
1098 }
1099
1100 static int loadstore(dstr *d)
1101 {
1102 octet *p = (octet *)d[0].buf;
1103 size_t sz = d[0].len;
1104 char bufa[1024], bufb[1024];
1105 octet bufc[64];
1106 mp x;
1107 int i;
1108 int ok = 1;
1109
1110 mp_wmax = 0x7ful;
1111 for (i = 8; ok && i <= 32; i++) {
1112 mpw_bits = i;
1113 mpw_max = (mp_wmax << 1) + 1;
1114
1115 mp_create(&x); mp_loadl(&x, p, sz);
1116 mpbindump(&x, bufa); bindumpl(p, sz, bufb);
1117 if (lscheck("mp_loadl", i, d, &x, bufa, bufb)) ok = 0;
1118 mp_storeb(&x, bufc, sizeof(bufc));
1119 bindumpb(bufc, sizeof(bufc), bufa);
1120 if (lscheck("mp_storeb", i, d, &x, bufa, bufb)) ok = 0;
1121 mp_destroy(&x);
1122
1123 mp_create(&x); mp_loadb(&x, p, sz);
1124 mpbindump(&x, bufa); bindumpb(p, sz, bufb);
1125 if (lscheck("mp_loadb", i, d, &x, bufa, bufb)) ok = 0;
1126 mp_storel(&x, bufc, sizeof(bufc));
1127 bindumpl(bufc, sizeof(bufc), bufa);
1128 if (lscheck("mp_storel", i, d, &x, bufa, bufb)) ok = 0;
1129 mp_destroy(&x);
1130 }
1131
1132 mpw_max = 0xffff;
1133 mpw_bits = 16;
1134 return (ok);
1135 }
1136
1137 /* --- Comparison test --- */
1138
1139 static int compare(dstr *d)
1140 {
1141 mp x, y;
1142 int r, s;
1143 int ok = 1;
1144
1145 mp_create(&x); mp_create(&y);
1146 mp_loadb(&x, (octet *)d[0].buf, d[0].len);
1147 mp_loadb(&y, (octet *)d[1].buf, d[1].len);
1148 r = *(int *)d[2].buf;
1149 s = mp_ucmp(&x, &y);
1150
1151 if (r != s) {
1152 ok = 0;
1153 printf("\nfailed compare:"
1154 "\n\tx = ");
1155 type_hex.dump(&d[0], stdout);
1156 printf("\n\ty = ");
1157 type_hex.dump(&d[1], stdout);
1158 printf("\n\texpected %i; found %i\n", r, s);
1159 }
1160
1161 return (ok);
1162 }
1163
1164 /* --- Addition and subtraction test --- */
1165
1166 static int addsub(dstr *d)
1167 {
1168 mp t, x, y, z;
1169 int ok = 1;
1170
1171 mp_create(&t); mp_create(&x); mp_create(&y); mp_create(&z);
1172 mp_loadb(&x, (octet *)d[0].buf, d[0].len);
1173 mp_loadb(&y, (octet *)d[1].buf, d[1].len);
1174 mp_loadb(&z, (octet *)d[2].buf, d[2].len);
1175
1176 mp_uadd(&t, &x, &y);
1177 if (mp_ucmp(&t, &z)) {
1178 ok = 0;
1179 printf("\nfailed add:");
1180 printf("\n\tx = "); type_hex.dump(&d[0], stdout);
1181 printf("\n\ty = "); type_hex.dump(&d[1], stdout);
1182 printf("\n\tz = "); type_hex.dump(&d[2], stdout);
1183 fputc('\n', stdout);
1184 }
1185
1186 mp_usub(&t, &z, &x);
1187 if (mp_ucmp(&t, &y)) {
1188 ok = 0;
1189 printf("\nfailed subtract:");
1190 printf("\n\tz = "); type_hex.dump(&d[2], stdout);
1191 printf("\n\tx = "); type_hex.dump(&d[0], stdout);
1192 printf("\n\ty = "); type_hex.dump(&d[1], stdout);
1193 fputc('\n', stdout);
1194 }
1195
1196 return (ok);
1197 }
1198
1199 /* --- Shifting --- */
1200
1201 static int shift(dstr *d)
1202 {
1203 char bufa[1024], bufb[1024];
1204 mp x, y;
1205 int n = *(int *)d[1].buf;
1206 char *p;
1207 int i;
1208 int ok = 1;
1209
1210 mp_create(&x);
1211 mp_create(&y);
1212
1213 mp_loadb(&x, (octet *)d[0].buf, d[0].len);
1214 p = bufa;
1215 for (i = 0; i < n; i++)
1216 *p++ = '0';
1217 mpbindump(&x, p);
1218
1219 mp_lsl(&y, &x, n);
1220 mpbindump(&y, bufb);
1221 if (lscheck("lsl", n, d, &x, bufb, bufa))
1222 ok = 0;
1223
1224 mp_lsr(&y, &x, n);
1225 mpbindump(&y, bufb);
1226 if (lscheck("lsr", n, d, &x, bufb, bufa + 2 * n))
1227 ok = 0;
1228
1229 return (ok);
1230 }
1231
1232 /* --- Test driver stub --- */
1233
1234 static test_chunk mp__defs[] = {
1235 { "mp-loadstore", loadstore, { &type_hex, 0 } },
1236 { "mp-cmp", compare, { &type_hex, &type_hex, &type_int, 0 } },
1237 { "mp-addsub", addsub, { &type_hex, &type_hex, &type_hex, 0 } },
1238 { "mp-shift", shift, { &type_hex, &type_int, 0 } },
1239 { 0, 0, { 0 } }
1240 };
1241
1242 int main(int argc, char *argv[])
1243 {
1244 test_run(argc, argv, mp__defs, SRCDIR"/tests/mp");
1245 return (0);
1246 }
1247
1248 #endif
1249
1250 /*----- That's all, folks -------------------------------------------------*/
1251
1252 #endif
Catacomb cryptographic library (downstream of http://git.distorted.org.uk/~mdw/catacomb/)