[u/mdw/catacomb] / mpx.h

/* -*-c-*-
 *
 * $Id: mpx.h,v 1.9 1999/12/22 15:49:07 mdw Exp $
 *
 * Low level multiprecision arithmetic
 *
 * (c) 1999 Straylight/Edgeware
 */

/*----- Licensing notice --------------------------------------------------* 
 *
 * This file is part of Catacomb.
 *
 * Catacomb is free software; you can redistribute it and/or modify
 * 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.9  1999/12/22 15:49:07  mdw
 * New function for division by a small integer.
 *
 * Revision 1.8  1999/12/11 10:57:43  mdw
 * Karatsuba squaring algorithm.
 *
 * Revision 1.7  1999/12/11 01:51:28  mdw
 * Change Karatsuba parameters slightly.
 *
 * Revision 1.6  1999/12/10 23:23:51  mdw
 * Karatsuba-Ofman multiplication algorithm.
 *
 * Revision 1.5  1999/11/20 22:23:27  mdw
 * Add function versions of some low-level macros with wider use.
 *
 * Revision 1.4  1999/11/17 18:04:43  mdw
 * Add two's complement support.  Fix a bug in MPX_UMLAN.
 *
 * 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 CATACOMB_MPX_H
#define CATACOMB_MPX_H

#ifdef __cplusplus
  extern "C" {
#endif

/*----- The idea ----------------------------------------------------------*
 *
 * This file provides functions and macros which work on vectors of words as
 * unsigned multiprecision integers.  The interface works in terms of base
 * and limit pointers (i.e., a pointer to the start of a vector, and a
 * pointer just past its end) rather than base pointer and length, because
 * that requires more arithmetic and state to work on.
 *
 * The interfaces are slightly bizarre in other ways.  Try to use the
 * higher-level functions where you can: they're rather better designed to
 * actually be friendly and useful.
 */

/*----- Header files ------------------------------------------------------*/

#include <string.h>

#ifndef CATACOMB_MPW_H
#  include "mpw.h"
#endif

/*----- General manipulation ----------------------------------------------*/

/* --- @MPX_SHRINK@ --- *
 *
 * Arguments:	@const mpw *v@ = pointer to vector of words
 *		@const mpw *vl@ = (updated) current limit of vector
 *
 * Use:		Shrinks down the limit of a multiprecision integer vector.
 */

#define MPX_SHRINK(v, vl) do {						\
  const mpw *_vv = (v), *_vvl = (vl);					\
  while (_vvl > _vv && !_vvl[-1])					\
    _vvl--;								\
  (vl) = (mpw *)_vvl;							\
} while (0)

/* --- @MPX_BITS@ --- *
 *
 * Arguments:	@unsigned long b@ = result variable
 *		@const mpw *v@ = pointer to array of words
 *		@const mpw *vl@ = limit of vector (from @MPX_SHRINK@)
 *
 * Use:		Calculates the number of bits in a multiprecision value.
 */

#define MPX_BITS(b, v, vl) do {						\
  const mpw *_v = (v), *_vl = (vl);					\
  MPX_SHRINK(_v, _vl);							\
  if (_v == _vl)							\
    (b) = 0;								\
  else {								\
    unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1;	       		\
    mpw _w = _vl[-1];							\
    unsigned _k = MPW_BITS / 2;						\
    while (_k) {							\
      if (_w >> _k) {							\
	_w >>= _k;							\
	_b += _k;							\
      }									\
      _k >>= 1;								\
    }									\
    (b) = _b;								\
  }									\
} while (0)

/* --- @MPX_OCTETS@ --- *
 *
 * 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.
 */

#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;								\
  }									\
} while (0)

/* --- @MPX_COPY@ --- *
 *
 * Arguments:	@dv, dvl@ = destination vector base and limit
 *		@av, avl@ = source vector base and limit
 *
 * Use:		Copies a multiprecision integer.
 */

#define MPX_COPY(dv, dvl, av, avl) do {					\
  mpw *_dv = (dv), *_dvl = (dvl);					\
  size_t _dn = _dvl - _dv;						\
  const mpw *_av = (av), *_avl = (avl);					\
  size_t _an = _avl - _av;						\
  if (_av == _dv) {							\
    if (_dvl > _avl)							\
      memset(_dv, 0, MPWS(_dn - _an));					\
  } else if (_an >= _dn)						\
    memmove(_dv, _av, MPWS(_dn));					\
  else {								\
    memmove(_dv, _av, MPWS(_an));					\
    memset(_dv + _an, 0, MPWS(_dn - _an));				\
  }									\
} while (0)

/* --- @MPX_ZERO@ --- *
 *
 * Arguments:	@v, vl@ = base and limit of vector to clear
 *
 * Use:		Zeroes the area between the two vector pointers.
 */

#define MPX_ZERO(v, vl) do {						\
  mpw *_v = (v), *_vl = (vl);						\
  if (_v < _vl)								\
    memset(_v, 0, MPWS(_vl - _v));					\
} while (0)

/*----- Loading and storing -----------------------------------------------*/

