3 * $Id: mpx.h,v 1.18 2004/04/08 01:36:15 mdw Exp $
5 * Low level multiprecision arithmetic
7 * (c) 1999 Straylight/Edgeware
10 /*----- Licensing notice --------------------------------------------------*
12 * This file is part of Catacomb.
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.
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.
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,
30 #ifndef CATACOMB_MPX_H
31 #define CATACOMB_MPX_H
37 /*----- The idea ----------------------------------------------------------*
39 * This file provides functions and macros which work on vectors of words as
40 * unsigned multiprecision integers. The interface works in terms of base
41 * and limit pointers (i.e., a pointer to the start of a vector, and a
42 * pointer just past its end) rather than base pointer and length, because
43 * that requires more arithmetic and state to work on.
45 * The interfaces are slightly bizarre in other ways. Try to use the
46 * higher-level functions where you can: they're rather better designed to
47 * actually be friendly and useful.
50 /*----- Header files ------------------------------------------------------*/
54 #ifndef CATACOMB_MPW_H
58 /*----- General manipulation ----------------------------------------------*/
60 /* --- @MPX_SHRINK@ --- *
62 * Arguments: @const mpw *v@ = pointer to vector of words
63 * @const mpw *vl@ = (updated) current limit of vector
65 * Use: Shrinks down the limit of a multiprecision integer vector.
68 #define MPX_SHRINK(v, vl) do { \
69 const mpw *_vv = (v), *_vvl = (vl); \
70 while (_vvl > _vv && !_vvl[-1]) \
75 /* --- @MPX_BITS@ --- *
77 * Arguments: @unsigned long b@ = result variable
78 * @const mpw *v@ = pointer to array of words
79 * @const mpw *vl@ = limit of vector (from @MPX_SHRINK@)
81 * Use: Calculates the number of bits in a multiprecision value.
84 #define MPX_BITS(b, v, vl) do { \
85 const mpw *_v = (v), *_vl = (vl); \
86 MPX_SHRINK(_v, _vl); \
90 unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1; \
92 unsigned _k = MPW_BITS / 2; \
104 /* --- @MPX_OCTETS@ --- *
106 * Arguments: @size_t o@ = result variable
107 * @const mpw *v, *vl@ = pointer to array of words
109 * Use: Calculates the number of octets in a multiprecision value.
112 #define MPX_OCTETS(o, v, vl) do { \
114 MPX_BITS(_bb, (v), (vl)); \
115 (o) = (_bb + 7) >> 3; \
118 /* --- @MPX_OCTETS2C@ --- *
120 * Arguments: @size_t o@ = result variable
121 * @const mpw *v, *vl@ = pointer to array of words
123 * Use: Calculates the number of octets in a multiprecision value, if
124 * you represent it as two's complement.
127 #define MPX_OCTETS2C(o, v, vl) do { \
129 MPX_BITS(_bb, (v), (vl)); \
130 (o) = (_bb >> 3) + 1; \
133 /* --- @MPX_COPY@ --- *
135 * Arguments: @dv, dvl@ = destination vector base and limit
136 * @av, avl@ = source vector base and limit
138 * Use: Copies a multiprecision integer.
141 #define MPX_COPY(dv, dvl, av, avl) do { \
142 mpw *_dv = (dv), *_dvl = (dvl); \
143 size_t _dn = _dvl - _dv; \
144 const mpw *_av = (av), *_avl = (avl); \
145 size_t _an = _avl - _av; \
148 memset(_dv, 0, MPWS(_dn - _an)); \
149 } else if (_an >= _dn) \
150 memmove(_dv, _av, MPWS(_dn)); \
152 memmove(_dv, _av, MPWS(_an)); \
153 memset(_dv + _an, 0, MPWS(_dn - _an)); \
157 /* --- @MPX_ZERO@ --- *
159 * Arguments: @v, vl@ = base and limit of vector to clear
161 * Use: Zeroes the area between the two vector pointers.
164 #define MPX_ZERO(v, vl) do { \
165 mpw *_v = (v), *_vl = (vl); \
167 memset(_v, 0, MPWS(_vl - _v)); \
170 /* --- @MPX_ONE@ --- *
172 * Arguments: @v, vl@ = base and limit of vector to clear
174 * Use: Fills the area between the two vector pointers with ones.
177 #define MPX_ONE(v, vl) do { \
179 const mpw *_vl = (vl); \
184 /*----- Loading and storing -----------------------------------------------*/
186 /* --- @mpx_storel@ --- *
188 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
189 * @void *p@ = pointer to octet array
190 * @size_t sz@ = size of octet array
194 * Use: Stores an MP in an octet array, least significant octet
195 * first. High-end octets are silently discarded if there
196 * isn't enough space for them.
199 extern void mpx_storel(const mpw */
*v*/
, const mpw */
*vl*/
,
200 void */
*p*/
, size_t /*sz*/);
202 /* --- @mpx_loadl@ --- *
204 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
205 * @const void *p@ = pointer to octet array
206 * @size_t sz@ = size of octet array
210 * Use: Loads an MP in an octet array, least significant octet
211 * first. High-end octets are ignored if there isn't enough
215 extern void mpx_loadl(mpw */
*v*/
, mpw */
*vl*/
,
216 const void */
*p*/
, size_t /*sz*/);
218 /* --- @mpx_storeb@ --- *
220 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
221 * @void *p@ = pointer to octet array
222 * @size_t sz@ = size of octet array
226 * Use: Stores an MP in an octet array, most significant octet
227 * first. High-end octets are silently discarded if there
228 * isn't enough space for them.
231 extern void mpx_storeb(const mpw */
*v*/
, const mpw */
*vl*/
,
232 void */
*p*/
, size_t /*sz*/);
234 /* --- @mpx_loadb@ --- *
236 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
237 * @const void *p@ = pointer to octet array
238 * @size_t sz@ = size of octet array
242 * Use: Loads an MP in an octet array, most significant octet
243 * first. High-end octets are ignored if there isn't enough
247 extern void mpx_loadb(mpw */
*v*/
, mpw */
*vl*/
,
248 const void */
*p*/
, size_t /*sz*/);
250 /* --- @mpx_storel2cn@ --- *
252 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
253 * @void *pp@ = pointer to octet array
254 * @size_t sz@ = size of octet array
258 * Use: Stores a negative MP in an octet array, least significant
259 * octet first, as two's complement. High-end octets are
260 * silently discarded if there isn't enough space for them.
261 * This obviously makes the output bad.
264 extern void mpx_storel2cn(const mpw */
*v*/
, const mpw */
*vl*/
,
265 void */
*p*/
, size_t /*sz*/);
267 /* --- @mpx_loadl2cn@ --- *
269 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
270 * @const void *pp@ = pointer to octet array
271 * @size_t sz@ = size of octet array
275 * Use: Loads a negative MP in an octet array, least significant
276 * octet first, as two's complement. High-end octets are
277 * ignored if there isn't enough space for them. This probably
278 * means you made the wrong choice coming here.
281 extern void mpx_loadl2cn(mpw */
*v*/
, mpw */
*vl*/
,
282 const void */
*p*/
, size_t /*sz*/);
284 /* --- @mpx_storeb2cn@ --- *
286 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
287 * @void *pp@ = pointer to octet array
288 * @size_t sz@ = size of octet array
292 * Use: Stores a negative MP in an octet array, most significant
293 * octet first, as two's complement. High-end octets are
294 * silently discarded if there isn't enough space for them,
295 * which probably isn't what you meant.
298 extern void mpx_storeb2cn(const mpw */
*v*/
, const mpw */
*vl*/
,
299 void */
*p*/
, size_t /*sz*/);
301 /* --- @mpx_loadb2cn@ --- *
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
309 * Use: Loads a negative MP in an octet array, most significant octet
310 * first as two's complement. High-end octets are ignored if
311 * there isn't enough space for them. This probably means you
312 * chose this function wrongly.
315 extern void mpx_loadb2cn(mpw */
*v*/
, mpw */
*vl*/
,
316 const void */
*p*/
, size_t /*sz*/);
319 /*----- Logical shifting --------------------------------------------------*/
321 /* --- @mpx_lsl@ --- *
323 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
324 * @const mpw *av, *avl@ = source vector base and limit
325 * @size_t n@ = number of bit positions to shift by
329 * Use: Performs a logical shift left operation on an integer.
332 extern void mpx_lsl(mpw */
*dv*/
, mpw */
*dvl*/
,
333 const mpw */
*av*/
, const mpw */
*avl*/
,
336 /* --- @mpx_lslc@ --- *
338 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
339 * @const mpw *av, *avl@ = source vector base and limit
340 * @size_t n@ = number of bit positions to shift by
344 * Use: Performs a logical shift left operation on an integer, only
345 * it fills in the bits with ones instead of zeroes.
348 extern void mpx_lslc(mpw */
*dv*/
, mpw */
*dvl*/
,
349 const mpw */
*av*/
, const mpw */
*avl*/
,
352 /* --- @mpx_lsr@ --- *
354 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
355 * @const mpw *av, *avl@ = source vector base and limit
356 * @size_t n@ = number of bit positions to shift by
360 * Use: Performs a logical shift right operation on an integer.
363 extern void mpx_lsr(mpw */
*dv*/
, mpw */
*dvl*/
,
364 const mpw */
*av*/
, const mpw */
*avl*/
,
367 /*----- Bitwise operations ------------------------------------------------*/
369 /* --- @mpx_bitop@ --- *
371 * Arguments: @mpw *dv, *dvl@ = destination vector
372 * @const mpw *av, *avl@ = first source vector
373 * @const mpw *bv, *bvl@ = second source vector
377 * Use: Provide the dyadic boolean functions. The functions are
378 * named after the truth table they generate:
385 #define MPX_DOBIN(what) \
386 what(0000) what(0001) what(0010) what(0011) \
387 what(0100) what(0101) what(0110) what(0111) \
388 what(1000) what(1001) what(1010) what(1011) \
389 what(1100) what(1101) what(1110) what(1111)
391 #define MPX_BITDECL(string) \
392 extern void mpx_bit##string(mpw */*dv*/, mpw */*dvl*/, \
393 const mpw */*av*/, const mpw */*avl*/, \
394 const mpw */*bv*/, const mpw */*bvl*/);
395 MPX_DOBIN(MPX_BITDECL
)
397 /* --- @mpx_[n]and@, @mpx_[n]or@, @mpx_xor@ --- *
399 * Synonyms for the commonly-used functions above.
402 #define mpx_and mpx_bit0001
403 #define mpx_or mpx_bit0111
404 #define mpx_nand mpx_bit1110
405 #define mpx_nor mpx_bit1000
406 #define mpx_xor mpx_bit0110
408 /* --- @mpx_not@ --- *
410 * Arguments: @mpw *dv, *dvl@ = destination vector
411 * @const mpw *av, *avl@ = first source vector
418 extern void mpx_not(mpw */
*dv*/
, mpw */
*dvl*/
,
419 const mpw */
*av*/
, const mpw */
*avl*/
);
421 /*----- Unsigned arithmetic -----------------------------------------------*/
423 /* --- @mpx_2c@ --- *
425 * Arguments: @mpw *dv, *dvl@ = destination vector
426 * @const mpw *v, *vl@ = source vector
430 * Use: Calculates the two's complement of @v@.
433 extern void mpx_2c(mpw */
*dv*/
, mpw */
*dvl*/
,
434 const mpw */
*v*/
, const mpw */
*vl*/
);
436 /* --- @mpx_ueq@ --- *
438 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
439 * @const mpw *bv, *bvl@ = second argument vector base and limit
441 * Returns: Nonzero if the two vectors are equal.
443 * Use: Performs an unsigned integer test for equality.
446 extern int mpx_ueq(const mpw */
*av*/
, const mpw */
*avl*/
,
447 const mpw */
*bv*/
, const mpw */
*bvl*/
);
449 /* --- @mpx_ucmp@ --- *
451 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
452 * @const mpw *bv, *bvl@ = second argument vector base and limit
454 * Returns: Less than, equal to, or greater than zero depending on
455 * whether @a@ is less than, equal to or greater than @b@,
458 * Use: Performs an unsigned integer comparison.
461 #define MPX_UCMP(av, avl, op, dv, dvl) \
462 (mpx_ucmp((av), (avl), (dv), (dvl)) op 0)
464 extern int mpx_ucmp(const mpw */
*av*/
, const mpw */
*avl*/
,
465 const mpw */
*bv*/
, const mpw */
*bvl*/
);
467 /* --- @mpx_uadd@ --- *
469 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
470 * @const mpw *av, *avl@ = first addend vector base and limit
471 * @const mpw *bv, *bvl@ = second addend vector base and limit
475 * Use: Performs unsigned integer addition. If the result overflows
476 * the destination vector, high-order bits are discarded. This
477 * means that two's complement addition happens more or less for
478 * free, although that's more a side-effect than anything else.
479 * The result vector may be equal to either or both source
480 * vectors, but may not otherwise overlap them.
483 extern void mpx_uadd(mpw */
*dv*/
, mpw */
*dvl*/
,
484 const mpw */
*av*/
, const mpw */
*avl*/
,
485 const mpw */
*bv*/
, const mpw */
*bvl*/
);
487 /* --- @mpx_uaddn@ --- *
489 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
490 * @mpw n@ = other addend
494 * Use: Adds a small integer to a multiprecision number.
497 #define MPX_UADDN(dv, dvl, n) do { \
498 mpw *_ddv = (dv), *_ddvl = (dvl); \
501 while (_c && _ddv < _ddvl) { \
502 mpd _x = (mpd)*_ddv + (mpd)_c; \
504 _c = _x >> MPW_BITS; \
508 extern void mpx_uaddn(mpw */
*dv*/
, mpw */
*dvl*/
, mpw
/*n*/);
510 /* --- @mpx_uaddnlsl@ --- *
512 * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
513 * @mpw a@ = second argument
514 * @unsigned o@ = offset in bits
518 * Use: Computes %$d + 2^o a$%. If the result overflows then
519 * high-order bits are discarded, as usual. We must have
520 * @0 < o < MPW_BITS@.
523 extern void mpx_uaddnlsl(mpw */
*dv*/
, mpw */
*dvl*/
,
524 mpw
/*a*/, unsigned /*o*/);
526 /* --- @mpx_usub@ --- *
528 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
529 * @const mpw *av, *avl@ = first argument vector base and limit
530 * @const mpw *bv, *bvl@ = second argument vector base and limit
534 * Use: Performs unsigned integer subtraction. If the result
535 * overflows the destination vector, high-order bits are
536 * discarded. This means that two's complement subtraction
537 * happens more or less for free, although that's more a side-
538 * effect than anything else. The result vector may be equal to
539 * either or both source vectors, but may not otherwise overlap
543 extern void mpx_usub(mpw */
*dv*/
, mpw */
*dvl*/
,
544 const mpw */
*av*/
, const mpw */
*avl*/
,
545 const mpw */
*bv*/
, const mpw */
*bvl*/
);
547 /* --- @mpx_usubn@ --- *
549 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
554 * Use: Subtracts a small integer from a multiprecision number.
557 #define MPX_USUBN(dv, dvl, n) do { \
558 mpw *_ddv = (dv), *_ddvl = (dvl); \
561 while (_ddv < _ddvl) { \
562 mpd _x = (mpd)*_ddv - (mpd)_c; \
564 if (_x >> MPW_BITS) \
571 extern void mpx_usubn(mpw */
*dv*/
, mpw */
*dvl*/
, mpw
/*n*/);
573 /* --- @mpx_usubnlsl@ --- *
575 * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
576 * @mpw a@ = second argument
577 * @unsigned o@ = offset in bits
581 * Use: Computes %$d - 2^o a$%. If the result overflows then
582 * high-order bits are discarded, as usual, so you get two's
583 * complement. Which might be what you wanted... We must have
584 * @0 < o < MPW_BITS@.
587 extern void mpx_usubnlsl(mpw */
*dv*/
, mpw */
*dvl*/
,
588 mpw
/*a*/, unsigned /*o*/);
590 /* --- @mpx_umul@ --- *
592 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
593 * @const mpw *av, *avl@ = multiplicand vector base and limit
594 * @const mpw *bv, *bvl@ = multiplier vector base and limit
598 * Use: Performs unsigned integer multiplication. If the result
599 * overflows the desination vector, high-order bits are
600 * discarded. The result vector may not overlap the argument
601 * vectors in any way.
604 extern void mpx_umul(mpw */
*dv*/
, mpw */
*dvl*/
,
605 const mpw */
*av*/
, const mpw */
*avl*/
,
606 const mpw */
*bv*/
, const mpw */
*bvl*/
);
608 /* --- @mpx_umuln@ --- *
610 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
611 * @const mpw *av, *avl@ = multiplicand vector base and limit
612 * @mpw m@ = multiplier
616 * Use: Multiplies a multiprecision integer by a single-word value.
617 * The destination and source may be equal. The destination
618 * is completely cleared after use.
621 #define MPX_UMULN(dv, dvl, av, avl, m) do { \
622 mpw *_dv = (dv), *_dvl = (dvl); \
623 const mpw *_av = (av), *_avl = (avl); \
627 while (_av < _avl) { \
631 _x = (mpd)_m * (mpd)*_av++ + _c; \
633 _c = _x >> MPW_BITS; \
637 MPX_ZERO(_dv, _dvl); \
641 extern void mpx_umuln(mpw */
*dv*/
, mpw */
*dvl*/
,
642 const mpw */
*av*/
, const mpw */
*avl*/
, mpw m
);
644 /* --- @mpx_umlan@ --- *
646 * Arguments: @mpw *dv, *dvl@ = destination/accumulator base and limit
647 * @const mpw *av, *avl@ = multiplicand vector base and limit
648 * @mpw m@ = multiplier
652 * Use: Multiplies a multiprecision integer by a single-word value
653 * and adds the result to an accumulator.
656 #define MPX_UMLAN(dv, dvl, av, avl, m) do { \
657 mpw *_dv = (dv), *_dvl = (dvl); \
658 const mpw *_av = (av), *_avl = (avl); \
662 while (_dv < _dvl && _av < _avl) { \
664 _x = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _cc; \
666 _cc = _x >> MPW_BITS; \
668 MPX_UADDN(_dv, _dvl, _cc); \
671 extern void mpx_umlan(mpw */
*dv*/
, mpw */
*dvl*/
,
672 const mpw */
*av*/
, const mpw */
*avl*/
, mpw m
);
674 /* --- @mpx_usqr@ --- *
676 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
677 * @const mpw *av, *av@ = source vector base and limit
681 * Use: Performs unsigned integer squaring. The result vector must
682 * not overlap the source vector in any way.
685 extern void mpx_usqr(mpw */
*dv*/
, mpw */
*dvl*/
,
686 const mpw */
*av*/
, const mpw */
*avl*/
);
688 /* --- @mpx_udiv@ --- *
690 * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit
691 * @mpw *rv, *rvl@ = dividend/remainder vector base and limit
692 * @const mpw *dv, *dvl@ = divisor vector base and limit
693 * @mpw *sv, *svl@ = scratch workspace
697 * Use: Performs unsigned integer division. If the result overflows
698 * the quotient vector, high-order bits are discarded. (Clearly
699 * the remainder vector can't overflow.) The various vectors
700 * may not overlap in any way. Yes, I know it's a bit odd
701 * requiring the dividend to be in the result position but it
702 * does make some sense really. The remainder must have
703 * headroom for at least two extra words. The scratch space
704 * must be at least one word larger than the divisor.
707 extern void mpx_udiv(mpw */
*qv*/
, mpw */
*qvl*/
, mpw */
*rv*/
, mpw */
*rvl*/
,
708 const mpw */
*dv*/
, const mpw */
*dvl*/
,
709 mpw */
*sv*/
, mpw */
*svl*/
);
711 /* --- @mpx_udivn@ --- *
713 * Arguments: @mpw *qv, *qvl@ = storage for the quotient (may overlap
715 * @const mpw *rv, *rvl@ = dividend
716 * @mpw d@ = single-precision divisor
718 * Returns: Remainder after divison.
720 * Use: Performs a single-precision division operation.
723 extern mpw
mpx_udivn(mpw */
*qv*/
, mpw */
*qvl*/
,
724 const mpw */
*rv*/
, const mpw */
*rvl*/
, mpw
/*d*/);
726 /*----- Karatsuba multiplication algorithms -------------------------------*/
728 /* --- @MPK_THRESH@ --- *
730 * This is the limiting length for using Karatsuba algorithms. It's best to
731 * use the simpler classical multiplication method on numbers smaller than
732 * this. It is unsafe to make this constant less than four (i.e., the
733 * algorithms will fail).
736 #define MPK_THRESH 16
738 /* --- @mpx_kmul@ --- *
740 * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer
741 * @const mpw *av, *avl@ = pointer to first argument
742 * @const mpw *bv, *bvl@ = pointer to second argument
743 * @mpw *sv, *svl@ = pointer to scratch workspace
747 * Use: Multiplies two multiprecision integers using Karatsuba's
748 * algorithm. This is rather faster than traditional long
749 * multiplication (e.g., @mpx_umul@) on large numbers, although
750 * more expensive on small ones.
752 * The destination must be three times as large as the larger
753 * argument. The scratch space must be five times as large as
754 * the larger argument.
757 extern void mpx_kmul(mpw */
*dv*/
, mpw */
*dvl*/
,
758 const mpw */
*av*/
, const mpw */
*avl*/
,
759 const mpw */
*bv*/
, const mpw */
*bvl*/
,
760 mpw */
*sv*/
, mpw */
*svl*/
);
762 /* --- @mpx_ksqr@ --- *
764 * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer
765 * @const mpw *av, *avl@ = pointer to first argument
766 * @mpw *sv, *svl@ = pointer to scratch workspace
770 * Use: Squares a multiprecision integers using something similar to
771 * Karatsuba's multiplication algorithm. This is rather faster
772 * than traditional long multiplication (e.g., @mpx_umul@) on
773 * large numbers, although more expensive on small ones, and
774 * rather simpler than full-blown Karatsuba multiplication.
776 * The destination must be three times as large as the larger
777 * argument. The scratch space must be five times as large as
778 * the larger argument.
781 extern void mpx_ksqr(mpw */
*dv*/
, mpw */
*dvl*/
,
782 const mpw */
*av*/
, const mpw */
*avl*/
,
783 mpw */
*sv*/
, mpw */
*svl*/
);
785 /*----- That's all, folks -------------------------------------------------*/