X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/3c9ede175eabf2276f9ff99e8b886562f8b661a3..6ec3a4cf4aaa7cd375e1aa18f85861986259b8e5:/mpx.h diff --git a/mpx.h b/mpx.h index 9429e91..19f5cc7 100644 --- a/mpx.h +++ b/mpx.h @@ -1,13 +1,13 @@ /* -*-c-*- * - * $Id: mpx.h,v 1.3 1999/11/13 01:51:29 mdw Exp $ + * $Id: mpx.h,v 1.18 2004/04/08 01:36:15 mdw Exp $ * * Low level multiprecision arithmetic * * (c) 1999 Straylight/Edgeware */ -/*----- Licensing notice --------------------------------------------------* +/*----- Licensing notice --------------------------------------------------* * * This file is part of Catacomb. * @@ -15,34 +15,20 @@ * it under the terms of the GNU Library General Public License as * published by the Free Software Foundation; either version 2 of the * License, or (at your option) any later version. - * + * * Catacomb is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Library General Public License for more details. - * + * * You should have received a copy of the GNU Library General Public * License along with Catacomb; if not, write to the Free * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, * MA 02111-1307, USA. */ -/*----- Revision history --------------------------------------------------* - * - * $Log: mpx.h,v $ - * Revision 1.3 1999/11/13 01:51:29 mdw - * Minor interface changes. Should be stable now. - * - * Revision 1.2 1999/11/11 17:47:55 mdw - * Minor changes for different `mptypes.h' format. - * - * Revision 1.1 1999/09/03 08:41:12 mdw - * Initial import. - * - */ - -#ifndef MPX_H -#define MPX_H +#ifndef CATACOMB_MPX_H +#define CATACOMB_MPX_H #ifdef __cplusplus extern "C" { @@ -65,7 +51,7 @@ #include -#ifndef MPW_H +#ifndef CATACOMB_MPW_H # include "mpw.h" #endif @@ -101,9 +87,9 @@ if (_v == _vl) \ (b) = 0; \ else { \ - unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1; \ + unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1; \ mpw _w = _vl[-1]; \ - unsigned _k = MPW_BITS / 2; \ + unsigned _k = MPW_P2; \ while (_k) { \ if (_w >> _k) { \ _w >>= _k; \ @@ -124,25 +110,24 @@ */ #define MPX_OCTETS(o, v, vl) do { \ - const mpw *_v = (v), *_vl = (vl); \ - MPX_SHRINK(_v, _vl); \ - if (_v == _vl) \ - (o) = 0; \ - else { \ - size_t _o = (MPW_BITS / 8) * (_vl - _v - 1); \ - mpw _w = _vl[-1]; \ - unsigned _k = MPW_BITS / 2; \ - while (_k >= 8) { \ - if (_w >> _k) { \ - _w >>= _k; \ - _o += _k >> 3; \ - } \ - _k >>= 1; \ - } \ - if (_w) \ - _o++; \ - (o) = _o; \ - } \ + unsigned long _bb; \ + MPX_BITS(_bb, (v), (vl)); \ + (o) = (_bb + 7) >> 3; \ +} while (0) + +/* --- @MPX_OCTETS2C@ --- * + * + * Arguments: @size_t o@ = result variable + * @const mpw *v, *vl@ = pointer to array of words + * + * Use: Calculates the number of octets in a multiprecision value, if + * you represent it as two's complement. + */ + +#define MPX_OCTETS2C(o, v, vl) do { \ + unsigned long _bb; \ + MPX_BITS(_bb, (v), (vl)); \ + (o) = (_bb >> 3) + 1; \ } while (0) /* --- @MPX_COPY@ --- * @@ -182,6 +167,20 @@ memset(_v, 0, MPWS(_vl - _v)); \ } while (0) +/* --- @MPX_ONE@ --- * + * + * Arguments: @v, vl@ = base and limit of vector to clear + * + * Use: Fills the area between the two vector pointers with ones. + */ + +#define MPX_ONE(v, vl) do { \ + mpw * _v = (v); \ + const mpw *_vl = (vl); \ + while (_v < _vl) \ + *_v++ = MPW_MAX; \ +} while (0) + /*----- Loading and storing -----------------------------------------------*/ /* --- @mpx_storel@ --- * @@ -248,6 +247,75 @@ extern void mpx_storeb(const mpw */*v*/, const mpw */*vl*/, extern void mpx_loadb(mpw */*v*/, mpw */*vl*/, const void */*p*/, size_t /*sz*/); +/* --- @mpx_storel2cn@ --- * + * + * Arguments: @const mpw *v, *vl@ = base and limit of source vector + * @void *pp@ = pointer to octet array + * @size_t sz@ = size of octet array + * + * Returns: --- + * + * Use: Stores a negative MP in an octet array, least significant + * octet first, as two's complement. High-end octets are + * silently discarded if there isn't enough space for them. + * This obviously makes the output bad. + */ + +extern void mpx_storel2cn(const mpw */*v*/, const mpw */*vl*/, + void */*p*/, size_t /*sz*/); + +/* --- @mpx_loadl2cn@ --- * + * + * Arguments: @mpw *v, *vl@ = base and limit of destination vector + * @const void *pp@ = pointer to octet array + * @size_t sz@ = size of octet array + * + * Returns: --- + * + * Use: Loads a negative MP in an octet array, least significant + * octet first, as two's complement. High-end octets are + * ignored if there isn't enough space for them. This probably + * means you made the wrong choice coming here. + */ + +extern void mpx_loadl2cn(mpw */*v*/, mpw */*vl*/, + const void */*p*/, size_t /*sz*/); + +/* --- @mpx_storeb2cn@ --- * + * + * Arguments: @const mpw *v, *vl@ = base and limit of source vector + * @void *pp@ = pointer to octet array + * @size_t sz@ = size of octet array + * + * Returns: --- + * + * Use: Stores a negative MP in an octet array, most significant + * octet first, as two's complement. High-end octets are + * silently discarded if there isn't enough space for them, + * which probably isn't what you meant. + */ + +extern void mpx_storeb2cn(const mpw */*v*/, const mpw */*vl*/, + void */*p*/, size_t /*sz*/); + +/* --- @mpx_loadb2cn@ --- * + * + * Arguments: @mpw *v, *vl@ = base and limit of destination vector + * @const void *pp@ = pointer to octet array + * @size_t sz@ = size of octet array + * + * Returns: --- + * + * Use: Loads a negative MP in an octet array, most significant octet + * first as two's complement. High-end octets are ignored if + * there isn't enough space for them. This probably means you + * chose this function wrongly. + */ + +extern void mpx_loadb2cn(mpw */*v*/, mpw */*vl*/, + const void */*p*/, size_t /*sz*/); + + /*----- Logical shifting --------------------------------------------------*/ /* --- @mpx_lsl@ --- * @@ -265,6 +333,22 @@ extern void mpx_lsl(mpw */*dv*/, mpw */*dvl*/, const mpw */*av*/, const mpw */*avl*/, size_t /*n*/); +/* --- @mpx_lslc@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector base and limit + * @const mpw *av, *avl@ = source vector base and limit + * @size_t n@ = number of bit positions to shift by + * + * Returns: --- + * + * Use: Performs a logical shift left operation on an integer, only + * it fills in the bits with ones instead of zeroes. + */ + +extern void mpx_lslc(mpw */*dv*/, mpw */*dvl*/, + const mpw */*av*/, const mpw */*avl*/, + size_t /*n*/); + /* --- @mpx_lsr@ --- * * * Arguments: @mpw *dv, *dvl@ = destination vector base and limit @@ -280,8 +364,88 @@ extern void mpx_lsr(mpw */*dv*/, mpw */*dvl*/, const mpw */*av*/, const mpw */*avl*/, size_t /*n*/); +/*----- Bitwise operations ------------------------------------------------*/ + +/* --- @mpx_bitop@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector + * @const mpw *av, *avl@ = first source vector + * @const mpw *bv, *bvl@ = second source vector + * + * Returns: --- + * + * Use: Provide the dyadic boolean functions. The functions are + * named after the truth table they generate: + * + * a: 0011 + * b: 0101 + * @mpx_bitXXXX@ + */ + +#define MPX_DOBIN(what) \ + what(0000) what(0001) what(0010) what(0011) \ + what(0100) what(0101) what(0110) what(0111) \ + what(1000) what(1001) what(1010) what(1011) \ + what(1100) what(1101) what(1110) what(1111) + +#define MPX_BITDECL(string) \ + extern void mpx_bit##string(mpw */*dv*/, mpw */*dvl*/, \ + const mpw */*av*/, const mpw */*avl*/, \ + const mpw */*bv*/, const mpw */*bvl*/); +MPX_DOBIN(MPX_BITDECL) + +/* --- @mpx_[n]and@, @mpx_[n]or@, @mpx_xor@ --- * + * + * Synonyms for the commonly-used functions above. + */ + +#define mpx_and mpx_bit0001 +#define mpx_or mpx_bit0111 +#define mpx_nand mpx_bit1110 +#define mpx_nor mpx_bit1000 +#define mpx_xor mpx_bit0110 + +/* --- @mpx_not@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector + * @const mpw *av, *avl@ = first source vector + * + * Returns: --- + * + * Use; Bitwise NOT. + */ + +extern void mpx_not(mpw */*dv*/, mpw */*dvl*/, + const mpw */*av*/, const mpw */*avl*/); + /*----- Unsigned arithmetic -----------------------------------------------*/ +/* --- @mpx_2c@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector + * @const mpw *v, *vl@ = source vector + * + * Returns: --- + * + * Use: Calculates the two's complement of @v@. + */ + +extern void mpx_2c(mpw */*dv*/, mpw */*dvl*/, + const mpw */*v*/, const mpw */*vl*/); + +/* --- @mpx_ueq@ --- * + * + * Arguments: @const mpw *av, *avl@ = first argument vector base and limit + * @const mpw *bv, *bvl@ = second argument vector base and limit + * + * Returns: Nonzero if the two vectors are equal. + * + * Use: Performs an unsigned integer test for equality. + */ + +extern int mpx_ueq(const mpw */*av*/, const mpw */*avl*/, + const mpw */*bv*/, const mpw */*bvl*/); + /* --- @mpx_ucmp@ --- * * * Arguments: @const mpw *av, *avl@ = first argument vector base and limit @@ -320,10 +484,12 @@ extern void mpx_uadd(mpw */*dv*/, mpw */*dvl*/, const mpw */*av*/, const mpw */*avl*/, const mpw */*bv*/, const mpw */*bvl*/); -/* --- @MPX_UADDN@ --- * +/* --- @mpx_uaddn@ --- * * - * Arguments: @dv, dvl@ = source and destination vector base and limit - * @n@ = other addend + * Arguments: @mpw *dv, *dvl@ = source and destination base and limit + * @mpw n@ = other addend + * + * Returns: --- * * Use: Adds a small integer to a multiprecision number. */ @@ -339,6 +505,24 @@ extern void mpx_uadd(mpw */*dv*/, mpw */*dvl*/, } \ } while (0) +extern void mpx_uaddn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/); + +/* --- @mpx_uaddnlsl@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination and first argument vector + * @mpw a@ = second argument + * @unsigned o@ = offset in bits + * + * Returns: --- + * + * Use: Computes %$d + 2^o a$%. If the result overflows then + * high-order bits are discarded, as usual. We must have + * @0 < o < MPW_BITS@. + */ + +extern void mpx_uaddnlsl(mpw */*dv*/, mpw */*dvl*/, + mpw /*a*/, unsigned /*o*/); + /* --- @mpx_usub@ --- * * * Arguments: @mpw *dv, *dvl@ = destination vector base and limit @@ -360,10 +544,12 @@ extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/, const mpw */*av*/, const mpw */*avl*/, const mpw */*bv*/, const mpw */*bvl*/); -/* --- @MPX_USUBN@ --- * +/* --- @mpx_usubn@ --- * + * + * Arguments: @mpw *dv, *dvl@ = source and destination base and limit + * @n@ = subtrahend * - * Arguments: @@dv, dvl@ = destination vector base and limit - * @n@ = other addend + * Returns: --- * * Use: Subtracts a small integer from a multiprecision number. */ @@ -382,6 +568,25 @@ extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/, } \ } while (0) +extern void mpx_usubn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/); + +/* --- @mpx_usubnlsl@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination and first argument vector + * @mpw a@ = second argument + * @unsigned o@ = offset in bits + * + * Returns: --- + * + * Use: Computes %$d - 2^o a$%. If the result overflows then + * high-order bits are discarded, as usual, so you get two's + * complement. Which might be what you wanted... We must have + * @0 < o < MPW_BITS@. + */ + +extern void mpx_usubnlsl(mpw */*dv*/, mpw */*dvl*/, + mpw /*a*/, unsigned /*o*/); + /* --- @mpx_umul@ --- * * * Arguments: @mpw *dv, *dvl@ = destination vector base and limit @@ -400,11 +605,13 @@ extern void mpx_umul(mpw */*dv*/, mpw */*dvl*/, const mpw */*av*/, const mpw */*avl*/, const mpw */*bv*/, const mpw */*bvl*/); -/* --- @MPX_UMULN@ --- * +/* --- @mpx_umuln@ --- * * - * Arguments: @dv, dvl@ = destination vector base and limit - * @av, avl@ = multiplicand vector base and limit - * @m@ = multiplier + * Arguments: @mpw *dv, *dvl@ = destination vector base and limit + * @const mpw *av, *avl@ = multiplicand vector base and limit + * @mpw m@ = multiplier + * + * Returns: --- * * Use: Multiplies a multiprecision integer by a single-word value. * The destination and source may be equal. The destination @@ -431,11 +638,16 @@ extern void mpx_umul(mpw */*dv*/, mpw */*dvl*/, } \ } while (0) -/* --- @MPX_UMLAN@ --- * +extern void mpx_umuln(mpw */*dv*/, mpw */*dvl*/, + const mpw */*av*/, const mpw */*avl*/, mpw m); + +/* --- @mpx_umlan@ --- * * - * Arguments: @dv, dvl@ = destination/accumulator vector base and limit - * @av, avl@ = multiplicand vector base and limit - * @m@ = multiplier + * Arguments: @mpw *dv, *dvl@ = destination/accumulator base and limit + * @const mpw *av, *avl@ = multiplicand vector base and limit + * @mpw m@ = multiplier + * + * Returns: --- * * Use: Multiplies a multiprecision integer by a single-word value * and adds the result to an accumulator. @@ -444,20 +656,21 @@ extern void mpx_umul(mpw */*dv*/, mpw */*dvl*/, #define MPX_UMLAN(dv, dvl, av, avl, m) do { \ mpw *_dv = (dv), *_dvl = (dvl); \ const mpw *_av = (av), *_avl = (avl); \ - mpw _c = 0; \ + mpw _cc = 0; \ mpd _m = (m); \ \ - while (_av < _avl) { \ + while (_dv < _dvl && _av < _avl) { \ mpd _x; \ - if (_dv >= _dvl) \ - break; \ - _x = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _c; \ + _x = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _cc; \ *_dv++ = MPW(_x); \ - _c = _x >> MPW_BITS; \ + _cc = _x >> MPW_BITS; \ } \ - MPX_UADDN(_dv, _dvl, _c); \ + MPX_UADDN(_dv, _dvl, _cc); \ } while (0) +extern void mpx_umlan(mpw */*dv*/, mpw */*dvl*/, + const mpw */*av*/, const mpw */*avl*/, mpw m); + /* --- @mpx_usqr@ --- * * * Arguments: @mpw *dv, *dvl@ = destination vector base and limit @@ -477,7 +690,7 @@ extern void mpx_usqr(mpw */*dv*/, mpw */*dvl*/, * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit * @mpw *rv, *rvl@ = dividend/remainder vector base and limit * @const mpw *dv, *dvl@ = divisor vector base and limit - * @mpw *sv, *svl@ = scratch workspace + * @mpw *sv, *svl@ = scratch workspace * * Returns: --- * @@ -488,14 +701,87 @@ extern void mpx_usqr(mpw */*dv*/, mpw */*dvl*/, * requiring the dividend to be in the result position but it * does make some sense really. The remainder must have * headroom for at least two extra words. The scratch space - * must be at least two words larger than twice the size of the - * divisor. + * must be at least one word larger than the divisor. */ extern void mpx_udiv(mpw */*qv*/, mpw */*qvl*/, mpw */*rv*/, mpw */*rvl*/, const mpw */*dv*/, const mpw */*dvl*/, mpw */*sv*/, mpw */*svl*/); +/* --- @mpx_udivn@ --- * + * + * Arguments: @mpw *qv, *qvl@ = storage for the quotient (may overlap + * dividend) + * @const mpw *rv, *rvl@ = dividend + * @mpw d@ = single-precision divisor + * + * Returns: Remainder after divison. + * + * Use: Performs a single-precision division operation. + */ + +extern mpw mpx_udivn(mpw */*qv*/, mpw */*qvl*/, + const mpw */*rv*/, const mpw */*rvl*/, mpw /*d*/); + +/*----- Karatsuba multiplication algorithms -------------------------------*/ + +/* --- @MPK_THRESH@ --- * + * + * This is the limiting length for using Karatsuba algorithms. It's best to + * use the simpler classical multiplication method on numbers smaller than + * this. It is unsafe to make this constant less than four (i.e., the + * algorithms will fail). + */ + +#define MPK_THRESH 16 + +/* --- @mpx_kmul@ --- * + * + * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer + * @const mpw *av, *avl@ = pointer to first argument + * @const mpw *bv, *bvl@ = pointer to second argument + * @mpw *sv, *svl@ = pointer to scratch workspace + * + * Returns: --- + * + * Use: Multiplies two multiprecision integers using Karatsuba's + * algorithm. This is rather faster than traditional long + * multiplication (e.g., @mpx_umul@) on large numbers, although + * more expensive on small ones. + * + * The destination must be three times as large as the larger + * argument. The scratch space must be five times as large as + * the larger argument. + */ + +extern void mpx_kmul(mpw */*dv*/, mpw */*dvl*/, + const mpw */*av*/, const mpw */*avl*/, + const mpw */*bv*/, const mpw */*bvl*/, + mpw */*sv*/, mpw */*svl*/); + +/* --- @mpx_ksqr@ --- * + * + * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer + * @const mpw *av, *avl@ = pointer to first argument + * @mpw *sv, *svl@ = pointer to scratch workspace + * + * Returns: --- + * + * Use: Squares a multiprecision integers using something similar to + * Karatsuba's multiplication algorithm. This is rather faster + * than traditional long multiplication (e.g., @mpx_umul@) on + * large numbers, although more expensive on small ones, and + * rather simpler than full-blown Karatsuba multiplication. + * + * The destination must be three times as large as the larger + * argument. The scratch space must be five times as large as + * the larger argument. + */ + +extern void mpx_ksqr(mpw */*dv*/, mpw */*dvl*/, + const mpw */*av*/, const mpw */*avl*/, + mpw */*sv*/, mpw */*svl*/); + /*----- That's all, folks -------------------------------------------------*/ #ifdef __cplusplus