/* --- @mpx_storel@ --- *
 *
 * Arguments:	@const mpw *v, *vl@ = base and limit of source vector
 *		@void *p@ = pointer to octet array
 *		@size_t sz@ = size of octet array
 *
 * Returns:	---
 *
 * Use:		Stores an MP in an octet array, least significant octet
 *		first.  High-end octets are silently discarded if there
 *		isn't enough space for them.
 */

extern void mpx_storel(const mpw */*v*/, const mpw */*vl*/,
		       void */*p*/, size_t /*sz*/);

/* --- @mpx_loadl@ --- *
 *
 * Arguments:	@mpw *v, *vl@ = base and limit of destination vector
 *		@const void *p@ = pointer to octet array
 *		@size_t sz@ = size of octet array
 *
 * Returns:	---
 *
 * Use:		Loads an MP in an octet array, least significant octet
 *		first.  High-end octets are ignored if there isn't enough
 *		space for them.
 */

extern void mpx_loadl(mpw */*v*/, mpw */*vl*/,
		      const void */*p*/, size_t /*sz*/);

/* --- @mpx_storeb@ --- *
 *
 * Arguments:	@const mpw *v, *vl@ = base and limit of source vector
 *		@void *p@ = pointer to octet array
 *		@size_t sz@ = size of octet array
 *
 * Returns:	---
 *
 * Use:		Stores an MP in an octet array, most significant octet
 *		first.  High-end octets are silently discarded if there
 *		isn't enough space for them.
 */

extern void mpx_storeb(const mpw */*v*/, const mpw */*vl*/,
		       void */*p*/, size_t /*sz*/);

/* --- @mpx_loadb@ --- *
 *
 * Arguments:	@mpw *v, *vl@ = base and limit of destination vector
 *		@const void *p@ = pointer to octet array
 *		@size_t sz@ = size of octet array
 *
 * Returns:	---
 *
 * Use:		Loads an MP in an octet array, most significant octet
 *		first.  High-end octets are ignored if there isn't enough
 *		space for them.
 */

extern void mpx_loadb(mpw */*v*/, mpw */*vl*/,
		      const void */*p*/, size_t /*sz*/);

/*----- Logical shifting --------------------------------------------------*/

/* --- @mpx_lsl@ --- *
 *
 * 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.
 */

extern void mpx_lsl(mpw */*dv*/, mpw */*dvl*/,
		    const mpw */*av*/, const mpw */*avl*/,
		    size_t /*n*/);

/* --- @mpx_lsr@ --- *
 *
 * 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 right operation on an integer.
 */

extern void mpx_lsr(mpw */*dv*/, mpw */*dvl*/,
		    const mpw */*av*/, const mpw */*avl*/,
		    size_t /*n*/);

/*----- 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_ucmp@ --- *
 *
 * Arguments:	@const mpw *av, *avl@ = first argument vector base and limit
 *		@const mpw *bv, *bvl@ = second argument vector base and limit
 *
 * Returns:	Less than, equal to, or greater than zero depending on
 *		whether @a@ is less than, equal to or greater than @b@,
 *		respectively.
 *
 * Use:		Performs an unsigned integer comparison.
 */

#define MPX_UCMP(av, avl, op, dv, dvl)					\
  (mpx_ucmp((av), (avl), (dv), (dvl)) op 0)

extern int mpx_ucmp(const mpw */*av*/, const mpw */*avl*/,
		    const mpw */*bv*/, const mpw */*bvl*/);

/* --- @mpx_uadd@ --- *
 *
 * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
 *		@const mpw *av, *avl@ = first addend vector base and limit
 *		@const mpw *bv, *bvl@ = second addend vector base and limit
 *
 * Returns:	---
 *
 * Use:		Performs unsigned integer addition.  If the result overflows
 *		the destination vector, high-order bits are discarded.  This
 *		means that two's complement addition happens more or less for
 *		free, although that's more a side-effect than anything else.
 *		The result vector may be equal to either or both source
 *		vectors, but may not otherwise overlap them.
 */

extern void mpx_uadd(mpw */*dv*/, mpw */*dvl*/,
		     const mpw */*av*/, const mpw */*avl*/,
		     const mpw */*bv*/, const mpw */*bvl*/);

/* --- @mpx_uaddn@ --- *
 *
 * Arguments:	@mpw *dv, *dvl@ = source and destination base and limit
 *		@mpw n@ = other addend
 *
 * Returns:	---
 *
 * Use:		Adds a small integer to a multiprecision number.
 */

#define MPX_UADDN(dv, dvl, n) do {					\
  mpw *_ddv = (dv), *_ddvl = (dvl);					\
  mpw _c = (n);								\
									\
  while (_c && _ddv < _ddvl) {						\
    mpd _x = (mpd)*_ddv + (mpd)_c;					\
    *_ddv++ = MPW(_x);							\
    _c = _x >> MPW_BITS;						\
  }									\
} while (0)

extern void mpx_uaddn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);

/* --- @mpx_usub@ --- *
 *
 * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
 *		@const mpw *av, *avl@ = first argument vector base and limit
 *		@const mpw *bv, *bvl@ = second argument vector base and limit
 *
 * Returns:	---
 *
 * Use:		Performs unsigned integer subtraction.  If the result
 *		overflows the destination vector, high-order bits are
 *		discarded.  This means that two's complement subtraction
 *		happens more or less for free, although that's more a side-
 *		effect than anything else.  The result vector may be equal to
 *		either or both source vectors, but may not otherwise overlap
 *		them.
 */

extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/,
		     const mpw */*av*/, const mpw */*avl*/,
		     const mpw */*bv*/, const mpw */*bvl*/);

/* --- @mpx_usubn@ --- *
 *
 * Arguments:	@mpw *dv, *dvl@ = source and destination base and limit
 *		@n@ = subtrahend
 *
 * Returns:	---
 *
 * Use:		Subtracts a small integer from a multiprecision number.
 */

#define MPX_USUBN(dv, dvl, n) do {					\
  mpw *_ddv = (dv), *_ddvl = (dvl);					\
  mpw _c = (n);								\
									\
  while (_ddv < _ddvl) {						\
    mpd _x = (mpd)*_ddv - (mpd)_c;					\
    *_ddv++ = MPW(_x);							\
    if (_x >> MPW_BITS)							\
      _c = 1;								\
    else								\
      break;								\
  }									\
} while (0)

extern void mpx_usubn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);

/* --- @mpx_umul@ --- *
 *
 * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
 *		@const mpw *av, *avl@ = multiplicand vector base and limit
 *		@const mpw *bv, *bvl@ = multiplier vector base and limit
 *
 * Returns:	---
 *
 * Use:		Performs unsigned integer multiplication.  If the result
 *		overflows the desination vector, high-order bits are
 *		discarded.  The result vector may not overlap the argument
 *		vectors in any way.
 */

extern void mpx_umul(mpw */*dv*/, mpw */*dvl*/,
		     const mpw */*av*/, const mpw */*avl*/,
		     const mpw */*bv*/, const mpw */*bvl*/);

/* --- @mpx_umuln@ --- *
 *
 * 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
 *		is completely cleared after use.
 */

#define MPX_UMULN(dv, dvl, av, avl, m) do {				\
  mpw *_dv = (dv), *_dvl = (dvl);					\
  const mpw *_av = (av), *_avl = (avl);					\
  mpw _c = 0;								\
  mpd _m = (m);								\
									\
  while (_av < _avl) {							\
    mpd _x;								\
    if (_dv >= _dvl)							\
      break;								\
    _x = (mpd)_m * (mpd)*_av++ + _c;					\
    *_dv++ = MPW(_x);							\
    _c = _x >> MPW_BITS;						\
  }									\
  if (_dv < _dvl) {							\
    *_dv++ = MPW(_c);							\
    MPX_ZERO(_dv, _dvl);						\
  }									\
} while (0)

extern void mpx_umuln(mpw */*dv*/, mpw */*dvl*/,
		      const mpw */*av*/, const mpw */*avl*/, mpw m);

/* --- @mpx_umlan@ --- *
 *
 * 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.
 */

#define MPX_UMLAN(dv, dvl, av, avl, m) do {				\
  mpw *_dv = (dv), *_dvl = (dvl);					\
  const mpw *_av = (av), *_avl = (avl);					\
  mpw _cc = 0;								\
  mpd _m = (m);								\
									\
  while (_dv < _dvl && _av < _avl) {					\
    mpd _x;								\
    _x = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _cc;			\
    *_dv++ = MPW(_x);							\
    _cc = _x >> MPW_BITS;						\
  }									\
  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
 *		@const mpw *av, *av@ = source vector base and limit
 *
 * Returns:	---
 *
 * Use:		Performs unsigned integer squaring.  The result vector must
 *		not overlap the source vector in any way.
 */

extern void mpx_usqr(mpw */*dv*/, mpw */*dvl*/,
		     const mpw */*av*/, const mpw */*avl*/);

/* --- @mpx_udiv@ --- *
 *
 * 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 
 *
 * Returns:	---
 *
 * Use:		Performs unsigned integer division.  If the result overflows
 *		the quotient vector, high-order bits are discarded.  (Clearly
 *		the remainder vector can't overflow.)  The various vectors
 *		may not overlap in any way.  Yes, I know it's a bit odd
 *		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 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 -------------------------------*/

/* --- @KARATSUBA_CUTOFF@ --- *
 *
 * This is the limiting length for using Karatsuba algorithms.  It's best to
 * use the simpler classical multiplication method on numbers smaller than
 * this.
 */

#define KARATSUBA_CUTOFF 16

/* --- @KARATSUBA_SLOP@ --- *
 *
 * The extra number of words required as scratch space by the Karatsuba
 * routines.  This is a (generous) guess, since the actual amount of space
 * required is proportional to the recursion depth.
 */

#define KARATSUBA_SLOP 64

/* --- @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 and scratch buffers must be twice as large as
 *		the larger argument.  The scratch space must be twice as
 *		large as the larger argument, plus the magic number
 *		@KARATSUBA_SLOP@.
 */

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 twice as large as the argument.  The
 *		scratch space must be twice as large as the argument, plus
 *		the magic number @KARATSUBA_SLOP@.
 */

extern void mpx_ksqr(mpw */*dv*/, mpw */*dvl*/,
		     const mpw */*av*/, const mpw */*avl*/,
		     mpw */*sv*/, mpw */*svl*/);

