a80a8e33726d6b68d981bfab2655ccd9d9a291b3
[u/mdw/catacomb] / mp.h
1 /* -*-c-*-
2 *
3 * $Id: mp.h,v 1.6 1999/12/10 23:19:46 mdw Exp $
4 *
5 * Simple 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: mp.h,v $
33 * Revision 1.6 1999/12/10 23:19:46 mdw
34 * Minor bugfixes. New interface for suggested destinations.
35 *
36 * Revision 1.5 1999/11/22 20:50:37 mdw
37 * Add support for computing Jacobi symbols.
38 *
39 * Revision 1.4 1999/11/21 22:13:02 mdw
40 * Add mp version of MPX_BITS.
41 *
42 * Revision 1.3 1999/11/19 13:19:14 mdw
43 * Fix const annotation.
44 *
45 * Revision 1.2 1999/11/17 18:02:16 mdw
46 * New multiprecision integer arithmetic suite.
47 *
48 */
49
50 #ifndef CATACOMB_MP_H
51 #define CATACOMB_MP_H
52
53 #ifdef __cplusplus
54 extern "C" {
55 #endif
56
57 /*----- Header files ------------------------------------------------------*/
58
59 #include <assert.h>
60 #include <string.h>
61
62 #include <mLib/sub.h>
63
64 #ifndef CATACOMB_MPW_H
65 # include "mpw.h"
66 #endif
67
68 #ifndef CATACOMB_MPX_H
69 # include "mpx.h"
70 #endif
71
72 /*----- Data structures ---------------------------------------------------*/
73
74 typedef struct mp {
75 mpw *v, *vl;
76 size_t sz;
77 unsigned f;
78 unsigned ref;
79 } mp;
80
81 #define MP_NEG 1u
82 #define MP_BURN 2u
83 #define MP_CONST 4u
84 #define MP_UNDEF 8u
85 #define MP_DESTROYED 16u
86
87 /*----- Useful constants --------------------------------------------------*/
88
89 extern mp mp_const[];
90
91 #define MP_ZERO (&mp_const[0])
92 #define MP_ONE (&mp_const[1])
93 #define MP_TWO (&mp_const[2])
94 #define MP_THREE (&mp_const[3])
95 #define MP_FOUR (&mp_const[4])
96 #define MP_FIVE (&mp_const[5])
97 #define MP_TEN (&mp_const[6])
98 #define MP_MONE (&mp_const[7])
99
100 #define MP_NEW ((mp *)0)
101
102 /*----- Memory allocation hooks -------------------------------------------*/
103
104 #ifndef CATACOMB_MPARENA_H
105 # include "mparena.h"
106 #endif
107
108 /* --- @MP_ARENA@ --- *
109 *
110 * This selects where memory is allocated from. Tweak to use more fancy
111 * things like custom arenas.
112 */
113
114 #ifndef MP_ARENA
115 # define MP_ARENA MPARENA_GLOBAL
116 #endif
117
118 /* --- @MP_ALLOC@ --- *
119 *
120 * Arguments: @size_t sz@ = size required
121 *
122 * Returns: Pointer to an allocated vector of the requested size.
123 *
124 * Use: Hook for vector allocation.
125 */
126
127 #ifndef MP_ALLOC
128 # define MP_ALLOC(sz) mpalloc(MP_ARENA, (sz))
129 #endif
130
131 /* --- @MP_FREE@ --- *
132 *
133 * Arguments: @mpw *v@ = pointer to vector
134 *
135 * Returns: ---
136 *
137 * Use: Hook for vector deallocation.
138 */
139
140 #ifndef MP_FREE
141 # define MP_FREE(v) mpfree(MP_ARENA, (v))
142 #endif
143
144 /*----- Paranoia management -----------------------------------------------*/
145
146 /* --- @mp_burn@ --- *
147 *
148 * Arguments: @mp *m@ = pointer to a multiprecision integer
149 *
150 * Returns: ---
151 *
152 * Use: Marks the integer as `burn-after-use'. When the integer's
153 * memory is deallocated, it is deleted so that traces can't
154 * remain in the swap file. In theory.
155 */
156
157 extern void mp_burn(mp */*m*/);
158
159 /*----- Trivial macros ----------------------------------------------------*/
160
161 /* --- @MP_LEN@ --- *
162 *
163 * Arguments: @mp *m@ = pointer to a multiprecision integer
164 *
165 * Returns: Length of the integer, in words.
166 */
167
168 #define MP_LEN(m) ((m)->vl - ((m)->v))
169
170 /*----- Memory management and reference counting --------------------------*/
171
172 /* --- @mp_create@ --- *
173 *
174 * Arguments: @size_t sz@ = size of vector required
175 *
176 * Returns: Pointer to pristine new MP structure with enough memory
177 * bolted onto it.
178 *
179 * Use: Creates a new multiprecision integer with indeterminate
180 * contents. The integer has a single reference.
181 */
182
183 extern mp *mp_create(size_t /*sz*/);
184
185 /* --- @mp_build@ --- *
186 *
187 * Arguments: @mp *m@ = pointer to an MP block to fill in
188 * @mpw *v@ = pointer to a word array
189 * @mpw *vl@ = pointer just past end of array
190 *
191 * Returns: ---
192 *
193 * Use: Creates a multiprecision integer representing some smallish
194 * number. You must provide storage for the number and dispose
195 * of it when you've finished with it. The number is marked as
196 * constant while it exists.
197 */
198
199 extern void mp_build(mp */*m*/, mpw */*v*/, mpw */*vl*/);
200
201 /* --- @mp_destroy@ --- *
202 *
203 * Arguments: @mp *m@ = pointer to a multiprecision integer
204 *
205 * Returns: ---
206 *
207 * Use: Destroys a multiprecision integer. The reference count isn't
208 * checked. Don't use this function if you don't know what
209 * you're doing: use @mp_drop@ instead.
210 */
211
212 extern void mp_destroy(mp */*m*/);
213
214 /* --- @mp_copy@ --- *
215 *
216 * Arguments: @mp *m@ = pointer to a multiprecision integer
217 *
218 * Returns: A copy of the given multiprecision integer.
219 *
220 * Use: Copies the given integer. In fact you just get another
221 * reference to the same old one again.
222 */
223
224 extern mp *mp_copy(mp */*m*/);
225
226 #define MP_COPY(m) ((m)->ref++, (m))
227
228 /* --- @mp_drop@ --- *
229 *
230 * Arguments: @mp *m@ = pointer to a multiprecision integer
231 *
232 * Returns: ---
233 *
234 * Use: Drops a reference to an integer which isn't wanted any more.
235 * If there are no more references, the integer is destroyed.
236 */
237
238 extern void mp_drop(mp */*m*/);
239
240 #define MP_DROP(m) do { \
241 mp *_mm = (m); \
242 if (_mm->ref > 1) \
243 _mm->ref--; \
244 else if (!(_mm->f & MP_CONST)) \
245 mp_destroy(_mm); \
246 } while (0)
247
248 /* --- @mp_split@ --- *
249 *
250 * Arguments: @mp *m@ = pointer to a multiprecision integer
251 *
252 * Returns: A reference to the same integer, possibly with a different
253 * address.
254 *
255 * Use: Splits off a modifiable version of the integer referred to.
256 */
257
258 extern mp *mp_split(mp */*m*/);
259
260 #define MP_SPLIT(m) do { \
261 mp *_mm = (m); \
262 if ((_mm->f & MP_CONST) || _mm->ref != 1) { \
263 mp *_dd = mp_create(_mm->sz); \
264 _dd->vl = _dd->v + MP_LEN(_mm); \
265 _dd->f = _mm->f & (MP_NEG | MP_BURN); \
266 memcpy(_dd->v, _mm->v, MPWS(MP_LEN(_mm))); \
267 _dd->ref = 1; \
268 _mm->ref--; \
269 (m) = _dd; \
270 } \
271 } while (0)
272
273 /* --- @mp_resize@ --- *
274 *
275 * Arguments: @mp *m@ = pointer to a multiprecision integer
276 * @size_t sz@ = new size
277 *
278 * Returns: ---
279 *
280 * Use: Resizes the vector containing the integer's digits. The new
281 * size must be at least as large as the current integer's
282 * length. The integer's length is increased and new digits are
283 * filled with zeroes. This isn't really intended for client
284 * use.
285 */
286
287 extern void mp_resize(mp */*m*/, size_t /*sz*/);
288
289 #define MP_RESIZE(m, ssz) do { \
290 mp *_m = (m); \
291 size_t _sz = (ssz); \
292 size_t _len = MP_LEN(_m); \
293 mpw *_v = MP_ALLOC(_sz); \
294 if (!(_m->f & MP_UNDEF)) \
295 memcpy(_v, _m->v, MPWS(_len)); \
296 if (_m->f & MP_BURN) \
297 memset(_m->v, 0, MPWS(_m->sz)); \
298 MP_FREE(_m->v); \
299 _m->v = _v; \
300 _m->vl = _v + _len; \
301 _m->sz = _sz; \
302 } while (0)
303
304 /* --- @mp_ensure@ --- *
305 *
306 * Arguments: @mp *m@ = pointer to a multiprecision integer
307 * @size_t sz@ = required size
308 *
309 * Returns: ---
310 *
311 * Use: Ensures that the integer has enough space for @sz@ digits.
312 * The value is not changed.
313 */
314
315 extern void mp_ensure(mp */*m*/, size_t /*sz*/);
316
317 #define MP_ENSURE(m, ssz) do { \
318 mp *_mm = (m); \
319 size_t _ssz = (ssz); \
320 size_t _len = MP_LEN(_mm); \
321 if (_ssz > _mm->sz) \
322 MP_RESIZE(_mm, _ssz); \
323 if (!(_mm->f & MP_UNDEF) && _ssz > _len) \
324 memset(_mm->vl, 0, MPWS(_ssz - _len)); \
325 _mm->vl = _mm->v + _ssz; \
326 } while (0)
327
328 /* --- @mp_modify@ --- *
329 *
330 * Arguments: @mp *m@ = pointer to a multiprecision integer
331 * @size_t sz@ = size required
332 *
333 * Returns: Pointer to the integer (possibly different).
334 *
335 * Use: Prepares an integer to be overwritten. It's split off from
336 * other references to the same integer, and sufficient space is
337 * allocated.
338 */
339
340 extern mp *mp_modify(mp */*m*/, size_t /*sz*/);
341
342 #define MP_MODIFY(m, sz) do { \
343 size_t _rq = (sz); \
344 mp *_m = (m); \
345 if (_m == MP_NEW || m->ref > 1 || (_m->f & MP_CONST)) { \
346 if (_m) \
347 MP_DROP(_m); \
348 _m = mp_create(_rq); \
349 } else \
350 MP_ENSURE(_m, _rq); \
351 (m) = _m; \
352 } while (0)
353
354 /*----- Size manipulation -------------------------------------------------*/
355
356 /* --- @mp_shrink@ --- *
357 *
358 * Arguments: @mp *m@ = pointer to a multiprecision integer
359 *
360 * Returns: ---
361 *
362 * Use: Reduces the recorded length of an integer. This doesn't
363 * reduce the amount of memory used, although it can improve
364 * performance a bit. To reduce memory, use @mp_minimize@
365 * instead. This can't change the value of an integer, and is
366 * therefore safe to use even when there are multiple
367 * references.
368 */
369
370 extern void mp_shrink(mp */*m*/);
371
372 #define MP_SHRINK(m) do { \
373 mp *_mm = (m); \
374 MPX_SHRINK(_mm->v, _mm->vl); \
375 if (!MP_LEN(_mm)) \
376 _mm->f &= ~MP_NEG; \
377 } while (0)
378
379 /* --- @mp_minimize@ --- *
380 *
381 * Arguments: @mp *m@ = pointer to a multiprecision integer
382 *
383 * Returns: ---
384 *
385 * Use: Reduces the amount of memory an integer uses. It's best to
386 * do this to numbers which aren't going to change in the
387 * future.
388 */
389
390 extern void mp_minimize(mp */*m*/);
391
392 /*----- Bit scanning ------------------------------------------------------*/
393
394 #ifndef CATACOMB_MPSCAN_H
395 # include "mpscan.h"
396 #endif
397
398 /* --- @mp_scan@ --- *
399 *
400 * Arguments: @mpscan *sc@ = pointer to bitscanner block
401 * @const mp *m@ = pointer to a multiprecision integer
402 *
403 * Returns: ---
404 *
405 * Use: Initializes a bitscanner on a multiprecision integer.
406 */
407
408 extern void mp_scan(mpscan */*sc*/, const mp */*m*/);
409
410 #define MP_SCAN(sc, m) do { \
411 const mp *_mm = (m); \
412 mpscan *_sc = (sc); \
413 MPSCAN_INITX(_sc, _mm->v, _mm->vl); \
414 } while (0)
415
416 /* --- Other bitscanning aliases --- */
417
418 #define mp_step mpscan_step
419 #define mp_bit mpscan_bit
420
421 #define MP_STEP MPSCAN_STEP
422 #define MP_BIT MPSCAN_BIT
423
424 /*----- Loading and storing -----------------------------------------------*/
425
426 /* --- @mp_octets@ --- *
427 *
428 * Arguments: @const mp *m@ = a multiprecision integer
429 *
430 * Returns: The number of octets required to represent @m@.
431 *
432 * Use: Calculates the external storage required for a multiprecision
433 * integer.
434 */
435
436 extern size_t mp_octets(const mp */*m*/);
437
438 /* --- @mp_bits@ --- *
439 *
440 * Arguments: @const mp *m@ = a multiprecision integer
441 *
442 * Returns: The number of bits required to represent @m@.
443 *
444 * Use: Calculates the external storage required for a multiprecision
445 * integer.
446 */
447
448 extern unsigned long mp_bits(const mp */*m*/);
449
450 /* --- @mp_loadl@ --- *
451 *
452 * Arguments: @mp *d@ = destination
453 * @const void *pv@ = pointer to source data
454 * @size_t sz@ = size of the source data
455 *
456 * Returns: Resulting multiprecision number.
457 *
458 * Use: Loads a multiprecision number from an array of octets. The
459 * first byte in the array is the least significant. More
460 * formally, if the bytes are %$b_0, b_1, \ldots, b_{n-1}$%
461 * then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%.
462 */
463
464 extern mp *mp_loadl(mp */*d*/, const void */*pv*/, size_t /*sz*/);
465
466 /* --- @mp_storel@ --- *
467 *
468 * Arguments: @const mp *m@ = source
469 * @void *pv@ = pointer to output array
470 * @size_t sz@ = size of the output array
471 *
472 * Returns: ---
473 *
474 * Use: Stores a multiprecision number in an array of octets. The
475 * first byte in the array is the least significant. If the
476 * array is too small to represent the number, high-order bits
477 * are truncated; if the array is too large, high order bytes
478 * are filled with zeros. More formally, if the number is
479 * %$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%,
480 * then the array is %$b_0, b_1, \ldots, b_{n-1}$%.
481 */
482
483 extern void mp_storel(const mp */*m*/, void */*pv*/, size_t /*sz*/);
484
485 /* --- @mp_loadb@ --- *
486 *
487 * Arguments: @mp *d@ = destination
488 * @const void *pv@ = pointer to source data
489 * @size_t sz@ = size of the source data
490 *
491 * Returns: Resulting multiprecision number.
492 *
493 * Use: Loads a multiprecision number from an array of octets. The
494 * last byte in the array is the least significant. More
495 * formally, if the bytes are %$b_{n-1}, b_{n-2}, \ldots, b_0$%
496 * then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%.
497 */
498
499 extern mp *mp_loadb(mp */*d*/, const void */*pv*/, size_t /*sz*/);
500
501 /* --- @mp_storeb@ --- *
502 *
503 * Arguments: @const mp *m@ = source
504 * @void *pv@ = pointer to output array
505 * @size_t sz@ = size of the output array
506 *
507 * Returns: ---
508 *
509 * Use: Stores a multiprecision number in an array of octets. The
510 * last byte in the array is the least significant. If the
511 * array is too small to represent the number, high-order bits
512 * are truncated; if the array is too large, high order bytes
513 * are filled with zeros. More formally, if the number is
514 * %$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%,
515 * then the array is %$b_{n-1}, b_{n-2}, \ldots, b_0$%.
516 */
517
518 extern void mp_storeb(const mp */*m*/, void */*pv*/, size_t /*sz*/);
519
520 /*----- Simple arithmetic -------------------------------------------------*/
521
522 /* --- @mp_2c@ --- *
523 *
524 * Arguments: @mp *d@ = destination
525 * @mp *a@ = source
526 *
527 * Returns: Result, @a@ converted to two's complement notation.
528 */
529
530 extern mp *mp_2c(mp */*d*/, mp */*a*/);
531
532 /* --- @mp_sm@ --- *
533 *
534 * Arguments: @mp *d@ = destination
535 * @mp *a@ = source
536 *
537 * Returns: Result, @a@ converted to the native signed-magnitude
538 * notation.
539 */
540
541 extern mp *mp_sm(mp */*d*/, mp */*a*/);
542
543 /* --- @mp_lsl@ --- *
544 *
545 * Arguments: @mp *d@ = destination
546 * @mp *a@ = source
547 * @size_t n@ = number of bits to move
548 *
549 * Returns: Result, @a@ shifted left by @n@.
550 */
551
552 extern mp *mp_lsl(mp */*d*/, mp */*a*/, size_t /*n*/);
553
554 /* --- @mp_lsr@ --- *
555 *
556 * Arguments: @mp *d@ = destination
557 * @mp *a@ = source
558 * @size_t n@ = number of bits to move
559 *
560 * Returns: Result, @a@ shifted left by @n@.
561 */
562
563 extern mp *mp_lsr(mp */*d*/, mp */*a*/, size_t /*n*/);
564
565 /* --- @mp_cmp@ --- *
566 *
567 * Arguments: @const mp *a, *b@ = two numbers
568 *
569 * Returns: Less than, equal to or greater than zero, according to
570 * whether @a@ is less than, equal to or greater than @b@.
571 */
572
573 extern int mp_cmp(const mp */*a*/, const mp */*b*/);
574
575 #define MP_CMP(a, op, b) (mp_cmp((a), (b)) op 0)
576
577 /* --- @mp_add@ --- *
578 *
579 * Arguments: @mp *d@ = destination
580 * @mp *a, *b@ = sources
581 *
582 * Returns: Result, @a@ added to @b@.
583 */
584
585 extern mp *mp_add(mp */*d*/, mp */*a*/, mp */*b*/);
586
587 /* --- @mp_sub@ --- *
588 *
589 * Arguments: @mp *d@ = destination
590 * @mp *a, *b@ = sources
591 *
592 * Returns: Result, @b@ subtracted from @a@.
593 */
594
595 extern mp *mp_sub(mp */*d*/, mp */*a*/, mp */*b*/);
596
597 /* --- @mp_mul@ --- *
598 *
599 * Arguments: @mp *d@ = destination
600 * @mp *a, *b@ = sources
601 *
602 * Returns: Result, @a@ multiplied by @b@.
603 */
604
605 extern mp *mp_mul(mp */*d*/, mp */*a*/, mp */*b*/);
606
607 /* --- @mp_sqr@ --- *
608 *
609 * Arguments: @mp *d@ = destination
610 * @mp *a@ = source
611 *
612 * Returns: Result, @a@ squared.
613 */
614
615 extern mp *mp_sqr(mp */*d*/, mp */*a*/);
616
617 /* --- @mp_div@ --- *
618 *
619 * Arguments: @mp **qq, **rr@ = destination, quotient and remainder
620 * @mp *a, *b@ = sources
621 *
622 * Use: Calculates the quotient and remainder when @a@ is divided by
623 * @b@.
624 */
625
626 extern void mp_div(mp **/*qq*/, mp **/*rr*/, mp */*a*/, mp */*b*/);
627
628 /*----- More advanced algorithms ------------------------------------------*/
629
630 /* --- @mp_gcd@ --- *
631 *
632 * Arguments: @mp **gcd, **xx, **yy@ = where to write the results
633 * @mp *a, *b@ = sources (must be nonzero)
634 *
635 * Returns: ---
636 *
637 * Use: Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that
638 * @ax + by = gcd(a, b)@. This is useful for computing modular
639 * inverses. Neither @a@ nor @b@ may be zero.
640 */
641
642 extern void mp_gcd(mp **/*gcd*/, mp **/*xx*/, mp **/*yy*/,
643 mp */*a*/, mp */*b*/);
644
645 /* --- @mp_jacobi@ --- *
646 *
647 * Arguments: @mp *a@ = an integer less than @n@
648 * @mp *n@ = an odd integer
649 *
650 * Returns: @-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
651 *
652 * Use: Computes the Jacobi symbol. If @n@ is prime, this is the
653 * Legendre symbol and is equal to 1 if and only if @a@ is a
654 * quadratic residue mod @n@. The result is zero if and only if
655 * @a@ and @n@ have a common factor greater than one.
656 */
657
658 int mp_jacobi(mp */*a*/, mp */*n*/);
659
660 /*----- Test harness support ----------------------------------------------*/
661
662 #include <mLib/testrig.h>
663
664 #ifndef CATACOMB_MPTEXT_H
665 # include "mptext.h"
666 #endif
667
668 extern const test_type type_mp;
669
670 /*----- That's all, folks -------------------------------------------------*/
671
672 #ifdef __cplusplus
673 }
674 #endif
675
676 #endif