3 * Low-level multiprecision arithmetic
5 * (c) 1999 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of Catacomb.
12 * Catacomb is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU Library General Public License as
14 * published by the Free Software Foundation; either version 2 of the
15 * License, or (at your option) any later version.
17 * Catacomb is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Library General Public License for more details.
22 * You should have received a copy of the GNU Library General Public
23 * License along with Catacomb; if not, write to the Free
24 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
28 /*----- Header files ------------------------------------------------------*/
35 #include <mLib/bits.h>
36 #include <mLib/macros.h>
42 /*----- Loading and storing -----------------------------------------------*/
44 /* --- These are all variations on a theme --- *
46 * Essentially we want to feed bits into a shift register, @ibits@ bits at a
47 * time, and extract them @obits@ bits at a time whenever there are enough.
48 * Of course, @i@ and @o@ will, in general, be different sizes, and we don't
49 * necessarily know which is larger.
51 * During an operation, we have a shift register @w@ and a most-recent input
52 * @t@. Together, these hold @bits@ significant bits of input. We arrange
53 * that @bits < ibits + obits <= 2*MPW_BITS@, so we can get away with using
54 * an @mpw@ for both of these quantitities.
57 /* --- @MPX_GETBITS@ --- *
59 * Arguments: @ibits@ = width of input units, in bits
60 * @obits@ = width of output units, in bits
61 * @iavail@ = condition expression: is input data available?
62 * @getbits@ = function or macro: set argument to next input
64 * Use: Read an input unit into @t@ and update the necessary
67 * It is assumed on entry that @bits < obits@. On exit, we have
68 * @bits < ibits + obits@, and @t@ is live.
71 #define MPX_GETBITS(ibits, obits, iavail, getbits) do { \
72 if (!iavail) goto flush; \
73 if (bits >= ibits) w |= t << (bits - ibits); \
78 /* --- @MPX_PUTBITS@ --- *
80 * Arguments: @ibits@ = width of input units, in bits
81 * @obits@ = width of output units, in bits
82 * @oavail@ = condition expression: is output space available?
83 * @putbits@ = function or macro: write its argument to output
85 * Use: Emit an output unit, and update the necessary variables. If
86 * the output buffer is full, then force an immediate return.
88 * We assume that @bits < ibits + obits@, and that @t@ is only
89 * relevant if @bits >= ibits@. (The @MPX_GETBITS@ macro
90 * ensures that this is true.)
93 #define SHRW(w, b) ((b) < MPW_BITS ? (w) >> (b) : 0)
95 #define MPX_PUTBITS(ibits, obits, oavail, putbits) do { \
96 if (!oavail) return; \
100 w = SHRW(w, obits); \
102 putbits(w | (t << (bits - ibits))); \
104 if (bits >= ibits) w = SHRW(w, obits) | (t << (bits - ibits)); \
105 else w = SHRW(w, obits) | (t >> (ibits - bits)); \
110 /* --- @MPX_LOADSTORE@ --- *
112 * Arguments: @name@ = name of function to create, without @mpx_@ prefix
113 * @wconst@ = qualifiers for @mpw *@ arguments
114 * @oconst@ = qualifiers for octet pointers
115 * @decls@ = additional declarations needed
116 * @ibits@ = width of input units, in bits
117 * @iavail@ = condition expression: is input data available?
118 * @getbits@ = function or macro: set argument to next input
119 * @obits@ = width of output units, in bits
120 * @oavail@ = condition expression: is output space available?
121 * @putbits@ = function or macro: write its argument to output
122 * @clear@ = statements to clear remainder of output
124 * Use: Generates a function to convert between a sequence of
125 * multiprecision words and a vector of octets.
127 * The arguments @ibits@, @iavail@ and @getbits@ are passed on
128 * to @MPX_GETBITS@; similarly, @obits@, @oavail@, and @putbits@
129 * are passed on to @MPX_PUTBITS@.
131 * The following variables are in scope: @v@ and @vl are the
132 * current base and limit of the word vector; @p@ and @q@ are
133 * the base and limit of the octet vector; @w@ and @t@ form the
134 * shift register used during the conversion (see commentary
135 * above); and @bits@ tracks the number of live bits in the
139 #define MPX_LOADSTORE(name, wconst, oconst, decls, \
140 ibits, iavail, getbits, obits, oavail, putbits, \
143 void mpx_##name(wconst mpw *v, wconst mpw *vl, \
144 oconst void *pp, size_t sz) \
147 oconst octet *p = pp, *q = p + sz; \
152 while (bits < obits) MPX_GETBITS(ibits, obits, iavail, getbits); \
153 while (bits >= obits) MPX_PUTBITS(ibits, obits, oavail, putbits); \
157 while (bits > 0) MPX_PUTBITS(ibits, obits, oavail, putbits); \
163 /* --- Macros for @getbits@ and @putbits@ --- */
165 #define GETMPW(t) do { t = *v++; } while (0)
166 #define PUTMPW(x) do { *v++ = MPW(x); } while (0)
168 #define GETOCTETI(t) do { t = *p++; } while (0)
169 #define PUTOCTETD(x) do { *--q = U8(x); } while (0)
171 #define PUTOCTETI(x) do { *p++ = U8(x); } while (0)
172 #define GETOCTETD(t) do { t = *--q; } while (0)
174 /* --- Machinery for two's complement I/O --- */
179 #define GETMPW_2CN(t) do { \
180 t = MPW(~*v++ + c); \
184 #define PUTMPW_2CN(t) do { \
185 mpw _t = MPW(~(t) + c); \
190 #define FLUSHW_2CN do { \
191 if (c) MPX_ONE(v, vl); \
192 else MPX_ZERO(v, vl); \
195 #define FLUSHO_2CN do { \
196 memset(p, c ? 0xff : 0, q - p); \
199 /* --- @mpx_storel@ --- *
201 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
202 * @void *pp@ = pointer to octet array
203 * @size_t sz@ = size of octet array
207 * Use: Stores an MP in an octet array, least significant octet
208 * first. High-end octets are silently discarded if there
209 * isn't enough space for them.
212 MPX_LOADSTORE(storel, const, EMPTY, EMPTY,
213 MPW_BITS, (v < vl), GETMPW,
214 8, (p < q), PUTOCTETI,
215 { memset(p, 0, q - p); })
217 /* --- @mpx_loadl@ --- *
219 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
220 * @const void *pp@ = pointer to octet array
221 * @size_t sz@ = size of octet array
225 * Use: Loads an MP in an octet array, least significant octet
226 * first. High-end octets are ignored if there isn't enough
230 MPX_LOADSTORE(loadl, EMPTY, const, EMPTY,
231 8, (p < q), GETOCTETI,
232 MPW_BITS, (v < vl), PUTMPW,
233 { MPX_ZERO(v, vl); })
236 /* --- @mpx_storeb@ --- *
238 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
239 * @void *pp@ = pointer to octet array
240 * @size_t sz@ = size of octet array
244 * Use: Stores an MP in an octet array, most significant octet
245 * first. High-end octets are silently discarded if there
246 * isn't enough space for them.
249 MPX_LOADSTORE(storeb, const, EMPTY, EMPTY,
250 MPW_BITS, (v < vl), GETMPW,
251 8, (p < q), PUTOCTETD,
252 { memset(p, 0, q - p); })
254 /* --- @mpx_loadb@ --- *
256 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
257 * @const void *pp@ = pointer to octet array
258 * @size_t sz@ = size of octet array
262 * Use: Loads an MP in an octet array, most significant octet
263 * first. High-end octets are ignored if there isn't enough
267 MPX_LOADSTORE(loadb, EMPTY, const, EMPTY,
268 8, (p < q), GETOCTETD,
269 MPW_BITS, (v < vl), PUTMPW,
270 { MPX_ZERO(v, vl); })
272 /* --- @mpx_storel2cn@ --- *
274 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
275 * @void *pp@ = pointer to octet array
276 * @size_t sz@ = size of octet array
280 * Use: Stores a negative MP in an octet array, least significant
281 * octet first, as two's complement. High-end octets are
282 * silently discarded if there isn't enough space for them.
283 * This obviously makes the output bad.
286 MPX_LOADSTORE(storel2cn, const, EMPTY, DECL_2CN,
287 MPW_BITS, (v < vl), GETMPW_2CN,
288 8, (p < q), PUTOCTETI,
291 /* --- @mpx_loadl2cn@ --- *
293 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
294 * @const void *pp@ = pointer to octet array
295 * @size_t sz@ = size of octet array
299 * Use: Loads a negative MP in an octet array, least significant
300 * octet first, as two's complement. High-end octets are
301 * ignored if there isn't enough space for them. This probably
302 * means you made the wrong choice coming here.
305 MPX_LOADSTORE(loadl2cn, EMPTY, const, DECL_2CN,
306 8, (p < q), GETOCTETI,
307 MPW_BITS, (v < vl), PUTMPW_2CN,
310 /* --- @mpx_storeb2cn@ --- *
312 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
313 * @void *pp@ = pointer to octet array
314 * @size_t sz@ = size of octet array
318 * Use: Stores a negative MP in an octet array, most significant
319 * octet first, as two's complement. High-end octets are
320 * silently discarded if there isn't enough space for them,
321 * which probably isn't what you meant.
324 MPX_LOADSTORE(storeb2cn, const, EMPTY, DECL_2CN,
325 MPW_BITS, (v < vl), GETMPW_2CN,
326 8, (p < q), PUTOCTETD,
329 /* --- @mpx_loadb2cn@ --- *
331 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
332 * @const void *pp@ = pointer to octet array
333 * @size_t sz@ = size of octet array
337 * Use: Loads a negative MP in an octet array, most significant octet
338 * first as two's complement. High-end octets are ignored if
339 * there isn't enough space for them. This probably means you
340 * chose this function wrongly.
343 MPX_LOADSTORE(loadb2cn, EMPTY, const, DECL_2CN,
344 8, (p < q), GETOCTETD,
345 MPW_BITS, (v < vl), PUTMPW_2CN,
348 /*----- Logical shifting --------------------------------------------------*/
350 /* --- @MPX_SHIFT1@ --- *
352 * Arguments: @init@ = initial accumulator value
353 * @out@ = expression to store in each output word
354 * @next@ = expression for next accumulator value
356 * Use: Performs a single-position shift. The input is scanned
357 * right-to-left. In the expressions @out@ and @next@, the
358 * accumulator is available in @w@ and the current input word is
361 * This macro is intended to be used in the @shift1@ argument of
362 * @MPX_SHIFTOP@, and expects variables describing the operation
363 * to be set up accordingly.
366 #define MPX_SHIFT1(init, out, next) do { \
369 if (dv >= dvl) break; \
374 if (dv < dvl) { *dv++ = MPW(w); MPX_ZERO(dv, dvl); } \
377 /* --- @MPX_SHIFTW@ --- *
379 * Arguments: @max@ = the maximum shift (in words) which is nontrivial
380 * @clear@ = function (or macro) to clear low-order output words
381 * @copy@ = statement to copy words from input to output
383 * Use: Performs a shift by a whole number of words. If the shift
384 * amount is @max@ or more words, then the destination is
385 * @clear@ed entirely; otherwise, @copy@ is executed.
387 * This macro is intended to be used in the @shiftw@ argument of
388 * @MPX_SHIFTOP@, and expects variables describing the operation
389 * to be set up accordingly.
392 #define MPX_SHIFTW(max, clear, copy) do { \
393 if (nw >= (max)) clear(dv, dvl); \
397 /* --- @MPX_SHIFTOP@ --- *
399 * Arguments: @name@ = name of function to define (without `@mpx_@' prefix)
400 * @shift1@ = statement to shift by a single bit
401 * @shiftw@ = statement to shift by a whole number of words
402 * @shift@ = statement to perform a general shift
404 * Use: Emits a shift operation. The input is @av@..@avl@; the
405 * output is @dv@..@dvl@; and the shift amount (in bits) is
406 * @n@. In @shiftw@ and @shift@, @nw@ and @nb@ are set up such
407 * that @n = nw*MPW_BITS + nb@ and @nb < MPW_BITS@.
410 #define MPX_SHIFTOP(name, shift1, shiftw, shift) \
412 void mpx_##name(mpw *dv, mpw *dvl, \
413 const mpw *av, const mpw *avl, \
418 MPX_COPY(dv, dvl, av, avl); \
420 do shift1 while (0); \
422 size_t nw = n/MPW_BITS; \
423 unsigned nb = n%MPW_BITS; \
424 if (!nb) do shiftw while (0); \
425 else do shift while (0); \
429 /* --- @MPX_SHIFT_LEFT@ --- *
431 * Arguments: @name@ = name of function to define (without `@mpx_@' prefix)
432 * @init1@ = initializer for single-bit shift accumulator
433 * @clear@ = function (or macro) to clear low-order output words
434 * @flush@ = expression for low-order nontrivial output word
436 * Use: Emits a left-shift operation. This expands to a call on
437 * @MPX_SHIFTOP@, but implements the complicated @shift@
440 * The @init1@ argument is as for @MPX_SHIFT1@, and @clear@ is
441 * as for @MPX_SHIFTW@ (though is used elsewhere). In a general
442 * shift, @nw@ whole low-order output words are set using
443 * @clear@; high-order words are zeroed; and the remaining words
444 * set with a left-to-right pass across the input; at the end of
445 * the operation, the least significant output word above those
446 * @clear@ed is set using @flush@, which may use the accumulator
447 * @w@ = @av[0] << nb@.
450 #define MPX_SHIFT_LEFT(name, init1, clear, flush) \
451 MPX_SHIFTOP(name, { \
454 t >> (MPW_BITS - 1)); \
456 MPX_SHIFTW(dvl - dv, clear, { \
457 MPX_COPY(dv + nw, dvl, av, avl); \
458 clear(dv, dv + nw); \
461 size_t nr = MPW_BITS - nb; \
462 size_t dvn = dvl - dv; \
463 size_t avn = avl - av; \
471 if (dvn <= avn + nw) { \
472 avl = av + dvn - nw; \
475 size_t off = avn + nw + 1; \
476 MPX_ZERO(dv + off, dvl); \
483 *--dvl = MPW(w | (t >> nr)); \
487 *--dvl = MPW(flush); \
491 /* --- @mpx_lsl@ --- *
493 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
494 * @const mpw *av, *avl@ = source vector base and limit
495 * @size_t n@ = number of bit positions to shift by
499 * Use: Performs a logical shift left operation on an integer.
502 MPX_SHIFT_LEFT(lsl, 0, MPX_ZERO, w)
504 /* --- @mpx_lslc@ --- *
506 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
507 * @const mpw *av, *avl@ = source vector base and limit
508 * @size_t n@ = number of bit positions to shift by
512 * Use: Performs a logical shift left operation on an integer, only
513 * it fills in the bits with ones instead of zeroes.
516 MPX_SHIFT_LEFT(lslc, 1, MPX_ONE, w | (MPW_MAX >> nr))
518 /* --- @mpx_lsr@ --- *
520 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
521 * @const mpw *av, *avl@ = source vector base and limit
522 * @size_t n@ = number of bit positions to shift by
526 * Use: Performs a logical shift right operation on an integer.
530 MPX_SHIFT1(av < avl ? *av++ >> 1 : 0,
531 w | (t << (MPW_BITS - 1)),
534 MPX_SHIFTW(avl - av, MPX_ZERO,
535 { MPX_COPY(dv, dvl, av + nw, avl); });
537 size_t nr = MPW_BITS - nb;
541 w = av < avl ? *av++ : 0;
544 if (dv >= dvl) goto done;
546 *dv++ = MPW((w >> nb) | (t << nr));
550 *dv++ = MPW(w >> nb);
556 /*----- Bitwise operations ------------------------------------------------*/
558 /* --- @mpx_bitop@ --- *
560 * Arguments: @mpw *dv, *dvl@ = destination vector
561 * @const mpw *av, *avl@ = first source vector
562 * @const mpw *bv, *bvl@ = second source vector
566 * Use; Provides the dyadic boolean functions.
569 #define MPX_BITBINOP(string) \
571 void mpx_bit##string(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, \
572 const mpw *bv, const mpw *bvl) \
574 MPX_SHRINK(av, avl); \
575 MPX_SHRINK(bv, bvl); \
579 a = (av < avl) ? *av++ : 0; \
580 b = (bv < bvl) ? *bv++ : 0; \
581 *dv++ = B##string(a, b); \
582 IGNORE(a); IGNORE(b); \
586 MPX_DOBIN(MPX_BITBINOP)
588 void mpx_not(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
594 a = (av < avl) ? *av++ : 0;
599 /*----- Unsigned arithmetic -----------------------------------------------*/
601 /* --- @mpx_2c@ --- *
603 * Arguments: @mpw *dv, *dvl@ = destination vector
604 * @const mpw *v, *vl@ = source vector
608 * Use: Calculates the two's complement of @v@.
611 void mpx_2c(mpw *dv, mpw *dvl, const mpw *v, const mpw *vl)
614 while (dv < dvl && v < vl)
615 *dv++ = c = MPW(~*v++);
622 MPX_UADDN(dv, dvl, 1);
625 /* --- @mpx_ueq@ --- *
627 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
628 * @const mpw *bv, *bvl@ = second argument vector base and limit
630 * Returns: Nonzero if the two vectors are equal.
632 * Use: Performs an unsigned integer test for equality.
635 int mpx_ueq(const mpw *av, const mpw *avl, const mpw *bv, const mpw *bvl)
639 if (avl - av != bvl - bv)
648 /* --- @mpx_ucmp@ --- *
650 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
651 * @const mpw *bv, *bvl@ = second argument vector base and limit
653 * Returns: Less than, equal to, or greater than zero depending on
654 * whether @a@ is less than, equal to or greater than @b@,
657 * Use: Performs an unsigned integer comparison.
660 int mpx_ucmp(const mpw *av, const mpw *avl, const mpw *bv, const mpw *bvl)
665 if (avl - av > bvl - bv)
667 else if (avl - av < bvl - bv)
669 else while (avl > av) {
670 mpw a = *--avl, b = *--bvl;
679 /* --- @mpx_uadd@ --- *
681 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
682 * @const mpw *av, *avl@ = first addend vector base and limit
683 * @const mpw *bv, *bvl@ = second addend vector base and limit
687 * Use: Performs unsigned integer addition. If the result overflows
688 * the destination vector, high-order bits are discarded. This
689 * means that two's complement addition happens more or less for
690 * free, although that's more a side-effect than anything else.
691 * The result vector may be equal to either or both source
692 * vectors, but may not otherwise overlap them.
695 void mpx_uadd(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
696 const mpw *bv, const mpw *bvl)
700 while (av < avl || bv < bvl) {
705 a = (av < avl) ? *av++ : 0;
706 b = (bv < bvl) ? *bv++ : 0;
707 x = (mpd)a + (mpd)b + c;
717 /* --- @mpx_uaddn@ --- *
719 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
720 * @mpw n@ = other addend
724 * Use: Adds a small integer to a multiprecision number.
727 void mpx_uaddn(mpw *dv, mpw *dvl, mpw n) { MPX_UADDN(dv, dvl, n); }
729 /* --- @mpx_uaddnlsl@ --- *
731 * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
732 * @mpw a@ = second argument
733 * @unsigned o@ = offset in bits
737 * Use: Computes %$d + 2^o a$%. If the result overflows then
738 * high-order bits are discarded, as usual. We must have
739 * @0 < o < MPW_BITS@.
742 void mpx_uaddnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
746 while (x && dv < dvl) {
753 /* --- @mpx_usub@ --- *
755 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
756 * @const mpw *av, *avl@ = first argument vector base and limit
757 * @const mpw *bv, *bvl@ = second argument vector base and limit
761 * Use: Performs unsigned integer subtraction. If the result
762 * overflows the destination vector, high-order bits are
763 * discarded. This means that two's complement subtraction
764 * happens more or less for free, althuogh that's more a side-
765 * effect than anything else. The result vector may be equal to
766 * either or both source vectors, but may not otherwise overlap
770 void mpx_usub(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
771 const mpw *bv, const mpw *bvl)
775 while (av < avl || bv < bvl) {
780 a = (av < avl) ? *av++ : 0;
781 b = (bv < bvl) ? *bv++ : 0;
782 x = (mpd)a - (mpd)b - c;
795 /* --- @mpx_usubn@ --- *
797 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
802 * Use: Subtracts a small integer from a multiprecision number.
805 void mpx_usubn(mpw *dv, mpw *dvl, mpw n) { MPX_USUBN(dv, dvl, n); }
807 /* --- @mpx_uaddnlsl@ --- *
809 * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
810 * @mpw a@ = second argument
811 * @unsigned o@ = offset in bits
815 * Use: Computes %$d + 2^o a$%. If the result overflows then
816 * high-order bits are discarded, as usual. We must have
817 * @0 < o < MPW_BITS@.
820 void mpx_usubnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
822 mpw b = a >> (MPW_BITS - o);
826 mpd x = (mpd)*dv - MPW(a);
830 MPX_USUBN(dv, dvl, b);
834 /* --- @mpx_umul@ --- *
836 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
837 * @const mpw *av, *avl@ = multiplicand vector base and limit
838 * @const mpw *bv, *bvl@ = multiplier vector base and limit
842 * Use: Performs unsigned integer multiplication. If the result
843 * overflows the desination vector, high-order bits are
844 * discarded. The result vector may not overlap the argument
845 * vectors in any way.
848 void mpx_umul(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
849 const mpw *bv, const mpw *bvl)
851 /* --- This is probably worthwhile on a multiply --- */
856 /* --- Deal with a multiply by zero --- */
863 /* --- Do the initial multiply and initialize the accumulator --- */
865 MPX_UMULN(dv, dvl, av, avl, *bv++);
867 /* --- Do the remaining multiply/accumulates --- */
869 while (dv < dvl && bv < bvl) {
879 x = (mpd)*dvv + (mpd)m * (mpd)*avv++ + c;
883 MPX_UADDN(dvv, dvl, c);
888 /* --- @mpx_umuln@ --- *
890 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
891 * @const mpw *av, *avl@ = multiplicand vector base and limit
892 * @mpw m@ = multiplier
896 * Use: Multiplies a multiprecision integer by a single-word value.
897 * The destination and source may be equal. The destination
898 * is completely cleared after use.
901 void mpx_umuln(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
902 { MPX_UMULN(dv, dvl, av, avl, m); }
904 /* --- @mpx_umlan@ --- *
906 * Arguments: @mpw *dv, *dvl@ = destination/accumulator base and limit
907 * @const mpw *av, *avl@ = multiplicand vector base and limit
908 * @mpw m@ = multiplier
912 * Use: Multiplies a multiprecision integer by a single-word value
913 * and adds the result to an accumulator.
916 void mpx_umlan(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
917 { MPX_UMLAN(dv, dvl, av, avl, m); }
919 /* --- @mpx_usqr@ --- *
921 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
922 * @const mpw *av, *av@ = source vector base and limit
926 * Use: Performs unsigned integer squaring. The result vector must
927 * not overlap the source vector in any way.
930 void mpx_usqr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
934 /* --- Main loop --- */
942 /* --- Stop if I've run out of destination --- */
947 /* --- Work out the square at this point in the proceedings --- */
950 mpd x = (mpd)a * (mpd)a + *dvv;
952 c = MPW(x >> MPW_BITS);
955 /* --- Now fix up the rest of the vector upwards --- */
958 while (dvv < dvl && avv < avl) {
959 mpd x = (mpd)a * (mpd)*avv++;
960 mpd y = ((x << 1) & MPW_MAX) + c + *dvv;
961 c = (x >> (MPW_BITS - 1)) + (y >> MPW_BITS);
964 while (dvv < dvl && c) {
970 /* --- Get ready for the next round --- */
977 /* --- @mpx_udiv@ --- *
979 * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit
980 * @mpw *rv, *rvl@ = dividend/remainder vector base and limit
981 * @const mpw *dv, *dvl@ = divisor vector base and limit
982 * @mpw *sv, *svl@ = scratch workspace
986 * Use: Performs unsigned integer division. If the result overflows
987 * the quotient vector, high-order bits are discarded. (Clearly
988 * the remainder vector can't overflow.) The various vectors
989 * may not overlap in any way. Yes, I know it's a bit odd
990 * requiring the dividend to be in the result position but it
991 * does make some sense really. The remainder must have
992 * headroom for at least two extra words. The scratch space
993 * must be at least one word larger than the divisor.
996 void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
997 const mpw *dv, const mpw *dvl,
1004 /* --- Initialize the quotient --- */
1008 /* --- Perform some sanity checks --- */
1010 MPX_SHRINK(dv, dvl);
1011 assert(((void)"division by zero in mpx_udiv", dv < dvl));
1013 /* --- Normalize the divisor --- *
1015 * The algorithm requires that the divisor be at least two digits long.
1016 * This is easy to fix.
1023 for (b = MPW_P2; b; b >>= 1) {
1024 if (d <= (MPW_MAX >> b)) {
1033 /* --- Normalize the dividend/remainder to match --- */
1036 mpx_lsl(rv, rvl, rv, rvl, norm);
1037 mpx_lsl(sv, svl, dv, dvl, norm);
1040 MPX_SHRINK(dv, dvl);
1043 MPX_SHRINK(rv, rvl);
1047 /* --- Work out the relative scales --- */
1050 size_t rvn = rvl - rv;
1051 size_t dvn = dvl - dv;
1053 /* --- If the divisor is clearly larger, notice this --- */
1056 mpx_lsr(rv, rvl, rv, rvl, norm);
1063 /* --- Calculate the most significant quotient digit --- *
1065 * Because the divisor has its top bit set, this can only happen once. The
1066 * pointer arithmetic is a little contorted, to make sure that the
1067 * behaviour is defined.
1070 if (MPX_UCMP(rv + scale, rvl, >=, dv, dvl)) {
1071 mpx_usub(rv + scale, rvl, rv + scale, rvl, dv, dvl);
1072 if (qvl - qv > scale)
1076 /* --- Now for the main loop --- */
1085 /* --- Get an estimate for the next quotient digit --- */
1092 rh = ((mpd)r << MPW_BITS) | rr;
1098 /* --- Refine the estimate --- */
1101 mpd yh = (mpd)d * q;
1102 mpd yy = (mpd)dd * q;
1106 yh += yy >> MPW_BITS;
1109 while (yh > rh || (yh == rh && yl > rrr)) {
1118 /* --- Remove a chunk from the dividend --- */
1125 /* --- Calculate the size of the chunk --- *
1127 * This does the whole job of calculating @r >> scale - qd@.
1130 for (svv = rv + scale, dvv = dv;
1131 dvv < dvl && svv < rvl;
1133 mpd x = (mpd)*dvv * (mpd)q + mc;
1135 x = (mpd)*svv - MPW(x) - sc;
1144 mpd x = (mpd)*svv - mc - sc;
1154 /* --- Fix if the quotient was too large --- *
1156 * This doesn't seem to happen very often.
1159 if (rvl[-1] > MPW_MAX / 2) {
1160 mpx_uadd(rv + scale, rvl, rv + scale, rvl, dv, dvl);
1165 /* --- Done for another iteration --- */
1167 if (qvl - qv > scale)
1174 /* --- Now fiddle with unnormalizing and things --- */
1176 mpx_lsr(rv, rvl, rv, rvl, norm);
1179 /* --- @mpx_udivn@ --- *
1181 * Arguments: @mpw *qv, *qvl@ = storage for the quotient (may overlap
1183 * @const mpw *rv, *rvl@ = dividend
1184 * @mpw d@ = single-precision divisor
1186 * Returns: Remainder after divison.
1188 * Use: Performs a single-precision division operation.
1191 mpw mpx_udivn(mpw *qv, mpw *qvl, const mpw *rv, const mpw *rvl, mpw d)
1194 size_t ql = qvl - qv;
1200 r = (r << MPW_BITS) | rv[i];
1208 /*----- Test rig ----------------------------------------------------------*/
1212 #include <mLib/alloc.h>
1213 #include <mLib/dstr.h>
1214 #include <mLib/quis.h>
1215 #include <mLib/testrig.h>
1219 #define ALLOC(v, vl, sz) do { \
1220 size_t _sz = (sz); \
1221 mpw *_vv = xmalloc(MPWS(_sz)); \
1222 mpw *_vvl = _vv + _sz; \
1227 #define LOAD(v, vl, d) do { \
1228 const dstr *_d = (d); \
1230 ALLOC(_v, _vl, MPW_RQ(_d->len)); \
1231 mpx_loadb(_v, _vl, _d->buf, _d->len); \
1236 #define MAX(x, y) ((x) > (y) ? (x) : (y))
1238 static void dumpbits(const char *msg, const void *pp, size_t sz)
1240 const octet *p = pp;
1243 fprintf(stderr, " %02x", *p++);
1244 fputc('\n', stderr);
1247 static void dumpmp(const char *msg, const mpw *v, const mpw *vl)
1252 fprintf(stderr, " %08lx", (unsigned long)*--vl);
1253 fputc('\n', stderr);
1256 static int chkscan(const mpw *v, const mpw *vl,
1257 const void *pp, size_t sz, int step)
1260 const octet *p = pp;
1264 mpscan_initx(&mps, v, vl);
1269 for (i = 0; i < 8 && MPSCAN_STEP(&mps); i++) {
1270 if (MPSCAN_BIT(&mps) != (x & 1)) {
1272 "\n*** error, step %i, bit %u, expected %u, found %u\n",
1273 step, bit, x & 1, MPSCAN_BIT(&mps));
1285 static int loadstore(dstr *v)
1288 size_t sz = MPW_RQ(v->len) * 2, diff;
1292 dstr_ensure(&d, v->len);
1293 m = xmalloc(MPWS(sz));
1295 for (diff = 0; diff < sz; diff += 5) {
1300 mpx_loadl(m, ml, v->buf, v->len);
1301 if (!chkscan(m, ml, v->buf, v->len, +1))
1303 MPX_OCTETS(oct, m, ml);
1304 mpx_storel(m, ml, d.buf, d.sz);
1305 if (memcmp(d.buf, v->buf, oct) != 0) {
1306 dumpbits("\n*** storel failed", d.buf, d.sz);
1310 mpx_loadb(m, ml, v->buf, v->len);
1311 if (!chkscan(m, ml, v->buf + v->len - 1, v->len, -1))
1313 MPX_OCTETS(oct, m, ml);
1314 mpx_storeb(m, ml, d.buf, d.sz);
1315 if (memcmp(d.buf + d.sz - oct, v->buf + v->len - oct, oct) != 0) {
1316 dumpbits("\n*** storeb failed", d.buf, d.sz);
1322 dumpbits("input data", v->buf, v->len);
1329 static int twocl(dstr *v)
1336 sz = v[0].len; if (v[1].len > sz) sz = v[1].len;
1337 dstr_ensure(&d, sz);
1340 m = xmalloc(MPWS(sz));
1343 mpx_loadl(m, ml, v[0].buf, v[0].len);
1344 mpx_storel2cn(m, ml, d.buf, v[1].len);
1345 if (memcmp(d.buf, v[1].buf, v[1].len)) {
1346 dumpbits("\n*** storel2cn failed", d.buf, v[1].len);
1350 mpx_loadl2cn(m, ml, v[1].buf, v[1].len);
1351 mpx_storel(m, ml, d.buf, v[0].len);
1352 if (memcmp(d.buf, v[0].buf, v[0].len)) {
1353 dumpbits("\n*** loadl2cn failed", d.buf, v[0].len);
1358 dumpbits("pos", v[0].buf, v[0].len);
1359 dumpbits("neg", v[1].buf, v[1].len);
1368 static int twocb(dstr *v)
1375 sz = v[0].len; if (v[1].len > sz) sz = v[1].len;
1376 dstr_ensure(&d, sz);
1379 m = xmalloc(MPWS(sz));
1382 mpx_loadb(m, ml, v[0].buf, v[0].len);
1383 mpx_storeb2cn(m, ml, d.buf, v[1].len);
1384 if (memcmp(d.buf, v[1].buf, v[1].len)) {
1385 dumpbits("\n*** storeb2cn failed", d.buf, v[1].len);
1389 mpx_loadb2cn(m, ml, v[1].buf, v[1].len);
1390 mpx_storeb(m, ml, d.buf, v[0].len);
1391 if (memcmp(d.buf, v[0].buf, v[0].len)) {
1392 dumpbits("\n*** loadb2cn failed", d.buf, v[0].len);
1397 dumpbits("pos", v[0].buf, v[0].len);
1398 dumpbits("neg", v[1].buf, v[1].len);
1407 static int lsl(dstr *v)
1410 int n = *(int *)v[1].buf;
1417 ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS);
1419 mpx_lsl(d, dl, a, al, n);
1420 if (!mpx_ueq(d, dl, c, cl)) {
1421 fprintf(stderr, "\n*** lsl(%i) failed\n", n);
1422 dumpmp(" a", a, al);
1423 dumpmp("expected", c, cl);
1424 dumpmp(" result", d, dl);
1428 xfree(a); xfree(c); xfree(d);
1432 static int lslc(dstr *v)
1435 int n = *(int *)v[1].buf;
1442 ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS);
1444 mpx_lslc(d, dl, a, al, n);
1445 if (!mpx_ueq(d, dl, c, cl)) {
1446 fprintf(stderr, "\n*** lslc(%i) failed\n", n);
1447 dumpmp(" a", a, al);
1448 dumpmp("expected", c, cl);
1449 dumpmp(" result", d, dl);
1453 xfree(a); xfree(c); xfree(d);
1457 static int lsr(dstr *v)
1460 int n = *(int *)v[1].buf;
1467 ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS + 1);
1469 mpx_lsr(d, dl, a, al, n);
1470 if (!mpx_ueq(d, dl, c, cl)) {
1471 fprintf(stderr, "\n*** lsr(%i) failed\n", n);
1472 dumpmp(" a", a, al);
1473 dumpmp("expected", c, cl);
1474 dumpmp(" result", d, dl);
1478 xfree(a); xfree(c); xfree(d);
1482 static int uadd(dstr *v)
1493 ALLOC(d, dl, MAX(al - a, bl - b) + 1);
1495 mpx_uadd(d, dl, a, al, b, bl);
1496 if (!mpx_ueq(d, dl, c, cl)) {
1497 fprintf(stderr, "\n*** uadd failed\n");
1498 dumpmp(" a", a, al);
1499 dumpmp(" b", b, bl);
1500 dumpmp("expected", c, cl);
1501 dumpmp(" result", d, dl);
1505 xfree(a); xfree(b); xfree(c); xfree(d);
1509 static int usub(dstr *v)
1520 ALLOC(d, dl, al - a);
1522 mpx_usub(d, dl, a, al, b, bl);
1523 if (!mpx_ueq(d, dl, c, cl)) {
1524 fprintf(stderr, "\n*** usub failed\n");
1525 dumpmp(" a", a, al);
1526 dumpmp(" b", b, bl);
1527 dumpmp("expected", c, cl);
1528 dumpmp(" result", d, dl);
1532 xfree(a); xfree(b); xfree(c); xfree(d);
1536 static int umul(dstr *v)
1547 ALLOC(d, dl, (al - a) + (bl - b));
1549 mpx_umul(d, dl, a, al, b, bl);
1550 if (!mpx_ueq(d, dl, c, cl)) {
1551 fprintf(stderr, "\n*** umul failed\n");
1552 dumpmp(" a", a, al);
1553 dumpmp(" b", b, bl);
1554 dumpmp("expected", c, cl);
1555 dumpmp(" result", d, dl);
1559 xfree(a); xfree(b); xfree(c); xfree(d);
1563 static int usqr(dstr *v)
1572 ALLOC(d, dl, 2 * (al - a));
1574 mpx_usqr(d, dl, a, al);
1575 if (!mpx_ueq(d, dl, c, cl)) {
1576 fprintf(stderr, "\n*** usqr failed\n");
1577 dumpmp(" a", a, al);
1578 dumpmp("expected", c, cl);
1579 dumpmp(" result", d, dl);
1583 xfree(a); xfree(c); xfree(d);
1587 static int udiv(dstr *v)
1597 ALLOC(a, al, MPW_RQ(v[0].len) + 2); mpx_loadb(a, al, v[0].buf, v[0].len);
1601 ALLOC(qq, qql, al - a);
1602 ALLOC(s, sl, (bl - b) + 1);
1604 mpx_udiv(qq, qql, a, al, b, bl, s, sl);
1605 if (!mpx_ueq(qq, qql, q, ql) ||
1606 !mpx_ueq(a, al, r, rl)) {
1607 fprintf(stderr, "\n*** udiv failed\n");
1608 dumpmp(" divisor", b, bl);
1609 dumpmp("expect r", r, rl);
1610 dumpmp("result r", a, al);
1611 dumpmp("expect q", q, ql);
1612 dumpmp("result q", qq, qql);
1616 xfree(a); xfree(b); xfree(r); xfree(q); xfree(s); xfree(qq);
1620 static test_chunk defs[] = {
1621 { "load-store", loadstore, { &type_hex, 0 } },
1622 { "2cl", twocl, { &type_hex, &type_hex, } },
1623 { "2cb", twocb, { &type_hex, &type_hex, } },
1624 { "lsl", lsl, { &type_hex, &type_int, &type_hex, 0 } },
1625 { "lslc", lslc, { &type_hex, &type_int, &type_hex, 0 } },
1626 { "lsr", lsr, { &type_hex, &type_int, &type_hex, 0 } },
1627 { "uadd", uadd, { &type_hex, &type_hex, &type_hex, 0 } },
1628 { "usub", usub, { &type_hex, &type_hex, &type_hex, 0 } },
1629 { "umul", umul, { &type_hex, &type_hex, &type_hex, 0 } },
1630 { "usqr", usqr, { &type_hex, &type_hex, 0 } },
1631 { "udiv", udiv, { &type_hex, &type_hex, &type_hex, &type_hex, 0 } },
1635 int main(int argc, char *argv[])
1637 test_run(argc, argv, defs, SRCDIR"/t/mpx");
1643 /*----- That's all, folks -------------------------------------------------*/