/*----- That's all, folks -------------------------------------------------*/

#ifdef __cplusplus
  }
#endif

#endif
Commit	Line	Data
d03ab969	1	/* --c--
d03ab969	2	*
698bd937	3	* $Id: mpx.h,v 1.9 1999/12/22 15:49:07 mdw Exp $
d03ab969	4	*
	5	* Low level 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: mpx.h,v $
698bd937	33	* Revision 1.9 1999/12/22 15:49:07 mdw
	34	* New function for division by a small integer.
	35	*
5bf74dea	36	* Revision 1.8 1999/12/11 10:57:43 mdw
	37	* Karatsuba squaring algorithm.
	38	*
652a6acf	39	* Revision 1.7 1999/12/11 01:51:28 mdw
	40	* Change Karatsuba parameters slightly.
	41	*
a86e33af	42	* Revision 1.6 1999/12/10 23:23:51 mdw
	43	* Karatsuba-Ofman multiplication algorithm.
	44	*
dd517851	45	* Revision 1.5 1999/11/20 22:23:27 mdw
	46	* Add function versions of some low-level macros with wider use.
	47	*
7c13f461	48	* Revision 1.4 1999/11/17 18:04:43 mdw
	49	* Add two's complement support. Fix a bug in MPX_UMLAN.
	50	*
3c9ede17	51	* Revision 1.3 1999/11/13 01:51:29 mdw
	52	* Minor interface changes. Should be stable now.
	53	*
b9b1c853	54	* Revision 1.2 1999/11/11 17:47:55 mdw
	55	* Minor changes for different `mptypes.h' format.
	56	*
d03ab969	57	* Revision 1.1 1999/09/03 08:41:12 mdw
	58	* Initial import.
	59	*
	60	*/
	61
a86e33af	62	#ifndef CATACOMB_MPX_H
a86e33af	63	#define CATACOMB_MPX_H
d03ab969	64
	65	#ifdef __cplusplus
	66	extern "C" {
	67	#endif
	68
	69	/----- The idea ----------------------------------------------------------
	70	*
	71	* This file provides functions and macros which work on vectors of words as
	72	* unsigned multiprecision integers. The interface works in terms of base
	73	* and limit pointers (i.e., a pointer to the start of a vector, and a
	74	* pointer just past its end) rather than base pointer and length, because
	75	* that requires more arithmetic and state to work on.
	76	*
	77	* The interfaces are slightly bizarre in other ways. Try to use the
	78	* higher-level functions where you can: they're rather better designed to
	79	* actually be friendly and useful.
	80	*/
	81
	82	/----- Header files ------------------------------------------------------/
	83
	84	#include <string.h>
	85
a86e33af	86	#ifndef CATACOMB_MPW_H
3c9ede17	87	# include "mpw.h"
d03ab969	88	#endif
	89
	90	/----- General manipulation ----------------------------------------------/
	91
	92	/* --- @MPX_SHRINK@ --- *
	93	*
	94	* Arguments: @const mpw *v@ = pointer to vector of words
	95	* @const mpw *vl@ = (updated) current limit of vector
	96	*
	97	* Use: Shrinks down the limit of a multiprecision integer vector.
	98	*/
	99
	100	#define MPX_SHRINK(v, vl) do { \
3c9ede17	101	const mpw _vv = (v), _vvl = (vl); \
	102	while (_vvl > _vv && !_vvl[-1]) \
	103	_vvl--; \
	104	(vl) = (mpw *)_vvl; \
d03ab969	105	} while (0)
	106
	107	/* --- @MPX_BITS@ --- *
	108	*
	109	* Arguments: @unsigned long b@ = result variable
	110	* @const mpw *v@ = pointer to array of words
	111	* @const mpw *vl@ = limit of vector (from @MPX_SHRINK@)
	112	*
	113	* Use: Calculates the number of bits in a multiprecision value.
	114	*/
	115
	116	#define MPX_BITS(b, v, vl) do { \
	117	const mpw _v = (v), _vl = (vl); \
3c9ede17	118	MPX_SHRINK(_v, _vl); \
d03ab969	119	if (_v == _vl) \
	120	(b) = 0; \
	121	else { \
	122	unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1; \
	123	mpw _w = _vl[-1]; \
	124	unsigned _k = MPW_BITS / 2; \
	125	while (_k) { \
	126	if (_w >> _k) { \
	127	_w >>= _k; \
	128	_b += _k; \
	129	} \
	130	_k >>= 1; \
	131	} \
	132	(b) = _b; \
	133	} \
	134	} while (0)
	135
	136	/* --- @MPX_OCTETS@ --- *
	137	*
	138	* Arguments: @size_t o@ = result variable
3c9ede17	139	* @const mpw v, vl@ = pointer to array of words
d03ab969	140	*
	141	* Use: Calculates the number of octets in a multiprecision value.
	142	*/
	143
