3 * $Id: mpx.h,v 1.7 1999/12/11 01:51:28 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.7 1999/12/11 01:51:28 mdw
34 * Change Karatsuba parameters slightly.
36 * Revision 1.6 1999/12/10 23:23:51 mdw
37 * Karatsuba-Ofman multiplication algorithm.
39 * Revision 1.5 1999/11/20 22:23:27 mdw
40 * Add function versions of some low-level macros with wider use.
42 * Revision 1.4 1999/11/17 18:04:43 mdw
43 * Add two's complement support. Fix a bug in MPX_UMLAN.
45 * Revision 1.3 1999/11/13 01:51:29 mdw
46 * Minor interface changes. Should be stable now.
48 * Revision 1.2 1999/11/11 17:47:55 mdw
49 * Minor changes for different `mptypes.h' format.
51 * Revision 1.1 1999/09/03 08:41:12 mdw
56 #ifndef CATACOMB_MPX_H
57 #define CATACOMB_MPX_H
63 /*----- The idea ----------------------------------------------------------*
65 * This file provides functions and macros which work on vectors of words as
66 * unsigned multiprecision integers. The interface works in terms of base
67 * and limit pointers (i.e., a pointer to the start of a vector, and a
68 * pointer just past its end) rather than base pointer and length, because
69 * that requires more arithmetic and state to work on.
71 * The interfaces are slightly bizarre in other ways. Try to use the
72 * higher-level functions where you can: they're rather better designed to
73 * actually be friendly and useful.
76 /*----- Header files ------------------------------------------------------*/
80 #ifndef CATACOMB_MPW_H
84 /*----- General manipulation ----------------------------------------------*/
86 /* --- @MPX_SHRINK@ --- *
88 * Arguments: @const mpw *v@ = pointer to vector of words
89 * @const mpw *vl@ = (updated) current limit of vector
91 * Use: Shrinks down the limit of a multiprecision integer vector.
94 #define MPX_SHRINK(v, vl) do { \
95 const mpw *_vv = (v), *_vvl = (vl); \
96 while (_vvl > _vv && !_vvl[-1]) \
101 /* --- @MPX_BITS@ --- *
103 * Arguments: @unsigned long b@ = result variable
104 * @const mpw *v@ = pointer to array of words
105 * @const mpw *vl@ = limit of vector (from @MPX_SHRINK@)
107 * Use: Calculates the number of bits in a multiprecision value.
110 #define MPX_BITS(b, v, vl) do { \
111 const mpw *_v = (v), *_vl = (vl); \
112 MPX_SHRINK(_v, _vl); \
116 unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1; \
118 unsigned _k = MPW_BITS / 2; \
130 /* --- @MPX_OCTETS@ --- *
132 * Arguments: @size_t o@ = result variable
133 * @const mpw *v, *vl@ = pointer to array of words
135 * Use: Calculates the number of octets in a multiprecision value.
138 #define MPX_OCTETS(o, v, vl) do { \
139 const mpw *_v = (v), *_vl = (vl); \
140 MPX_SHRINK(_v, _vl); \
144 size_t _o = (MPW_BITS / 8) * (_vl - _v - 1); \
146 unsigned _k = MPW_BITS / 2; \
160 /* --- @MPX_COPY@ --- *
162 * Arguments: @dv, dvl@ = destination vector base and limit
163 * @av, avl@ = source vector base and limit
165 * Use: Copies a multiprecision integer.
168 #define MPX_COPY(dv, dvl, av, avl) do { \
169 mpw *_dv = (dv), *_dvl = (dvl); \
170 size_t _dn = _dvl - _dv; \
171 const mpw *_av = (av), *_avl = (avl); \
172 size_t _an = _avl - _av; \
175 memset(_dv, 0, MPWS(_dn - _an)); \
176 } else if (_an >= _dn) \
177 memmove(_dv, _av, MPWS(_dn)); \
179 memmove(_dv, _av, MPWS(_an)); \
180 memset(_dv + _an, 0, MPWS(_dn - _an)); \
184 /* --- @MPX_ZERO@ --- *
186 * Arguments: @v, vl@ = base and limit of vector to clear
188 * Use: Zeroes the area between the two vector pointers.
191 #define MPX_ZERO(v, vl) do { \
192 mpw *_v = (v), *_vl = (vl); \
194 memset(_v, 0, MPWS(_vl - _v)); \
197 /*----- Loading and storing -----------------------------------------------*/
199 /* --- @mpx_storel@ --- *
201 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
202 * @void *p@ = 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 extern void mpx_storel(const mpw */
*v*/
, const mpw */
*vl*/
,
213 void */
*p*/
, size_t /*sz*/);
215 /* --- @mpx_loadl@ --- *
217 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
218 * @const void *p@ = pointer to octet array
219 * @size_t sz@ = size of octet array
223 * Use: Loads an MP in an octet array, least significant octet
224 * first. High-end octets are ignored if there isn't enough
228 extern void mpx_loadl(mpw */
*v*/
, mpw */
*vl*/
,
229 const void */
*p*/
, size_t /*sz*/);
231 /* --- @mpx_storeb@ --- *
233 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
234 * @void *p@ = pointer to octet array
235 * @size_t sz@ = size of octet array
239 * Use: Stores an MP in an octet array, most significant octet
240 * first. High-end octets are silently discarded if there
241 * isn't enough space for them.
244 extern void mpx_storeb(const mpw */
*v*/
, const mpw */
*vl*/
,
245 void */
*p*/
, size_t /*sz*/);
247 /* --- @mpx_loadb@ --- *
249 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
250 * @const void *p@ = pointer to octet array
251 * @size_t sz@ = size of octet array
255 * Use: Loads an MP in an octet array, most significant octet
256 * first. High-end octets are ignored if there isn't enough
260 extern void mpx_loadb(mpw */
*v*/
, mpw */
*vl*/
,
261 const void */
*p*/
, size_t /*sz*/);
263 /*----- Logical shifting --------------------------------------------------*/
265 /* --- @mpx_lsl@ --- *
267 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
268 * @const mpw *av, *avl@ = source vector base and limit
269 * @size_t n@ = number of bit positions to shift by
273 * Use: Performs a logical shift left operation on an integer.
276 extern void mpx_lsl(mpw */
*dv*/
, mpw */
*dvl*/
,
277 const mpw */
*av*/
, const mpw */
*avl*/
,
280 /* --- @mpx_lsr@ --- *
282 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
283 * @const mpw *av, *avl@ = source vector base and limit
284 * @size_t n@ = number of bit positions to shift by
288 * Use: Performs a logical shift right operation on an integer.
291 extern void mpx_lsr(mpw */
*dv*/
, mpw */
*dvl*/
,
292 const mpw */
*av*/
, const mpw */
*avl*/
,
295 /*----- Unsigned arithmetic -----------------------------------------------*/
297 /* --- @mpx_2c@ --- *
299 * Arguments: @mpw *dv, *dvl@ = destination vector
300 * @const mpw *v, *vl@ = source vector
304 * Use: Calculates the two's complement of @v@.
307 extern void mpx_2c(mpw */
*dv*/
, mpw */
*dvl*/
,
308 const mpw */
*v*/
, const mpw */
*vl*/
);
310 /* --- @mpx_ucmp@ --- *
312 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
313 * @const mpw *bv, *bvl@ = second argument vector base and limit
315 * Returns: Less than, equal to, or greater than zero depending on
316 * whether @a@ is less than, equal to or greater than @b@,
319 * Use: Performs an unsigned integer comparison.
322 #define MPX_UCMP(av, avl, op, dv, dvl) \
323 (mpx_ucmp((av), (avl), (dv), (dvl)) op 0)
325 extern int mpx_ucmp(const mpw */
*av*/
, const mpw */
*avl*/
,
326 const mpw */
*bv*/
, const mpw */
*bvl*/
);
328 /* --- @mpx_uadd@ --- *
330 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
331 * @const mpw *av, *avl@ = first addend vector base and limit
332 * @const mpw *bv, *bvl@ = second addend vector base and limit
336 * Use: Performs unsigned integer addition. If the result overflows
337 * the destination vector, high-order bits are discarded. This
338 * means that two's complement addition happens more or less for
339 * free, although that's more a side-effect than anything else.
340 * The result vector may be equal to either or both source
341 * vectors, but may not otherwise overlap them.
344 extern void mpx_uadd(mpw */
*dv*/
, mpw */
*dvl*/
,
345 const mpw */
*av*/
, const mpw */
*avl*/
,
346 const mpw */
*bv*/
, const mpw */
*bvl*/
);
348 /* --- @mpx_uaddn@ --- *
350 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
351 * @mpw n@ = other addend
355 * Use: Adds a small integer to a multiprecision number.
358 #define MPX_UADDN(dv, dvl, n) do { \
359 mpw *_ddv = (dv), *_ddvl = (dvl); \
362 while (_c && _ddv < _ddvl) { \
363 mpd _x = (mpd)*_ddv + (mpd)_c; \
365 _c = _x >> MPW_BITS; \
369 extern void mpx_uaddn(mpw */
*dv*/
, mpw */
*dvl*/
, mpw
/*n*/);
371 /* --- @mpx_usub@ --- *
373 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
374 * @const mpw *av, *avl@ = first argument vector base and limit
375 * @const mpw *bv, *bvl@ = second argument vector base and limit
379 * Use: Performs unsigned integer subtraction. If the result
380 * overflows the destination vector, high-order bits are
381 * discarded. This means that two's complement subtraction
382 * happens more or less for free, although that's more a side-
383 * effect than anything else. The result vector may be equal to
384 * either or both source vectors, but may not otherwise overlap
388 extern void mpx_usub(mpw */
*dv*/
, mpw */
*dvl*/
,
389 const mpw */
*av*/
, const mpw */
*avl*/
,
390 const mpw */
*bv*/
, const mpw */
*bvl*/
);
392 /* --- @mpx_usubn@ --- *
394 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
399 * Use: Subtracts a small integer from a multiprecision number.
402 #define MPX_USUBN(dv, dvl, n) do { \
403 mpw *_ddv = (dv), *_ddvl = (dvl); \
406 while (_ddv < _ddvl) { \
407 mpd _x = (mpd)*_ddv - (mpd)_c; \
409 if (_x >> MPW_BITS) \
416 extern void mpx_usubn(mpw */
*dv*/
, mpw */
*dvl*/
, mpw
/*n*/);
418 /* --- @mpx_umul@ --- *
420 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
421 * @const mpw *av, *avl@ = multiplicand vector base and limit
422 * @const mpw *bv, *bvl@ = multiplier vector base and limit
426 * Use: Performs unsigned integer multiplication. If the result
427 * overflows the desination vector, high-order bits are
428 * discarded. The result vector may not overlap the argument
429 * vectors in any way.
432 extern void mpx_umul(mpw */
*dv*/
, mpw */
*dvl*/
,
433 const mpw */
*av*/
, const mpw */
*avl*/
,
434 const mpw */
*bv*/
, const mpw */
*bvl*/
);
436 /* --- @mpx_umuln@ --- *
438 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
439 * @const mpw *av, *avl@ = multiplicand vector base and limit
440 * @mpw m@ = multiplier
444 * Use: Multiplies a multiprecision integer by a single-word value.
445 * The destination and source may be equal. The destination
446 * is completely cleared after use.
449 #define MPX_UMULN(dv, dvl, av, avl, m) do { \
450 mpw *_dv = (dv), *_dvl = (dvl); \
451 const mpw *_av = (av), *_avl = (avl); \
455 while (_av < _avl) { \
459 _x = (mpd)_m * (mpd)*_av++ + _c; \
461 _c = _x >> MPW_BITS; \
465 MPX_ZERO(_dv, _dvl); \
469 extern void mpx_umuln(mpw */
*dv*/
, mpw */
*dvl*/
,
470 const mpw */
*av*/
, const mpw */
*avl*/
, mpw m
);
472 /* --- @mpx_umlan@ --- *
474 * Arguments: @mpw *dv, *dvl@ = destination/accumulator base and limit
475 * @const mpw *av, *avl@ = multiplicand vector base and limit
476 * @mpw m@ = multiplier
480 * Use: Multiplies a multiprecision integer by a single-word value
481 * and adds the result to an accumulator.
484 #define MPX_UMLAN(dv, dvl, av, avl, m) do { \
485 mpw *_dv = (dv), *_dvl = (dvl); \
486 const mpw *_av = (av), *_avl = (avl); \
490 while (_av < _avl) { \
494 _x = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _cc; \
496 _cc = _x >> MPW_BITS; \
498 MPX_UADDN(_dv, _dvl, _cc); \
501 extern void mpx_umlan(mpw */
*dv*/
, mpw */
*dvl*/
,
502 const mpw */
*av*/
, const mpw */
*avl*/
, mpw m
);
504 /* --- @mpx_usqr@ --- *
506 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
507 * @const mpw *av, *av@ = source vector base and limit
511 * Use: Performs unsigned integer squaring. The result vector must
512 * not overlap the source vector in any way.
515 extern void mpx_usqr(mpw */
*dv*/
, mpw */
*dvl*/
,
516 const mpw */
*av*/
, const mpw */
*avl*/
);
518 /* --- @mpx_kmul@ --- *
520 * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer
521 * @const mpw *av, *avl@ = pointer to first argument
522 * @const mpw *bv, *bvl@ = pointer to second argument
523 * @mpw *sv, *svl@ = pointer to scratch workspace
527 * Use: Multiplies two multiprecision integers using Karatsuba's
528 * algorithm. This is rather faster than traditional long
529 * multiplication (e.g., @mpx_umul@) on large numbers, although
530 * more expensive on small ones.
532 * The destination and scratch buffers must be twice as large as
533 * the larger argument.
536 #define KARATSUBA_CUTOFF 20
537 #define KARATSUBA_SLOP 32
539 extern void mpx_kmul(mpw */
*dv*/
, mpw */
*dvl*/
,
540 const mpw */
*av*/
, const mpw */
*avl*/
,
541 const mpw */
*bv*/
, const mpw */
*bvl*/
,
542 mpw */
*sv*/
, mpw */
*svl*/
);
544 /* --- @mpx_udiv@ --- *
546 * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit
547 * @mpw *rv, *rvl@ = dividend/remainder vector base and limit
548 * @const mpw *dv, *dvl@ = divisor vector base and limit
549 * @mpw *sv, *svl@ = scratch workspace
553 * Use: Performs unsigned integer division. If the result overflows
554 * the quotient vector, high-order bits are discarded. (Clearly
555 * the remainder vector can't overflow.) The various vectors
556 * may not overlap in any way. Yes, I know it's a bit odd
557 * requiring the dividend to be in the result position but it
558 * does make some sense really. The remainder must have
559 * headroom for at least two extra words. The scratch space
560 * must be at least one word larger than the divisor.
563 extern void mpx_udiv(mpw */
*qv*/
, mpw */
*qvl*/
, mpw */
*rv*/
, mpw */
*rvl*/
,
564 const mpw */
*dv*/
, const mpw */
*dvl*/
,
565 mpw */
*sv*/
, mpw */
*svl*/
);
567 /*----- That's all, folks -------------------------------------------------*/