d03ab969 |
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 |