3c9ede17	144	#define MPX_OCTETS(o, v, vl) do { \
d03ab969	145	const mpw _v = (v), _vl = (vl); \
3c9ede17	146	MPX_SHRINK(_v, _vl); \
d03ab969	147	if (_v == _vl) \
	148	(o) = 0; \
	149	else { \
3c9ede17	150	size_t _o = (MPW_BITS / 8) * (_vl - _v - 1); \
d03ab969	151	mpw _w = _vl[-1]; \
d03ab969	152	unsigned _k = MPW_BITS / 2; \
3c9ede17	153	while (_k >= 8) { \
d03ab969	154	if (_w >> _k) { \
d03ab969	155	_w >>= _k; \
3c9ede17	156	_o += _k >> 3; \
d03ab969	157	} \
	158	_k >>= 1; \
	159	} \
	160	if (_w) \
	161	_o++; \
	162	(o) = _o; \
	163	} \
	164	} while (0)
	165
	166	/* --- @MPX_COPY@ --- *
	167	*
	168	* Arguments: @dv, dvl@ = destination vector base and limit
	169	* @av, avl@ = source vector base and limit
	170	*
	171	* Use: Copies a multiprecision integer.
	172	*/
	173
3c9ede17	174	#define MPX_COPY(dv, dvl, av, avl) do { \
	175	mpw _dv = (dv), _dvl = (dvl); \
	176	size_t _dn = _dvl - _dv; \
	177	const mpw _av = (av), _avl = (avl); \
	178	size_t _an = _avl - _av; \
d03ab969	179	if (_av == _dv) { \
d03ab969	180	if (_dvl > _avl) \
3c9ede17	181	memset(_dv, 0, MPWS(_dn - _an)); \
d03ab969	182	} else if (_an >= _dn) \
	183	memmove(_dv, _av, MPWS(_dn)); \
	184	else { \
	185	memmove(_dv, _av, MPWS(_an)); \
	186	memset(_dv + _an, 0, MPWS(_dn - _an)); \
	187	} \
	188	} while (0)
	189
	190	/* --- @MPX_ZERO@ --- *
	191	*
	192	* Arguments: @v, vl@ = base and limit of vector to clear
	193	*
	194	* Use: Zeroes the area between the two vector pointers.
	195	*/
	196
3c9ede17	197	#define MPX_ZERO(v, vl) do { \
d03ab969	198	mpw _v = (v), _vl = (vl); \
3c9ede17	199	if (_v < _vl) \
3c9ede17	200	memset(_v, 0, MPWS(_vl - _v)); \
d03ab969	201	} while (0)
	202
	203	/----- Loading and storing -----------------------------------------------/
	204
	205	/* --- @mpx_storel@ --- *
	206	*
	207	* Arguments: @const mpw v, vl@ = base and limit of source vector
3c9ede17	208	* @void *p@ = pointer to octet array
d03ab969	209	* @size_t sz@ = size of octet array
	210	*
	211	* Returns: ---
	212	*
	213	* Use: Stores an MP in an octet array, least significant octet
	214	* first. High-end octets are silently discarded if there
	215	* isn't enough space for them.
	216	*/
	217
	218	extern void mpx_storel(const mpw /v/, const mpw /vl/,
3c9ede17	219	void /p/, size_t /sz*/);
d03ab969	220
	221	/* --- @mpx_loadl@ --- *
	222	*
	223	* Arguments: @mpw v, vl@ = base and limit of destination vector
3c9ede17	224	* @const void *p@ = pointer to octet array
d03ab969	225	* @size_t sz@ = size of octet array
	226	*
	227	* Returns: ---
	228	*
	229	* Use: Loads an MP in an octet array, least significant octet
	230	* first. High-end octets are ignored if there isn't enough
	231	* space for them.
	232	*/
	233
	234	extern void mpx_loadl(mpw /v/, mpw /vl/,
3c9ede17	235	const void /p/, size_t /sz*/);
d03ab969	236
	237	/* --- @mpx_storeb@ --- *
	238	*
	239	* Arguments: @const mpw v, vl@ = base and limit of source vector
3c9ede17	240	* @void *p@ = pointer to octet array
d03ab969	241	* @size_t sz@ = size of octet array
	242	*
	243	* Returns: ---
	244	*
	245	* Use: Stores an MP in an octet array, most significant octet
	246	* first. High-end octets are silently discarded if there
	247	* isn't enough space for them.
	248	*/
	249
	250	extern void mpx_storeb(const mpw /v/, const mpw /vl/,
3c9ede17	251	void /p/, size_t /sz*/);
d03ab969	252
	253	/* --- @mpx_loadb@ --- *
	254	*
	255	* Arguments: @mpw v, vl@ = base and limit of destination vector
3c9ede17	256	* @const void *p@ = pointer to octet array
d03ab969	257	* @size_t sz@ = size of octet array
	258	*
	259	* Returns: ---
	260	*
	261	* Use: Loads an MP in an octet array, most significant octet
	262	* first. High-end octets are ignored if there isn't enough
	263	* space for them.
	264	*/
	265
	266	extern void mpx_loadb(mpw /v/, mpw /vl/,
3c9ede17	267	const void /p/, size_t /sz*/);
d03ab969	268
	269	/----- Logical shifting --------------------------------------------------/
	270
	271	/* --- @mpx_lsl@ --- *
	272	*
	273	* Arguments: @mpw dv, dvl@ = destination vector base and limit
	274	* @const mpw av, avl@ = source vector base and limit
	275	* @size_t n@ = number of bit positions to shift by
	276	*
	277	* Returns: ---
	278	*
	279	* Use: Performs a logical shift left operation on an integer.
	280	*/
	281
	282	extern void mpx_lsl(mpw /dv/, mpw /dvl/,
	283	const mpw /av/, const mpw /avl/,
	284	size_t /n/);
	285
	286	/* --- @mpx_lsr@ --- *
	287	*
	288	* Arguments: @mpw dv, dvl@ = destination vector base and limit
	289	* @const mpw av, avl@ = source vector base and limit
	290	* @size_t n@ = number of bit positions to shift by
	291	*
	292	* Returns: ---
	293	*
	294	* Use: Performs a logical shift right operation on an integer.
	295	*/
	296
	297	extern void mpx_lsr(mpw /dv/, mpw /dvl/,
	298	const mpw /av/, const mpw /avl/,
	299	size_t /n/);
	300
	301	/----- Unsigned arithmetic -----------------------------------------------/
	302
