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