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 ------------------------------------------------------*/
37 #include <mLib/bits.h>
38 #include <mLib/macros.h>
45 /*----- Loading and storing -----------------------------------------------*/
47 /* --- These are all variations on a theme --- *
49 * Essentially we want to feed bits into a shift register, @ibits@ bits at a
50 * time, and extract them @obits@ bits at a time whenever there are enough.
51 * Of course, @i@ and @o@ will, in general, be different sizes, and we don't
52 * necessarily know which is larger.
54 * During an operation, we have a shift register @w@ and a most-recent input
55 * @t@. Together, these hold @bits@ significant bits of input. We arrange
56 * that @bits < ibits + obits <= 2*MPW_BITS@, so we can get away with using
57 * an @mpw@ for both of these quantitities.
60 /* --- @MPX_GETBITS@ --- *
62 * Arguments: @ibits@ = width of input units, in bits
63 * @obits@ = width of output units, in bits
64 * @iavail@ = condition expression: is input data available?
65 * @getbits@ = function or macro: set argument to next input
67 * Use: Read an input unit into @t@ and update the necessary
70 * It is assumed on entry that @bits < obits@. On exit, we have
71 * @bits < ibits + obits@, and @t@ is live.
74 #define MPX_GETBITS(ibits, obits, iavail, getbits) do { \
75 if (!iavail) goto flush; \
76 if (bits >= ibits) w |= t << (bits - ibits); \
81 /* --- @MPX_PUTBITS@ --- *
83 * Arguments: @ibits@ = width of input units, in bits
84 * @obits@ = width of output units, in bits
85 * @oavail@ = condition expression: is output space available?
86 * @putbits@ = function or macro: write its argument to output
88 * Use: Emit an output unit, and update the necessary variables. If
89 * the output buffer is full, then force an immediate return.
91 * We assume that @bits < ibits + obits@, and that @t@ is only
92 * relevant if @bits >= ibits@. (The @MPX_GETBITS@ macro
93 * ensures that this is true.)
96 #define SHRW(w, b) ((b) < MPW_BITS ? (w) >> (b) : 0)
98 #define MPX_PUTBITS(ibits, obits, oavail, putbits) do { \
99 if (!oavail) return; \
100 if (bits < ibits) { \
103 w = SHRW(w, obits); \
105 putbits(w | (t << (bits - ibits))); \
107 if (bits >= ibits) w = SHRW(w, obits) | (t << (bits - ibits)); \
108 else w = SHRW(w, obits) | (t >> (ibits - bits)); \
113 /* --- @MPX_LOADSTORE@ --- *
115 * Arguments: @name@ = name of function to create, without @mpx_@ prefix
116 * @wconst@ = qualifiers for @mpw *@ arguments
117 * @oconst@ = qualifiers for octet pointers
118 * @decls@ = additional declarations needed
119 * @ibits@ = width of input units, in bits
120 * @iavail@ = condition expression: is input data available?
121 * @getbits@ = function or macro: set argument to next input
122 * @obits@ = width of output units, in bits
123 * @oavail@ = condition expression: is output space available?
124 * @putbits@ = function or macro: write its argument to output
125 * @fixfinal@ = statements to fix shift register at the end
126 * @clear@ = statements to clear remainder of output
128 * Use: Generates a function to convert between a sequence of
129 * multiprecision words and a vector of octets.
131 * The arguments @ibits@, @iavail@ and @getbits@ are passed on
132 * to @MPX_GETBITS@; similarly, @obits@, @oavail@, and @putbits@
133 * are passed on to @MPX_PUTBITS@.
135 * The following variables are in scope: @v@ and @vl are the
136 * current base and limit of the word vector; @p@ and @q@ are
137 * the base and limit of the octet vector; @w@ and @t@ form the
138 * shift register used during the conversion (see commentary
139 * above); and @bits@ tracks the number of live bits in the
143 #define MPX_LOADSTORE(name, wconst, oconst, decls, \
144 ibits, iavail, getbits, obits, oavail, putbits, \
147 void mpx_##name(wconst mpw *v, wconst mpw *vl, \
148 oconst void *pp, size_t sz) \
151 oconst octet *p = pp, *q = p + sz; \
156 while (bits < obits) MPX_GETBITS(ibits, obits, iavail, getbits); \
157 while (bits >= obits) MPX_PUTBITS(ibits, obits, oavail, putbits); \
163 while (bits > 0) MPX_PUTBITS(ibits, obits, oavail, putbits); \
170 /* --- Macros for @getbits@ and @putbits@ --- */
172 #define GETMPW(t) do { t = *v++; } while (0)
173 #define PUTMPW(x) do { *v++ = MPW(x); } while (0)
175 #define GETOCTETI(t) do { t = *p++; } while (0)
176 #define PUTOCTETD(x) do { *--q = U8(x); } while (0)
178 #define PUTOCTETI(x) do { *p++ = U8(x); } while (0)
179 #define GETOCTETD(t) do { t = *--q; } while (0)
181 /* --- Machinery for two's complement I/O --- */
186 #define GETMPW_2CN(t) do { \
187 t = MPW(~*v++ + c); \
191 #define PUTMPW_2CN(t) do { \
192 mpw _t = MPW(~(t) + c); \
197 #define FIXFINALW_2CN do { \
198 if (c && !w && !t); \
199 else if (bits == 8) t ^= ~(mpw)0xffu; \
200 else t ^= ((mpw)1 << (MPW_BITS - bits + 8)) - 256u; \
203 #define FLUSHO_2CN do { \
204 memset(p, c ? 0 : 0xff, q - p); \
207 /* --- @mpx_storel@ --- *
209 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
210 * @void *pp@ = pointer to octet array
211 * @size_t sz@ = size of octet array
215 * Use: Stores an MP in an octet array, least significant octet
216 * first. High-end octets are silently discarded if there
217 * isn't enough space for them.
220 MPX_LOADSTORE(storel
, const, EMPTY
, EMPTY
,
221 MPW_BITS
, (v
< vl
), GETMPW
,
222 8, (p
< q
), PUTOCTETI
,
223 EMPTY
, { memset(p
, 0, q
- p
); })
225 /* --- @mpx_loadl@ --- *
227 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
228 * @const void *pp@ = pointer to octet array
229 * @size_t sz@ = size of octet array
233 * Use: Loads an MP in an octet array, least significant octet
234 * first. High-end octets are ignored if there isn't enough
238 MPX_LOADSTORE(loadl
, EMPTY
, const, EMPTY
,
239 8, (p
< q
), GETOCTETI
,
240 MPW_BITS
, (v
< vl
), PUTMPW
,
241 EMPTY
, { MPX_ZERO(v
, vl
); })
244 /* --- @mpx_storeb@ --- *
246 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
247 * @void *pp@ = pointer to octet array
248 * @size_t sz@ = size of octet array
252 * Use: Stores an MP in an octet array, most significant octet
253 * first. High-end octets are silently discarded if there
254 * isn't enough space for them.
257 MPX_LOADSTORE(storeb
, const, EMPTY
, EMPTY
,
258 MPW_BITS
, (v
< vl
), GETMPW
,
259 8, (p
< q
), PUTOCTETD
,
260 EMPTY
, { memset(p
, 0, q
- p
); })
262 /* --- @mpx_loadb@ --- *
264 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
265 * @const void *pp@ = pointer to octet array
266 * @size_t sz@ = size of octet array
270 * Use: Loads an MP in an octet array, most significant octet
271 * first. High-end octets are ignored if there isn't enough
275 MPX_LOADSTORE(loadb
, EMPTY
, const, EMPTY
,
276 8, (p
< q
), GETOCTETD
,
277 MPW_BITS
, (v
< vl
), PUTMPW
,
278 EMPTY
, { MPX_ZERO(v
, vl
); })
280 /* --- @mpx_storel2cn@ --- *
282 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
283 * @void *pp@ = pointer to octet array
284 * @size_t sz@ = size of octet array
288 * Use: Stores a negative MP in an octet array, least significant
289 * octet first, as two's complement. High-end octets are
290 * silently discarded if there isn't enough space for them.
291 * This obviously makes the output bad.
294 MPX_LOADSTORE(storel2cn
, const, EMPTY
, DECL_2CN
,
295 MPW_BITS
, (v
< vl
), GETMPW_2CN
,
296 8, (p
< q
), PUTOCTETI
,
297 EMPTY
, { FLUSHO_2CN
; })
299 /* --- @mpx_loadl2cn@ --- *
301 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
302 * @const void *pp@ = pointer to octet array
303 * @size_t sz@ = size of octet array
307 * Use: Loads a negative MP in an octet array, least significant
308 * octet first, as two's complement. High-end octets are
309 * ignored if there isn't enough space for them. This probably
310 * means you made the wrong choice coming here.
313 MPX_LOADSTORE(loadl2cn
, EMPTY
, const, DECL_2CN
,
314 8, (p
< q
), GETOCTETI
,
315 MPW_BITS
, (v
< vl
), PUTMPW_2CN
,
316 { FIXFINALW_2CN
; }, { MPX_ZERO(v
, vl
); })
318 /* --- @mpx_storeb2cn@ --- *
320 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
321 * @void *pp@ = pointer to octet array
322 * @size_t sz@ = size of octet array
326 * Use: Stores a negative MP in an octet array, most significant
327 * octet first, as two's complement. High-end octets are
328 * silently discarded if there isn't enough space for them,
329 * which probably isn't what you meant.
332 MPX_LOADSTORE(storeb2cn
, const, EMPTY
, DECL_2CN
,
333 MPW_BITS
, (v
< vl
), GETMPW_2CN
,
334 8, (p
< q
), PUTOCTETD
,
335 EMPTY
, { FLUSHO_2CN
; })
337 /* --- @mpx_loadb2cn@ --- *
339 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
340 * @const void *pp@ = pointer to octet array
341 * @size_t sz@ = size of octet array
345 * Use: Loads a negative MP in an octet array, most significant octet
346 * first as two's complement. High-end octets are ignored if
347 * there isn't enough space for them. This probably means you
348 * chose this function wrongly.
351 MPX_LOADSTORE(loadb2cn
, EMPTY
, const, DECL_2CN
,
352 8, (p
< q
), GETOCTETD
,
353 MPW_BITS
, (v
< vl
), PUTMPW_2CN
,
354 { FIXFINALW_2CN
; }, { MPX_ZERO(v
, vl
); })
356 /*----- Logical shifting --------------------------------------------------*/
358 /* --- @MPX_SHIFT1@ --- *
360 * Arguments: @init@ = initial accumulator value
361 * @out@ = expression to store in each output word
362 * @next@ = expression for next accumulator value
364 * Use: Performs a single-position shift. The input is scanned
365 * right-to-left. In the expressions @out@ and @next@, the
366 * accumulator is available in @w@ and the current input word is
369 * This macro is intended to be used in the @shift1@ argument of
370 * @MPX_SHIFTOP@, and expects variables describing the operation
371 * to be set up accordingly.
374 #define MPX_SHIFT1(init, out, next) do { \
377 if (dv >= dvl) break; \
382 if (dv < dvl) { *dv++ = MPW(w); MPX_ZERO(dv, dvl); } \
385 /* --- @MPX_SHIFTW@ --- *
387 * Arguments: @max@ = the maximum shift (in words) which is nontrivial
388 * @clear@ = function (or macro) to clear low-order output words
389 * @copy@ = statement to copy words from input to output
391 * Use: Performs a shift by a whole number of words. If the shift
392 * amount is @max@ or more words, then the destination is
393 * @clear@ed entirely; otherwise, @copy@ is executed.
395 * This macro is intended to be used in the @shiftw@ argument of
396 * @MPX_SHIFTOP@, and expects variables describing the operation
397 * to be set up accordingly.
400 #define MPX_SHIFTW(max, clear, copy) do { \
401 if (nw >= (max)) clear(dv, dvl); \
405 /* --- @MPX_SHIFTOP@ --- *
407 * Arguments: @name@ = name of function to define (without `@mpx_@' prefix)
408 * @shift1@ = statement to shift by a single bit
409 * @shiftw@ = statement to shift by a whole number of words
410 * @shift@ = statement to perform a general shift
412 * Use: Emits a shift operation. The input is @av@..@avl@; the
413 * output is @dv@..@dvl@; and the shift amount (in bits) is
414 * @n@. In @shiftw@ and @shift@, @nw@ and @nb@ are set up such
415 * that @n = nw*MPW_BITS + nb@ and @nb < MPW_BITS@.
418 #define MPX_SHIFTOP(name, shift1, shiftw, shift) \
420 void mpx_##name(mpw *dv, mpw *dvl, \
421 const mpw *av, const mpw *avl, \
426 MPX_COPY(dv, dvl, av, avl); \
428 do shift1 while (0); \
430 size_t nw = n/MPW_BITS; \
431 unsigned nb = n%MPW_BITS; \
432 if (!nb) do shiftw while (0); \
433 else do shift while (0); \
437 /* --- @MPX_SHIFT_LEFT@ --- *
439 * Arguments: @name@ = name of function to define (without `@mpx_@' prefix)
440 * @init1@ = initializer for single-bit shift accumulator
441 * @clear@ = function (or macro) to clear low-order output words
442 * @flush@ = expression for low-order nontrivial output word
444 * Use: Emits a left-shift operation. This expands to a call on
445 * @MPX_SHIFTOP@, but implements the complicated @shift@
448 * The @init1@ argument is as for @MPX_SHIFT1@, and @clear@ is
449 * as for @MPX_SHIFTW@ (though is used elsewhere). In a general
450 * shift, @nw@ whole low-order output words are set using
451 * @clear@; high-order words are zeroed; and the remaining words
452 * set with a left-to-right pass across the input; at the end of
453 * the operation, the least significant output word above those
454 * @clear@ed is set using @flush@, which may use the accumulator
455 * @w@ = @av[0] << nb@.
458 #define MPX_SHIFT_LEFT(name, init1, clear, flush) \
459 MPX_SHIFTOP(name, { \
462 t >> (MPW_BITS - 1)); \
464 MPX_SHIFTW(dvl - dv, clear, { \
465 MPX_COPY(dv + nw, dvl, av, avl); \
466 clear(dv, dv + nw); \
469 size_t nr = MPW_BITS - nb; \
470 size_t dvn = dvl - dv; \
471 size_t avn = avl - av; \
479 if (dvn <= avn + nw) { \
480 avl = av + dvn - nw; \
483 size_t off = avn + nw + 1; \
484 MPX_ZERO(dv + off, dvl); \
491 *--dvl = MPW(w | (t >> nr)); \
495 *--dvl = MPW(flush); \
499 /* --- @mpx_lsl@ --- *
501 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
502 * @const mpw *av, *avl@ = source vector base and limit
503 * @size_t n@ = number of bit positions to shift by
507 * Use: Performs a logical shift left operation on an integer.
510 MPX_SHIFT_LEFT(lsl
, 0, MPX_ZERO
, w
)
512 /* --- @mpx_lslc@ --- *
514 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
515 * @const mpw *av, *avl@ = source vector base and limit
516 * @size_t n@ = number of bit positions to shift by
520 * Use: Performs a logical shift left operation on an integer, only
521 * it fills in the bits with ones instead of zeroes.
524 MPX_SHIFT_LEFT(lslc
, 1, MPX_ONE
, w
| (MPW_MAX
>> nr
))
526 /* --- @mpx_lsr@ --- *
528 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
529 * @const mpw *av, *avl@ = source vector base and limit
530 * @size_t n@ = number of bit positions to shift by
534 * Use: Performs a logical shift right operation on an integer.
538 MPX_SHIFT1(av
< avl ?
*av
++ >> 1 : 0,
539 w
| (t
<< (MPW_BITS
- 1)),
542 MPX_SHIFTW(avl
- av
, MPX_ZERO
,
543 { MPX_COPY(dv
, dvl
, av
+ nw
, avl
); });
545 size_t nr
= MPW_BITS
- nb
;
556 if (dv
>= dvl
) goto done
;
558 *dv
++ = MPW((w
>> nb
) | (t
<< nr
));
564 *dv
++ = MPW(w
>> nb
);
570 /*----- Bitwise operations ------------------------------------------------*/
572 /* --- @mpx_bitop@ --- *
574 * Arguments: @mpw *dv, *dvl@ = destination vector
575 * @const mpw *av, *avl@ = first source vector
576 * @const mpw *bv, *bvl@ = second source vector
580 * Use; Provides the dyadic boolean functions.
583 #define MPX_BITBINOP(string) \
585 void mpx_bit##string(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, \
586 const mpw *bv, const mpw *bvl) \
588 MPX_SHRINK(av, avl); \
589 MPX_SHRINK(bv, bvl); \
593 a = (av < avl) ? *av++ : 0; \
594 b = (bv < bvl) ? *bv++ : 0; \
595 *dv++ = B##string(a, b); \
596 IGNORE(a); IGNORE(b); \
600 MPX_DOBIN(MPX_BITBINOP
)
602 void mpx_not(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
)
608 a
= (av
< avl
) ?
*av
++ : 0;
613 /*----- Unsigned arithmetic -----------------------------------------------*/
615 /* --- @mpx_2c@ --- *
617 * Arguments: @mpw *dv, *dvl@ = destination vector
618 * @const mpw *v, *vl@ = source vector
622 * Use: Calculates the two's complement of @v@.
625 void mpx_2c(mpw
*dv
, mpw
*dvl
, const mpw
*v
, const mpw
*vl
)
628 while (dv
< dvl
&& v
< vl
)
629 *dv
++ = c
= MPW(~*v
++);
636 MPX_UADDN(dv
, dvl
, 1);
639 /* --- @mpx_ueq@ --- *
641 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
642 * @const mpw *bv, *bvl@ = second argument vector base and limit
644 * Returns: Nonzero if the two vectors are equal.
646 * Use: Performs an unsigned integer test for equality.
649 int mpx_ueq(const mpw
*av
, const mpw
*avl
, const mpw
*bv
, const mpw
*bvl
)
653 if (avl
- av
!= bvl
- bv
)
662 /* --- @mpx_ucmp@ --- *
664 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
665 * @const mpw *bv, *bvl@ = second argument vector base and limit
667 * Returns: Less than, equal to, or greater than zero depending on
668 * whether @a@ is less than, equal to or greater than @b@,
671 * Use: Performs an unsigned integer comparison.
674 int mpx_ucmp(const mpw
*av
, const mpw
*avl
, const mpw
*bv
, const mpw
*bvl
)
679 if (avl
- av
> bvl
- bv
)
681 else if (avl
- av
< bvl
- bv
)
683 else while (avl
> av
) {
684 mpw a
= *--avl
, b
= *--bvl
;
693 /* --- @mpx_uadd@ --- *
695 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
696 * @const mpw *av, *avl@ = first addend vector base and limit
697 * @const mpw *bv, *bvl@ = second addend vector base and limit
701 * Use: Performs unsigned integer addition. If the result overflows
702 * the destination vector, high-order bits are discarded. This
703 * means that two's complement addition happens more or less for
704 * free, although that's more a side-effect than anything else.
705 * The result vector may be equal to either or both source
706 * vectors, but may not otherwise overlap them.
709 void mpx_uadd(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
,
710 const mpw
*bv
, const mpw
*bvl
)
714 while (av
< avl
|| bv
< bvl
) {
719 a
= (av
< avl
) ?
*av
++ : 0;
720 b
= (bv
< bvl
) ?
*bv
++ : 0;
721 x
= (mpd
)a
+ (mpd
)b
+ c
;
731 /* --- @mpx_uaddn@ --- *
733 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
734 * @mpw n@ = other addend
738 * Use: Adds a small integer to a multiprecision number.
741 void mpx_uaddn(mpw
*dv
, mpw
*dvl
, mpw n
) { MPX_UADDN(dv
, dvl
, n
); }
743 /* --- @mpx_uaddnlsl@ --- *
745 * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
746 * @mpw a@ = second argument
747 * @unsigned o@ = offset in bits
751 * Use: Computes %$d + 2^o a$%. If the result overflows then
752 * high-order bits are discarded, as usual. We must have
753 * @0 < o < MPW_BITS@.
756 void mpx_uaddnlsl(mpw
*dv
, mpw
*dvl
, mpw a
, unsigned o
)
760 while (x
&& dv
< dvl
) {
767 /* --- @mpx_usub@ --- *
769 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
770 * @const mpw *av, *avl@ = first argument vector base and limit
771 * @const mpw *bv, *bvl@ = second argument vector base and limit
775 * Use: Performs unsigned integer subtraction. If the result
776 * overflows the destination vector, high-order bits are
777 * discarded. This means that two's complement subtraction
778 * happens more or less for free, althuogh that's more a side-
779 * effect than anything else. The result vector may be equal to
780 * either or both source vectors, but may not otherwise overlap
784 void mpx_usub(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
,
785 const mpw
*bv
, const mpw
*bvl
)
789 while (av
< avl
|| bv
< bvl
) {
794 a
= (av
< avl
) ?
*av
++ : 0;
795 b
= (bv
< bvl
) ?
*bv
++ : 0;
796 x
= (mpd
)a
- (mpd
)b
- c
;
809 /* --- @mpx_usubn@ --- *
811 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
816 * Use: Subtracts a small integer from a multiprecision number.
819 void mpx_usubn(mpw
*dv
, mpw
*dvl
, mpw n
) { MPX_USUBN(dv
, dvl
, n
); }
821 /* --- @mpx_usubnlsl@ --- *
823 * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
824 * @mpw a@ = second argument
825 * @unsigned o@ = offset in bits
829 * Use: Computes %$d + 2^o a$%. If the result overflows then
830 * high-order bits are discarded, as usual. We must have
831 * @0 < o < MPW_BITS@.
834 void mpx_usubnlsl(mpw
*dv
, mpw
*dvl
, mpw a
, unsigned o
)
836 mpw b
= a
>> (MPW_BITS
- o
);
840 mpd x
= (mpd
)*dv
- MPW(a
);
844 MPX_USUBN(dv
, dvl
, b
);
848 /* --- @mpx_umul@ --- *
850 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
851 * @const mpw *av, *avl@ = multiplicand vector base and limit
852 * @const mpw *bv, *bvl@ = multiplier vector base and limit
856 * Use: Performs unsigned integer multiplication. If the result
857 * overflows the desination vector, high-order bits are
858 * discarded. The result vector may not overlap the argument
859 * vectors in any way.
862 CPU_DISPATCH(EMPTY
, (void), void, mpx_umul
,
863 (mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
,
864 const mpw
*bv
, const mpw
*bvl
),
865 (dv
, dvl
, av
, avl
, bv
, bvl
), pick_umul
, simple_umul
);
867 static void simple_umul(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
,
868 const mpw
*bv
, const mpw
*bvl
)
870 /* --- This is probably worthwhile on a multiply --- */
875 /* --- Deal with a multiply by zero --- */
882 /* --- Do the initial multiply and initialize the accumulator --- */
884 MPX_UMULN(dv
, dvl
, av
, avl
, *bv
++);
886 /* --- Do the remaining multiply/accumulates --- */
888 while (dv
< dvl
&& bv
< bvl
) {
898 x
= (mpd
)*dvv
+ (mpd
)m
* (mpd
)*avv
++ + c
;
902 MPX_UADDN(dvv
, dvl
, c
);
907 #define MAYBE_UMUL4(impl) \
908 extern void mpx_umul4_##impl(mpw */*dv*/, \
909 const mpw */*av*/, const mpw */*avl*/, \
910 const mpw */*bv*/, const mpw */*bvl*/); \
911 static void maybe_umul4_##impl(mpw *dv, mpw *dvl, \
912 const mpw *av, const mpw *avl, \
913 const mpw *bv, const mpw *bvl) \
915 size_t an = avl - av, bn = bvl - bv, dn = dvl - dv; \
916 if (!an || an%4 != 0 || !bn || bn%4 != 0 || dn < an + bn) \
917 simple_umul(dv, dvl, av, avl, bv, bvl); \
919 mpx_umul4_##impl(dv, av, avl, bv, bvl); \
920 MPX_ZERO(dv + an + bn, dvl); \
925 MAYBE_UMUL4(x86_sse2
)
930 MAYBE_UMUL4(amd64_sse2
)
931 MAYBE_UMUL4(amd64_avx
)
935 MAYBE_UMUL4(arm_neon
)
939 MAYBE_UMUL4(arm64_simd
)
942 static mpx_umul__functype
*pick_umul(void)
945 DISPATCH_PICK_COND(mpx_umul
, maybe_umul4_x86_avx
,
946 cpu_feature_p(CPUFEAT_X86_AVX
));
947 DISPATCH_PICK_COND(mpx_umul
, maybe_umul4_x86_sse2
,
948 cpu_feature_p(CPUFEAT_X86_SSE2
));
951 DISPATCH_PICK_COND(mpx_umul
, maybe_umul4_amd64_avx
,
952 cpu_feature_p(CPUFEAT_X86_AVX
));
953 DISPATCH_PICK_COND(mpx_umul
, maybe_umul4_amd64_sse2
,
954 cpu_feature_p(CPUFEAT_X86_SSE2
));
957 DISPATCH_PICK_COND(mpx_umul
, maybe_umul4_arm_neon
,
958 cpu_feature_p(CPUFEAT_ARM_NEON
));
961 DISPATCH_PICK_COND(mpx_umul
, maybe_umul4_arm64_simd
,
962 cpu_feature_p(CPUFEAT_ARM_NEON
));
964 DISPATCH_PICK_FALLBACK(mpx_umul
, simple_umul
);
967 /* --- @mpx_umuln@ --- *
969 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
970 * @const mpw *av, *avl@ = multiplicand vector base and limit
971 * @mpw m@ = multiplier
975 * Use: Multiplies a multiprecision integer by a single-word value.
976 * The destination and source may be equal. The destination
977 * is completely cleared after use.
980 void mpx_umuln(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
, mpw m
)
981 { MPX_UMULN(dv
, dvl
, av
, avl
, m
); }
983 /* --- @mpx_umlan@ --- *
985 * Arguments: @mpw *dv, *dvl@ = destination/accumulator base and limit
986 * @const mpw *av, *avl@ = multiplicand vector base and limit
987 * @mpw m@ = multiplier
991 * Use: Multiplies a multiprecision integer by a single-word value
992 * and adds the result to an accumulator.
995 void mpx_umlan(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
, mpw m
)
996 { MPX_UMLAN(dv
, dvl
, av
, avl
, m
); }
998 /* --- @mpx_usqr@ --- *
1000 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
1001 * @const mpw *av, *av@ = source vector base and limit
1005 * Use: Performs unsigned integer squaring. The result vector must
1006 * not overlap the source vector in any way.
1009 void mpx_usqr(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
)
1013 /* --- Main loop --- */
1016 const mpw
*avv
= av
;
1021 /* --- Stop if I've run out of destination --- */
1026 /* --- Work out the square at this point in the proceedings --- */
1029 mpd x
= (mpd
)a
* (mpd
)a
+ *dvv
;
1031 c
= MPW(x
>> MPW_BITS
);
1034 /* --- Now fix up the rest of the vector upwards --- */
1037 while (dvv
< dvl
&& avv
< avl
) {
1038 mpd x
= (mpd
)a
* (mpd
)*avv
++;
1039 mpd y
= ((x
<< 1) & MPW_MAX
) + c
+ *dvv
;
1040 c
= (x
>> (MPW_BITS
- 1)) + (y
>> MPW_BITS
);
1043 while (dvv
< dvl
&& c
) {
1049 /* --- Get ready for the next round --- */
1056 /* --- @mpx_udiv@ --- *
1058 * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit
1059 * @mpw *rv, *rvl@ = dividend/remainder vector base and limit
1060 * @const mpw *dv, *dvl@ = divisor vector base and limit
1061 * @mpw *sv, *svl@ = scratch workspace
1065 * Use: Performs unsigned integer division. If the result overflows
1066 * the quotient vector, high-order bits are discarded. (Clearly
1067 * the remainder vector can't overflow.) The various vectors
1068 * may not overlap in any way. Yes, I know it's a bit odd
1069 * requiring the dividend to be in the result position but it
1070 * does make some sense really. The remainder must have
1071 * headroom for at least two extra words. The scratch space
1072 * must be at least one word larger than the divisor.
1075 void mpx_udiv(mpw
*qv
, mpw
*qvl
, mpw
*rv
, mpw
*rvl
,
1076 const mpw
*dv
, const mpw
*dvl
,
1083 /* --- Initialize the quotient --- */
1087 /* --- Perform some sanity checks --- */
1089 MPX_SHRINK(dv
, dvl
);
1090 assert(((void)"division by zero in mpx_udiv", dv
< dvl
));
1092 /* --- Normalize the divisor --- *
1094 * The algorithm requires that the divisor be at least two digits long.
1095 * This is easy to fix.
1102 for (b
= MPW_P2
; b
; b
>>= 1) {
1103 if (d
<= (MPW_MAX
>> b
)) {
1112 /* --- Normalize the dividend/remainder to match --- */
1115 mpx_lsl(rv
, rvl
, rv
, rvl
, norm
);
1116 mpx_lsl(sv
, svl
, dv
, dvl
, norm
);
1119 MPX_SHRINK(dv
, dvl
);
1122 MPX_SHRINK(rv
, rvl
);
1126 /* --- Work out the relative scales --- */
1129 size_t rvn
= rvl
- rv
;
1130 size_t dvn
= dvl
- dv
;
1132 /* --- If the divisor is clearly larger, notice this --- */
1135 mpx_lsr(rv
, rvl
, rv
, rvl
, norm
);
1142 /* --- Calculate the most significant quotient digit --- *
1144 * Because the divisor has its top bit set, this can only happen once. The
1145 * pointer arithmetic is a little contorted, to make sure that the
1146 * behaviour is defined.
1149 if (MPX_UCMP(rv
+ scale
, rvl
, >=, dv
, dvl
)) {
1150 mpx_usub(rv
+ scale
, rvl
, rv
+ scale
, rvl
, dv
, dvl
);
1151 if (qvl
- qv
> scale
)
1155 /* --- Now for the main loop --- */
1164 /* --- Get an estimate for the next quotient digit --- */
1171 rh
= ((mpd
)r
<< MPW_BITS
) | rr
;
1177 /* --- Refine the estimate --- */
1180 mpd yh
= (mpd
)d
* q
;
1181 mpd yy
= (mpd
)dd
* q
;
1185 yh
+= yy
>> MPW_BITS
;
1188 while (yh
> rh
|| (yh
== rh
&& yl
> rrr
)) {
1197 /* --- Remove a chunk from the dividend --- */
1204 /* --- Calculate the size of the chunk --- *
1206 * This does the whole job of calculating @r >> scale - qd@.
1209 for (svv
= rv
+ scale
, dvv
= dv
;
1210 dvv
< dvl
&& svv
< rvl
;
1212 mpd x
= (mpd
)*dvv
* (mpd
)q
+ mc
;
1214 x
= (mpd
)*svv
- MPW(x
) - sc
;
1223 mpd x
= (mpd
)*svv
- mc
- sc
;
1233 /* --- Fix if the quotient was too large --- *
1235 * This doesn't seem to happen very often.
1238 if (rvl
[-1] > MPW_MAX
/ 2) {
1239 mpx_uadd(rv
+ scale
, rvl
, rv
+ scale
, rvl
, dv
, dvl
);
1244 /* --- Done for another iteration --- */
1246 if (qvl
- qv
> scale
)
1253 /* --- Now fiddle with unnormalizing and things --- */
1255 mpx_lsr(rv
, rvl
, rv
, rvl
, norm
);
1258 /* --- @mpx_udivn@ --- *
1260 * Arguments: @mpw *qv, *qvl@ = storage for the quotient (may overlap
1262 * @const mpw *rv, *rvl@ = dividend
1263 * @mpw d@ = single-precision divisor
1265 * Returns: Remainder after divison.
1267 * Use: Performs a single-precision division operation.
1270 mpw
mpx_udivn(mpw
*qv
, mpw
*qvl
, const mpw
*rv
, const mpw
*rvl
, mpw d
)
1273 size_t ql
= qvl
- qv
;
1279 r
= (r
<< MPW_BITS
) | rv
[i
];
1287 /*----- Test rig ----------------------------------------------------------*/
1291 #include <mLib/alloc.h>
1292 #include <mLib/dstr.h>
1293 #include <mLib/macros.h>
1294 #include <mLib/quis.h>
1295 #include <mLib/testrig.h>
1297 #ifdef ENABLE_ASM_DEBUG
1298 # include "regdump.h"
1303 #define ALLOC(v, vl, sz) do { \
1304 size_t _sz = (sz); \
1305 mpw *_vv = xmalloc(MPWS(_sz)); \
1306 mpw *_vvl = _vv + _sz; \
1307 memset(_vv, 0xa5, MPWS(_sz)); \
1312 #define LOAD(v, vl, d) do { \
1313 const dstr *_d = (d); \
1315 ALLOC(_v, _vl, MPW_RQ(_d->len)); \
1316 mpx_loadb(_v, _vl, _d->buf, _d->len); \
1321 #define MAX(x, y) ((x) > (y) ? (x) : (y))
1323 static void dumpbits(const char *msg
, const void *pp
, size_t sz
)
1325 const octet
*p
= pp
;
1328 fprintf(stderr
, " %02x", *p
++);
1329 fputc('\n', stderr
);
1332 static void dumpmp(const char *msg
, const mpw
*v
, const mpw
*vl
)
1337 fprintf(stderr
, " %08lx", (unsigned long)*--vl
);
1338 fputc('\n', stderr
);
1341 static int chkscan(const mpw
*v
, const mpw
*vl
,
1342 const void *pp
, size_t sz
, int step
)
1345 const octet
*p
= pp
;
1349 mpscan_initx(&mps
, v
, vl
);
1354 for (i
= 0; i
< 8 && MPSCAN_STEP(&mps
); i
++) {
1355 if (MPSCAN_BIT(&mps
) != (x
& 1)) {
1357 "\n*** error, step %i, bit %u, expected %u, found %u\n",
1358 step
, bit
, x
& 1, MPSCAN_BIT(&mps
));
1370 static int loadstore(dstr
*v
)
1373 size_t sz
= MPW_RQ(v
->len
) * 2, diff
;
1377 dstr_ensure(&d
, v
->len
);
1378 m
= xmalloc(MPWS(sz
));
1380 for (diff
= 0; diff
< sz
; diff
+= 5) {
1385 mpx_loadl(m
, ml
, v
->buf
, v
->len
);
1386 if (!chkscan(m
, ml
, v
->buf
, v
->len
, +1))
1388 MPX_OCTETS(oct
, m
, ml
);
1389 mpx_storel(m
, ml
, d
.buf
, d
.sz
);
1390 if (MEMCMP(d
.buf
, !=, v
->buf
, oct
)) {
1391 dumpbits("\n*** storel failed", d
.buf
, d
.sz
);
1395 mpx_loadb(m
, ml
, v
->buf
, v
->len
);
1396 if (!chkscan(m
, ml
, v
->buf
+ v
->len
- 1, v
->len
, -1))
1398 MPX_OCTETS(oct
, m
, ml
);
1399 mpx_storeb(m
, ml
, d
.buf
, d
.sz
);
1400 if (MEMCMP(d
.buf
+ d
.sz
- oct
, !=, v
->buf
+ v
->len
- oct
, oct
)) {
1401 dumpbits("\n*** storeb failed", d
.buf
, d
.sz
);
1407 dumpbits("input data", v
->buf
, v
->len
);
1414 static int twocl(dstr
*v
)
1418 size_t sz0
, sz1
, szmax
;
1422 sz0
= MPW_RQ(v
[0].len
); sz1
= MPW_RQ(v
[1].len
);
1423 dstr_ensure(&d
, v
[0].len
> v
[1].len ? v
[0].len
: v
[1].len
);
1425 szmax
= sz0
> sz1 ? sz0
: sz1
;
1426 m
= xmalloc(MPWS(szmax
));
1427 ml0
= m
+ sz0
; ml1
= m
+ sz1
;
1429 for (i
= 0; i
< 2; i
++) {
1430 if (i
) ml0
= ml1
= m
+ szmax
;
1432 mpx_loadl(m
, ml0
, v
[0].buf
, v
[0].len
);
1433 mpx_storel2cn(m
, ml0
, d
.buf
, v
[1].len
);
1434 if (MEMCMP(d
.buf
, !=, v
[1].buf
, v
[1].len
)) {
1435 dumpbits("\n*** storel2cn failed", d
.buf
, v
[1].len
);
1439 mpx_loadl2cn(m
, ml1
, v
[1].buf
, v
[1].len
);
1440 mpx_storel(m
, ml1
, d
.buf
, v
[0].len
);
1441 if (MEMCMP(d
.buf
, !=, v
[0].buf
, v
[0].len
)) {
1442 dumpbits("\n*** loadl2cn failed", d
.buf
, v
[0].len
);
1448 dumpbits("pos", v
[0].buf
, v
[0].len
);
1449 dumpbits("neg", v
[1].buf
, v
[1].len
);
1458 static int twocb(dstr
*v
)
1462 size_t sz0
, sz1
, szmax
;
1466 sz0
= MPW_RQ(v
[0].len
); sz1
= MPW_RQ(v
[1].len
);
1467 dstr_ensure(&d
, v
[0].len
> v
[1].len ? v
[0].len
: v
[1].len
);
1469 szmax
= sz0
> sz1 ? sz0
: sz1
;
1470 m
= xmalloc(MPWS(szmax
));
1471 ml0
= m
+ sz0
; ml1
= m
+ sz1
;
1473 for (i
= 0; i
< 2; i
++) {
1474 if (i
) ml0
= ml1
= m
+ szmax
;
1476 mpx_loadb(m
, ml0
, v
[0].buf
, v
[0].len
);
1477 mpx_storeb2cn(m
, ml0
, d
.buf
, v
[1].len
);
1478 if (MEMCMP(d
.buf
, !=, v
[1].buf
, v
[1].len
)) {
1479 dumpbits("\n*** storeb2cn failed", d
.buf
, v
[1].len
);
1483 mpx_loadb2cn(m
, ml1
, v
[1].buf
, v
[1].len
);
1484 mpx_storeb(m
, ml1
, d
.buf
, v
[0].len
);
1485 if (MEMCMP(d
.buf
, !=, v
[0].buf
, v
[0].len
)) {
1486 dumpbits("\n*** loadb2cn failed", d
.buf
, v
[0].len
);
1492 dumpbits("pos", v
[0].buf
, v
[0].len
);
1493 dumpbits("neg", v
[1].buf
, v
[1].len
);
1502 static int lsl(dstr
*v
)
1505 int n
= *(int *)v
[1].buf
;
1512 ALLOC(d
, dl
, al
- a
+ (n
+ MPW_BITS
- 1) / MPW_BITS
);
1514 mpx_lsl(d
, dl
, a
, al
, n
);
1515 if (!mpx_ueq(d
, dl
, c
, cl
)) {
1516 fprintf(stderr
, "\n*** lsl(%i) failed\n", n
);
1517 dumpmp(" a", a
, al
);
1518 dumpmp("expected", c
, cl
);
1519 dumpmp(" result", d
, dl
);
1523 xfree(a
); xfree(c
); xfree(d
);
1527 static int lslc(dstr
*v
)
1530 int n
= *(int *)v
[1].buf
;
1537 ALLOC(d
, dl
, al
- a
+ (n
+ MPW_BITS
- 1) / MPW_BITS
);
1539 mpx_lslc(d
, dl
, a
, al
, n
);
1540 if (!mpx_ueq(d
, dl
, c
, cl
)) {
1541 fprintf(stderr
, "\n*** lslc(%i) failed\n", n
);
1542 dumpmp(" a", a
, al
);
1543 dumpmp("expected", c
, cl
);
1544 dumpmp(" result", d
, dl
);
1548 xfree(a
); xfree(c
); xfree(d
);
1552 static int lsr(dstr
*v
)
1555 int n
= *(int *)v
[1].buf
;
1562 ALLOC(d
, dl
, al
- a
+ (n
+ MPW_BITS
- 1) / MPW_BITS
+ 1);
1564 mpx_lsr(d
, dl
, a
, al
, n
);
1565 if (!mpx_ueq(d
, dl
, c
, cl
)) {
1566 fprintf(stderr
, "\n*** lsr(%i) failed\n", n
);
1567 dumpmp(" a", a
, al
);
1568 dumpmp("expected", c
, cl
);
1569 dumpmp(" result", d
, dl
);
1573 xfree(a
); xfree(c
); xfree(d
);
1577 static int uadd(dstr
*v
)
1588 ALLOC(d
, dl
, MAX(al
- a
, bl
- b
) + 1);
1590 mpx_uadd(d
, dl
, a
, al
, b
, bl
);
1591 if (!mpx_ueq(d
, dl
, c
, cl
)) {
1592 fprintf(stderr
, "\n*** uadd failed\n");
1593 dumpmp(" a", a
, al
);
1594 dumpmp(" b", b
, bl
);
1595 dumpmp("expected", c
, cl
);
1596 dumpmp(" result", d
, dl
);
1600 xfree(a
); xfree(b
); xfree(c
); xfree(d
);
1604 static int usub(dstr
*v
)
1615 ALLOC(d
, dl
, al
- a
);
1617 mpx_usub(d
, dl
, a
, al
, b
, bl
);
1618 if (!mpx_ueq(d
, dl
, c
, cl
)) {
1619 fprintf(stderr
, "\n*** usub failed\n");
1620 dumpmp(" a", a
, al
);
1621 dumpmp(" b", b
, bl
);
1622 dumpmp("expected", c
, cl
);
1623 dumpmp(" result", d
, dl
);
1627 xfree(a
); xfree(b
); xfree(c
); xfree(d
);
1631 static int umul(dstr
*v
)
1642 ALLOC(d
, dl
, (al
- a
) + (bl
- b
));
1644 mpx_umul(d
, dl
, a
, al
, b
, bl
);
1645 if (!mpx_ueq(d
, dl
, c
, cl
)) {
1646 fprintf(stderr
, "\n*** umul failed\n");
1647 dumpmp(" a", a
, al
);
1648 dumpmp(" b", b
, bl
);
1649 dumpmp("expected", c
, cl
);
1650 dumpmp(" result", d
, dl
);
1654 xfree(a
); xfree(b
); xfree(c
); xfree(d
);
1658 static int usqr(dstr
*v
)
1667 ALLOC(d
, dl
, 2 * (al
- a
));
1669 mpx_usqr(d
, dl
, a
, al
);
1670 if (!mpx_ueq(d
, dl
, c
, cl
)) {
1671 fprintf(stderr
, "\n*** usqr failed\n");
1672 dumpmp(" a", a
, al
);
1673 dumpmp("expected", c
, cl
);
1674 dumpmp(" result", d
, dl
);
1678 xfree(a
); xfree(c
); xfree(d
);
1682 static int udiv(dstr
*v
)
1692 ALLOC(a
, al
, MPW_RQ(v
[0].len
) + 2); mpx_loadb(a
, al
, v
[0].buf
, v
[0].len
);
1696 ALLOC(qq
, qql
, al
- a
);
1697 ALLOC(s
, sl
, (bl
- b
) + 1);
1699 mpx_udiv(qq
, qql
, a
, al
, b
, bl
, s
, sl
);
1700 if (!mpx_ueq(qq
, qql
, q
, ql
) ||
1701 !mpx_ueq(a
, al
, r
, rl
)) {
1702 fprintf(stderr
, "\n*** udiv failed\n");
1703 dumpmp(" divisor", b
, bl
);
1704 dumpmp("expect r", r
, rl
);
1705 dumpmp("result r", a
, al
);
1706 dumpmp("expect q", q
, ql
);
1707 dumpmp("result q", qq
, qql
);
1711 xfree(a
); xfree(b
); xfree(r
); xfree(q
); xfree(s
); xfree(qq
);
1715 static test_chunk defs
[] = {
1716 { "load-store", loadstore
, { &type_hex
, 0 } },
1717 { "2cl", twocl
, { &type_hex
, &type_hex
, } },
1718 { "2cb", twocb
, { &type_hex
, &type_hex
, } },
1719 { "lsl", lsl
, { &type_hex
, &type_int
, &type_hex
, 0 } },
1720 { "lslc", lslc
, { &type_hex
, &type_int
, &type_hex
, 0 } },
1721 { "lsr", lsr
, { &type_hex
, &type_int
, &type_hex
, 0 } },
1722 { "uadd", uadd
, { &type_hex
, &type_hex
, &type_hex
, 0 } },
1723 { "usub", usub
, { &type_hex
, &type_hex
, &type_hex
, 0 } },
1724 { "umul", umul
, { &type_hex
, &type_hex
, &type_hex
, 0 } },
1725 { "usqr", usqr
, { &type_hex
, &type_hex
, 0 } },
1726 { "udiv", udiv
, { &type_hex
, &type_hex
, &type_hex
, &type_hex
, 0 } },
1730 int main(int argc
, char *argv
[])
1732 #ifdef ENABLE_ASM_DEBUG
1735 test_run(argc
, argv
, defs
, SRCDIR
"/t/mpx");
1741 /*----- That's all, folks -------------------------------------------------*/