7c13f461	303	/* --- @mpx_2c@ --- *
	304	*
	305	* Arguments: @mpw dv, dvl@ = destination vector
	306	* @const mpw v, vl@ = source vector
	307	*
	308	* Returns: ---
	309	*
	310	* Use: Calculates the two's complement of @v@.
	311	*/
	312
	313	extern void mpx_2c(mpw /dv/, mpw /dvl/,
	314	const mpw /v/, const mpw /vl/);
	315
d03ab969	316	/* --- @mpx_ucmp@ --- *
	317	*
	318	* Arguments: @const mpw av, avl@ = first argument vector base and limit
	319	* @const mpw bv, bvl@ = second argument vector base and limit
	320	*
	321	* Returns: Less than, equal to, or greater than zero depending on
	322	* whether @a@ is less than, equal to or greater than @b@,
	323	* respectively.
	324	*
	325	* Use: Performs an unsigned integer comparison.
	326	*/
	327
	328	#define MPX_UCMP(av, avl, op, dv, dvl) \
	329	(mpx_ucmp((av), (avl), (dv), (dvl)) op 0)
	330
	331	extern int mpx_ucmp(const mpw /av/, const mpw /avl/,
	332	const mpw /bv/, const mpw /bvl/);
	333
	334	/* --- @mpx_uadd@ --- *
	335	*
	336	* Arguments: @mpw dv, dvl@ = destination vector base and limit
	337	* @const mpw av, avl@ = first addend vector base and limit
	338	* @const mpw bv, bvl@ = second addend vector base and limit
	339	*
	340	* Returns: ---
	341	*
	342	* Use: Performs unsigned integer addition. If the result overflows
	343	* the destination vector, high-order bits are discarded. This
	344	* means that two's complement addition happens more or less for
	345	* free, although that's more a side-effect than anything else.
	346	* The result vector may be equal to either or both source
	347	* vectors, but may not otherwise overlap them.
	348	*/
	349
	350	extern void mpx_uadd(mpw /dv/, mpw /dvl/,
	351	const mpw /av/, const mpw /avl/,
	352	const mpw /bv/, const mpw /bvl/);
	353
dd517851	354	/* --- @mpx_uaddn@ --- *
	355	*
	356	* Arguments: @mpw dv, dvl@ = source and destination base and limit
	357	* @mpw n@ = other addend
3c9ede17	358	*
dd517851	359	* Returns: ---
3c9ede17	360	*
	361	* Use: Adds a small integer to a multiprecision number.
	362	*/
	363
	364	#define MPX_UADDN(dv, dvl, n) do { \
	365	mpw _ddv = (dv), _ddvl = (dvl); \
	366	mpw _c = (n); \
	367	\
	368	while (_c && _ddv < _ddvl) { \
	369	mpd _x = (mpd)*_ddv + (mpd)_c; \
	370	*_ddv++ = MPW(_x); \
	371	_c = _x >> MPW_BITS; \
	372	} \
	373	} while (0)
	374
dd517851	375	extern void mpx_uaddn(mpw /dv/, mpw /dvl/, mpw /n/);
dd517851	376
d03ab969	377	/* --- @mpx_usub@ --- *
	378	*
	379	* Arguments: @mpw dv, dvl@ = destination vector base and limit
	380	* @const mpw av, avl@ = first argument vector base and limit
	381	* @const mpw bv, bvl@ = second argument vector base and limit
	382	*
	383	* Returns: ---
	384	*
	385	* Use: Performs unsigned integer subtraction. If the result
	386	* overflows the destination vector, high-order bits are
	387	* discarded. This means that two's complement subtraction
3c9ede17	388	* happens more or less for free, although that's more a side-
d03ab969	389	* effect than anything else. The result vector may be equal to
	390	* either or both source vectors, but may not otherwise overlap
	391	* them.
	392	*/
	393
	394	extern void mpx_usub(mpw /dv/, mpw /dvl/,
	395	const mpw /av/, const mpw /avl/,
	396	const mpw /bv/, const mpw /bvl/);
	397
