3 * Simple multiprecision arithmetic
5 * (c) 1999 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of Catacomb.
12 * Catacomb is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU Library General Public License as
14 * published by the Free Software Foundation; either version 2 of the
15 * License, or (at your option) any later version.
17 * Catacomb is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Library General Public License for more details.
22 * You should have received a copy of the GNU Library General Public
23 * License along with Catacomb; if not, write to the Free
24 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
35 /*----- Header files ------------------------------------------------------*/
40 #include <mLib/macros.h>
43 #ifndef CATACOMB_MPW_H
47 #ifndef CATACOMB_ARENA_H
51 #ifndef CATACOMB_MPARENA_H
55 #ifndef CATACOMB_MPX_H
59 /*----- Data structures ---------------------------------------------------*/
61 /* --- A multiprecision integer --- */
64 mpw *v, *vl; /* Vector of digits, current limit */
65 size_t sz; /* Size of digit buffer in words */
66 mparena *a; /* Arena for buffer allocation */
67 unsigned f; /* Flags (see below) */
68 unsigned ref; /* Reference counter */
71 #define MP_NEG 1u /* Negative (signed magnitude) */
72 #define MP_BURN 2u /* Secret (viral flag) */
73 #define MP_CONST 4u /* Uses strange memory allocation */
74 #define MP_UNDEF 8u /* Contains nothing interesting */
75 #define MP_DESTROYED 16u /* Has been destroyed */
77 /* --- A factor for simultaneous exponentation --- *
79 * Used by the Montgomery and Barrett exponentiators.
82 typedef struct mp_expfactor {
87 /*----- Useful constants --------------------------------------------------*/
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_256 (&mp_const[7])
99 #define MP_MONE (&mp_const[8])
101 #define MP_NEW ((mp *)0)
102 #define MP_NEWSEC (&mp_const[9])
104 /*----- Trivial macros ----------------------------------------------------*/
106 /* --- @MP_LEN@ --- *
108 * Arguments: @mp *m@ = pointer to a multiprecision integer
110 * Returns: Length of the integer, in words.
113 #define MP_LEN(m) ((m)->vl - ((m)->v))
115 /*----- Memory management and reference counting --------------------------*/
117 /* --- @mp_new@ --- *
119 * Arguments: @size_t sz@ = size of vector required
120 * @unsigned f@ = flags to set
122 * Returns: Pointer to a new MP structure.
124 * Use: Allocates a new multiprecision integer. The data space is
125 * allocated from either the standard global or secret arena,
126 * depending on the initial flags requested.
129 extern mp *mp_new(size_t /*sz*/, unsigned /*f*/);
131 /* --- @mp_create@ --- *
133 * Arguments: @size_t sz@ = size of vector required
135 * Returns: Pointer to pristine new MP structure with enough memory
138 * Use: Creates a new multiprecision integer with indeterminate
139 * contents. The integer has a single reference.
142 extern mp *mp_create(size_t /*sz*/);
144 /* --- @mp_createsecure@ --- *
146 * Arguments: @size_t sz@ = size of vector required
148 * Returns: Pointer to pristine new MP structure with enough memory
151 * Use: Creates a new multiprecision integer with indeterminate
152 * contents. The integer has a single reference. The integer's
153 * data space is allocated from the secure arena. Its burn flag
157 extern mp *mp_createsecure(size_t /*sz*/);
159 /* --- @mp_build@ --- *
161 * Arguments: @mp *m@ = pointer to an MP block to fill in
162 * @mpw *v@ = pointer to a word array
163 * @mpw *vl@ = pointer just past end of array
167 * Use: Creates a multiprecision integer representing some smallish
168 * number. You must provide storage for the number and dispose
169 * of it when you've finished with it. The number is marked as
170 * constant while it exists.
173 extern void mp_build(mp */*m*/, mpw */*v*/, mpw */*vl*/);
175 /* --- @mp_destroy@ --- *
177 * Arguments: @mp *m@ = pointer to a multiprecision integer
181 * Use: Destroys a multiprecision integer. The reference count isn't
182 * checked. Don't use this function if you don't know what
183 * you're doing: use @mp_drop@ instead.
186 extern void mp_destroy(mp */*m*/);
188 /* --- @mp_copy@ --- *
190 * Arguments: @mp *m@ = pointer to a multiprecision integer
192 * Returns: A copy of the given multiprecision integer.
194 * Use: Copies the given integer. In fact you just get another
195 * reference to the same old one again.
198 extern mp *mp_copy(mp */*m*/);
200 #define MP_COPY(m) MUFFLE_WARNINGS_EXPR(GCC_WARNING("-Wunused-value"), \
203 /* --- @mp_drop@ --- *
205 * Arguments: @mp *m@ = pointer to a multiprecision integer
209 * Use: Drops a reference to an integer which isn't wanted any more.
210 * If there are no more references, the integer is destroyed.
213 extern void mp_drop(mp */*m*/);
215 #define MP_DROP(m) do { \
218 if (_mm->ref == 0 && !(_mm->f & MP_CONST)) \
222 /* --- @mp_split@ --- *
224 * Arguments: @mp *m@ = pointer to a multiprecision integer
226 * Returns: A reference to the same integer, possibly with a different
229 * Use: Splits off a modifiable version of the integer referred to.
232 extern mp *mp_split(mp */*m*/);
234 #define MP_SPLIT(m) do { \
236 if ((_m->f & MP_CONST) || _m->ref > 1) { \
237 size_t _len = MP_LEN(_m); \
238 mp *_mm = mp_new(_len, _m->f); \
239 if (!(_m->f & MP_UNDEF)) \
240 memcpy(_mm->v, _m->v, MPWS(_len)); \
247 /* --- @mp_resize@ --- *
249 * Arguments: @mp *m@ = pointer to a multiprecision integer
250 * @size_t sz@ = new size
254 * Use: Resizes the vector containing the integer's digits. The new
255 * size must be at least as large as the current integer's
256 * length. This isn't really intended for client use.
259 extern void mp_resize(mp */*m*/, size_t /*sz*/);
261 #define MP_RESIZE(m, ssz) do { \
263 size_t _sz = (ssz); \
264 mparena *_a = (_m->f & MP_BURN) ? MPARENA_SECURE : MPARENA_GLOBAL; \
266 size_t _len = MP_LEN(_m); \
267 assert(((void)"can't make size less than length", _sz >= _len)); \
268 _v = mpalloc(_a, _sz); \
269 if (!(_m->f & MP_UNDEF)) \
270 memcpy(_v, _m->v, MPWS(_len)); \
271 if (_m->f & MP_BURN) \
272 memset(_m->v, 0, MPWS(_m->sz)); \
273 mpfree(_m->a, _m->v); \
276 _m->vl = _v + _len; \
279 /* --- @mp_ensure@ --- *
281 * Arguments: @mp *m@ = pointer to a multiprecision integer
282 * @size_t sz@ = required size
286 * Use: Ensures that the integer has enough space for @sz@ digits.
287 * The value is not changed.
290 extern void mp_ensure(mp */*m*/, size_t /*sz*/);
292 #define MP_ENSURE(m, ssz) do { \
294 size_t _ssz = (ssz); \
295 size_t _len = MP_LEN(_m); \
296 if (_ssz >= _len) { \
298 mp_resize(_m, _ssz); \
299 if (!(_m->f & MP_UNDEF) && _ssz > _len) \
300 memset(_m->vl, 0, MPWS(_ssz - _len)); \
301 _m->vl = _m->v + _ssz; \
305 /* --- @mp_dest@ --- *
307 * Arguments: @mp *m@ = a suggested destination integer
308 * @size_t sz@ = size required for result, in digits
309 * @unsigned f@ = various flags
311 * Returns: A pointer to an appropriate destination.
313 * Use: Converts a suggested destination into a real destination with
314 * the required properties. If the real destination is @d@,
315 * then the following properties will hold:
317 * * @d@ will have exactly one reference.
319 * * If @m@ is not @MP_NEW@, then the contents of @m@ will not
320 * change, unless @f@ has the @MP_UNDEF@ flag set.
322 * * If @m@ is not @MP_NEW@, then he reference count of @m@ on
323 * entry is equal to the sum of the counts of @d@ and @m@ on
326 * * The size of @d@ will be at least @sz@.
328 * * If @f@ has the @MP_BURN@ flag set, then @d@ will be
329 * allocated from @MPARENA_SECURE@.
331 * Understanding this function is crucial to using Catacomb's
332 * multiprecision integer library effectively.
335 extern mp *mp_dest(mp */*m*/, size_t /*sz*/, unsigned /*f*/);
337 #define MP_DEST(m, ssz, f) do { \
339 size_t _ssz = (ssz); \
341 _m = mp_dest(_m, _ssz, _f); \
345 /*----- Size manipulation -------------------------------------------------*/
347 /* --- @mp_shrink@ --- *
349 * Arguments: @mp *m@ = pointer to a multiprecision integer
353 * Use: Reduces the recorded length of an integer. This doesn't
354 * reduce the amount of memory used, although it can improve
355 * performance a bit. To reduce memory, use @mp_minimize@
356 * instead. This can't change the value of an integer, and is
357 * therefore safe to use even when there are multiple
361 extern void mp_shrink(mp */*m*/);
363 #define MP_SHRINK(m) do { \
365 MPX_SHRINK(_mm->v, _mm->vl); \
370 /* --- @mp_minimize@ --- *
372 * Arguments: @mp *m@ = pointer to a multiprecision integer
376 * Use: Reduces the amount of memory an integer uses. It's best to
377 * do this to numbers which aren't going to change in the
381 extern void mp_minimize(mp */*m*/);
383 /*----- Bit scanning ------------------------------------------------------*/
385 #ifndef CATACOMB_MPSCAN_H
389 /* --- @mp_scan@ --- *
391 * Arguments: @mpscan *sc@ = pointer to bitscanner block
392 * @const mp *m@ = pointer to a multiprecision integer
396 * Use: Initializes a bitscanner on a multiprecision integer.
399 extern void mp_scan(mpscan */*sc*/, const mp */*m*/);
401 #define MP_SCAN(sc, m) do { \
402 const mp *_mm = (m); \
403 mpscan *_sc = (sc); \
404 MPSCAN_INITX(_sc, _mm->v, _mm->vl); \
407 /* --- @mp_rscan@ --- *
409 * Arguments: @mpscan *sc@ = pointer to bitscanner block
410 * @const mp *m@ = pointer to a multiprecision integer
414 * Use: Initializes a reverse bitscanner on a multiprecision
418 extern void mp_rscan(mpscan */*sc*/, const mp */*m*/);
420 #define MP_RSCAN(sc, m) do { \
421 const mp *_mm = (m); \
422 mpscan *_sc = (sc); \
423 MPSCAN_RINITX(_sc, _mm->v, _mm->vl); \
426 /* --- Other bitscanning aliases --- */
428 #define mp_step mpscan_step
429 #define mp_bit mpscan_bit
430 #define mp_rstep mpscan_rstep
431 #define mp_rbit mpscan_rbit
433 #define MP_STEP MPSCAN_STEP
434 #define MP_BIT MPSCAN_BIT
435 #define MP_RSTEP MPSCAN_RSTEP
436 #define MP_RBIT MPSCAN_RBIT
438 /*----- Loading and storing -----------------------------------------------*/
440 /* --- @mp_octets@ --- *
442 * Arguments: @const mp *m@ = a multiprecision integer
444 * Returns: The number of octets required to represent @m@.
446 * Use: Calculates the external storage required for a multiprecision
450 extern size_t mp_octets(const mp */*m*/);
452 /* --- @mp_octets2c@ --- *
454 * Arguments: @const mp *m@ = a multiprecision integer
456 * Returns: The number of octets required to represent @m@.
458 * Use: Calculates the external storage required for a multiprecision
459 * integer represented as two's complement.
462 extern size_t mp_octets2c(const mp */*m*/);
464 /* --- @mp_bits@ --- *
466 * Arguments: @const mp *m@ = a multiprecision integer
468 * Returns: The number of bits required to represent @m@.
470 * Use: Calculates the external storage required for a multiprecision
474 extern unsigned long mp_bits(const mp */*m*/);
476 /* --- @mp_loadl@ --- *
478 * Arguments: @mp *d@ = destination
479 * @const void *pv@ = pointer to source data
480 * @size_t sz@ = size of the source data
482 * Returns: Resulting multiprecision number.
484 * Use: Loads a multiprecision number from an array of octets. The
485 * first byte in the array is the least significant. More
486 * formally, if the bytes are %$b_0, b_1, \ldots, b_{n-1}$%
487 * then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%.
490 extern mp *mp_loadl(mp */*d*/, const void */*pv*/, size_t /*sz*/);
492 /* --- @mp_storel@ --- *
494 * Arguments: @const mp *m@ = source
495 * @void *pv@ = pointer to output array
496 * @size_t sz@ = size of the output array
500 * Use: Stores a multiprecision number in an array of octets. The
501 * first byte in the array is the least significant. If the
502 * array is too small to represent the number, high-order bits
503 * are truncated; if the array is too large, high order bytes
504 * are filled with zeros. More formally, if the number is
505 * %$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%,
506 * then the array is %$b_0, b_1, \ldots, b_{n-1}$%.
509 extern void mp_storel(const mp */*m*/, void */*pv*/, size_t /*sz*/);
511 /* --- @mp_loadb@ --- *
513 * Arguments: @mp *d@ = destination
514 * @const void *pv@ = pointer to source data
515 * @size_t sz@ = size of the source data
517 * Returns: Resulting multiprecision number.
519 * Use: Loads a multiprecision number from an array of octets. The
520 * last byte in the array is the least significant. More
521 * formally, if the bytes are %$b_{n-1}, b_{n-2}, \ldots, b_0$%
522 * then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%.
525 extern mp *mp_loadb(mp */*d*/, const void */*pv*/, size_t /*sz*/);
527 /* --- @mp_storeb@ --- *
529 * Arguments: @const mp *m@ = source
530 * @void *pv@ = pointer to output array
531 * @size_t sz@ = size of the output array
535 * Use: Stores a multiprecision number in an array of octets. The
536 * last byte in the array is the least significant. If the
537 * array is too small to represent the number, high-order bits
538 * are truncated; if the array is too large, high order bytes
539 * are filled with zeros. More formally, if the number is
540 * %$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%,
541 * then the array is %$b_{n-1}, b_{n-2}, \ldots, b_0$%.
544 extern void mp_storeb(const mp */*m*/, void */*pv*/, size_t /*sz*/);
546 /* --- @mp_loadl2c@ --- *
548 * Arguments: @mp *d@ = destination
549 * @const void *pv@ = pointer to source data
550 * @size_t sz@ = size of the source data
552 * Returns: Resulting multiprecision number.
554 * Use: Loads a multiprecision number from an array of octets as
555 * two's complement. The first byte in the array is the least
559 extern mp *mp_loadl2c(mp */*d*/, const void */*pv*/, size_t /*sz*/);
561 /* --- @mp_storel2c@ --- *
563 * Arguments: @const mp *m@ = source
564 * @void *pv@ = pointer to output array
565 * @size_t sz@ = size of the output array
569 * Use: Stores a multiprecision number in an array of octets as two's
570 * complement. The first byte in the array is the least
571 * significant. If the array is too small to represent the
572 * number, high-order bits are truncated; if the array is too
573 * large, high order bytes are sign-extended.
576 extern void mp_storel2c(const mp */*m*/, void */*pv*/, size_t /*sz*/);
578 /* --- @mp_loadb2c@ --- *
580 * Arguments: @mp *d@ = destination
581 * @const void *pv@ = pointer to source data
582 * @size_t sz@ = size of the source data
584 * Returns: Resulting multiprecision number.
586 * Use: Loads a multiprecision number from an array of octets as
587 * two's complement. The last byte in the array is the least
591 extern mp *mp_loadb2c(mp */*d*/, const void */*pv*/, size_t /*sz*/);
593 /* --- @mp_storeb2c@ --- *
595 * Arguments: @const mp *m@ = source
596 * @void *pv@ = pointer to output array
597 * @size_t sz@ = size of the output array
601 * Use: Stores a multiprecision number in an array of octets, as
602 * two's complement. The last byte in the array is the least
603 * significant. If the array is too small to represent the
604 * number, high-order bits are truncated; if the array is too
605 * large, high order bytes are sign-extended.
608 extern void mp_storeb2c(const mp */*m*/, void */*pv*/, size_t /*sz*/);
610 /*----- Bit operations ----------------------------------------------------*/
612 /* --- @mp_not@ --- *
614 * Arguments: @mp *d@ = destination
617 * Returns: The bitwise complement of the source.
620 extern mp *mp_not(mp */*d*/, mp */*a*/);
622 /* --- @mp_bitop@ --- *
624 * Arguments: @mp *d@ = destination
625 * @mp *a, *b@ = sources
627 * Returns: The result of the given bitwise operation. These functions
628 * don't handle negative numbers at all sensibly. For that, use
629 * the @...2c@ variants. The functions are named after the
630 * truth tables they generate:
637 #define MP_BITDECL(string) \
638 extern mp *mp_bit##string(mp */*d*/, mp */*a*/, mp */*b*/);
639 MPX_DOBIN(MP_BITDECL)
641 /* --- @mp_[n]and@, @mp_[n]or@, @mp_[n]xor@, @mp_not@ --- *
643 * Synonyms for the commonly-used functions.
646 #define mp_and mp_bit0001
647 #define mp_or mp_bit0111
648 #define mp_nand mp_bit1110
649 #define mp_nor mp_bit1000
650 #define mp_xor mp_bit0110
652 /* --- @mp_testbit@ --- *
654 * Arguments: @mp *x@ = a large integer
655 * @unsigned long n@ = which bit to test
657 * Returns: Nonzero if the bit is set, zero if not.
660 extern int mp_testbit(mp */*x*/, unsigned long /*n*/);
662 /* --- @mp_setbit@, @mp_clearbit@ --- *
664 * Arguments: @mp *d@ = a destination
665 * @mp *x@ = a large integer
666 * @unsigned long n@ = which bit to modify
668 * Returns: The argument @x@, with the appropriate bit set or cleared.
671 extern mp *mp_setbit(mp */*d*/, mp */*x*/, unsigned long /*n*/);
672 extern mp *mp_clearbit(mp */*d*/, mp */*x*/, unsigned long /*n*/);
674 /* --- @mp_lsl@, @mp_lslc@, @mp_lsr@ --- *
676 * Arguments: @mp *d@ = destination
678 * @size_t n@ = number of bits to move
680 * Returns: Result, @a@ shifted left or right by @n@.
682 * Use: Bitwise shift operators. @mp_lslc@ fills the bits introduced
683 * on the right with ones instead of zeroes: it's used
684 * internally by @mp_lsl2c@, though it may be useful on its
688 extern mp *mp_lsl(mp */*d*/, mp */*a*/, size_t /*n*/);
689 extern mp *mp_lslc(mp */*d*/, mp */*a*/, size_t /*n*/);
690 extern mp *mp_lsr(mp */*d*/, mp */*a*/, size_t /*n*/);
692 /* --- @mp_not2c@ --- *
694 * Arguments: @mp *d@ = destination
697 * Returns: The sign-extended complement of the argument.
700 extern mp *mp_not2c(mp */*d*/, mp */*a*/);
702 /* --- @mp_bitop2c@ --- *
704 * Arguments: @mp *d@ = destination
705 * @mp *a, *b@ = sources
707 * Returns: The result of the given bitwise operation. Negative numbers
708 * are treated as two's complement, sign-extended infinitely to
709 * the left. The functions are named after the truth tables
717 #define MP_BIT2CDECL(string) \
718 extern mp *mp_bit##string##2c(mp */*d*/, mp */*a*/, mp */*b*/);
719 MPX_DOBIN(MP_BIT2CDECL)
721 /* --- @mp_[n]and@, @mp_[n]or@, @mp_[n]xor@, @mp_not@ --- *
723 * Synonyms for the commonly-used functions.
726 #define mp_and2c mp_bit00012c
727 #define mp_or2c mp_bit01112c
728 #define mp_nand2c mp_bit11102c
729 #define mp_nor2c mp_bit10002c
730 #define mp_xor2c mp_bit01102c
732 /* --- @mp_lsl2c@, @mp_lsr2c@ --- *
734 * Arguments: @mp *d@ = destination
736 * @size_t n@ = number of bits to move
738 * Returns: Result, @a@ shifted left or right by @n@. Handles the
739 * pretence of sign-extension for negative numbers.
742 extern mp *mp_lsl2c(mp */*d*/, mp */*a*/, size_t /*n*/);
743 extern mp *mp_lsr2c(mp */*d*/, mp */*a*/, size_t /*n*/);
745 /* --- @mp_testbit2c@ --- *
747 * Arguments: @mp *x@ = a large integer
748 * @unsigned long n@ = which bit to test
750 * Returns: Nonzero if the bit is set, zero if not. Fakes up two's
751 * complement representation.
754 extern int mp_testbit2c(mp */*x*/, unsigned long /*n*/);
756 /* --- @mp_setbit2c@, @mp_clearbit2c@ --- *
758 * Arguments: @mp *d@ = a destination
759 * @mp *x@ = a large integer
760 * @unsigned long n@ = which bit to modify
762 * Returns: The argument @x@, with the appropriate bit set or cleared.
763 * Fakes up two's complement representation.
766 extern mp *mp_setbit2c(mp */*d*/, mp */*x*/, unsigned long /*n*/);
767 extern mp *mp_clearbit2c(mp */*d*/, mp */*x*/, unsigned long /*n*/);
769 /*----- Comparisons -------------------------------------------------------*/
773 * Arguments: @const mp *a, *b@ = two numbers
775 * Returns: Nonzero if the numbers are equal.
778 extern int mp_eq(const mp */*a*/, const mp */*b*/);
780 #define MP_EQ(a, b) \
781 ((((a)->f ^ (b)->f) & MP_NEG) == 0 && \
782 mpx_ueq((a)->v, (a)->vl, (b)->v, (b)->vl))
784 /* --- @mp_cmp@ --- *
786 * Arguments: @const mp *a, *b@ = two numbers
788 * Returns: Less than, equal to or greater than zero, according to
789 * whether @a@ is less than, equal to or greater than @b@.
792 extern int mp_cmp(const mp */*a*/, const mp */*b*/);
794 #define MP_CMP(a, op, b) (mp_cmp((a), (b)) op 0)
796 /* --- Other handy macros --- */
798 #define MP_NEGP(x) ((x)->f & MP_NEG)
799 #define MP_ZEROP(x) (!MP_LEN(x))
800 #define MP_POSP(x) (!MP_NEGP(x) && !MP_ZEROP(x))
801 #define MP_ODDP(x) (!MP_ZEROP(x) && ((x)->v[0] & 1u))
802 #define MP_EVENP(x) (!MP_ODDP(x))
804 /*----- Arithmetic operations ---------------------------------------------*/
806 /* --- @mp_neg@ --- *
808 * Arguments: @mp *d@ = destination
811 * Returns: The negation of the argument.
813 * Use: Negates its argument.
816 extern mp *mp_neg(mp */*d*/, mp */*a*/);
818 /* --- @mp_add@ --- *
820 * Arguments: @mp *d@ = destination
821 * @mp *a, *b@ = sources
823 * Returns: Result, @a@ added to @b@.
826 extern mp *mp_add(mp */*d*/, mp */*a*/, mp */*b*/);
828 /* --- @mp_sub@ --- *
830 * Arguments: @mp *d@ = destination
831 * @mp *a, *b@ = sources
833 * Returns: Result, @b@ subtracted from @a@.
836 extern mp *mp_sub(mp */*d*/, mp */*a*/, mp */*b*/);
838 /* --- @mp_mul@ --- *
840 * Arguments: @mp *d@ = destination
841 * @mp *a, *b@ = sources
843 * Returns: Result, @a@ multiplied by @b@.
846 extern mp *mp_mul(mp */*d*/, mp */*a*/, mp */*b*/);
848 /* --- @mp_sqr@ --- *
850 * Arguments: @mp *d@ = destination
853 * Returns: Result, @a@ squared.
856 extern mp *mp_sqr(mp */*d*/, mp */*a*/);
858 /* --- @mp_div@ --- *
860 * Arguments: @mp **qq, **rr@ = destination, quotient and remainder
861 * @mp *a, *b@ = sources
863 * Use: Calculates the quotient and remainder when @a@ is divided by
867 extern void mp_div(mp **/*qq*/, mp **/*rr*/, mp */*a*/, mp */*b*/);
869 /* --- @mp_exp@ --- *
871 * Arguments: @mp *d@ = fake destination
875 * Returns: Result, %$a^e$%.
878 extern mp *mp_exp(mp */*d*/, mp */*a*/, mp */*e*/);
880 /* --- @mp_odd@ --- *
882 * Arguments: @mp *d@ = pointer to destination integer
883 * @mp *m@ = pointer to source integer
884 * @size_t *s@ = where to store the power of 2
886 * Returns: An odd integer integer %$t$% such that %$m = 2^s t$%.
888 * Use: Computes a power of two and an odd integer which, when
889 * multiplied, give a specified result. This sort of thing is
890 * useful in number theory quite often.
893 extern mp *mp_odd(mp */*d*/, mp */*m*/, size_t */*s*/);
895 /* --- @mp_leastcongruent@ --- *
897 * Arguments: @mp *d@ = pointer to destination
898 * @mp *b@ = lower bound
899 * @mp *r@ = representative
902 * Returns: The smallest integer %$x \equiv r \pmod{m}$% such that
906 extern mp *mp_leastcongruent(mp */*d*/, mp */*b*/, mp */*r*/, mp */*m*/);
908 /*----- More advanced algorithms ------------------------------------------*/
910 /* --- @mp_sqrt@ --- *
912 * Arguments: @mp *d@ = pointer to destination integer
913 * @mp *a@ = (nonnegative) integer to take square root of
915 * Returns: The largest integer %$x$% such that %$x^2 \le a$%.
917 * Use: Computes integer square roots.
919 * The current implementation isn't very good: it uses the
920 * Newton-Raphson method to find an approximation to %$a$%. If
921 * there's any demand for a better version, I'll write one.
924 extern mp *mp_sqrt(mp */*d*/, mp */*a*/);
926 /* --- @mp_gcd@ --- *
928 * Arguments: @mp **gcd, **xx, **yy@ = where to write the results
929 * @mp *a, *b@ = sources (must be nonzero)
933 * Use: Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that
934 * @ax + by = gcd(a, b)@. This is useful for computing modular
935 * inverses. Neither @a@ nor @b@ may be zero.
938 extern void mp_gcd(mp **/*gcd*/, mp **/*xx*/, mp **/*yy*/,
939 mp */*a*/, mp */*b*/);
941 /* -- @mp_modinv@ --- *
943 * Arguments: @mp *d@ = destination
947 * Returns: The inverse %$x^{-1} \bmod p$%.
949 * Use: Computes a modular inverse. An assertion fails if %$p$%
953 extern mp *mp_modinv(mp */*d*/, mp */*x*/, mp */*p*/);
955 /* --- @mp_jacobi@ --- *
957 * Arguments: @mp *a@ = an integer
958 * @mp *n@ = another integer
960 * Returns: @-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
962 * Use: Computes the Kronecker symbol %$\jacobi{a}{n}$%. If @n@ is
963 * prime, this is the Legendre symbol and is equal to 1 if and
964 * only if @a@ is a quadratic residue mod @n@. The result is
965 * zero if and only if @a@ and @n@ have a common factor greater
968 * If @n@ is composite, then this computes the Kronecker symbol
970 * %$\jacobi{a}{n}=\jacobi{a}{u}\prod_i\jacobi{a}{p_i}^{e_i}$%
972 * where %$n = u p_0^{e_0} \ldots p_{n-1}^{e_{n-1}}$% is the
973 * prime factorization of %$n$%. The missing bits are:
975 * * %$\jacobi{a}{1} = 1$%;
976 * * %$\jacobi{a}{-1} = 1$% if @a@ is negative, or 1 if
978 * * %$\jacobi{a}{0} = 0$%;
979 * * %$\jacobi{a}{2}$ is 0 if @a@ is even, 1 if @a@ is
980 * congruent to 1 or 7 (mod 8), or %$-1$% otherwise.
982 * If %$n$% is positive and odd, then this is the Jacobi
983 * symbol. (The Kronecker symbol is a consistent domain
984 * extension; the Jacobi symbol was implemented first, and the
988 extern int mp_jacobi(mp */*a*/, mp */*n*/);
990 /* --- @mp_modsqrt@ --- *
992 * Arguments: @mp *d@ = destination integer
993 * @mp *a@ = source integer
994 * @mp *p@ = modulus (must be prime)
996 * Returns: If %$a$% is a quadratic residue, a square root of %$a$%; else
999 * Use: Returns an integer %$x$% such that %$x^2 \equiv a \pmod{p}$%,
1000 * if one exists; else a null pointer. This function will not
1001 * work if %$p$% is composite: you must factor the modulus, take
1002 * a square root mod each factor, and recombine the results
1003 * using the Chinese Remainder Theorem.
1005 * We guarantee that the square root returned is the smallest
1006 * one (i.e., the `positive' square root).
1009 extern mp *mp_modsqrt(mp */*d*/, mp */*a*/, mp */*p*/);
1011 /* --- @mp_modexp@ --- *
1013 * Arguments: @mp *d@ = fake destination
1014 * @mp *x@ = base of exponentiation
1015 * @mp *e@ = exponent
1016 * @mp *n@ = modulus (must be positive)
1018 * Returns: The value %$x^e \bmod n$%.
1021 extern mp *mp_modexp(mp */*d*/, mp */*x*/, mp */*e*/, mp */*n*/);
1023 /*----- Test harness support ----------------------------------------------*/
1025 #include <mLib/testrig.h>
1027 #ifndef CATACOMB_MPTEXT_H
1028 # include "mptext.h"
1031 extern const test_type type_mp;
1033 /*----- That's all, folks -------------------------------------------------*/