projects / u / mdw / catacomb / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Keep quiet about expected errors on incoming connections.
[u/mdw/catacomb] / mpx.c
CommitLineData
d03ab969 1/* -*-c-*-
2 *
81578196 3 * $Id: mpx.c,v 1.16 2003/05/16 09:09:24 mdw Exp $
d03ab969 4 *
5 * Low-level multiprecision arithmetic
6 *
7 * (c) 1999 Straylight/Edgeware
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: mpx.c,v $
81578196 33 * Revision 1.16 2003/05/16 09:09:24 mdw
34 * Fix @mp_lsl2c@. Turns out to be surprisingly tricky.
35 *
2bd53494 36 * Revision 1.15 2002/10/20 01:12:31 mdw
37 * Two's complement I/O fixes.
38 *
4f29a732 39 * Revision 1.14 2002/10/19 18:55:08 mdw
40 * Fix overflows in shift primitives.
41 *
75263f25 42 * Revision 1.13 2002/10/19 17:56:50 mdw
43 * Fix bit operations. Test them (a bit) better.
44 *
f09e814a 45 * Revision 1.12 2002/10/06 22:52:50 mdw
46 * Pile of changes for supporting two's complement properly.
47 *
0f32e0f8 48 * Revision 1.11 2001/04/03 19:36:05 mdw
49 * Add some simple bitwise operations so that Perl can use them.
50 *
1a05a8ef 51 * Revision 1.10 2000/10/08 12:06:12 mdw
52 * Provide @mpx_ueq@ for rapidly testing equality of two integers.
53 *
99b30c23 54 * Revision 1.9 2000/06/26 07:52:50 mdw
55 * Portability fix for the bug fix.
56 *
ce76ff16 57 * Revision 1.8 2000/06/25 12:59:02 mdw
58 * (mpx_udiv): Fix bug in quotient digit estimation.
59 *
698bd937 60 * Revision 1.7 1999/12/22 15:49:07 mdw
61 * New function for division by a small integer.
62 *
42684bdb 63 * Revision 1.6 1999/11/20 22:43:44 mdw
64 * Integrate testing for MPX routines.
65 *
dd517851 66 * Revision 1.5 1999/11/20 22:23:27 mdw
67 * Add function versions of some low-level macros with wider use.
68 *
f45a00c6 69 * Revision 1.4 1999/11/17 18:04:09 mdw
70 * Add two's-complement functionality. Improve mpx_udiv a little by
71 * performing the multiplication of the divisor by q with the subtraction
72 * from r.
73 *
501da53c 74 * Revision 1.3 1999/11/13 01:57:31 mdw
75 * Remove stray debugging code.
76 *
c8a2f9ef 77 * Revision 1.2 1999/11/13 01:50:59 mdw
78 * Multiprecision routines finished and tested.
79 *
d03ab969 80 * Revision 1.1 1999/09/03 08:41:12 mdw
81 * Initial import.
82 *
83 */
84
85/*----- Header files ------------------------------------------------------*/
86
c8a2f9ef 87#include <assert.h>
d03ab969 88#include <stdio.h>
89#include <stdlib.h>
90#include <string.h>
91
92#include <mLib/bits.h>
93
94#include "mptypes.h"
95#include "mpx.h"
75263f25 96#include "bitops.h"
d03ab969 97
98/*----- Loading and storing -----------------------------------------------*/
99
100/* --- @mpx_storel@ --- *
101 *
102 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
c8a2f9ef 103 * @void *pp@ = pointer to octet array
d03ab969 104 * @size_t sz@ = size of octet array
105 *
106 * Returns: ---
107 *
108 * Use: Stores an MP in an octet array, least significant octet
109 * first. High-end octets are silently discarded if there
110 * isn't enough space for them.
111 */
112
c8a2f9ef 113void mpx_storel(const mpw *v, const mpw *vl, void *pp, size_t sz)
d03ab969 114{
115 mpw n, w = 0;
c8a2f9ef 116 octet *p = pp, *q = p + sz;
d03ab969 117 unsigned bits = 0;
118
119 while (p < q) {
120 if (bits < 8) {
121 if (v >= vl) {
122 *p++ = U8(w);
123 break;
124 }
125 n = *v++;
126 *p++ = U8(w | n << bits);
127 w = n >> (8 - bits);
128 bits += MPW_BITS - 8;
129 } else {
130 *p++ = U8(w);
131 w >>= 8;
132 bits -= 8;
133 }
134 }
135 memset(p, 0, q - p);
136}
137
138/* --- @mpx_loadl@ --- *
139 *
140 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
c8a2f9ef 141 * @const void *pp@ = pointer to octet array
d03ab969 142 * @size_t sz@ = size of octet array
143 *
144 * Returns: ---
145 *
146 * Use: Loads an MP in an octet array, least significant octet
147 * first. High-end octets are ignored if there isn't enough
148 * space for them.
149 */
150
c8a2f9ef 151void mpx_loadl(mpw *v, mpw *vl, const void *pp, size_t sz)
d03ab969 152{
153 unsigned n;
c8a2f9ef 154 mpw w = 0;
155 const octet *p = pp, *q = p + sz;
d03ab969 156 unsigned bits = 0;
157
158 if (v >= vl)
159 return;
160 while (p < q) {
161 n = U8(*p++);
162 w |= n << bits;
163 bits += 8;
164 if (bits >= MPW_BITS) {
165 *v++ = MPW(w);
166 w = n >> (MPW_BITS - bits + 8);
167 bits -= MPW_BITS;
168 if (v >= vl)
169 return;
170 }
171 }
172 *v++ = w;
173 MPX_ZERO(v, vl);
174}
175
176/* --- @mpx_storeb@ --- *
177 *
178 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
c8a2f9ef 179 * @void *pp@ = pointer to octet array
d03ab969 180 * @size_t sz@ = size of octet array
181 *
182 * Returns: ---
183 *
184 * Use: Stores an MP in an octet array, most significant octet
185 * first. High-end octets are silently discarded if there
186 * isn't enough space for them.
187 */
188
c8a2f9ef 189void mpx_storeb(const mpw *v, const mpw *vl, void *pp, size_t sz)
d03ab969 190{
191 mpw n, w = 0;
c8a2f9ef 192 octet *p = pp, *q = p + sz;
d03ab969 193 unsigned bits = 0;
194
195 while (q > p) {
196 if (bits < 8) {
197 if (v >= vl) {
198 *--q = U8(w);
199 break;
200 }
201 n = *v++;
202 *--q = U8(w | n << bits);
203 w = n >> (8 - bits);
204 bits += MPW_BITS - 8;
205 } else {
206 *--q = U8(w);
207 w >>= 8;
208 bits -= 8;
209 }
210 }
211 memset(p, 0, q - p);
212}
213
214/* --- @mpx_loadb@ --- *
215 *
216 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
c8a2f9ef 217 * @const void *pp@ = pointer to octet array
d03ab969 218 * @size_t sz@ = size of octet array
219 *
220 * Returns: ---
221 *
222 * Use: Loads an MP in an octet array, most significant octet
223 * first. High-end octets are ignored if there isn't enough
224 * space for them.
225 */
226
c8a2f9ef 227void mpx_loadb(mpw *v, mpw *vl, const void *pp, size_t sz)
d03ab969 228{
229 unsigned n;
c8a2f9ef 230 mpw w = 0;
231 const octet *p = pp, *q = p + sz;
d03ab969 232 unsigned bits = 0;
233
234 if (v >= vl)
235 return;
236 while (q > p) {
237 n = U8(*--q);
238 w |= n << bits;
239 bits += 8;
240 if (bits >= MPW_BITS) {
241 *v++ = MPW(w);
242 w = n >> (MPW_BITS - bits + 8);
243 bits -= MPW_BITS;
244 if (v >= vl)
245 return;
246 }
247 }
248 *v++ = w;
249 MPX_ZERO(v, vl);
250}
251
f09e814a 252/* --- @mpx_storel2cn@ --- *
253 *
254 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
255 * @void *pp@ = pointer to octet array
256 * @size_t sz@ = size of octet array
257 *
258 * Returns: ---
259 *
260 * Use: Stores a negative MP in an octet array, least significant
261 * octet first, as two's complement. High-end octets are
262 * silently discarded if there isn't enough space for them.
263 * This obviously makes the output bad.
264 */
265
266void mpx_storel2cn(const mpw *v, const mpw *vl, void *pp, size_t sz)
267{
268 unsigned c = 1;
269 unsigned b = 0;
270 mpw n, w = 0;
271 octet *p = pp, *q = p + sz;
272 unsigned bits = 0;
273
274 while (p < q) {
275 if (bits < 8) {
276 if (v >= vl) {
277 b = w;
278 break;
279 }
280 n = *v++;
281 b = w | n << bits;
282 w = n >> (8 - bits);
283 bits += MPW_BITS - 8;
284 } else {
285 b = w;
286 w >>= 8;
287 bits -= 8;
288 }
289 b = U8(~b + c);
2bd53494 290 c = c && !b;
f09e814a 291 *p++ = b;
292 }
293 while (p < q) {
294 b = U8(~b + c);
2bd53494 295 c = c && !b;
f09e814a 296 *p++ = b;
297 b = 0;
298 }
299}
300
301/* --- @mpx_loadl2cn@ --- *
302 *
303 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
304 * @const void *pp@ = pointer to octet array
305 * @size_t sz@ = size of octet array
306 *
307 * Returns: ---
308 *
309 * Use: Loads a negative MP in an octet array, least significant
310 * octet first, as two's complement. High-end octets are
311 * ignored if there isn't enough space for them. This probably
312 * means you made the wrong choice coming here.
313 */
314
315void mpx_loadl2cn(mpw *v, mpw *vl, const void *pp, size_t sz)
316{
317 unsigned n;
318 unsigned c = 1;
319 mpw w = 0;
320 const octet *p = pp, *q = p + sz;
321 unsigned bits = 0;
322
323 if (v >= vl)
324 return;
325 while (p < q) {
326 n = U8(~(*p++) + c);
2bd53494 327 c = c && !n;
f09e814a 328 w |= n << bits;
329 bits += 8;
330 if (bits >= MPW_BITS) {
331 *v++ = MPW(w);
332 w = n >> (MPW_BITS - bits + 8);
333 bits -= MPW_BITS;
334 if (v >= vl)
335 return;
336 }
337 }
338 *v++ = w;
339 MPX_ZERO(v, vl);
340}
341
342/* --- @mpx_storeb2cn@ --- *
343 *
344 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
345 * @void *pp@ = pointer to octet array
346 * @size_t sz@ = size of octet array
347 *
348 * Returns: ---
349 *
350 * Use: Stores a negative MP in an octet array, most significant
351 * octet first, as two's complement. High-end octets are
352 * silently discarded if there isn't enough space for them,
353 * which probably isn't what you meant.
354 */
355
356void mpx_storeb2cn(const mpw *v, const mpw *vl, void *pp, size_t sz)
357{
358 mpw n, w = 0;
359 unsigned b = 0;
360 unsigned c = 1;
361 octet *p = pp, *q = p + sz;
362 unsigned bits = 0;
363
364 while (q > p) {
365 if (bits < 8) {
366 if (v >= vl) {
367 b = w;
368 break;
369 }
370 n = *v++;
371 b = w | n << bits;
372 w = n >> (8 - bits);
373 bits += MPW_BITS - 8;
374 } else {
375 b = w;
376 w >>= 8;
377 bits -= 8;
378 }
379 b = U8(~b + c);
2bd53494 380 c = c && !b;
f09e814a 381 *--q = b;
382 }
383 while (q > p) {
384 b = ~b + c;
2bd53494 385 c = c && !(b & 0xff);
f09e814a 386 *--q = b;
387 b = 0;
388 }
389}
390
391/* --- @mpx_loadb2cn@ --- *
392 *
393 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
394 * @const void *pp@ = pointer to octet array
395 * @size_t sz@ = size of octet array
396 *
397 * Returns: ---
398 *
399 * Use: Loads a negative MP in an octet array, most significant octet
400 * first as two's complement. High-end octets are ignored if
401 * there isn't enough space for them. This probably means you
402 * chose this function wrongly.
403 */
404
405void mpx_loadb2cn(mpw *v, mpw *vl, const void *pp, size_t sz)
406{
407 unsigned n;
408 unsigned c = 1;
409 mpw w = 0;
410 const octet *p = pp, *q = p + sz;
411 unsigned bits = 0;
412
413 if (v >= vl)
414 return;
415 while (q > p) {
416 n = U8(~(*--q) + c);
2bd53494 417 c = c && !n;
f09e814a 418 w |= n << bits;
419 bits += 8;
420 if (bits >= MPW_BITS) {
421 *v++ = MPW(w);
422 w = n >> (MPW_BITS - bits + 8);
423 bits -= MPW_BITS;
424 if (v >= vl)
425 return;
426 }
427 }
428 *v++ = w;
429 MPX_ZERO(v, vl);
430}
431
d03ab969 432/*----- Logical shifting --------------------------------------------------*/
433
434/* --- @mpx_lsl@ --- *
435 *
436 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
437 * @const mpw *av, *avl@ = source vector base and limit
438 * @size_t n@ = number of bit positions to shift by
439 *
440 * Returns: ---
441 *
442 * Use: Performs a logical shift left operation on an integer.
443 */
444
445void mpx_lsl(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
446{
447 size_t nw;
448 unsigned nb;
449
450 /* --- Trivial special case --- */
451
452 if (n == 0)
453 MPX_COPY(dv, dvl, av, avl);
454
455 /* --- Single bit shifting --- */
456
457 else if (n == 1) {
458 mpw w = 0;
459 while (av < avl) {
460 mpw t;
461 if (dv >= dvl)
462 goto done;
463 t = *av++;
464 *dv++ = MPW((t << 1) | w);
465 w = t >> (MPW_BITS - 1);
466 }
467 if (dv >= dvl)
468 goto done;
469 *dv++ = MPW(w);
470 MPX_ZERO(dv, dvl);
c8a2f9ef 471 goto done;
d03ab969 472 }
473
474 /* --- Break out word and bit shifts for more sophisticated work --- */
475
476 nw = n / MPW_BITS;
477 nb = n % MPW_BITS;
478
479 /* --- Handle a shift by a multiple of the word size --- */
480
481 if (nb == 0) {
4f29a732 482 if (nw >= dvl - dv)
483 MPX_ZERO(dv, dvl);
484 else {
485 MPX_COPY(dv + nw, dvl, av, avl);
486 memset(dv, 0, MPWS(nw));
487 }
d03ab969 488 }
489
c8a2f9ef 490 /* --- And finally the difficult case --- *
491 *
492 * This is a little convoluted, because I have to start from the end and
493 * work backwards to avoid overwriting the source, if they're both the same
494 * block of memory.
495 */
d03ab969 496
497 else {
498 mpw w;
499 size_t nr = MPW_BITS - nb;
c8a2f9ef 500 size_t dvn = dvl - dv;
501 size_t avn = avl - av;
d03ab969 502
c8a2f9ef 503 if (dvn <= nw) {
d03ab969 504 MPX_ZERO(dv, dvl);
505 goto done;
506 }
d03ab969 507
c8a2f9ef 508 if (dvn > avn + nw) {
509 size_t off = avn + nw + 1;
510 MPX_ZERO(dv + off, dvl);
511 dvl = dv + off;
512 w = 0;
513 } else {
514 avl = av + dvn - nw;
515 w = *--avl << nb;
d03ab969 516 }
517
c8a2f9ef 518 while (avl > av) {
519 mpw t = *--avl;
520 *--dvl = (t >> nr) | w;
521 w = t << nb;
d03ab969 522 }
c8a2f9ef 523
524 *--dvl = w;
525 MPX_ZERO(dv, dvl);
d03ab969 526 }
527
528done:;
529}
530
81578196 531/* --- @mpx_lslc@ --- *
532 *
533 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
534 * @const mpw *av, *avl@ = source vector base and limit
535 * @size_t n@ = number of bit positions to shift by
536 *
537 * Returns: ---
538 *
539 * Use: Performs a logical shift left operation on an integer, only
540 * it fills in the bits with ones instead of zeroes.
541 */
542
543void mpx_lslc(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
544{
545 size_t nw;
546 unsigned nb;
547
548 /* --- Trivial special case --- */
549
550 if (n == 0)
551 MPX_COPY(dv, dvl, av, avl);
552
553 /* --- Single bit shifting --- */
554
555 else if (n == 1) {
556 mpw w = 1;
557 while (av < avl) {
558 mpw t;
559 if (dv >= dvl)
560 goto done;
561 t = *av++;
562 *dv++ = MPW((t << 1) | w);
563 w = t >> (MPW_BITS - 1);
564 }
565 if (dv >= dvl)
566 goto done;
567 *dv++ = MPW(w);
568 MPX_ZERO(dv, dvl);
569 goto done;
570 }
571
572 /* --- Break out word and bit shifts for more sophisticated work --- */
573
574 nw = n / MPW_BITS;
575 nb = n % MPW_BITS;
576
577 /* --- Handle a shift by a multiple of the word size --- */
578
579 if (nb == 0) {
580 if (nw >= dvl - dv)
581 MPX_ONE(dv, dvl);
582 else {
583 MPX_COPY(dv + nw, dvl, av, avl);
584 MPX_ONE(dv, dv + nw);
585 }
586 }
587
588 /* --- And finally the difficult case --- *
589 *
590 * This is a little convoluted, because I have to start from the end and
591 * work backwards to avoid overwriting the source, if they're both the same
592 * block of memory.
593 */
594
595 else {
596 mpw w;
597 size_t nr = MPW_BITS - nb;
598 size_t dvn = dvl - dv;
599 size_t avn = avl - av;
600
601 if (dvn <= nw) {
602 MPX_ONE(dv, dvl);
603 goto done;
604 }
605
606 if (dvn > avn + nw) {
607 size_t off = avn + nw + 1;
608 MPX_ZERO(dv + off, dvl);
609 dvl = dv + off;
610 w = 0;
611 } else {
612 avl = av + dvn - nw;
613 w = *--avl << nb;
614 }
615
616 while (avl > av) {
617 mpw t = *--avl;
618 *--dvl = (t >> nr) | w;
619 w = t << nb;
620 }
621
622 *--dvl = (MPW_MAX >> nr) | w;
623 MPX_ONE(dv, dvl);
624 }
625
626done:;
627}
628
d03ab969 629/* --- @mpx_lsr@ --- *
630 *
631 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
632 * @const mpw *av, *avl@ = source vector base and limit
633 * @size_t n@ = number of bit positions to shift by
634 *
635 * Returns: ---
636 *
637 * Use: Performs a logical shift right operation on an integer.
638 */
639
640void mpx_lsr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
641{
642 size_t nw;
643 unsigned nb;
644
645 /* --- Trivial special case --- */
646
647 if (n == 0)
648 MPX_COPY(dv, dvl, av, avl);
649
650 /* --- Single bit shifting --- */
651
652 else if (n == 1) {
653 mpw w = *av++ >> 1;
654 while (av < avl) {
655 mpw t;
656 if (dv >= dvl)
657 goto done;
658 t = *av++;
659 *dv++ = MPW((t << (MPW_BITS - 1)) | w);
660 w = t >> 1;
661 }
662 if (dv >= dvl)
663 goto done;
664 *dv++ = MPW(w);
665 MPX_ZERO(dv, dvl);
c8a2f9ef 666 goto done;
d03ab969 667 }
668
669 /* --- Break out word and bit shifts for more sophisticated work --- */
670
671 nw = n / MPW_BITS;
672 nb = n % MPW_BITS;
673
674 /* --- Handle a shift by a multiple of the word size --- */
675
4f29a732 676 if (nb == 0) {
677 if (nw >= avl - av)
678 MPX_ZERO(dv, dvl);
679 else
680 MPX_COPY(dv, dvl, av + nw, avl);
681 }
d03ab969 682
683 /* --- And finally the difficult case --- */
684
685 else {
686 mpw w;
687 size_t nr = MPW_BITS - nb;
688
689 av += nw;
4f29a732 690 w = av < avl ? *av++ : 0;
d03ab969 691 while (av < avl) {
692 mpw t;
693 if (dv >= dvl)
694 goto done;
695 t = *av++;
696 *dv++ = MPW((w >> nb) | (t << nr));
697 w = t;
698 }
699 if (dv < dvl) {
700 *dv++ = MPW(w >> nb);
701 MPX_ZERO(dv, dvl);
702 }
703 }
704
705done:;
706}
707
0f32e0f8 708/*----- Bitwise operations ------------------------------------------------*/
709
f09e814a 710/* --- @mpx_bitop@ --- *
0f32e0f8 711 *
712 * Arguments: @mpw *dv, *dvl@ = destination vector
713 * @const mpw *av, *avl@ = first source vector
714 * @const mpw *bv, *bvl@ = second source vector
715 *
716 * Returns: ---
717 *
f09e814a 718 * Use; Provides the dyadic boolean functions.
0f32e0f8 719 */
720
f09e814a 721#define MPX_BITBINOP(string) \
0f32e0f8 722 \
f09e814a 723void mpx_bit##string(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, \
724 const mpw *bv, const mpw *bvl) \
0f32e0f8 725{ \
726 MPX_SHRINK(av, avl); \
727 MPX_SHRINK(bv, bvl); \
728 \
729 while (dv < dvl) { \
730 mpw a, b; \
731 a = (av < avl) ? *av++ : 0; \
732 b = (bv < bvl) ? *bv++ : 0; \
75263f25 733 *dv++ = B##string(a, b); \
0f32e0f8 734 } \
735}
736
f09e814a 737MPX_DOBIN(MPX_BITBINOP)
0f32e0f8 738
739void mpx_not(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
740{
741 MPX_SHRINK(av, avl);
742
743 while (dv < dvl) {
744 mpw a;
745 a = (av < avl) ? *av++ : 0;
746 *dv++ = ~a;
747 }
748}
749
d03ab969 750/*----- Unsigned arithmetic -----------------------------------------------*/
751
f45a00c6 752/* --- @mpx_2c@ --- *
753 *
754 * Arguments: @mpw *dv, *dvl@ = destination vector
755 * @const mpw *v, *vl@ = source vector
756 *
757 * Returns: ---
758 *
759 * Use: Calculates the two's complement of @v@.
760 */
761
762void mpx_2c(mpw *dv, mpw *dvl, const mpw *v, const mpw *vl)
763{
764 mpw c = 0;
765 while (dv < dvl && v < vl)
766 *dv++ = c = MPW(~*v++);
767 if (dv < dvl) {
768 if (c > MPW_MAX / 2)
769 c = MPW(~0);
770 while (dv < dvl)
771 *dv++ = c;
772 }
773 MPX_UADDN(dv, dvl, 1);
774}
775
1a05a8ef 776/* --- @mpx_ueq@ --- *
777 *
778 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
779 * @const mpw *bv, *bvl@ = second argument vector base and limit
780 *
781 * Returns: Nonzero if the two vectors are equal.
782 *
783 * Use: Performs an unsigned integer test for equality.
784 */
785
786int mpx_ueq(const mpw *av, const mpw *avl, const mpw *bv, const mpw *bvl)
787{
788 MPX_SHRINK(av, avl);
789 MPX_SHRINK(bv, bvl);
790 if (avl - av != bvl - bv)
791 return (0);
792 while (av < avl) {
793 if (*av++ != *bv++)
794 return (0);
795 }
796 return (1);
797}
798
d03ab969 799/* --- @mpx_ucmp@ --- *
800 *
801 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
802 * @const mpw *bv, *bvl@ = second argument vector base and limit
803 *
804 * Returns: Less than, equal to, or greater than zero depending on
805 * whether @a@ is less than, equal to or greater than @b@,
806 * respectively.
807 *
808 * Use: Performs an unsigned integer comparison.
809 */
810
811int mpx_ucmp(const mpw *av, const mpw *avl, const mpw *bv, const mpw *bvl)
812{
813 MPX_SHRINK(av, avl);
814 MPX_SHRINK(bv, bvl);
815
816 if (avl - av > bvl - bv)
817 return (+1);
818 else if (avl - av < bvl - bv)
819 return (-1);
820 else while (avl > av) {
821 mpw a = *--avl, b = *--bvl;
822 if (a > b)
823 return (+1);
824 else if (a < b)
825 return (-1);
826 }
827 return (0);
828}
1a05a8ef 829
d03ab969 830/* --- @mpx_uadd@ --- *
831 *
832 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
833 * @const mpw *av, *avl@ = first addend vector base and limit
834 * @const mpw *bv, *bvl@ = second addend vector base and limit
835 *
836 * Returns: ---
837 *
838 * Use: Performs unsigned integer addition. If the result overflows
839 * the destination vector, high-order bits are discarded. This
840 * means that two's complement addition happens more or less for
841 * free, although that's more a side-effect than anything else.
842 * The result vector may be equal to either or both source
843 * vectors, but may not otherwise overlap them.
844 */
845
846void mpx_uadd(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
847 const mpw *bv, const mpw *bvl)
848{
849 mpw c = 0;
850
851 while (av < avl || bv < bvl) {
852 mpw a, b;
853 mpd x;
854 if (dv >= dvl)
855 return;
856 a = (av < avl) ? *av++ : 0;
857 b = (bv < bvl) ? *bv++ : 0;
858 x = (mpd)a + (mpd)b + c;
859 *dv++ = MPW(x);
860 c = x >> MPW_BITS;
861 }
862 if (dv < dvl) {
863 *dv++ = c;
864 MPX_ZERO(dv, dvl);
865 }
866}
867
dd517851 868/* --- @mpx_uaddn@ --- *
869 *
870 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
871 * @mpw n@ = other addend
872 *
873 * Returns: ---
874 *
875 * Use: Adds a small integer to a multiprecision number.
876 */
877
878void mpx_uaddn(mpw *dv, mpw *dvl, mpw n) { MPX_UADDN(dv, dvl, n); }
879
d03ab969 880/* --- @mpx_usub@ --- *
881 *
882 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
883 * @const mpw *av, *avl@ = first argument vector base and limit
884 * @const mpw *bv, *bvl@ = second argument vector base and limit
885 *
886 * Returns: ---
887 *
888 * Use: Performs unsigned integer subtraction. If the result
889 * overflows the destination vector, high-order bits are
890 * discarded. This means that two's complement subtraction
891 * happens more or less for free, althuogh that's more a side-
892 * effect than anything else. The result vector may be equal to
893 * either or both source vectors, but may not otherwise overlap
894 * them.
895 */
896
897void mpx_usub(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
898 const mpw *bv, const mpw *bvl)
899{
900 mpw c = 0;
901
902 while (av < avl || bv < bvl) {
903 mpw a, b;
904 mpd x;
905 if (dv >= dvl)
906 return;
907 a = (av < avl) ? *av++ : 0;
908 b = (bv < bvl) ? *bv++ : 0;
c8a2f9ef 909 x = (mpd)a - (mpd)b - c;
d03ab969 910 *dv++ = MPW(x);
c8a2f9ef 911 if (x >> MPW_BITS)
912 c = 1;
913 else
914 c = 0;
d03ab969 915 }
c8a2f9ef 916 if (c)
917 c = MPW_MAX;
d03ab969 918 while (dv < dvl)
c8a2f9ef 919 *dv++ = c;
d03ab969 920}
921
dd517851 922/* --- @mpx_usubn@ --- *
923 *
924 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
925 * @n@ = subtrahend
926 *
927 * Returns: ---
928 *
929 * Use: Subtracts a small integer from a multiprecision number.
930 */
931
932void mpx_usubn(mpw *dv, mpw *dvl, mpw n) { MPX_USUBN(dv, dvl, n); }
933
d03ab969 934/* --- @mpx_umul@ --- *
935 *
936 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
937 * @const mpw *av, *avl@ = multiplicand vector base and limit
938 * @const mpw *bv, *bvl@ = multiplier vector base and limit
939 *
940 * Returns: ---
941 *
942 * Use: Performs unsigned integer multiplication. If the result
943 * overflows the desination vector, high-order bits are
944 * discarded. The result vector may not overlap the argument
945 * vectors in any way.
946 */
947
948void mpx_umul(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
949 const mpw *bv, const mpw *bvl)
950{
951 /* --- This is probably worthwhile on a multiply --- */
952
953 MPX_SHRINK(av, avl);
954 MPX_SHRINK(bv, bvl);
955
956 /* --- Deal with a multiply by zero --- */
957
958 if (bv == bvl) {
c8a2f9ef 959 MPX_ZERO(dv, dvl);
d03ab969 960 return;
961 }
962
963 /* --- Do the initial multiply and initialize the accumulator --- */
964
965 MPX_UMULN(dv, dvl, av, avl, *bv++);
966
967 /* --- Do the remaining multiply/accumulates --- */
968
c8a2f9ef 969 while (dv < dvl && bv < bvl) {
d03ab969 970 mpw m = *bv++;
c8a2f9ef 971 mpw c = 0;
d03ab969 972 const mpw *avv = av;
973 mpw *dvv = ++dv;
974
975 while (avv < avl) {
976 mpd x;
977 if (dvv >= dvl)
978 goto next;
c8a2f9ef 979 x = (mpd)*dvv + (mpd)m * (mpd)*avv++ + c;
980 *dvv++ = MPW(x);
d03ab969 981 c = x >> MPW_BITS;
982 }
c8a2f9ef 983 MPX_UADDN(dvv, dvl, c);
d03ab969 984 next:;
985 }
986}
987
dd517851 988/* --- @mpx_umuln@ --- *
989 *
990 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
991 * @const mpw *av, *avl@ = multiplicand vector base and limit
992 * @mpw m@ = multiplier
993 *
994 * Returns: ---
995 *
996 * Use: Multiplies a multiprecision integer by a single-word value.
997 * The destination and source may be equal. The destination
998 * is completely cleared after use.
999 */
1000
1001void mpx_umuln(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
1002{
1003 MPX_UMULN(dv, dvl, av, avl, m);
1004}
1005
1006/* --- @mpx_umlan@ --- *
1007 *
1008 * Arguments: @mpw *dv, *dvl@ = destination/accumulator base and limit
1009 * @const mpw *av, *avl@ = multiplicand vector base and limit
1010 * @mpw m@ = multiplier
1011 *
1012 * Returns: ---
1013 *
1014 * Use: Multiplies a multiprecision integer by a single-word value
1015 * and adds the result to an accumulator.
1016 */
1017
1018void mpx_umlan(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
1019{
1020 MPX_UMLAN(dv, dvl, av, avl, m);
1021}
1022
c8a2f9ef 1023/* --- @mpx_usqr@ --- *
1024 *
1025 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
1026 * @const mpw *av, *av@ = source vector base and limit
1027 *
1028 * Returns: ---
1029 *
1030 * Use: Performs unsigned integer squaring. The result vector must
1031 * not overlap the source vector in any way.
1032 */
1033
1034void mpx_usqr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
1035{
1036 MPX_ZERO(dv, dvl);
1037
1038 /* --- Main loop --- */
1039
1040 while (av < avl) {
1041 const mpw *avv = av;
1042 mpw *dvv = dv;
1043 mpw a = *av;
1044 mpd c;
1045
1046 /* --- Stop if I've run out of destination --- */
1047
1048 if (dvv >= dvl)
1049 break;
1050
1051 /* --- Work out the square at this point in the proceedings --- */
1052
1053 {
c8a2f9ef 1054 mpd x = (mpd)a * (mpd)a + *dvv;
1055 *dvv++ = MPW(x);
1056 c = MPW(x >> MPW_BITS);
1057 }
1058
1059 /* --- Now fix up the rest of the vector upwards --- */
1060
1061 avv++;
1062 while (dvv < dvl && avv < avl) {
c8a2f9ef 1063 mpd x = (mpd)a * (mpd)*avv++;
1064 mpd y = ((x << 1) & MPW_MAX) + c + *dvv;
1065 c = (x >> (MPW_BITS - 1)) + (y >> MPW_BITS);
1066 *dvv++ = MPW(y);
1067 }
1068 while (dvv < dvl && c) {
1069 mpd x = c + *dvv;
1070 *dvv++ = MPW(x);
1071 c = x >> MPW_BITS;
1072 }
1073
1074 /* --- Get ready for the next round --- */
1075
1076 av++;
1077 dv += 2;
1078 }
1079}
1080
d03ab969 1081/* --- @mpx_udiv@ --- *
1082 *
1083 * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit
1084 * @mpw *rv, *rvl@ = dividend/remainder vector base and limit
1085 * @const mpw *dv, *dvl@ = divisor vector base and limit
c8a2f9ef 1086 * @mpw *sv, *svl@ = scratch workspace
d03ab969 1087 *
1088 * Returns: ---
1089 *
1090 * Use: Performs unsigned integer division. If the result overflows
1091 * the quotient vector, high-order bits are discarded. (Clearly
1092 * the remainder vector can't overflow.) The various vectors
1093 * may not overlap in any way. Yes, I know it's a bit odd
1094 * requiring the dividend to be in the result position but it
1095 * does make some sense really. The remainder must have
c8a2f9ef 1096 * headroom for at least two extra words. The scratch space
f45a00c6 1097 * must be at least one word larger than the divisor.
d03ab969 1098 */
1099
1100void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
c8a2f9ef 1101 const mpw *dv, const mpw *dvl,
1102 mpw *sv, mpw *svl)
d03ab969 1103{
d03ab969 1104 unsigned norm = 0;
1105 size_t scale;
1106 mpw d, dd;
1107
1108 /* --- Initialize the quotient --- */
1109
1110 MPX_ZERO(qv, qvl);
1111
c8a2f9ef 1112 /* --- Perform some sanity checks --- */
1113
1114 MPX_SHRINK(dv, dvl);
1115 assert(((void)"division by zero in mpx_udiv", dv < dvl));
1116
d03ab969 1117 /* --- Normalize the divisor --- *
1118 *
1119 * The algorithm requires that the divisor be at least two digits long.
1120 * This is easy to fix.
1121 */
1122
c8a2f9ef 1123 {
1124 unsigned b;
d03ab969 1125
c8a2f9ef 1126 d = dvl[-1];
1127 for (b = MPW_BITS / 2; b; b >>= 1) {
1128 if (d < (MPW_MAX >> b)) {
1129 d <<= b;
1130 norm += b;
1131 }
1132 }
1133 if (dv + 1 == dvl)
1134 norm += MPW_BITS;
d03ab969 1135 }
d03ab969 1136
1137 /* --- Normalize the dividend/remainder to match --- */
1138
c8a2f9ef 1139 if (norm) {
c8a2f9ef 1140 mpx_lsl(rv, rvl, rv, rvl, norm);
f45a00c6 1141 mpx_lsl(sv, svl, dv, dvl, norm);
c8a2f9ef 1142 dv = sv;
f45a00c6 1143 dvl = svl;
c8a2f9ef 1144 MPX_SHRINK(dv, dvl);
1145 }
1146
d03ab969 1147 MPX_SHRINK(rv, rvl);
c8a2f9ef 1148 d = dvl[-1];
1149 dd = dvl[-2];
d03ab969 1150
1151 /* --- Work out the relative scales --- */
1152
1153 {
1154 size_t rvn = rvl - rv;
c8a2f9ef 1155 size_t dvn = dvl - dv;
d03ab969 1156
1157 /* --- If the divisor is clearly larger, notice this --- */
1158
1159 if (dvn > rvn) {
1160 mpx_lsr(rv, rvl, rv, rvl, norm);
1161 return;
1162 }
1163
1164 scale = rvn - dvn;
1165 }
1166
1167 /* --- Calculate the most significant quotient digit --- *
1168 *
1169 * Because the divisor has its top bit set, this can only happen once. The
1170 * pointer arithmetic is a little contorted, to make sure that the
1171 * behaviour is defined.
1172 */
1173
1174 if (MPX_UCMP(rv + scale, rvl, >=, dv, dvl)) {
1175 mpx_usub(rv + scale, rvl, rv + scale, rvl, dv, dvl);
1176 if (qvl - qv > scale)
1177 qv[scale] = 1;
1178 }
1179
1180 /* --- Now for the main loop --- */
1181
1182 {
c8a2f9ef 1183 mpw *rvv = rvl - 2;
d03ab969 1184
1185 while (scale) {
c8a2f9ef 1186 mpw q;
1187 mpd rh;
d03ab969 1188
1189 /* --- Get an estimate for the next quotient digit --- */
1190
c8a2f9ef 1191 mpw r = rvv[1];
1192 mpw rr = rvv[0];
1193 mpw rrr = *--rvv;
1194
1195 scale--;
1196 rh = ((mpd)r << MPW_BITS) | rr;
d03ab969 1197 if (r == d)
1198 q = MPW_MAX;
c8a2f9ef 1199 else
1200 q = MPW(rh / d);
d03ab969 1201
1202 /* --- Refine the estimate --- */
1203
1204 {
1205 mpd yh = (mpd)d * q;
ce76ff16 1206 mpd yy = (mpd)dd * q;
1207 mpw yl;
c8a2f9ef 1208
ce76ff16 1209 if (yy > MPW_MAX)
1210 yh += yy >> MPW_BITS;
1211 yl = MPW(yy);
c8a2f9ef 1212
1213 while (yh > rh || (yh == rh && yl > rrr)) {
1214 q--;
1215 yh -= d;
ce76ff16 1216 if (yl < dd)
1217 yh--;
99b30c23 1218 yl = MPW(yl - dd);
c8a2f9ef 1219 }
1220 }
1221
1222 /* --- Remove a chunk from the dividend --- */
1223
1224 {
1225 mpw *svv;
1226 const mpw *dvv;
f45a00c6 1227 mpw mc = 0, sc = 0;
c8a2f9ef 1228
f45a00c6 1229 /* --- Calculate the size of the chunk --- *
1230 *
1231 * This does the whole job of calculating @r >> scale - qd@.
1232 */
c8a2f9ef 1233
f45a00c6 1234 for (svv = rv + scale, dvv = dv;
1235 dvv < dvl && svv < rvl;
1236 svv++, dvv++) {
1237 mpd x = (mpd)*dvv * (mpd)q + mc;
1238 mc = x >> MPW_BITS;
1239 x = (mpd)*svv - MPW(x) - sc;
c8a2f9ef 1240 *svv = MPW(x);
f45a00c6 1241 if (x >> MPW_BITS)
1242 sc = 1;
1243 else
1244 sc = 0;
1245 }
1246
1247 if (svv < rvl) {
1248 mpd x = (mpd)*svv - mc - sc;
1249 *svv++ = MPW(x);
1250 if (x >> MPW_BITS)
1251 sc = MPW_MAX;
1252 else
1253 sc = 0;
1254 while (svv < rvl)
1255 *svv++ = sc;
c8a2f9ef 1256 }
c8a2f9ef 1257
f45a00c6 1258 /* --- Fix if the quotient was too large --- *
c8a2f9ef 1259 *
f45a00c6 1260 * This doesn't seem to happen very often.
c8a2f9ef 1261 */
1262
c8a2f9ef 1263 if (rvl[-1] > MPW_MAX / 2) {
1264 mpx_uadd(rv + scale, rvl, rv + scale, rvl, dv, dvl);
1265 q--;
1266 }
1267 }
1268
1269 /* --- Done for another iteration --- */
1270
1271 if (qvl - qv > scale)
1272 qv[scale] = q;
1273 r = rr;
1274 rr = rrr;
1275 }
1276 }
1277
1278 /* --- Now fiddle with unnormalizing and things --- */
1279
1280 mpx_lsr(rv, rvl, rv, rvl, norm);
d03ab969 1281}
1282
698bd937 1283/* --- @mpx_udivn@ --- *
1284 *
1285 * Arguments: @mpw *qv, *qvl@ = storage for the quotient (may overlap
1286 * dividend)
1287 * @const mpw *rv, *rvl@ = dividend
1288 * @mpw d@ = single-precision divisor
1289 *
1290 * Returns: Remainder after divison.
1291 *
1292 * Use: Performs a single-precision division operation.
1293 */
1294
1295mpw mpx_udivn(mpw *qv, mpw *qvl, const mpw *rv, const mpw *rvl, mpw d)
1296{
1297 size_t i;
1298 size_t ql = qvl - qv;
1299 mpd r = 0;
1300
1301 i = rvl - rv;
1302 while (i > 0) {
1303 i--;
1304 r = (r << MPW_BITS) | rv[i];
1305 if (i < ql)
1306 qv[i] = r / d;
1307 r %= d;
1308 }
1309 return (MPW(r));
1310}
1311
42684bdb 1312/*----- Test rig ----------------------------------------------------------*/
1313
1314#ifdef TEST_RIG
1315
1316#include <mLib/alloc.h>
1317#include <mLib/dstr.h>
1318#include <mLib/quis.h>
1319#include <mLib/testrig.h>
1320
1321#include "mpscan.h"
1322
1323#define ALLOC(v, vl, sz) do { \
1324 size_t _sz = (sz); \
1325 mpw *_vv = xmalloc(MPWS(_sz)); \
1326 mpw *_vvl = _vv + _sz; \
1327 (v) = _vv; \
1328 (vl) = _vvl; \
1329} while (0)
1330
1331#define LOAD(v, vl, d) do { \
1332 const dstr *_d = (d); \
1333 mpw *_v, *_vl; \
1334 ALLOC(_v, _vl, MPW_RQ(_d->len)); \
1335 mpx_loadb(_v, _vl, _d->buf, _d->len); \
1336 (v) = _v; \
1337 (vl) = _vl; \
1338} while (0)
1339
1340#define MAX(x, y) ((x) > (y) ? (x) : (y))
1341
1342static void dumpbits(const char *msg, const void *pp, size_t sz)
1343{
1344 const octet *p = pp;
1345 fputs(msg, stderr);
1346 for (; sz; sz--)
1347 fprintf(stderr, " %02x", *p++);
1348 fputc('\n', stderr);
1349}
1350
1351static void dumpmp(const char *msg, const mpw *v, const mpw *vl)
1352{
1353 fputs(msg, stderr);
1354 MPX_SHRINK(v, vl);
1355 while (v < vl)
1356 fprintf(stderr, " %08lx", (unsigned long)*--vl);
1357 fputc('\n', stderr);
1358}
1359
1360static int chkscan(const mpw *v, const mpw *vl,
1361 const void *pp, size_t sz, int step)
1362{
1363 mpscan mps;
1364 const octet *p = pp;
1365 unsigned bit = 0;
1366 int ok = 1;
1367
1368 mpscan_initx(&mps, v, vl);
1369 while (sz) {
1370 unsigned x = *p;
1371 int i;
1372 p += step;
1373 for (i = 0; i < 8 && MPSCAN_STEP(&mps); i++) {
1374 if (MPSCAN_BIT(&mps) != (x & 1)) {
1375 fprintf(stderr,
1376 "\n*** error, step %i, bit %u, expected %u, found %u\n",
1377 step, bit, x & 1, MPSCAN_BIT(&mps));
1378 ok = 0;
1379 }
1380 x >>= 1;
1381 bit++;
1382 }
1383 sz--;
1384 }
1385
1386 return (ok);
1387}
1388
1389static int loadstore(dstr *v)
1390{
1391 dstr d = DSTR_INIT;
1392 size_t sz = MPW_RQ(v->len) * 2, diff;
1393 mpw *m, *ml;
1394 int ok = 1;
1395
1396 dstr_ensure(&d, v->len);
1397 m = xmalloc(MPWS(sz));
1398
1399 for (diff = 0; diff < sz; diff += 5) {
1400 size_t oct;
1401
1402 ml = m + sz - diff;
1403
1404 mpx_loadl(m, ml, v->buf, v->len);
1405 if (!chkscan(m, ml, v->buf, v->len, +1))
1406 ok = 0;
1407 MPX_OCTETS(oct, m, ml);
1408 mpx_storel(m, ml, d.buf, d.sz);
1409 if (memcmp(d.buf, v->buf, oct) != 0) {
1410 dumpbits("\n*** storel failed", d.buf, d.sz);
1411 ok = 0;
1412 }
1413
1414 mpx_loadb(m, ml, v->buf, v->len);
1415 if (!chkscan(m, ml, v->buf + v->len - 1, v->len, -1))
1416 ok = 0;
1417 MPX_OCTETS(oct, m, ml);
1418 mpx_storeb(m, ml, d.buf, d.sz);
1419 if (memcmp(d.buf + d.sz - oct, v->buf + v->len - oct, oct) != 0) {
1420 dumpbits("\n*** storeb failed", d.buf, d.sz);
1421 ok = 0;
1422 }
1423 }
1424
1425 if (!ok)
1426 dumpbits("input data", v->buf, v->len);
1427
1428 free(m);
1429 dstr_destroy(&d);
1430 return (ok);
1431}
1432
f09e814a 1433static int twocl(dstr *v)
1434{
1435 dstr d = DSTR_INIT;
1436 mpw *m, *ml;
1437 size_t sz;
1438 int ok = 1;
1439
1440 sz = v[0].len; if (v[1].len > sz) sz = v[1].len;
1441 dstr_ensure(&d, sz);
1442
1443 sz = MPW_RQ(sz);
1444 m = xmalloc(MPWS(sz));
1445 ml = m + sz;
1446
1447 mpx_loadl(m, ml, v[0].buf, v[0].len);
1448 mpx_storel2cn(m, ml, d.buf, v[1].len);
1449 if (memcmp(d.buf, v[1].buf, v[1].len)) {
1450 dumpbits("\n*** storel2cn failed", d.buf, v[1].len);
1451 ok = 0;
1452 }
1453
1454 mpx_loadl2cn(m, ml, v[1].buf, v[1].len);
1455 mpx_storel(m, ml, d.buf, v[0].len);
1456 if (memcmp(d.buf, v[0].buf, v[0].len)) {
1457 dumpbits("\n*** loadl2cn failed", d.buf, v[0].len);
1458 ok = 0;
1459 }
1460
1461 if (!ok) {
1462 dumpbits("pos", v[0].buf, v[0].len);
1463 dumpbits("neg", v[1].buf, v[1].len);
1464 }
1465
1466 free(m);
1467 dstr_destroy(&d);
1468
1469 return (ok);
1470}
1471
1472static int twocb(dstr *v)
1473{
1474 dstr d = DSTR_INIT;
1475 mpw *m, *ml;
1476 size_t sz;
1477 int ok = 1;
1478
1479 sz = v[0].len; if (v[1].len > sz) sz = v[1].len;
1480 dstr_ensure(&d, sz);
1481
1482 sz = MPW_RQ(sz);
1483 m = xmalloc(MPWS(sz));
1484 ml = m + sz;
1485
1486 mpx_loadb(m, ml, v[0].buf, v[0].len);
1487 mpx_storeb2cn(m, ml, d.buf, v[1].len);
1488 if (memcmp(d.buf, v[1].buf, v[1].len)) {
1489 dumpbits("\n*** storeb2cn failed", d.buf, v[1].len);
1490 ok = 0;
1491 }
1492
1493 mpx_loadb2cn(m, ml, v[1].buf, v[1].len);
1494 mpx_storeb(m, ml, d.buf, v[0].len);
1495 if (memcmp(d.buf, v[0].buf, v[0].len)) {
1496 dumpbits("\n*** loadb2cn failed", d.buf, v[0].len);
1497 ok = 0;
1498 }
1499
1500 if (!ok) {
1501 dumpbits("pos", v[0].buf, v[0].len);
1502 dumpbits("neg", v[1].buf, v[1].len);
1503 }
1504
1505 free(m);
1506 dstr_destroy(&d);
1507
1508 return (ok);
1509}
1510
42684bdb 1511static int lsl(dstr *v)
1512{
1513 mpw *a, *al;
1514 int n = *(int *)v[1].buf;
1515 mpw *c, *cl;
1516 mpw *d, *dl;
1517 int ok = 1;
1518
1519 LOAD(a, al, &v[0]);
1520 LOAD(c, cl, &v[2]);
1521 ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS);
1522
1523 mpx_lsl(d, dl, a, al, n);
1a05a8ef 1524 if (!mpx_ueq(d, dl, c, cl)) {
42684bdb 1525 fprintf(stderr, "\n*** lsl(%i) failed\n", n);
1526 dumpmp(" a", a, al);
1527 dumpmp("expected", c, cl);
1528 dumpmp(" result", d, dl);
1529 ok = 0;
1530 }
1531
1532 free(a); free(c); free(d);
1533 return (ok);
1534}
1535
81578196 1536static int lslc(dstr *v)
1537{
1538 mpw *a, *al;
1539 int n = *(int *)v[1].buf;
1540 mpw *c, *cl;
1541 mpw *d, *dl;
1542 int ok = 1;
1543
1544 LOAD(a, al, &v[0]);
1545 LOAD(c, cl, &v[2]);
1546 ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS);
1547
1548 mpx_lslc(d, dl, a, al, n);
1549 if (!mpx_ueq(d, dl, c, cl)) {
1550 fprintf(stderr, "\n*** lslc(%i) failed\n", n);
1551 dumpmp(" a", a, al);
1552 dumpmp("expected", c, cl);
1553 dumpmp(" result", d, dl);
1554 ok = 0;
1555 }
1556
1557 free(a); free(c); free(d);
1558 return (ok);
1559}
1560
42684bdb 1561static int lsr(dstr *v)
1562{
1563 mpw *a, *al;
1564 int n = *(int *)v[1].buf;
1565 mpw *c, *cl;
1566 mpw *d, *dl;
1567 int ok = 1;
1568
1569 LOAD(a, al, &v[0]);
1570 LOAD(c, cl, &v[2]);
1571 ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS + 1);
1572
1573 mpx_lsr(d, dl, a, al, n);
1a05a8ef 1574 if (!mpx_ueq(d, dl, c, cl)) {
42684bdb 1575 fprintf(stderr, "\n*** lsr(%i) failed\n", n);
1576 dumpmp(" a", a, al);
1577 dumpmp("expected", c, cl);
1578 dumpmp(" result", d, dl);
1579 ok = 0;
1580 }
1581
1582 free(a); free(c); free(d);
1583 return (ok);
1584}
1585
1586static int uadd(dstr *v)
1587{
1588 mpw *a, *al;
1589 mpw *b, *bl;
1590 mpw *c, *cl;
1591 mpw *d, *dl;
1592 int ok = 1;
1593
1594 LOAD(a, al, &v[0]);
1595 LOAD(b, bl, &v[1]);
1596 LOAD(c, cl, &v[2]);
1597 ALLOC(d, dl, MAX(al - a, bl - b) + 1);
1598
1599 mpx_uadd(d, dl, a, al, b, bl);
1a05a8ef 1600 if (!mpx_ueq(d, dl, c, cl)) {
42684bdb 1601 fprintf(stderr, "\n*** uadd failed\n");
1602 dumpmp(" a", a, al);
1603 dumpmp(" b", b, bl);
1604 dumpmp("expected", c, cl);
1605 dumpmp(" result", d, dl);
1606 ok = 0;
1607 }
1608
1609 free(a); free(b); free(c); free(d);
1610 return (ok);
1611}
1612
1613static int usub(dstr *v)
1614{
1615 mpw *a, *al;
1616 mpw *b, *bl;
1617 mpw *c, *cl;
1618 mpw *d, *dl;
1619 int ok = 1;
1620
1621 LOAD(a, al, &v[0]);
1622 LOAD(b, bl, &v[1]);
1623 LOAD(c, cl, &v[2]);
1624 ALLOC(d, dl, al - a);
1625
1626 mpx_usub(d, dl, a, al, b, bl);
1a05a8ef 1627 if (!mpx_ueq(d, dl, c, cl)) {
42684bdb 1628 fprintf(stderr, "\n*** usub failed\n");
1629 dumpmp(" a", a, al);
1630 dumpmp(" b", b, bl);
1631 dumpmp("expected", c, cl);
1632 dumpmp(" result", d, dl);
1633 ok = 0;
1634 }
1635
1636 free(a); free(b); free(c); free(d);
1637 return (ok);
1638}
1639
1640static int umul(dstr *v)
1641{
1642 mpw *a, *al;
1643 mpw *b, *bl;
1644 mpw *c, *cl;
1645 mpw *d, *dl;
1646 int ok = 1;
1647
1648 LOAD(a, al, &v[0]);
1649 LOAD(b, bl, &v[1]);
1650 LOAD(c, cl, &v[2]);
1651 ALLOC(d, dl, (al - a) + (bl - b));
1652
1653 mpx_umul(d, dl, a, al, b, bl);
1a05a8ef 1654 if (!mpx_ueq(d, dl, c, cl)) {
42684bdb 1655 fprintf(stderr, "\n*** umul failed\n");
1656 dumpmp(" a", a, al);
1657 dumpmp(" b", b, bl);
1658 dumpmp("expected", c, cl);
1659 dumpmp(" result", d, dl);
1660 ok = 0;
1661 }
1662
1663 free(a); free(b); free(c); free(d);
1664 return (ok);
1665}
1666
1667static int usqr(dstr *v)
1668{
1669 mpw *a, *al;
1670 mpw *c, *cl;
1671 mpw *d, *dl;
1672 int ok = 1;
1673
1674 LOAD(a, al, &v[0]);
1675 LOAD(c, cl, &v[1]);
1676 ALLOC(d, dl, 2 * (al - a));
1677
1678 mpx_usqr(d, dl, a, al);
1a05a8ef 1679 if (!mpx_ueq(d, dl, c, cl)) {
42684bdb 1680 fprintf(stderr, "\n*** usqr failed\n");
1681 dumpmp(" a", a, al);
1682 dumpmp("expected", c, cl);
1683 dumpmp(" result", d, dl);
1684 ok = 0;
1685 }
1686
1687 free(a); free(c); free(d);
1688 return (ok);
1689}
1690
1691static int udiv(dstr *v)
1692{
1693 mpw *a, *al;
1694 mpw *b, *bl;
1695 mpw *q, *ql;
1696 mpw *r, *rl;
1697 mpw *qq, *qql;
1698 mpw *s, *sl;
1699 int ok = 1;
1700
1701 ALLOC(a, al, MPW_RQ(v[0].len) + 2); mpx_loadb(a, al, v[0].buf, v[0].len);
1702 LOAD(b, bl, &v[1]);
1703 LOAD(q, ql, &v[2]);
1704 LOAD(r, rl, &v[3]);
1705 ALLOC(qq, qql, al - a);
1706 ALLOC(s, sl, (bl - b) + 1);
1707
1708 mpx_udiv(qq, qql, a, al, b, bl, s, sl);
1a05a8ef 1709 if (!mpx_ueq(qq, qql, q, ql) ||
1710 !mpx_ueq(a, al, r, rl)) {
42684bdb 1711 fprintf(stderr, "\n*** udiv failed\n");
1712 dumpmp(" divisor", b, bl);
1713 dumpmp("expect r", r, rl);
1714 dumpmp("result r", a, al);
1715 dumpmp("expect q", q, ql);
1716 dumpmp("result q", qq, qql);
1717 ok = 0;
1718 }
1719
1720 free(a); free(b); free(r); free(q); free(s); free(qq);
1721 return (ok);
1722}
1723
1724static test_chunk defs[] = {
1725 { "load-store", loadstore, { &type_hex, 0 } },
f09e814a 1726 { "2cl", twocl, { &type_hex, &type_hex, } },
1727 { "2cb", twocb, { &type_hex, &type_hex, } },
42684bdb 1728 { "lsl", lsl, { &type_hex, &type_int, &type_hex, 0 } },
81578196 1729 { "lslc", lslc, { &type_hex, &type_int, &type_hex, 0 } },
42684bdb 1730 { "lsr", lsr, { &type_hex, &type_int, &type_hex, 0 } },
1731 { "uadd", uadd, { &type_hex, &type_hex, &type_hex, 0 } },
1732 { "usub", usub, { &type_hex, &type_hex, &type_hex, 0 } },
1733 { "umul", umul, { &type_hex, &type_hex, &type_hex, 0 } },
1734 { "usqr", usqr, { &type_hex, &type_hex, 0 } },
1735 { "udiv", udiv, { &type_hex, &type_hex, &type_hex, &type_hex, 0 } },
1736 { 0, 0, { 0 } }
1737};
1738
1739int main(int argc, char *argv[])
1740{
1741 test_run(argc, argv, defs, SRCDIR"/tests/mpx");
1742 return (0);
1743}
1744
42684bdb 1745#endif
1746
d03ab969 1747/*----- That's all, folks -------------------------------------------------*/
Catacomb cryptographic library (downstream of http://git.distorted.org.uk/~mdw/catacomb/)