dd517851	398	/* --- @mpx_usubn@ --- *
3c9ede17	399	*
dd517851	400	* Arguments: @mpw dv, dvl@ = source and destination base and limit
	401	* @n@ = subtrahend
	402	*
	403	* Returns: ---
3c9ede17	404	*
	405	* Use: Subtracts a small integer from a multiprecision number.
	406	*/
	407
	408	#define MPX_USUBN(dv, dvl, n) do { \
	409	mpw _ddv = (dv), _ddvl = (dvl); \
	410	mpw _c = (n); \
	411	\
	412	while (_ddv < _ddvl) { \
	413	mpd _x = (mpd)*_ddv - (mpd)_c; \
	414	*_ddv++ = MPW(_x); \
	415	if (_x >> MPW_BITS) \
	416	_c = 1; \
	417	else \
	418	break; \
	419	} \
	420	} while (0)
	421
dd517851	422	extern void mpx_usubn(mpw /dv/, mpw /dvl/, mpw /n/);
dd517851	423
3c9ede17	424	/* --- @mpx_umul@ --- *
	425	*
	426	* Arguments: @mpw dv, dvl@ = destination vector base and limit
	427	* @const mpw av, avl@ = multiplicand vector base and limit
	428	* @const mpw bv, bvl@ = multiplier vector base and limit
	429	*
	430	* Returns: ---
	431	*
	432	* Use: Performs unsigned integer multiplication. If the result
	433	* overflows the desination vector, high-order bits are
	434	* discarded. The result vector may not overlap the argument
	435	* vectors in any way.
	436	*/
	437
	438	extern void mpx_umul(mpw /dv/, mpw /dvl/,
	439	const mpw /av/, const mpw /avl/,
	440	const mpw /bv/, const mpw /bvl/);
	441
dd517851	442	/* --- @mpx_umuln@ --- *
d03ab969	443	*
dd517851	444	* Arguments: @mpw dv, dvl@ = destination vector base and limit
	445	* @const mpw av, avl@ = multiplicand vector base and limit
	446	* @mpw m@ = multiplier
	447	*
	448	* Returns: ---
d03ab969	449	*
	450	* Use: Multiplies a multiprecision integer by a single-word value.
	451	* The destination and source may be equal. The destination
	452	* is completely cleared after use.
	453	*/
	454
	455	#define MPX_UMULN(dv, dvl, av, avl, m) do { \
	456	mpw _dv = (dv), _dvl = (dvl); \
	457	const mpw _av = (av), _avl = (avl); \
	458	mpw _c = 0; \
	459	mpd _m = (m); \
	460	\
	461	while (_av < _avl) { \
	462	mpd _x; \
	463	if (_dv >= _dvl) \
	464	break; \
3c9ede17	465	_x = (mpd)_m * (mpd)*_av++ + _c; \
d03ab969	466	*_dv++ = MPW(_x); \
	467	_c = _x >> MPW_BITS; \
	468	} \
	469	if (_dv < _dvl) { \
	470	*_dv++ = MPW(_c); \
	471	MPX_ZERO(_dv, _dvl); \
	472	} \
	473	} while (0)
	474
dd517851	475	extern void mpx_umuln(mpw /dv/, mpw /dvl/,
	476	const mpw /av/, const mpw /avl/, mpw m);
	477
	478	/* --- @mpx_umlan@ --- *
d03ab969	479	*
dd517851	480	* Arguments: @mpw dv, dvl@ = destination/accumulator base and limit
	481	* @const mpw av, avl@ = multiplicand vector base and limit
	482	* @mpw m@ = multiplier
	483	*
	484	* Returns: ---
d03ab969	485	*
	486	* Use: Multiplies a multiprecision integer by a single-word value
	487	* and adds the result to an accumulator.
	488	*/
	489
	490	#define MPX_UMLAN(dv, dvl, av, avl, m) do { \
	491	mpw _dv = (dv), _dvl = (dvl); \
	492	const mpw _av = (av), _avl = (avl); \
7c13f461	493	mpw _cc = 0; \
d03ab969	494	mpd _m = (m); \
d03ab969	495	\
5bf74dea	496	while (_dv < _dvl && _av < _avl) { \
d03ab969	497	mpd _x; \
7c13f461	498	_x = (mpd)_dv + (mpd)_m (mpd)*_av++ + _cc; \
d03ab969	499	*_dv++ = MPW(_x); \
7c13f461	500	_cc = _x >> MPW_BITS; \
d03ab969	501	} \
7c13f461	502	MPX_UADDN(_dv, _dvl, _cc); \
d03ab969	503	} while (0)
d03ab969	504
dd517851	505	extern void mpx_umlan(mpw /dv/, mpw /dvl/,
	506	const mpw /av/, const mpw /avl/, mpw m);
	507
