projects / u / mdw / catacomb / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
Provide @mpx_ueq@ for rapidly testing equality of two integers.
[u/mdw/catacomb] / mpx.h
1 /* -*-c-*-
2 *
3 * $Id: mpx.h,v 1.10 2000/10/08 12:06:12 mdw Exp $
4 *
5 * Low level multiprecision arithmetic
6 *
7 * (c) 1999 Straylight/Edgeware
8 */
9
10 /*----- Licensing notice --------------------------------------------------*
11 *
12 * This file is part of Catacomb.
13 *
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.
18 *
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.
23 *
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,
27 * MA 02111-1307, USA.
28 */
29
30 /*----- Revision history --------------------------------------------------*
31 *
32 * $Log: mpx.h,v $
33 * Revision 1.10 2000/10/08 12:06:12 mdw
34 * Provide @mpx_ueq@ for rapidly testing equality of two integers.
35 *
36 * Revision 1.9 1999/12/22 15:49:07 mdw
37 * New function for division by a small integer.
38 *
39 * Revision 1.8 1999/12/11 10:57:43 mdw
40 * Karatsuba squaring algorithm.
41 *
42 * Revision 1.7 1999/12/11 01:51:28 mdw
43 * Change Karatsuba parameters slightly.
44 *
45 * Revision 1.6 1999/12/10 23:23:51 mdw
46 * Karatsuba-Ofman multiplication algorithm.
47 *
48 * Revision 1.5 1999/11/20 22:23:27 mdw
49 * Add function versions of some low-level macros with wider use.
50 *
51 * Revision 1.4 1999/11/17 18:04:43 mdw
52 * Add two's complement support. Fix a bug in MPX_UMLAN.
53 *
54 * Revision 1.3 1999/11/13 01:51:29 mdw
55 * Minor interface changes. Should be stable now.
56 *
57 * Revision 1.2 1999/11/11 17:47:55 mdw
58 * Minor changes for different `mptypes.h' format.
59 *
60 * Revision 1.1 1999/09/03 08:41:12 mdw
61 * Initial import.
62 *
63 */
64
65 #ifndef CATACOMB_MPX_H
66 #define CATACOMB_MPX_H
67
68 #ifdef __cplusplus
69 extern "C" {
70 #endif
71
72 /*----- The idea ----------------------------------------------------------*
73 *
74 * This file provides functions and macros which work on vectors of words as
75 * unsigned multiprecision integers. The interface works in terms of base
76 * and limit pointers (i.e., a pointer to the start of a vector, and a
77 * pointer just past its end) rather than base pointer and length, because
78 * that requires more arithmetic and state to work on.
79 *
80 * The interfaces are slightly bizarre in other ways. Try to use the
81 * higher-level functions where you can: they're rather better designed to
82 * actually be friendly and useful.
83 */
84
85 /*----- Header files ------------------------------------------------------*/
86
87 #include <string.h>
88
89 #ifndef CATACOMB_MPW_H
90 # include "mpw.h"
91 #endif
92
93 /*----- General manipulation ----------------------------------------------*/
94
95 /* --- @MPX_SHRINK@ --- *
96 *
97 * Arguments: @const mpw *v@ = pointer to vector of words
98 * @const mpw *vl@ = (updated) current limit of vector
99 *
100 * Use: Shrinks down the limit of a multiprecision integer vector.
101 */
102
103 #define MPX_SHRINK(v, vl) do { \
104 const mpw *_vv = (v), *_vvl = (vl); \
105 while (_vvl > _vv && !_vvl[-1]) \
106 _vvl--; \
107 (vl) = (mpw *)_vvl; \
108 } while (0)
109
110 /* --- @MPX_BITS@ --- *
111 *
112 * Arguments: @unsigned long b@ = result variable
113 * @const mpw *v@ = pointer to array of words
114 * @const mpw *vl@ = limit of vector (from @MPX_SHRINK@)
115 *
116 * Use: Calculates the number of bits in a multiprecision value.
117 */
118
119 #define MPX_BITS(b, v, vl) do { \
120 const mpw *_v = (v), *_vl = (vl); \
121 MPX_SHRINK(_v, _vl); \
122 if (_v == _vl) \
123 (b) = 0; \
124 else { \
125 unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1; \
126 mpw _w = _vl[-1]; \
127 unsigned _k = MPW_BITS / 2; \
128 while (_k) { \
129 if (_w >> _k) { \
130 _w >>= _k; \
131 _b += _k; \
132 } \
133 _k >>= 1; \
134 } \
135 (b) = _b; \
136 } \
137 } while (0)
138
139 /* --- @MPX_OCTETS@ --- *
140 *
141 * Arguments: @size_t o@ = result variable
142 * @const mpw *v, *vl@ = pointer to array of words
143 *
144 * Use: Calculates the number of octets in a multiprecision value.
145 */
146
147 #define MPX_OCTETS(o, v, vl) do { \
148 const mpw *_v = (v), *_vl = (vl); \
149 MPX_SHRINK(_v, _vl); \
150 if (_v == _vl) \
151 (o) = 0; \
152 else { \
153 size_t _o = (MPW_BITS / 8) * (_vl - _v - 1); \
154 mpw _w = _vl[-1]; \
155 unsigned _k = MPW_BITS / 2; \
156 while (_k >= 8) { \
157 if (_w >> _k) { \
158 _w >>= _k; \
159 _o += _k >> 3; \
160 } \
161 _k >>= 1; \
162 } \
163 if (_w) \
164 _o++; \
165 (o) = _o; \
166 } \
167 } while (0)
168
169 /* --- @MPX_COPY@ --- *
170 *
171 * Arguments: @dv, dvl@ = destination vector base and limit
172 * @av, avl@ = source vector base and limit
173 *
174 * Use: Copies a multiprecision integer.
175 */
176
177 #define MPX_COPY(dv, dvl, av, avl) do { \
178 mpw *_dv = (dv), *_dvl = (dvl); \
179 size_t _dn = _dvl - _dv; \
180 const mpw *_av = (av), *_avl = (avl); \
181 size_t _an = _avl - _av; \
182 if (_av == _dv) { \
183 if (_dvl > _avl) \
184 memset(_dv, 0, MPWS(_dn - _an)); \
185 } else if (_an >= _dn) \
186 memmove(_dv, _av, MPWS(_dn)); \
187 else { \
188 memmove(_dv, _av, MPWS(_an)); \
189 memset(_dv + _an, 0, MPWS(_dn - _an)); \
190 } \
191 } while (0)
192
193 /* --- @MPX_ZERO@ --- *
194 *
195 * Arguments: @v, vl@ = base and limit of vector to clear
196 *
197 * Use: Zeroes the area between the two vector pointers.
198 */
199
200 #define MPX_ZERO(v, vl) do { \
201 mpw *_v = (v), *_vl = (vl); \
202 if (_v < _vl) \
203 memset(_v, 0, MPWS(_vl - _v)); \
204 } while (0)
205
206 /*----- Loading and storing -----------------------------------------------*/
207
208 /* --- @mpx_storel@ --- *
209 *
210 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
211 * @void *p@ = pointer to octet array
212 * @size_t sz@ = size of octet array
213 *
214 * Returns: ---
215 *
216 * Use: Stores an MP in an octet array, least significant octet
217 * first. High-end octets are silently discarded if there
218 * isn't enough space for them.
219 */
220
221 extern void mpx_storel(const mpw */*v*/, const mpw */*vl*/,
222 void */*p*/, size_t /*sz*/);
223
224 /* --- @mpx_loadl@ --- *
225 *
226 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
227 * @const void *p@ = pointer to octet array
228 * @size_t sz@ = size of octet array
229 *
230 * Returns: ---
231 *
232 * Use: Loads an MP in an octet array, least significant octet
233 * first. High-end octets are ignored if there isn't enough
234 * space for them.
235 */
236
237 extern void mpx_loadl(mpw */*v*/, mpw */*vl*/,
238 const void */*p*/, size_t /*sz*/);
239
240 /* --- @mpx_storeb@ --- *
241 *
242 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
243 * @void *p@ = pointer to octet array
244 * @size_t sz@ = size of octet array
245 *
246 * Returns: ---
247 *
248 * Use: Stores an MP in an octet array, most significant octet
249 * first. High-end octets are silently discarded if there
250 * isn't enough space for them.
251 */
252
253 extern void mpx_storeb(const mpw */*v*/, const mpw */*vl*/,
254 void */*p*/, size_t /*sz*/);
255
256 /* --- @mpx_loadb@ --- *
257 *
258 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
259 * @const void *p@ = pointer to octet array
260 * @size_t sz@ = size of octet array
261 *
262 * Returns: ---
263 *
264 * Use: Loads an MP in an octet array, most significant octet
265 * first. High-end octets are ignored if there isn't enough
266 * space for them.
267 */
268
269 extern void mpx_loadb(mpw */*v*/, mpw */*vl*/,
270 const void */*p*/, size_t /*sz*/);
271
272 /*----- Logical shifting --------------------------------------------------*/
273
274 /* --- @mpx_lsl@ --- *
275 *
276 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
277 * @const mpw *av, *avl@ = source vector base and limit
278 * @size_t n@ = number of bit positions to shift by
279 *
280 * Returns: ---
281 *
282 * Use: Performs a logical shift left operation on an integer.
283 */
284
285 extern void mpx_lsl(mpw */*dv*/, mpw */*dvl*/,
286 const mpw */*av*/, const mpw */*avl*/,
287 size_t /*n*/);
288
289 /* --- @mpx_lsr@ --- *
290 *
291 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
292 * @const mpw *av, *avl@ = source vector base and limit
293 * @size_t n@ = number of bit positions to shift by
294 *
295 * Returns: ---
296 *
297 * Use: Performs a logical shift right operation on an integer.
298 */
299
300 extern void mpx_lsr(mpw */*dv*/, mpw */*dvl*/,
301 const mpw */*av*/, const mpw */*avl*/,
302 size_t /*n*/);
303
304 /*----- Unsigned arithmetic -----------------------------------------------*/
305
306 /* --- @mpx_2c@ --- *
307 *
308 * Arguments: @mpw *dv, *dvl@ = destination vector
309 * @const mpw *v, *vl@ = source vector
310 *
311 * Returns: ---
312 *
313 * Use: Calculates the two's complement of @v@.
314 */
315
316 extern void mpx_2c(mpw */*dv*/, mpw */*dvl*/,
317 const mpw */*v*/, const mpw */*vl*/);
318
319 /* --- @mpx_ueq@ --- *
320 *
321 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
322 * @const mpw *bv, *bvl@ = second argument vector base and limit
323 *
324 * Returns: Nonzero if the two vectors are equal.
325 *
326 * Use: Performs an unsigned integer test for equality.
327 */
328
329 extern int mpx_ueq(const mpw */*av*/, const mpw */*avl*/,
330 const mpw */*bv*/, const mpw */*bvl*/);
331
332 /* --- @mpx_ucmp@ --- *
333 *
334 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
335 * @const mpw *bv, *bvl@ = second argument vector base and limit
336 *
337 * Returns: Less than, equal to, or greater than zero depending on
338 * whether @a@ is less than, equal to or greater than @b@,
339 * respectively.
340 *
341 * Use: Performs an unsigned integer comparison.
342 */
343
344 #define MPX_UCMP(av, avl, op, dv, dvl) \
345 (mpx_ucmp((av), (avl), (dv), (dvl)) op 0)
346
347 extern int mpx_ucmp(const mpw */*av*/, const mpw */*avl*/,
348 const mpw */*bv*/, const mpw */*bvl*/);
349
350 /* --- @mpx_uadd@ --- *
351 *
352 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
353 * @const mpw *av, *avl@ = first addend vector base and limit
354 * @const mpw *bv, *bvl@ = second addend vector base and limit
355 *
356 * Returns: ---
357 *
358 * Use: Performs unsigned integer addition. If the result overflows
359 * the destination vector, high-order bits are discarded. This
360 * means that two's complement addition happens more or less for
361 * free, although that's more a side-effect than anything else.
362 * The result vector may be equal to either or both source
363 * vectors, but may not otherwise overlap them.
364 */
365
366 extern void mpx_uadd(mpw */*dv*/, mpw */*dvl*/,
367 const mpw */*av*/, const mpw */*avl*/,
368 const mpw */*bv*/, const mpw */*bvl*/);
369
370 /* --- @mpx_uaddn@ --- *
371 *
372 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
373 * @mpw n@ = other addend
374 *
375 * Returns: ---
376 *
377 * Use: Adds a small integer to a multiprecision number.
378 */
379
380 #define MPX_UADDN(dv, dvl, n) do { \
381 mpw *_ddv = (dv), *_ddvl = (dvl); \
382 mpw _c = (n); \
383 \
384 while (_c && _ddv < _ddvl) { \
385 mpd _x = (mpd)*_ddv + (mpd)_c; \
386 *_ddv++ = MPW(_x); \
387 _c = _x >> MPW_BITS; \
388 } \
389 } while (0)
390
391 extern void mpx_uaddn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);
392
393 /* --- @mpx_usub@ --- *
394 *
395 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
396 * @const mpw *av, *avl@ = first argument vector base and limit
397 * @const mpw *bv, *bvl@ = second argument vector base and limit
398 *
399 * Returns: ---
400 *
401 * Use: Performs unsigned integer subtraction. If the result
402 * overflows the destination vector, high-order bits are
403 * discarded. This means that two's complement subtraction
404 * happens more or less for free, although that's more a side-
405 * effect than anything else. The result vector may be equal to
406 * either or both source vectors, but may not otherwise overlap
407 * them.
408 */
409
410 extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/,
411 const mpw */*av*/, const mpw */*avl*/,
412 const mpw */*bv*/, const mpw */*bvl*/);
413
414 /* --- @mpx_usubn@ --- *
415 *
416 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
417 * @n@ = subtrahend
418 *
419 * Returns: ---
420 *
421 * Use: Subtracts a small integer from a multiprecision number.
422 */
423
424 #define MPX_USUBN(dv, dvl, n) do { \
425 mpw *_ddv = (dv), *_ddvl = (dvl); \
426 mpw _c = (n); \
427 \
428 while (_ddv < _ddvl) { \
429 mpd _x = (mpd)*_ddv - (mpd)_c; \
430 *_ddv++ = MPW(_x); \
431 if (_x >> MPW_BITS) \
432 _c = 1; \
433 else \
434 break; \
435 } \
436 } while (0)
437
438 extern void mpx_usubn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);
439
440 /* --- @mpx_umul@ --- *
441 *
442 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
443 * @const mpw *av, *avl@ = multiplicand vector base and limit
444 * @const mpw *bv, *bvl@ = multiplier vector base and limit
445 *
446 * Returns: ---
447 *
448 * Use: Performs unsigned integer multiplication. If the result
449 * overflows the desination vector, high-order bits are
450 * discarded. The result vector may not overlap the argument
451 * vectors in any way.
452 */
453
454 extern void mpx_umul(mpw */*dv*/, mpw */*dvl*/,
455 const mpw */*av*/, const mpw */*avl*/,
456 const mpw */*bv*/, const mpw */*bvl*/);
457
458 /* --- @mpx_umuln@ --- *
459 *
460 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
461 * @const mpw *av, *avl@ = multiplicand vector base and limit
462 * @mpw m@ = multiplier
463 *
464 * Returns: ---
465 *
466 * Use: Multiplies a multiprecision integer by a single-word value.
467 * The destination and source may be equal. The destination
468 * is completely cleared after use.
469 */
470
471 #define MPX_UMULN(dv, dvl, av, avl, m) do { \
472 mpw *_dv = (dv), *_dvl = (dvl); \
473 const mpw *_av = (av), *_avl = (avl); \
474 mpw _c = 0; \
475 mpd _m = (m); \
476 \
477 while (_av < _avl) { \
478 mpd _x; \
479 if (_dv >= _dvl) \
480 break; \
481 _x = (mpd)_m * (mpd)*_av++ + _c; \
482 *_dv++ = MPW(_x); \
483 _c = _x >> MPW_BITS; \
484 } \
485 if (_dv < _dvl) { \
486 *_dv++ = MPW(_c); \
487 MPX_ZERO(_dv, _dvl); \
488 } \
489 } while (0)
490
491 extern void mpx_umuln(mpw */*dv*/, mpw */*dvl*/,
492 const mpw */*av*/, const mpw */*avl*/, mpw m);
493
494 /* --- @mpx_umlan@ --- *
495 *
496 * Arguments: @mpw *dv, *dvl@ = destination/accumulator base and limit
497 * @const mpw *av, *avl@ = multiplicand vector base and limit
498 * @mpw m@ = multiplier
499 *
500 * Returns: ---
501 *
502 * Use: Multiplies a multiprecision integer by a single-word value
503 * and adds the result to an accumulator.
504 */
505
506 #define MPX_UMLAN(dv, dvl, av, avl, m) do { \
507 mpw *_dv = (dv), *_dvl = (dvl); \
508 const mpw *_av = (av), *_avl = (avl); \
509 mpw _cc = 0; \
510 mpd _m = (m); \
511 \
512 while (_dv < _dvl && _av < _avl) { \
513 mpd _x; \
514 _x = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _cc; \
515 *_dv++ = MPW(_x); \
516 _cc = _x >> MPW_BITS; \
517 } \
518 MPX_UADDN(_dv, _dvl, _cc); \
519 } while (0)
520
521 extern void mpx_umlan(mpw */*dv*/, mpw */*dvl*/,
522 const mpw */*av*/, const mpw */*avl*/, mpw m);
523
524 /* --- @mpx_usqr@ --- *
525 *
526 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
527 * @const mpw *av, *av@ = source vector base and limit
528 *
529 * Returns: ---
530 *
531 * Use: Performs unsigned integer squaring. The result vector must
532 * not overlap the source vector in any way.
533 */
534
535 extern void mpx_usqr(mpw */*dv*/, mpw */*dvl*/,
536 const mpw */*av*/, const mpw */*avl*/);
537
538 /* --- @mpx_udiv@ --- *
539 *
540 * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit
541 * @mpw *rv, *rvl@ = dividend/remainder vector base and limit
542 * @const mpw *dv, *dvl@ = divisor vector base and limit
543 * @mpw *sv, *svl@ = scratch workspace
544 *
545 * Returns: ---
546 *
547 * Use: Performs unsigned integer division. If the result overflows
548 * the quotient vector, high-order bits are discarded. (Clearly
549 * the remainder vector can't overflow.) The various vectors
550 * may not overlap in any way. Yes, I know it's a bit odd
551 * requiring the dividend to be in the result position but it
552 * does make some sense really. The remainder must have
553 * headroom for at least two extra words. The scratch space
554 * must be at least one word larger than the divisor.
555 */
556
557 extern void mpx_udiv(mpw */*qv*/, mpw */*qvl*/, mpw */*rv*/, mpw */*rvl*/,
558 const mpw */*dv*/, const mpw */*dvl*/,
559 mpw */*sv*/, mpw */*svl*/);
560
561 /* --- @mpx_udivn@ --- *
562 *
563 * Arguments: @mpw *qv, *qvl@ = storage for the quotient (may overlap
564 * dividend)
565 * @const mpw *rv, *rvl@ = dividend
566 * @mpw d@ = single-precision divisor
567 *
568 * Returns: Remainder after divison.
569 *
570 * Use: Performs a single-precision division operation.
571 */
572
573 extern mpw mpx_udivn(mpw */*qv*/, mpw */*qvl*/,
574 const mpw */*rv*/, const mpw */*rvl*/, mpw /*d*/);
575
576 /*----- Karatsuba multiplication algorithms -------------------------------*/
577
578 /* --- @KARATSUBA_CUTOFF@ --- *
579 *
580 * This is the limiting length for using Karatsuba algorithms. It's best to
581 * use the simpler classical multiplication method on numbers smaller than
582 * this.
583 */
584
585 #define KARATSUBA_CUTOFF 16
586
587 /* --- @KARATSUBA_SLOP@ --- *
588 *
589 * The extra number of words required as scratch space by the Karatsuba
590 * routines. This is a (generous) guess, since the actual amount of space
591 * required is proportional to the recursion depth.
592 */
593
594 #define KARATSUBA_SLOP 64
595
596 /* --- @mpx_kmul@ --- *
597 *
598 * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer
599 * @const mpw *av, *avl@ = pointer to first argument
600 * @const mpw *bv, *bvl@ = pointer to second argument
601 * @mpw *sv, *svl@ = pointer to scratch workspace
602 *
603 * Returns: ---
604 *
605 * Use: Multiplies two multiprecision integers using Karatsuba's
606 * algorithm. This is rather faster than traditional long
607 * multiplication (e.g., @mpx_umul@) on large numbers, although
608 * more expensive on small ones.
609 *
610 * The destination and scratch buffers must be twice as large as
611 * the larger argument. The scratch space must be twice as
612 * large as the larger argument, plus the magic number
613 * @KARATSUBA_SLOP@.
614 */
615
616 extern void mpx_kmul(mpw */*dv*/, mpw */*dvl*/,
617 const mpw */*av*/, const mpw */*avl*/,
618 const mpw */*bv*/, const mpw */*bvl*/,
619 mpw */*sv*/, mpw */*svl*/);
620
621 /* --- @mpx_ksqr@ --- *
622 *
623 * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer
624 * @const mpw *av, *avl@ = pointer to first argument
625 * @mpw *sv, *svl@ = pointer to scratch workspace
626 *
627 * Returns: ---
628 *
629 * Use: Squares a multiprecision integers using something similar to
630 * Karatsuba's multiplication algorithm. This is rather faster
631 * than traditional long multiplication (e.g., @mpx_umul@) on
632 * large numbers, although more expensive on small ones, and
633 * rather simpler than full-blown Karatsuba multiplication.
634 *
635 * The destination must be twice as large as the argument. The
636 * scratch space must be twice as large as the argument, plus
637 * the magic number @KARATSUBA_SLOP@.
638 */
639
640 extern void mpx_ksqr(mpw */*dv*/, mpw */*dvl*/,
641 const mpw */*av*/, const mpw */*avl*/,
642 mpw */*sv*/, mpw */*svl*/);
643
644 /*----- That's all, folks -------------------------------------------------*/
645
646 #ifdef __cplusplus
647 }
648 #endif
649
650 #endif
Catacomb cryptographic library (downstream of http://git.distorted.org.uk/~mdw/catacomb/)