3 * $Id: mpx.h,v 1.8 1999/12/11 10:57:43 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 /*----- Revision history --------------------------------------------------*
33 * Revision 1.8 1999/12/11 10:57:43 mdw
34 * Karatsuba squaring algorithm.
36 * Revision 1.7 1999/12/11 01:51:28 mdw
37 * Change Karatsuba parameters slightly.
39 * Revision 1.6 1999/12/10 23:23:51 mdw
40 * Karatsuba-Ofman multiplication algorithm.
42 * Revision 1.5 1999/11/20 22:23:27 mdw
43 * Add function versions of some low-level macros with wider use.
45 * Revision 1.4 1999/11/17 18:04:43 mdw
46 * Add two's complement support. Fix a bug in MPX_UMLAN.
48 * Revision 1.3 1999/11/13 01:51:29 mdw
49 * Minor interface changes. Should be stable now.
51 * Revision 1.2 1999/11/11 17:47:55 mdw
52 * Minor changes for different `mptypes.h' format.
54 * Revision 1.1 1999/09/03 08:41:12 mdw
59 #ifndef CATACOMB_MPX_H
60 #define CATACOMB_MPX_H
66 /*----- The idea ----------------------------------------------------------*
68 * This file provides functions and macros which work on vectors of words as
69 * unsigned multiprecision integers. The interface works in terms of base
70 * and limit pointers (i.e., a pointer to the start of a vector, and a
71 * pointer just past its end) rather than base pointer and length, because
72 * that requires more arithmetic and state to work on.
74 * The interfaces are slightly bizarre in other ways. Try to use the
75 * higher-level functions where you can: they're rather better designed to
76 * actually be friendly and useful.
79 /*----- Header files ------------------------------------------------------*/
83 #ifndef CATACOMB_MPW_H
87 /*----- General manipulation ----------------------------------------------*/
89 /* --- @MPX_SHRINK@ --- *
91 * Arguments: @const mpw *v@ = pointer to vector of words
92 * @const mpw *vl@ = (updated) current limit of vector
94 * Use: Shrinks down the limit of a multiprecision integer vector.
97 #define MPX_SHRINK(v, vl) do { \
98 const mpw *_vv = (v), *_vvl = (vl); \
99 while (_vvl > _vv && !_vvl[-1]) \
101 (vl) = (mpw *)_vvl; \
104 /* --- @MPX_BITS@ --- *
106 * Arguments: @unsigned long b@ = result variable
107 * @const mpw *v@ = pointer to array of words
108 * @const mpw *vl@ = limit of vector (from @MPX_SHRINK@)
110 * Use: Calculates the number of bits in a multiprecision value.
113 #define MPX_BITS(b, v, vl) do { \
114 const mpw *_v = (v), *_vl = (vl); \
115 MPX_SHRINK(_v, _vl); \
119 unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1; \
121 unsigned _k = MPW_BITS / 2; \
133 /* --- @MPX_OCTETS@ --- *
135 * Arguments: @size_t o@ = result variable
136 * @const mpw *v, *vl@ = pointer to array of words
138 * Use: Calculates the number of octets in a multiprecision value.
141 #define MPX_OCTETS(o, v, vl) do { \
142 const mpw *_v = (v), *_vl = (vl); \
143 MPX_SHRINK(_v, _vl); \
147 size_t _o = (MPW_BITS / 8) * (_vl - _v - 1); \
149 unsigned _k = MPW_BITS / 2; \
163 /* --- @MPX_COPY@ --- *
165 * Arguments: @dv, dvl@ = destination vector base and limit
166 * @av, avl@ = source vector base and limit
168 * Use: Copies a multiprecision integer.
171 #define MPX_COPY(dv, dvl, av, avl) do { \
172 mpw *_dv = (dv), *_dvl = (dvl); \
173 size_t _dn = _dvl - _dv; \
174 const mpw *_av = (av), *_avl = (avl); \
175 size_t _an = _avl - _av; \
178 memset(_dv, 0, MPWS(_dn - _an)); \
179 } else if (_an >= _dn) \
180 memmove(_dv, _av, MPWS(_dn)); \
182 memmove(_dv, _av, MPWS(_an)); \
183 memset(_dv + _an, 0, MPWS(_dn - _an)); \
187 /* --- @MPX_ZERO@ --- *
189 * Arguments: @v, vl@ = base and limit of vector to clear
191 * Use: Zeroes the area between the two vector pointers.
194 #define MPX_ZERO(v, vl) do { \
195 mpw *_v = (v), *_vl = (vl); \
197 memset(_v, 0, MPWS(_vl - _v)); \
200 /*----- Loading and storing -----------------------------------------------*/
202 /* --- @mpx_storel@ --- *
204 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
205 * @void *p@ = pointer to octet array
206 * @size_t sz@ = size of octet array
210 * Use: Stores an MP in an octet array, least significant octet
211 * first. High-end octets are silently discarded if there
212 * isn't enough space for them.
215 extern void mpx_storel(const mpw */
*v*/
, const mpw */
*vl*/
,
216 void */
*p*/
, size_t /*sz*/);
218 /* --- @mpx_loadl@ --- *
220 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
221 * @const void *p@ = pointer to octet array
222 * @size_t sz@ = size of octet array
226 * Use: Loads an MP in an octet array, least significant octet
227 * first. High-end octets are ignored if there isn't enough
231 extern void mpx_loadl(mpw */
*v*/
, mpw */
*vl*/
,
232 const void */
*p*/
, size_t /*sz*/);
234 /* --- @mpx_storeb@ --- *
236 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
237 * @void *p@ = pointer to octet array
238 * @size_t sz@ = size of octet array
242 * Use: Stores an MP in an octet array, most significant octet
243 * first. High-end octets are silently discarded if there
244 * isn't enough space for them.
247 extern void mpx_storeb(const mpw */
*v*/
, const mpw */
*vl*/
,
248 void */
*p*/
, size_t /*sz*/);
250 /* --- @mpx_loadb@ --- *
252 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
253 * @const void *p@ = pointer to octet array
254 * @size_t sz@ = size of octet array
258 * Use: Loads an MP in an octet array, most significant octet
259 * first. High-end octets are ignored if there isn't enough
263 extern void mpx_loadb(mpw */
*v*/
, mpw */
*vl*/
,
264 const void */
*p*/
, size_t /*sz*/);
266 /*----- Logical shifting --------------------------------------------------*/
268 /* --- @mpx_lsl@ --- *
270 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
271 * @const mpw *av, *avl@ = source vector base and limit
272 * @size_t n@ = number of bit positions to shift by
276 * Use: Performs a logical shift left operation on an integer.
279 extern void mpx_lsl(mpw */
*dv*/
, mpw */
*dvl*/
,
280 const mpw */
*av*/
, const mpw */
*avl*/
,
283 /* --- @mpx_lsr@ --- *
285 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
286 * @const mpw *av, *avl@ = source vector base and limit
287 * @size_t n@ = number of bit positions to shift by
291 * Use: Performs a logical shift right operation on an integer.
294 extern void mpx_lsr(mpw */
*dv*/
, mpw */
*dvl*/
,
295 const mpw */
*av*/
, const mpw */
*avl*/
,
298 /*----- Unsigned arithmetic -----------------------------------------------*/
300 /* --- @mpx_2c@ --- *
302 * Arguments: @mpw *dv, *dvl@ = destination vector
303 * @const mpw *v, *vl@ = source vector
307 * Use: Calculates the two's complement of @v@.
310 extern void mpx_2c(mpw */
*dv*/
, mpw */
*dvl*/
,
311 const mpw */
*v*/
, const mpw */
*vl*/
);
313 /* --- @mpx_ucmp@ --- *
315 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
316 * @const mpw *bv, *bvl@ = second argument vector base and limit
318 * Returns: Less than, equal to, or greater than zero depending on
319 * whether @a@ is less than, equal to or greater than @b@,
322 * Use: Performs an unsigned integer comparison.
325 #define MPX_UCMP(av, avl, op, dv, dvl) \
326 (mpx_ucmp((av), (avl), (dv), (dvl)) op 0)
328 extern int mpx_ucmp(const mpw */
*av*/
, const mpw */
*avl*/
,
329 const mpw */
*bv*/
, const mpw */
*bvl*/
);
331 /* --- @mpx_uadd@ --- *
333 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
334 * @const mpw *av, *avl@ = first addend vector base and limit
335 * @const mpw *bv, *bvl@ = second addend vector base and limit
339 * Use: Performs unsigned integer addition. If the result overflows
340 * the destination vector, high-order bits are discarded. This
341 * means that two's complement addition happens more or less for
342 * free, although that's more a side-effect than anything else.
343 * The result vector may be equal to either or both source
344 * vectors, but may not otherwise overlap them.
347 extern void mpx_uadd(mpw */
*dv*/
, mpw */
*dvl*/
,
348 const mpw */
*av*/
, const mpw */
*avl*/
,
349 const mpw */
*bv*/
, const mpw */
*bvl*/
);
351 /* --- @mpx_uaddn@ --- *
353 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
354 * @mpw n@ = other addend
358 * Use: Adds a small integer to a multiprecision number.
361 #define MPX_UADDN(dv, dvl, n) do { \
362 mpw *_ddv = (dv), *_ddvl = (dvl); \
365 while (_c && _ddv < _ddvl) { \
366 mpd _x = (mpd)*_ddv + (mpd)_c; \
368 _c = _x >> MPW_BITS; \
372 extern void mpx_uaddn(mpw */
*dv*/
, mpw */
*dvl*/
, mpw
/*n*/);
374 /* --- @mpx_usub@ --- *
376 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
377 * @const mpw *av, *avl@ = first argument vector base and limit
378 * @const mpw *bv, *bvl@ = second argument vector base and limit
382 * Use: Performs unsigned integer subtraction. If the result
383 * overflows the destination vector, high-order bits are
384 * discarded. This means that two's complement subtraction
385 * happens more or less for free, although that's more a side-
386 * effect than anything else. The result vector may be equal to
387 * either or both source vectors, but may not otherwise overlap
391 extern void mpx_usub(mpw */
*dv*/
, mpw */
*dvl*/
,
392 const mpw */
*av*/
, const mpw */
*avl*/
,
393 const mpw */
*bv*/
, const mpw */
*bvl*/
);
395 /* --- @mpx_usubn@ --- *
397 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
402 * Use: Subtracts a small integer from a multiprecision number.
405 #define MPX_USUBN(dv, dvl, n) do { \
406 mpw *_ddv = (dv), *_ddvl = (dvl); \
409 while (_ddv < _ddvl) { \
410 mpd _x = (mpd)*_ddv - (mpd)_c; \
412 if (_x >> MPW_BITS) \
419 extern void mpx_usubn(mpw */
*dv*/
, mpw */
*dvl*/
, mpw
/*n*/);
421 /* --- @mpx_umul@ --- *
423 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
424 * @const mpw *av, *avl@ = multiplicand vector base and limit
425 * @const mpw *bv, *bvl@ = multiplier vector base and limit
429 * Use: Performs unsigned integer multiplication. If the result
430 * overflows the desination vector, high-order bits are
431 * discarded. The result vector may not overlap the argument
432 * vectors in any way.
435 extern void mpx_umul(mpw */
*dv*/
, mpw */
*dvl*/
,
436 const mpw */
*av*/
, const mpw */
*avl*/
,
437 const mpw */
*bv*/
, const mpw */
*bvl*/
);
439 /* --- @mpx_umuln@ --- *
441 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
442 * @const mpw *av, *avl@ = multiplicand vector base and limit
443 * @mpw m@ = multiplier
447 * Use: Multiplies a multiprecision integer by a single-word value.
448 * The destination and source may be equal. The destination
449 * is completely cleared after use.
452 #define MPX_UMULN(dv, dvl, av, avl, m) do { \
453 mpw *_dv = (dv), *_dvl = (dvl); \
454 const mpw *_av = (av), *_avl = (avl); \
458 while (_av < _avl) { \
462 _x = (mpd)_m * (mpd)*_av++ + _c; \
464 _c = _x >> MPW_BITS; \
468 MPX_ZERO(_dv, _dvl); \
472 extern void mpx_umuln(mpw */
*dv*/
, mpw */
*dvl*/
,
473 const mpw */
*av*/
, const mpw */
*avl*/
, mpw m
);
475 /* --- @mpx_umlan@ --- *
477 * Arguments: @mpw *dv, *dvl@ = destination/accumulator base and limit
478 * @const mpw *av, *avl@ = multiplicand vector base and limit
479 * @mpw m@ = multiplier
483 * Use: Multiplies a multiprecision integer by a single-word value
484 * and adds the result to an accumulator.
487 #define MPX_UMLAN(dv, dvl, av, avl, m) do { \
488 mpw *_dv = (dv), *_dvl = (dvl); \
489 const mpw *_av = (av), *_avl = (avl); \
493 while (_dv < _dvl && _av < _avl) { \
495 _x = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _cc; \
497 _cc = _x >> MPW_BITS; \
499 MPX_UADDN(_dv, _dvl, _cc); \
502 extern void mpx_umlan(mpw */
*dv*/
, mpw */
*dvl*/
,
503 const mpw */
*av*/
, const mpw */
*avl*/
, mpw m
);
505 /* --- @mpx_usqr@ --- *
507 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
508 * @const mpw *av, *av@ = source vector base and limit
512 * Use: Performs unsigned integer squaring. The result vector must
513 * not overlap the source vector in any way.
516 extern void mpx_usqr(mpw */
*dv*/
, mpw */
*dvl*/
,
517 const mpw */
*av*/
, const mpw */
*avl*/
);
519 /* --- @mpx_udiv@ --- *
521 * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit
522 * @mpw *rv, *rvl@ = dividend/remainder vector base and limit
523 * @const mpw *dv, *dvl@ = divisor vector base and limit
524 * @mpw *sv, *svl@ = scratch workspace
528 * Use: Performs unsigned integer division. If the result overflows
529 * the quotient vector, high-order bits are discarded. (Clearly
530 * the remainder vector can't overflow.) The various vectors
531 * may not overlap in any way. Yes, I know it's a bit odd
532 * requiring the dividend to be in the result position but it
533 * does make some sense really. The remainder must have
534 * headroom for at least two extra words. The scratch space
535 * must be at least one word larger than the divisor.
538 extern void mpx_udiv(mpw */
*qv*/
, mpw */
*qvl*/
, mpw */
*rv*/
, mpw */
*rvl*/
,
539 const mpw */
*dv*/
, const mpw */
*dvl*/
,
540 mpw */
*sv*/
, mpw */
*svl*/
);
542 /*----- Karatsuba multiplication algorithms -------------------------------*/
544 /* --- @KARATSUBA_CUTOFF@ --- *
546 * This is the limiting length for using Karatsuba algorithms. It's best to
547 * use the simpler classical multiplication method on numbers smaller than
551 #define KARATSUBA_CUTOFF 16
553 /* --- @KARATSUBA_SLOP@ --- *
555 * The extra number of words required as scratch space by the Karatsuba
556 * routines. This is a (generous) guess, since the actual amount of space
557 * required is proportional to the recursion depth.
560 #define KARATSUBA_SLOP 32
562 /* --- @mpx_kmul@ --- *
564 * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer
565 * @const mpw *av, *avl@ = pointer to first argument
566 * @const mpw *bv, *bvl@ = pointer to second argument
567 * @mpw *sv, *svl@ = pointer to scratch workspace
571 * Use: Multiplies two multiprecision integers using Karatsuba's
572 * algorithm. This is rather faster than traditional long
573 * multiplication (e.g., @mpx_umul@) on large numbers, although
574 * more expensive on small ones.
576 * The destination and scratch buffers must be twice as large as
577 * the larger argument. The scratch space must be twice as
578 * large as the larger argument, plus the magic number
582 extern void mpx_kmul(mpw */
*dv*/
, mpw */
*dvl*/
,
583 const mpw */
*av*/
, const mpw */
*avl*/
,
584 const mpw */
*bv*/
, const mpw */
*bvl*/
,
585 mpw */
*sv*/
, mpw */
*svl*/
);
587 /* --- @mpx_ksqr@ --- *
589 * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer
590 * @const mpw *av, *avl@ = pointer to first argument
591 * @mpw *sv, *svl@ = pointer to scratch workspace
595 * Use: Squares a multiprecision integers using something similar to
596 * Karatsuba's multiplication algorithm. This is rather faster
597 * than traditional long multiplication (e.g., @mpx_umul@) on
598 * large numbers, although more expensive on small ones, and
599 * rather simpler than full-blown Karatsuba multiplication.
601 * The destination must be twice as large as the argument. The
602 * scratch space must be twice as large as the argument, plus
603 * the magic number @KARATSUBA_SLOP@.
606 extern void mpx_ksqr(mpw */
*dv*/
, mpw */
*dvl*/
,
607 const mpw */
*av*/
, const mpw */
*avl*/
,
608 mpw */
*sv*/
, mpw */
*svl*/
);
610 /*----- That's all, folks -------------------------------------------------*/