3c9ede17	508	/* --- @mpx_usqr@ --- *
d03ab969	509	*
d03ab969	510	* Arguments: @mpw dv, dvl@ = destination vector base and limit
3c9ede17	511	* @const mpw av, av@ = source vector base and limit
d03ab969	512	*
	513	* Returns: ---
	514	*
3c9ede17	515	* Use: Performs unsigned integer squaring. The result vector must
3c9ede17	516	* not overlap the source vector in any way.
d03ab969	517	*/
d03ab969	518
3c9ede17	519	extern void mpx_usqr(mpw /dv/, mpw /dvl/,
3c9ede17	520	const mpw /av/, const mpw /avl/);
d03ab969	521
5bf74dea	522	/* --- @mpx_udiv@ --- *
	523	*
	524	* Arguments: @mpw qv, qvl@ = quotient vector base and limit
	525	* @mpw rv, rvl@ = dividend/remainder vector base and limit
	526	* @const mpw dv, dvl@ = divisor vector base and limit
	527	* @mpw sv, svl@ = scratch workspace
	528	*
	529	* Returns: ---
	530	*
	531	* Use: Performs unsigned integer division. If the result overflows
	532	* the quotient vector, high-order bits are discarded. (Clearly
	533	* the remainder vector can't overflow.) The various vectors
	534	* may not overlap in any way. Yes, I know it's a bit odd
	535	* requiring the dividend to be in the result position but it
	536	* does make some sense really. The remainder must have
	537	* headroom for at least two extra words. The scratch space
	538	* must be at least one word larger than the divisor.
	539	*/
	540
	541	extern void mpx_udiv(mpw /qv/, mpw /qvl/, mpw /rv/, mpw /rvl/,
	542	const mpw /dv/, const mpw /dvl/,
	543	mpw /sv/, mpw /svl/);
	544
698bd937	545	/* --- @mpx_udivn@ --- *
	546	*
	547	* Arguments: @mpw qv, qvl@ = storage for the quotient (may overlap
	548	* dividend)
	549	* @const mpw rv, rvl@ = dividend
	550	* @mpw d@ = single-precision divisor
	551	*
	552	* Returns: Remainder after divison.
	553	*
	554	* Use: Performs a single-precision division operation.
	555	*/
	556
	557	extern mpw mpx_udivn(mpw /qv/, mpw /qvl/,
	558	const mpw /rv/, const mpw /rvl/, mpw /d/);
	559
5bf74dea	560	/----- Karatsuba multiplication algorithms -------------------------------/
	561
	562	/* --- @KARATSUBA_CUTOFF@ --- *
	563	*
	564	* This is the limiting length for using Karatsuba algorithms. It's best to
	565	* use the simpler classical multiplication method on numbers smaller than
	566	* this.
	567	*/
	568
	569	#define KARATSUBA_CUTOFF 16
	570
	571	/* --- @KARATSUBA_SLOP@ --- *
	572	*
	573	* The extra number of words required as scratch space by the Karatsuba
	574	* routines. This is a (generous) guess, since the actual amount of space
	575	* required is proportional to the recursion depth.
	576	*/
	577
698bd937	578	#define KARATSUBA_SLOP 64
5bf74dea	579
a86e33af	580	/* --- @mpx_kmul@ --- *
	581	*
	582	* Arguments: @mpw dv, dvl@ = pointer to destination buffer
	583	* @const mpw av, avl@ = pointer to first argument
	584	* @const mpw bv, bvl@ = pointer to second argument
	585	* @mpw sv, svl@ = pointer to scratch workspace
	586	*
	587	* Returns: ---
	588	*
	589	* Use: Multiplies two multiprecision integers using Karatsuba's
	590	* algorithm. This is rather faster than traditional long
	591	* multiplication (e.g., @mpx_umul@) on large numbers, although
	592	* more expensive on small ones.
	593	*
	594	* The destination and scratch buffers must be twice as large as
5bf74dea	595	* the larger argument. The scratch space must be twice as
	596	* large as the larger argument, plus the magic number
	597	* @KARATSUBA_SLOP@.
a86e33af	598	*/
a86e33af	599
a86e33af	600	extern void mpx_kmul(mpw /dv/, mpw /dvl/,
	601	const mpw /av/, const mpw /avl/,
	602	const mpw /bv/, const mpw /bvl/,
	603	mpw /sv/, mpw /svl/);
	604
5bf74dea	605	/* --- @mpx_ksqr@ --- *
d03ab969	606	*
5bf74dea	607	* Arguments: @mpw dv, dvl@ = pointer to destination buffer
	608	* @const mpw av, avl@ = pointer to first argument
	609	* @mpw sv, svl@ = pointer to scratch workspace
d03ab969	610	*
	611	* Returns: ---
	612	*
5bf74dea	613	* Use: Squares a multiprecision integers using something similar to
	614	* Karatsuba's multiplication algorithm. This is rather faster
	615	* than traditional long multiplication (e.g., @mpx_umul@) on
	616	* large numbers, although more expensive on small ones, and
	617	* rather simpler than full-blown Karatsuba multiplication.
	618	*
	619	* The destination must be twice as large as the argument. The
	620	* scratch space must be twice as large as the argument, plus
	621	* the magic number @KARATSUBA_SLOP@.
d03ab969	622	*/
d03ab969	623
5bf74dea	624	extern void mpx_ksqr(mpw /dv/, mpw /dvl/,
5bf74dea	625	const mpw /av/, const mpw /avl/,
3c9ede17	626	mpw /sv/, mpw /svl/);
d03ab969	627
	628	/----- That's all, folks -------------------------------------------------/
	629
	630	#ifdef __cplusplus
	631	}
	632	#endif
	633
	634	#endif