3 * $Id: mpx.h,v 1.17 2004/03/27 00:04:46 mdw Exp $
5 * Low level 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.17 2004/03/27 00:04:46 mdw
34 * Implement efficient reduction for pleasant-looking primes.
36 * Revision 1.16 2003/05/16 09:09:24 mdw
37 * Fix @mp_lsl2c@. Turns out to be surprisingly tricky.
39 * Revision 1.15 2002/10/19 17:56:50 mdw
40 * Fix bit operations. Test them (a bit) better.
42 * Revision 1.14 2002/10/09 00:36:03 mdw
43 * Fix bounds on workspace for Karatsuba operations.
45 * Revision 1.13 2002/10/06 22:52:50 mdw
46 * Pile of changes for supporting two's complement properly.
48 * Revision 1.12 2001/04/03 19:36:05 mdw
49 * Add some simple bitwise operations so that Perl can use them.
51 * Revision 1.11 2000/10/08 15:48:35 mdw
52 * Rename Karatsuba constants now that we have @gfx_kmul@ too.
54 * Revision 1.10 2000/10/08 12:06:12 mdw
55 * Provide @mpx_ueq@ for rapidly testing equality of two integers.
57 * Revision 1.9 1999/12/22 15:49:07 mdw
58 * New function for division by a small integer.
60 * Revision 1.8 1999/12/11 10:57:43 mdw
61 * Karatsuba squaring algorithm.
63 * Revision 1.7 1999/12/11 01:51:28 mdw
64 * Change Karatsuba parameters slightly.
66 * Revision 1.6 1999/12/10 23:23:51 mdw
67 * Karatsuba-Ofman multiplication algorithm.
69 * Revision 1.5 1999/11/20 22:23:27 mdw
70 * Add function versions of some low-level macros with wider use.
72 * Revision 1.4 1999/11/17 18:04:43 mdw
73 * Add two's complement support. Fix a bug in MPX_UMLAN.
75 * Revision 1.3 1999/11/13 01:51:29 mdw
76 * Minor interface changes. Should be stable now.
78 * Revision 1.2 1999/11/11 17:47:55 mdw
79 * Minor changes for different `mptypes.h' format.
81 * Revision 1.1 1999/09/03 08:41:12 mdw
86 #ifndef CATACOMB_MPX_H
87 #define CATACOMB_MPX_H
93 /*----- The idea ----------------------------------------------------------*
95 * This file provides functions and macros which work on vectors of words as
96 * unsigned multiprecision integers. The interface works in terms of base
97 * and limit pointers (i.e., a pointer to the start of a vector, and a
98 * pointer just past its end) rather than base pointer and length, because
99 * that requires more arithmetic and state to work on.
101 * The interfaces are slightly bizarre in other ways. Try to use the
102 * higher-level functions where you can: they're rather better designed to
103 * actually be friendly and useful.
106 /*----- Header files ------------------------------------------------------*/
110 #ifndef CATACOMB_MPW_H
114 /*----- General manipulation ----------------------------------------------*/
116 /* --- @MPX_SHRINK@ --- *
118 * Arguments: @const mpw *v@ = pointer to vector of words
119 * @const mpw *vl@ = (updated) current limit of vector
121 * Use: Shrinks down the limit of a multiprecision integer vector.
124 #define MPX_SHRINK(v, vl) do { \
125 const mpw *_vv = (v), *_vvl = (vl); \
126 while (_vvl > _vv && !_vvl[-1]) \
128 (vl) = (mpw *)_vvl; \
131 /* --- @MPX_BITS@ --- *
133 * Arguments: @unsigned long b@ = result variable
134 * @const mpw *v@ = pointer to array of words
135 * @const mpw *vl@ = limit of vector (from @MPX_SHRINK@)
137 * Use: Calculates the number of bits in a multiprecision value.
140 #define MPX_BITS(b, v, vl) do { \
141 const mpw *_v = (v), *_vl = (vl); \
142 MPX_SHRINK(_v, _vl); \
146 unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1; \
148 unsigned _k = MPW_BITS / 2; \
160 /* --- @MPX_OCTETS@ --- *
162 * Arguments: @size_t o@ = result variable
163 * @const mpw *v, *vl@ = pointer to array of words
165 * Use: Calculates the number of octets in a multiprecision value.
168 #define MPX_OCTETS(o, v, vl) do { \
170 MPX_BITS(_bb, (v), (vl)); \
171 (o) = (_bb + 7) >> 3; \
174 /* --- @MPX_OCTETS2C@ --- *
176 * Arguments: @size_t o@ = result variable
177 * @const mpw *v, *vl@ = pointer to array of words
179 * Use: Calculates the number of octets in a multiprecision value, if
180 * you represent it as two's complement.
183 #define MPX_OCTETS2C(o, v, vl) do { \
185 MPX_BITS(_bb, (v), (vl)); \
186 (o) = (_bb >> 3) + 1; \
189 /* --- @MPX_COPY@ --- *
191 * Arguments: @dv, dvl@ = destination vector base and limit
192 * @av, avl@ = source vector base and limit
194 * Use: Copies a multiprecision integer.
197 #define MPX_COPY(dv, dvl, av, avl) do { \
198 mpw *_dv = (dv), *_dvl = (dvl); \
199 size_t _dn = _dvl - _dv; \
200 const mpw *_av = (av), *_avl = (avl); \
201 size_t _an = _avl - _av; \
204 memset(_dv, 0, MPWS(_dn - _an)); \
205 } else if (_an >= _dn) \
206 memmove(_dv, _av, MPWS(_dn)); \
208 memmove(_dv, _av, MPWS(_an)); \
209 memset(_dv + _an, 0, MPWS(_dn - _an)); \
213 /* --- @MPX_ZERO@ --- *
215 * Arguments: @v, vl@ = base and limit of vector to clear
217 * Use: Zeroes the area between the two vector pointers.
220 #define MPX_ZERO(v, vl) do { \
221 mpw *_v = (v), *_vl = (vl); \
223 memset(_v, 0, MPWS(_vl - _v)); \
226 /* --- @MPX_ONE@ --- *
228 * Arguments: @v, vl@ = base and limit of vector to clear
230 * Use: Fills the area between the two vector pointers with ones.
233 #define MPX_ONE(v, vl) do { \
235 const mpw *_vl = (vl); \
240 /*----- Loading and storing -----------------------------------------------*/
242 /* --- @mpx_storel@ --- *
244 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
245 * @void *p@ = pointer to octet array
246 * @size_t sz@ = size of octet array
250 * Use: Stores an MP in an octet array, least significant octet
251 * first. High-end octets are silently discarded if there
252 * isn't enough space for them.
255 extern void mpx_storel(const mpw */
*v*/
, const mpw */
*vl*/
,
256 void */
*p*/
, size_t /*sz*/);
258 /* --- @mpx_loadl@ --- *
260 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
261 * @const void *p@ = pointer to octet array
262 * @size_t sz@ = size of octet array
266 * Use: Loads an MP in an octet array, least significant octet
267 * first. High-end octets are ignored if there isn't enough
271 extern void mpx_loadl(mpw */
*v*/
, mpw */
*vl*/
,
272 const void */
*p*/
, size_t /*sz*/);
274 /* --- @mpx_storeb@ --- *
276 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
277 * @void *p@ = pointer to octet array
278 * @size_t sz@ = size of octet array
282 * Use: Stores an MP in an octet array, most significant octet
283 * first. High-end octets are silently discarded if there
284 * isn't enough space for them.
287 extern void mpx_storeb(const mpw */
*v*/
, const mpw */
*vl*/
,
288 void */
*p*/
, size_t /*sz*/);
290 /* --- @mpx_loadb@ --- *
292 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
293 * @const void *p@ = pointer to octet array
294 * @size_t sz@ = size of octet array
298 * Use: Loads an MP in an octet array, most significant octet
299 * first. High-end octets are ignored if there isn't enough
303 extern void mpx_loadb(mpw */
*v*/
, mpw */
*vl*/
,
304 const void */
*p*/
, size_t /*sz*/);
306 /* --- @mpx_storel2cn@ --- *
308 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
309 * @void *pp@ = pointer to octet array
310 * @size_t sz@ = size of octet array
314 * Use: Stores a negative MP in an octet array, least significant
315 * octet first, as two's complement. High-end octets are
316 * silently discarded if there isn't enough space for them.
317 * This obviously makes the output bad.
320 extern void mpx_storel2cn(const mpw */
*v*/
, const mpw */
*vl*/
,
321 void */
*p*/
, size_t /*sz*/);
323 /* --- @mpx_loadl2cn@ --- *
325 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
326 * @const void *pp@ = pointer to octet array
327 * @size_t sz@ = size of octet array
331 * Use: Loads a negative MP in an octet array, least significant
332 * octet first, as two's complement. High-end octets are
333 * ignored if there isn't enough space for them. This probably
334 * means you made the wrong choice coming here.
337 extern void mpx_loadl2cn(mpw */
*v*/
, mpw */
*vl*/
,
338 const void */
*p*/
, size_t /*sz*/);
340 /* --- @mpx_storeb2cn@ --- *
342 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
343 * @void *pp@ = pointer to octet array
344 * @size_t sz@ = size of octet array
348 * Use: Stores a negative MP in an octet array, most significant
349 * octet first, as two's complement. High-end octets are
350 * silently discarded if there isn't enough space for them,
351 * which probably isn't what you meant.
354 extern void mpx_storeb2cn(const mpw */
*v*/
, const mpw */
*vl*/
,
355 void */
*p*/
, size_t /*sz*/);
357 /* --- @mpx_loadb2cn@ --- *
359 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
360 * @const void *pp@ = pointer to octet array
361 * @size_t sz@ = size of octet array
365 * Use: Loads a negative MP in an octet array, most significant octet
366 * first as two's complement. High-end octets are ignored if
367 * there isn't enough space for them. This probably means you
368 * chose this function wrongly.
371 extern void mpx_loadb2cn(mpw */
*v*/
, mpw */
*vl*/
,
372 const void */
*p*/
, size_t /*sz*/);
375 /*----- Logical shifting --------------------------------------------------*/
377 /* --- @mpx_lsl@ --- *
379 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
380 * @const mpw *av, *avl@ = source vector base and limit
381 * @size_t n@ = number of bit positions to shift by
385 * Use: Performs a logical shift left operation on an integer.
388 extern void mpx_lsl(mpw */
*dv*/
, mpw */
*dvl*/
,
389 const mpw */
*av*/
, const mpw */
*avl*/
,
392 /* --- @mpx_lslc@ --- *
394 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
395 * @const mpw *av, *avl@ = source vector base and limit
396 * @size_t n@ = number of bit positions to shift by
400 * Use: Performs a logical shift left operation on an integer, only
401 * it fills in the bits with ones instead of zeroes.
404 extern void mpx_lslc(mpw */
*dv*/
, mpw */
*dvl*/
,
405 const mpw */
*av*/
, const mpw */
*avl*/
,
408 /* --- @mpx_lsr@ --- *
410 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
411 * @const mpw *av, *avl@ = source vector base and limit
412 * @size_t n@ = number of bit positions to shift by
416 * Use: Performs a logical shift right operation on an integer.
419 extern void mpx_lsr(mpw */
*dv*/
, mpw */
*dvl*/
,
420 const mpw */
*av*/
, const mpw */
*avl*/
,
423 /*----- Bitwise operations ------------------------------------------------*/
425 /* --- @mpx_bitop@ --- *
427 * Arguments: @mpw *dv, *dvl@ = destination vector
428 * @const mpw *av, *avl@ = first source vector
429 * @const mpw *bv, *bvl@ = second source vector
433 * Use: Provide the dyadic boolean functions. The functions are
434 * named after the truth table they generate:
441 #define MPX_DOBIN(what) \
442 what(0000) what(0001) what(0010) what(0011) \
443 what(0100) what(0101) what(0110) what(0111) \
444 what(1000) what(1001) what(1010) what(1011) \
445 what(1100) what(1101) what(1110) what(1111)
447 #define MPX_BITDECL(string) \
448 extern void mpx_bit##string(mpw */*dv*/, mpw */*dvl*/, \
449 const mpw */*av*/, const mpw */*avl*/, \
450 const mpw */*bv*/, const mpw */*bvl*/);
451 MPX_DOBIN(MPX_BITDECL
)
453 /* --- @mpx_[n]and@, @mpx_[n]or@, @mpx_xor@ --- *
455 * Synonyms for the commonly-used functions above.
458 #define mpx_and mpx_bit0001
459 #define mpx_or mpx_bit0111
460 #define mpx_nand mpx_bit1110
461 #define mpx_nor mpx_bit1000
462 #define mpx_xor mpx_bit0110
464 /* --- @mpx_not@ --- *
466 * Arguments: @mpw *dv, *dvl@ = destination vector
467 * @const mpw *av, *avl@ = first source vector
474 extern void mpx_not(mpw */
*dv*/
, mpw */
*dvl*/
,
475 const mpw */
*av*/
, const mpw */
*avl*/
);
477 /*----- Unsigned arithmetic -----------------------------------------------*/
479 /* --- @mpx_2c@ --- *
481 * Arguments: @mpw *dv, *dvl@ = destination vector
482 * @const mpw *v, *vl@ = source vector
486 * Use: Calculates the two's complement of @v@.
489 extern void mpx_2c(mpw */
*dv*/
, mpw */
*dvl*/
,
490 const mpw */
*v*/
, const mpw */
*vl*/
);
492 /* --- @mpx_ueq@ --- *
494 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
495 * @const mpw *bv, *bvl@ = second argument vector base and limit
497 * Returns: Nonzero if the two vectors are equal.
499 * Use: Performs an unsigned integer test for equality.
502 extern int mpx_ueq(const mpw */
*av*/
, const mpw */
*avl*/
,
503 const mpw */
*bv*/
, const mpw */
*bvl*/
);
505 /* --- @mpx_ucmp@ --- *
507 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
508 * @const mpw *bv, *bvl@ = second argument vector base and limit
510 * Returns: Less than, equal to, or greater than zero depending on
511 * whether @a@ is less than, equal to or greater than @b@,
514 * Use: Performs an unsigned integer comparison.
517 #define MPX_UCMP(av, avl, op, dv, dvl) \
518 (mpx_ucmp((av), (avl), (dv), (dvl)) op 0)
520 extern int mpx_ucmp(const mpw */
*av*/
, const mpw */
*avl*/
,
521 const mpw */
*bv*/
, const mpw */
*bvl*/
);
523 /* --- @mpx_uadd@ --- *
525 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
526 * @const mpw *av, *avl@ = first addend vector base and limit
527 * @const mpw *bv, *bvl@ = second addend vector base and limit
531 * Use: Performs unsigned integer addition. If the result overflows
532 * the destination vector, high-order bits are discarded. This
533 * means that two's complement addition happens more or less for
534 * free, although that's more a side-effect than anything else.
535 * The result vector may be equal to either or both source
536 * vectors, but may not otherwise overlap them.
539 extern void mpx_uadd(mpw */
*dv*/
, mpw */
*dvl*/
,
540 const mpw */
*av*/
, const mpw */
*avl*/
,
541 const mpw */
*bv*/
, const mpw */
*bvl*/
);
543 /* --- @mpx_uaddn@ --- *
545 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
546 * @mpw n@ = other addend
550 * Use: Adds a small integer to a multiprecision number.
553 #define MPX_UADDN(dv, dvl, n) do { \
554 mpw *_ddv = (dv), *_ddvl = (dvl); \
557 while (_c && _ddv < _ddvl) { \
558 mpd _x = (mpd)*_ddv + (mpd)_c; \
560 _c = _x >> MPW_BITS; \
564 extern void mpx_uaddn(mpw */
*dv*/
, mpw */
*dvl*/
, mpw
/*n*/);
566 /* --- @mpx_uaddnlsl@ --- *
568 * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
569 * @mpw a@ = second argument
570 * @unsigned o@ = offset in bits
574 * Use: Computes %$d + 2^o a$%. If the result overflows then
575 * high-order bits are discarded, as usual. We must have
576 * @0 < o < MPW_BITS@.
579 extern void mpx_uaddnlsl(mpw */
*dv*/
, mpw */
*dvl*/
,
580 mpw
/*a*/, unsigned /*o*/);
582 /* --- @mpx_usub@ --- *
584 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
585 * @const mpw *av, *avl@ = first argument vector base and limit
586 * @const mpw *bv, *bvl@ = second argument vector base and limit
590 * Use: Performs unsigned integer subtraction. If the result
591 * overflows the destination vector, high-order bits are
592 * discarded. This means that two's complement subtraction
593 * happens more or less for free, although that's more a side-
594 * effect than anything else. The result vector may be equal to
595 * either or both source vectors, but may not otherwise overlap
599 extern void mpx_usub(mpw */
*dv*/
, mpw */
*dvl*/
,
600 const mpw */
*av*/
, const mpw */
*avl*/
,
601 const mpw */
*bv*/
, const mpw */
*bvl*/
);
603 /* --- @mpx_usubn@ --- *
605 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
610 * Use: Subtracts a small integer from a multiprecision number.
613 #define MPX_USUBN(dv, dvl, n) do { \
614 mpw *_ddv = (dv), *_ddvl = (dvl); \
617 while (_ddv < _ddvl) { \
618 mpd _x = (mpd)*_ddv - (mpd)_c; \
620 if (_x >> MPW_BITS) \
627 extern void mpx_usubn(mpw */
*dv*/
, mpw */
*dvl*/
, mpw
/*n*/);
629 /* --- @mpx_usubnlsl@ --- *
631 * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
632 * @mpw a@ = second argument
633 * @unsigned o@ = offset in bits
637 * Use: Computes %$d - 2^o a$%. If the result overflows then
638 * high-order bits are discarded, as usual, so you get two's
639 * complement. Which might be what you wanted... We must have
640 * @0 < o < MPW_BITS@.
643 extern void mpx_usubnlsl(mpw */
*dv*/
, mpw */
*dvl*/
,
644 mpw
/*a*/, unsigned /*o*/);
646 /* --- @mpx_umul@ --- *
648 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
649 * @const mpw *av, *avl@ = multiplicand vector base and limit
650 * @const mpw *bv, *bvl@ = multiplier vector base and limit
654 * Use: Performs unsigned integer multiplication. If the result
655 * overflows the desination vector, high-order bits are
656 * discarded. The result vector may not overlap the argument
657 * vectors in any way.
660 extern void mpx_umul(mpw */
*dv*/
, mpw */
*dvl*/
,
661 const mpw */
*av*/
, const mpw */
*avl*/
,
662 const mpw */
*bv*/
, const mpw */
*bvl*/
);
664 /* --- @mpx_umuln@ --- *
666 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
667 * @const mpw *av, *avl@ = multiplicand vector base and limit
668 * @mpw m@ = multiplier
672 * Use: Multiplies a multiprecision integer by a single-word value.
673 * The destination and source may be equal. The destination
674 * is completely cleared after use.
677 #define MPX_UMULN(dv, dvl, av, avl, m) do { \
678 mpw *_dv = (dv), *_dvl = (dvl); \
679 const mpw *_av = (av), *_avl = (avl); \
683 while (_av < _avl) { \
687 _x = (mpd)_m * (mpd)*_av++ + _c; \
689 _c = _x >> MPW_BITS; \
693 MPX_ZERO(_dv, _dvl); \
697 extern void mpx_umuln(mpw */
*dv*/
, mpw */
*dvl*/
,
698 const mpw */
*av*/
, const mpw */
*avl*/
, mpw m
);
700 /* --- @mpx_umlan@ --- *
702 * Arguments: @mpw *dv, *dvl@ = destination/accumulator base and limit
703 * @const mpw *av, *avl@ = multiplicand vector base and limit
704 * @mpw m@ = multiplier
708 * Use: Multiplies a multiprecision integer by a single-word value
709 * and adds the result to an accumulator.
712 #define MPX_UMLAN(dv, dvl, av, avl, m) do { \
713 mpw *_dv = (dv), *_dvl = (dvl); \
714 const mpw *_av = (av), *_avl = (avl); \
718 while (_dv < _dvl && _av < _avl) { \
720 _x = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _cc; \
722 _cc = _x >> MPW_BITS; \
724 MPX_UADDN(_dv, _dvl, _cc); \
727 extern void mpx_umlan(mpw */
*dv*/
, mpw */
*dvl*/
,
728 const mpw */
*av*/
, const mpw */
*avl*/
, mpw m
);
730 /* --- @mpx_usqr@ --- *
732 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
733 * @const mpw *av, *av@ = source vector base and limit
737 * Use: Performs unsigned integer squaring. The result vector must
738 * not overlap the source vector in any way.
741 extern void mpx_usqr(mpw */
*dv*/
, mpw */
*dvl*/
,
742 const mpw */
*av*/
, const mpw */
*avl*/
);
744 /* --- @mpx_udiv@ --- *
746 * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit
747 * @mpw *rv, *rvl@ = dividend/remainder vector base and limit
748 * @const mpw *dv, *dvl@ = divisor vector base and limit
749 * @mpw *sv, *svl@ = scratch workspace
753 * Use: Performs unsigned integer division. If the result overflows
754 * the quotient vector, high-order bits are discarded. (Clearly
755 * the remainder vector can't overflow.) The various vectors
756 * may not overlap in any way. Yes, I know it's a bit odd
757 * requiring the dividend to be in the result position but it
758 * does make some sense really. The remainder must have
759 * headroom for at least two extra words. The scratch space
760 * must be at least one word larger than the divisor.
763 extern void mpx_udiv(mpw */
*qv*/
, mpw */
*qvl*/
, mpw */
*rv*/
, mpw */
*rvl*/
,
764 const mpw */
*dv*/
, const mpw */
*dvl*/
,
765 mpw */
*sv*/
, mpw */
*svl*/
);
767 /* --- @mpx_udivn@ --- *
769 * Arguments: @mpw *qv, *qvl@ = storage for the quotient (may overlap
771 * @const mpw *rv, *rvl@ = dividend
772 * @mpw d@ = single-precision divisor
774 * Returns: Remainder after divison.
776 * Use: Performs a single-precision division operation.
779 extern mpw
mpx_udivn(mpw */
*qv*/
, mpw */
*qvl*/
,
780 const mpw */
*rv*/
, const mpw */
*rvl*/
, mpw
/*d*/);
782 /*----- Karatsuba multiplication algorithms -------------------------------*/
784 /* --- @MPK_THRESH@ --- *
786 * This is the limiting length for using Karatsuba algorithms. It's best to
787 * use the simpler classical multiplication method on numbers smaller than
788 * this. It is unsafe to make this constant less than four (i.e., the
789 * algorithms will fail).
792 #define MPK_THRESH 16
794 /* --- @mpx_kmul@ --- *
796 * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer
797 * @const mpw *av, *avl@ = pointer to first argument
798 * @const mpw *bv, *bvl@ = pointer to second argument
799 * @mpw *sv, *svl@ = pointer to scratch workspace
803 * Use: Multiplies two multiprecision integers using Karatsuba's
804 * algorithm. This is rather faster than traditional long
805 * multiplication (e.g., @mpx_umul@) on large numbers, although
806 * more expensive on small ones.
808 * The destination must be three times as large as the larger
809 * argument. The scratch space must be five times as large as
810 * the larger argument.
813 extern void mpx_kmul(mpw */
*dv*/
, mpw */
*dvl*/
,
814 const mpw */
*av*/
, const mpw */
*avl*/
,
815 const mpw */
*bv*/
, const mpw */
*bvl*/
,
816 mpw */
*sv*/
, mpw */
*svl*/
);
818 /* --- @mpx_ksqr@ --- *
820 * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer
821 * @const mpw *av, *avl@ = pointer to first argument
822 * @mpw *sv, *svl@ = pointer to scratch workspace
826 * Use: Squares a multiprecision integers using something similar to
827 * Karatsuba's multiplication algorithm. This is rather faster
828 * than traditional long multiplication (e.g., @mpx_umul@) on
829 * large numbers, although more expensive on small ones, and
830 * rather simpler than full-blown Karatsuba multiplication.
832 * The destination must be three times as large as the larger
833 * argument. The scratch space must be five times as large as
834 * the larger argument.
837 extern void mpx_ksqr(mpw */
*dv*/
, mpw */
*dvl*/
,
838 const mpw */
*av*/
, const mpw */
*avl*/
,
839 mpw */
*sv*/
, mpw */
*svl*/
);
841 /*----- That's all, folks -------------------------------------------------*/