3 * $Id: mp.h,v 1.10 2000/10/08 12:03:16 mdw Exp $
5 * Simple 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.10 2000/10/08 12:03:16 mdw
34 * Provide @mp_eq@ and @MP_EQ@ for rapidly testing equality of two
37 * Revision 1.9 2000/07/29 17:03:31 mdw
38 * Add support for left-to-right bitscanning, for use in modular
41 * Revision 1.8 2000/06/22 19:02:01 mdw
44 * Revision 1.7 2000/06/17 11:45:09 mdw
45 * Major memory management overhaul. Added arena support. Use the secure
46 * arena for secret integers. Replace and improve the MP management macros
47 * (e.g., replace MP_MODIFY by MP_DEST).
49 * Revision 1.6 1999/12/10 23:19:46 mdw
50 * Minor bugfixes. New interface for suggested destinations.
52 * Revision 1.5 1999/11/22 20:50:37 mdw
53 * Add support for computing Jacobi symbols.
55 * Revision 1.4 1999/11/21 22:13:02 mdw
56 * Add mp version of MPX_BITS.
58 * Revision 1.3 1999/11/19 13:19:14 mdw
59 * Fix const annotation.
61 * Revision 1.2 1999/11/17 18:02:16 mdw
62 * New multiprecision integer arithmetic suite.
73 /*----- Header files ------------------------------------------------------*/
80 #ifndef CATACOMB_MPW_H
84 #ifndef CATACOMB_ARENA_H
88 #ifndef CATACOMB_MPARENA_H
92 #ifndef CATACOMB_MPX_H
96 /*----- Data structures ---------------------------------------------------*/
110 #define MP_DESTROYED 16u
112 /*----- Useful constants --------------------------------------------------*/
114 extern mp mp_const
[];
116 #define MP_ZERO (&mp_const[0])
117 #define MP_ONE (&mp_const[1])
118 #define MP_TWO (&mp_const[2])
119 #define MP_THREE (&mp_const[3])
120 #define MP_FOUR (&mp_const[4])
121 #define MP_FIVE (&mp_const[5])
122 #define MP_TEN (&mp_const[6])
123 #define MP_256 (&mp_const[7])
124 #define MP_MONE (&mp_const[8])
126 #define MP_NEW ((mp *)0)
127 #define MP_NEWSEC (&mp_const[9])
129 /*----- Trivial macros ----------------------------------------------------*/
131 /* --- @MP_LEN@ --- *
133 * Arguments: @mp *m@ = pointer to a multiprecision integer
135 * Returns: Length of the integer, in words.
138 #define MP_LEN(m) ((m)->vl - ((m)->v))
140 /*----- Memory management and reference counting --------------------------*/
142 /* --- @mp_new@ --- *
144 * Arguments: @size_t sz@ = size of vector required
145 * @unsigned f@ = flags to set
147 * Returns: Pointer to a new MP structure.
149 * Use: Allocates a new multiprecision integer. The data space is
150 * allocated from either the standard global or secret arena,
151 * depending on the initial flags requested.
154 extern mp
*mp_new(size_t /*sz*/, unsigned /*f*/);
156 /* --- @mp_create@ --- *
158 * Arguments: @size_t sz@ = size of vector required
160 * Returns: Pointer to pristine new MP structure with enough memory
163 * Use: Creates a new multiprecision integer with indeterminate
164 * contents. The integer has a single reference.
167 extern mp
*mp_create(size_t /*sz*/);
169 /* --- @mp_createsecure@ --- *
171 * Arguments: @size_t sz@ = size of vector required
173 * Returns: Pointer to pristine new MP structure with enough memory
176 * Use: Creates a new multiprecision integer with indeterminate
177 * contents. The integer has a single reference. The integer's
178 * data space is allocated from the secure arena. Its burn flag
182 extern mp
*mp_createsecure(size_t /*sz*/);
184 /* --- @mp_build@ --- *
186 * Arguments: @mp *m@ = pointer to an MP block to fill in
187 * @mpw *v@ = pointer to a word array
188 * @mpw *vl@ = pointer just past end of array
192 * Use: Creates a multiprecision integer representing some smallish
193 * number. You must provide storage for the number and dispose
194 * of it when you've finished with it. The number is marked as
195 * constant while it exists.
198 extern void mp_build(mp */
*m*/
, mpw */
*v*/
, mpw */
*vl*/
);
200 /* --- @mp_destroy@ --- *
202 * Arguments: @mp *m@ = pointer to a multiprecision integer
206 * Use: Destroys a multiprecision integer. The reference count isn't
207 * checked. Don't use this function if you don't know what
208 * you're doing: use @mp_drop@ instead.
211 extern void mp_destroy(mp */
*m*/
);
213 /* --- @mp_copy@ --- *
215 * Arguments: @mp *m@ = pointer to a multiprecision integer
217 * Returns: A copy of the given multiprecision integer.
219 * Use: Copies the given integer. In fact you just get another
220 * reference to the same old one again.
223 extern mp
*mp_copy(mp */
*m*/
);
225 #define MP_COPY(m) ((m)->ref++, (m))
227 /* --- @mp_drop@ --- *
229 * Arguments: @mp *m@ = pointer to a multiprecision integer
233 * Use: Drops a reference to an integer which isn't wanted any more.
234 * If there are no more references, the integer is destroyed.
237 extern void mp_drop(mp */
*m*/
);
239 #define MP_DROP(m) do { \
242 if (_mm->ref == 0 && !(_mm->f & MP_CONST)) \
246 /* --- @mp_split@ --- *
248 * Arguments: @mp *m@ = pointer to a multiprecision integer
250 * Returns: A reference to the same integer, possibly with a different
253 * Use: Splits off a modifiable version of the integer referred to.
256 extern mp
*mp_split(mp */
*m*/
);
258 #define MP_SPLIT(m) do { \
260 if ((_m->f & MP_CONST) || _m->ref > 1) { \
261 size_t _len = MP_LEN(_m); \
262 mp *_mm = mp_new(_len, _m->f); \
263 if (!(_m->f & MP_UNDEF)) \
264 memcpy(_mm->v, _m->v, MPWS(_len)); \
271 /* --- @mp_resize@ --- *
273 * Arguments: @mp *m@ = pointer to a multiprecision integer
274 * @size_t sz@ = new size
278 * Use: Resizes the vector containing the integer's digits. The new
279 * size must be at least as large as the current integer's
280 * length. This isn't really intended for client use.
283 extern void mp_resize(mp */
*m*/
, size_t /*sz*/);
285 #define MP_RESIZE(m, ssz) do { \
287 size_t _sz = (ssz); \
288 mparena *_a = (_m->f & MP_BURN) ? MPARENA_SECURE : MPARENA_GLOBAL; \
290 size_t _len = MP_LEN(_m); \
291 assert(((void)"can't make size less than length", _sz >= _len)); \
292 _v = mpalloc(_a, _sz); \
293 if (!(_m->f & MP_UNDEF)) \
294 memcpy(_v, _m->v, MPWS(_len)); \
295 if (_m->f & MP_BURN) \
296 memset(_m->v, 0, MPWS(_m->sz)); \
297 mpfree(_m->a, _m->v); \
300 _m->vl = _v + _len; \
303 /* --- @mp_ensure@ --- *
305 * Arguments: @mp *m@ = pointer to a multiprecision integer
306 * @size_t sz@ = required size
310 * Use: Ensures that the integer has enough space for @sz@ digits.
311 * The value is not changed.
314 extern void mp_ensure(mp */
*m*/
, size_t /*sz*/);
316 #define MP_ENSURE(m, ssz) do { \
318 size_t _ssz = (ssz); \
319 size_t _len = MP_LEN(_m); \
320 if (_ssz >= _len) { \
322 mp_resize(_m, _ssz); \
323 if (!(_m->f & MP_UNDEF) && _ssz > _len) \
324 memset(_m->vl, 0, MPWS(_ssz - _len)); \
325 _m->vl = _m->v + _ssz; \
329 /* --- @mp_dest@ --- *
331 * Arguments: @mp *m@ = a suggested destination integer
332 * @size_t sz@ = size required for result, in digits
333 * @unsigned f@ = various flags
335 * Returns: A pointer to an appropriate destination.
337 * Use: Converts a suggested destination into a real destination with
338 * the required properties. If the real destination is @d@,
339 * then the following properties will hold:
341 * * @d@ will have exactly one reference.
343 * * If @m@ is not @MP_NEW@, then the contents of @m@ will not
344 * change, unless @f@ has the @MP_UNDEF@ flag set.
346 * * If @m@ is not @MP_NEW@, then he reference count of @m@ on
347 * entry is equal to the sum of the counts of @d@ and @m@ on
350 * * The size of @d@ will be at least @sz@.
352 * * If @f@ has the @MP_BURN@ flag set, then @d@ will be
353 * allocated from @MPARENA_SECURE@.
355 * Understanding this function is crucial to using Catacomb's
356 * multiprecision integer library effectively.
359 extern mp
*mp_dest(mp */
*m*/
, size_t /*sz*/, unsigned /*f*/);
361 #define MP_DEST(m, ssz, f) do { \
363 size_t _ssz = (ssz); \
365 _m = mp_dest(_m, _ssz, _f); \
369 /*----- Size manipulation -------------------------------------------------*/
371 /* --- @mp_shrink@ --- *
373 * Arguments: @mp *m@ = pointer to a multiprecision integer
377 * Use: Reduces the recorded length of an integer. This doesn't
378 * reduce the amount of memory used, although it can improve
379 * performance a bit. To reduce memory, use @mp_minimize@
380 * instead. This can't change the value of an integer, and is
381 * therefore safe to use even when there are multiple
385 extern void mp_shrink(mp */
*m*/
);
387 #define MP_SHRINK(m) do { \
389 MPX_SHRINK(_mm->v, _mm->vl); \
394 /* --- @mp_minimize@ --- *
396 * Arguments: @mp *m@ = pointer to a multiprecision integer
400 * Use: Reduces the amount of memory an integer uses. It's best to
401 * do this to numbers which aren't going to change in the
405 extern void mp_minimize(mp */
*m*/
);
407 /*----- Bit scanning ------------------------------------------------------*/
409 #ifndef CATACOMB_MPSCAN_H
413 /* --- @mp_scan@ --- *
415 * Arguments: @mpscan *sc@ = pointer to bitscanner block
416 * @const mp *m@ = pointer to a multiprecision integer
420 * Use: Initializes a bitscanner on a multiprecision integer.
423 extern void mp_scan(mpscan */
*sc*/
, const mp */
*m*/
);
425 #define MP_SCAN(sc, m) do { \
426 const mp *_mm = (m); \
427 mpscan *_sc = (sc); \
428 MPSCAN_INITX(_sc, _mm->v, _mm->vl); \
431 /* --- @mp_rscan@ --- *
433 * Arguments: @mpscan *sc@ = pointer to bitscanner block
434 * @const mp *m@ = pointer to a multiprecision integer
438 * Use: Initializes a reverse bitscanner on a multiprecision
442 extern void mp_rscan(mpscan */
*sc*/
, const mp */
*m*/
);
444 #define MP_RSCAN(sc, m) do { \
445 const mp *_mm = (m); \
446 mpscan *_sc = (sc); \
447 MPSCAN_RINITX(_sc, _mm->v, _mm->vl); \
450 /* --- Other bitscanning aliases --- */
452 #define mp_step mpscan_step
453 #define mp_bit mpscan_bit
454 #define mp_rstep mpscan_rstep
455 #define mp_rbit mpscan_rbit
457 #define MP_STEP MPSCAN_STEP
458 #define MP_BIT MPSCAN_BIT
459 #define MP_RSTEP MPSCAN_RSTEP
460 #define MP_RBIT MPSCAN_RBIT
462 /*----- Loading and storing -----------------------------------------------*/
464 /* --- @mp_octets@ --- *
466 * Arguments: @const mp *m@ = a multiprecision integer
468 * Returns: The number of octets required to represent @m@.
470 * Use: Calculates the external storage required for a multiprecision
474 extern size_t mp_octets(const mp */
*m*/
);
476 /* --- @mp_bits@ --- *
478 * Arguments: @const mp *m@ = a multiprecision integer
480 * Returns: The number of bits required to represent @m@.
482 * Use: Calculates the external storage required for a multiprecision
486 extern unsigned long mp_bits(const mp */
*m*/
);
488 /* --- @mp_loadl@ --- *
490 * Arguments: @mp *d@ = destination
491 * @const void *pv@ = pointer to source data
492 * @size_t sz@ = size of the source data
494 * Returns: Resulting multiprecision number.
496 * Use: Loads a multiprecision number from an array of octets. The
497 * first byte in the array is the least significant. More
498 * formally, if the bytes are %$b_0, b_1, \ldots, b_{n-1}$%
499 * then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%.
502 extern mp
*mp_loadl(mp */
*d*/
, const void */
*pv*/
, size_t /*sz*/);
504 /* --- @mp_storel@ --- *
506 * Arguments: @const mp *m@ = source
507 * @void *pv@ = pointer to output array
508 * @size_t sz@ = size of the output array
512 * Use: Stores a multiprecision number in an array of octets. The
513 * first byte in the array is the least significant. If the
514 * array is too small to represent the number, high-order bits
515 * are truncated; if the array is too large, high order bytes
516 * are filled with zeros. More formally, if the number is
517 * %$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%,
518 * then the array is %$b_0, b_1, \ldots, b_{n-1}$%.
521 extern void mp_storel(const mp */
*m*/
, void */
*pv*/
, size_t /*sz*/);
523 /* --- @mp_loadb@ --- *
525 * Arguments: @mp *d@ = destination
526 * @const void *pv@ = pointer to source data
527 * @size_t sz@ = size of the source data
529 * Returns: Resulting multiprecision number.
531 * Use: Loads a multiprecision number from an array of octets. The
532 * last byte in the array is the least significant. More
533 * formally, if the bytes are %$b_{n-1}, b_{n-2}, \ldots, b_0$%
534 * then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%.
537 extern mp
*mp_loadb(mp */
*d*/
, const void */
*pv*/
, size_t /*sz*/);
539 /* --- @mp_storeb@ --- *
541 * Arguments: @const mp *m@ = source
542 * @void *pv@ = pointer to output array
543 * @size_t sz@ = size of the output array
547 * Use: Stores a multiprecision number in an array of octets. The
548 * last byte in the array is the least significant. If the
549 * array is too small to represent the number, high-order bits
550 * are truncated; if the array is too large, high order bytes
551 * are filled with zeros. More formally, if the number is
552 * %$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%,
553 * then the array is %$b_{n-1}, b_{n-2}, \ldots, b_0$%.
556 extern void mp_storeb(const mp */
*m*/
, void */
*pv*/
, size_t /*sz*/);
558 /*----- Simple arithmetic -------------------------------------------------*/
562 * Arguments: @mp *d@ = destination
565 * Returns: Result, @a@ converted to two's complement notation.
568 extern mp
*mp_2c(mp */
*d*/
, mp */
*a*/
);
572 * Arguments: @mp *d@ = destination
575 * Returns: Result, @a@ converted to the native signed-magnitude
579 extern mp
*mp_sm(mp */
*d*/
, mp */
*a*/
);
581 /* --- @mp_lsl@ --- *
583 * Arguments: @mp *d@ = destination
585 * @size_t n@ = number of bits to move
587 * Returns: Result, @a@ shifted left by @n@.
590 extern mp
*mp_lsl(mp */
*d*/
, mp */
*a*/
, size_t /*n*/);
592 /* --- @mp_lsr@ --- *
594 * Arguments: @mp *d@ = destination
596 * @size_t n@ = number of bits to move
598 * Returns: Result, @a@ shifted left by @n@.
601 extern mp
*mp_lsr(mp */
*d*/
, mp */
*a*/
, size_t /*n*/);
605 * Arguments: @const mp *a, *b@ = two numbers
607 * Returns: Nonzero if the numbers are equal.
610 extern int mp_eq(const mp */
*a*/
, const mp */
*b*/
);
612 #define MP_EQ(a, b) \
613 ((((a)->f ^ (b)->f) & MP_NEG) == 0 && \
614 mpx_ueq((a)->v, (a)->vl, (b)->v, (b)->vl))
616 /* --- @mp_cmp@ --- *
618 * Arguments: @const mp *a, *b@ = two numbers
620 * Returns: Less than, equal to or greater than zero, according to
621 * whether @a@ is less than, equal to or greater than @b@.
624 extern int mp_cmp(const mp */
*a*/
, const mp */
*b*/
);
626 #define MP_CMP(a, op, b) (mp_cmp((a), (b)) op 0)
628 /* --- @mp_add@ --- *
630 * Arguments: @mp *d@ = destination
631 * @mp *a, *b@ = sources
633 * Returns: Result, @a@ added to @b@.
636 extern mp
*mp_add(mp */
*d*/
, mp */
*a*/
, mp */
*b*/
);
638 /* --- @mp_sub@ --- *
640 * Arguments: @mp *d@ = destination
641 * @mp *a, *b@ = sources
643 * Returns: Result, @b@ subtracted from @a@.
646 extern mp
*mp_sub(mp */
*d*/
, mp */
*a*/
, mp */
*b*/
);
648 /* --- @mp_mul@ --- *
650 * Arguments: @mp *d@ = destination
651 * @mp *a, *b@ = sources
653 * Returns: Result, @a@ multiplied by @b@.
656 extern mp
*mp_mul(mp */
*d*/
, mp */
*a*/
, mp */
*b*/
);
658 /* --- @mp_sqr@ --- *
660 * Arguments: @mp *d@ = destination
663 * Returns: Result, @a@ squared.
666 extern mp
*mp_sqr(mp */
*d*/
, mp */
*a*/
);
668 /* --- @mp_div@ --- *
670 * Arguments: @mp **qq, **rr@ = destination, quotient and remainder
671 * @mp *a, *b@ = sources
673 * Use: Calculates the quotient and remainder when @a@ is divided by
677 extern void mp_div(mp
**/
*qq*/
, mp
**/
*rr*/
, mp */
*a*/
, mp */
*b*/
);
679 /* --- @mp_odd@ --- *
681 * Arguments: @mp *d@ = pointer to destination integer
682 * @mp *m@ = pointer to source integer
683 * @size_t *s@ = where to store the power of 2
685 * Returns: An odd integer integer %$t$% such that %$m = 2^s t$%.
687 * Use: Computes a power of two and an odd integer which, when
688 * multiplied, give a specified result. This sort of thing is
689 * useful in number theory quite often.
692 extern mp
*mp_odd(mp */
*d*/
, mp */
*m*/
, size_t */
*s*/
);
694 /*----- More advanced algorithms ------------------------------------------*/
696 /* --- @mp_sqrt@ --- *
698 * Arguments: @mp *d@ = pointer to destination integer
699 * @mp *a@ = (nonnegative) integer to take square root of
701 * Returns: The largest integer %$x$% such that %$x^2 \le a$%.
703 * Use: Computes integer square roots.
705 * The current implementation isn't very good: it uses the
706 * Newton-Raphson method to find an approximation to %$a$%. If
707 * there's any demand for a better version, I'll write one.
710 extern mp
*mp_sqrt(mp */
*d*/
, mp */
*a*/
);
712 /* --- @mp_gcd@ --- *
714 * Arguments: @mp **gcd, **xx, **yy@ = where to write the results
715 * @mp *a, *b@ = sources (must be nonzero)
719 * Use: Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that
720 * @ax + by = gcd(a, b)@. This is useful for computing modular
721 * inverses. Neither @a@ nor @b@ may be zero.
724 extern void mp_gcd(mp
**/
*gcd*/
, mp
**/
*xx*/
, mp
**/
*yy*/
,
725 mp */
*a*/
, mp */
*b*/
);
727 /* --- @mp_jacobi@ --- *
729 * Arguments: @mp *a@ = an integer less than @n@
730 * @mp *n@ = an odd integer
732 * Returns: @-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
734 * Use: Computes the Jacobi symbol. If @n@ is prime, this is the
735 * Legendre symbol and is equal to 1 if and only if @a@ is a
736 * quadratic residue mod @n@. The result is zero if and only if
737 * @a@ and @n@ have a common factor greater than one.
740 extern int mp_jacobi(mp */
*a*/
, mp */
*n*/
);
742 /* --- @mp_modsqrt@ --- *
744 * Arguments: @mp *d@ = destination integer
745 * @mp *a@ = source integer
746 * @mp *p@ = modulus (must be prime)
748 * Returns: If %$a$% is a quadratic residue, a square root of %$a$%; else
751 * Use: Returns an integer %$x$% such that %$x^2 \equiv a \pmod{p}$%,
752 * if one exists; else a null pointer. This function will not
753 * work if %$p$% is composite: you must factor the modulus, take
754 * a square root mod each factor, and recombine the results
755 * using the Chinese Remainder Theorem.
758 extern mp
*mp_modsqrt(mp */
*d*/
, mp */
*a*/
, mp */
*p*/
);
760 /*----- Test harness support ----------------------------------------------*/
762 #include <mLib/testrig.h>
764 #ifndef CATACOMB_MPTEXT_H
768 extern const test_type type_mp
;
770 /*----- That's all, folks -------------------------------------------------*/