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