Moved the Karatsuba macros into a separate file for better sharing.

[u/mdw/catacomb] / mp.h

/* -*-c-*-
 *
 * $Id: mp.h,v 1.6 1999/12/10 23:19:46 mdw Exp $
 *
 * Simple 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: mp.h,v $
 * Revision 1.6  1999/12/10 23:19:46  mdw
 * Minor bugfixes.  New interface for suggested destinations.
 *
 * Revision 1.5  1999/11/22 20:50:37  mdw
 * Add support for computing Jacobi symbols.
 *
 * Revision 1.4  1999/11/21 22:13:02  mdw
 * Add mp version of MPX_BITS.
 *
 * Revision 1.3  1999/11/19 13:19:14  mdw
 * Fix const annotation.
 *
 * Revision 1.2  1999/11/17 18:02:16  mdw
 * New multiprecision integer arithmetic suite.
 *
 */

#ifndef CATACOMB_MP_H
#define CATACOMB_MP_H

#ifdef __cplusplus
  extern "C" {
#endif

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

#include <assert.h>
#include <string.h>

#include <mLib/sub.h>

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

#ifndef CATACOMB_MPX_H
#  include "mpx.h"
#endif

/*----- Data structures ---------------------------------------------------*/

typedef struct mp {
  mpw *v, *vl;
  size_t sz;
  unsigned f;
  unsigned ref;
} mp;

#define MP_NEG 1u
#define MP_BURN 2u
#define MP_CONST 4u
#define MP_UNDEF 8u
#define MP_DESTROYED 16u

/*----- Useful constants --------------------------------------------------*/

extern mp mp_const[];

#define MP_ZERO	 (&mp_const[0])
#define MP_ONE   (&mp_const[1])
#define MP_TWO   (&mp_const[2])
#define MP_THREE (&mp_const[3])
#define MP_FOUR  (&mp_const[4])
#define MP_FIVE  (&mp_const[5])
#define MP_TEN   (&mp_const[6])
#define MP_MONE  (&mp_const[7])

#define MP_NEW ((mp *)0)

/*----- Memory allocation hooks -------------------------------------------*/

#ifndef CATACOMB_MPARENA_H
#  include "mparena.h"
#endif

/* --- @MP_ARENA@ --- *
 *
 * This selects where memory is allocated from.  Tweak to use more fancy
 * things like custom arenas.
 */

#ifndef MP_ARENA
#  define MP_ARENA MPARENA_GLOBAL
#endif

/* --- @MP_ALLOC@ --- *
 *
 * Arguments:	@size_t sz@ = size required
 *
 * Returns:	Pointer to an allocated vector of the requested size.
 *
 * Use:		Hook for vector allocation.
 */

#ifndef MP_ALLOC
#  define MP_ALLOC(sz) mpalloc(MP_ARENA, (sz))
#endif

/* --- @MP_FREE@ --- *
 *
 * Arguments:	@mpw *v@ = pointer to vector
 *
 * Returns:	---
 *
 * Use:		Hook for vector deallocation.
 */

#ifndef MP_FREE
#  define MP_FREE(v) mpfree(MP_ARENA, (v))
#endif

/*----- Paranoia management -----------------------------------------------*/

/* --- @mp_burn@ --- *
 *
 * Arguments:	@mp *m@ = pointer to a multiprecision integer
 *
 * Returns:	---
 *
 * Use:		Marks the integer as `burn-after-use'.  When the integer's
 *		memory is deallocated, it is deleted so that traces can't
 *		remain in the swap file.  In theory.
 */

extern void mp_burn(mp */*m*/);

/*----- Trivial macros ----------------------------------------------------*/

/* --- @MP_LEN@ --- *
 *
 * Arguments:	@mp *m@ = pointer to a multiprecision integer
 *
 * Returns:	Length of the integer, in words.
 */

#define MP_LEN(m) ((m)->vl - ((m)->v))

/*----- Memory management and reference counting --------------------------*/

/* --- @mp_create@ --- *
 *
 * Arguments:	@size_t sz@ = size of vector required
 *
 * Returns:	Pointer to pristine new MP structure with enough memory
 *		bolted onto it.
 *
 * Use:		Creates a new multiprecision integer with indeterminate
 *		contents.  The integer has a single reference.
 */

extern mp *mp_create(size_t /*sz*/);

/* --- @mp_build@ --- *
 *
 * Arguments:	@mp *m@ = pointer to an MP block to fill in
 *		@mpw *v@ = pointer to a word array
 *		@mpw *vl@ = pointer just past end of array
 *
 * Returns:	---
 *
 * Use:		Creates a multiprecision integer representing some smallish
 *		number.  You must provide storage for the number and dispose
 *		of it when you've finished with it.  The number is marked as
 *		constant while it exists.
 */

extern void mp_build(mp */*m*/, mpw */*v*/, mpw */*vl*/);

/* --- @mp_destroy@ --- *
 *
 * Arguments:	@mp *m@ = pointer to a multiprecision integer
 *
 * Returns:	---
 *
 * Use:		Destroys a multiprecision integer. The reference count isn't
 *		checked.  Don't use this function if you don't know what
 *		you're doing: use @mp_drop@ instead.
 */

extern void mp_destroy(mp */*m*/);

/* --- @mp_copy@ --- *
 *
 * Arguments:	@mp *m@ = pointer to a multiprecision integer
 *
 * Returns:	A copy of the given multiprecision integer.
 *
 * Use:		Copies the given integer.  In fact you just get another
 *		reference to the same old one again.
 */

extern mp *mp_copy(mp */*m*/);

#define MP_COPY(m) ((m)->ref++, (m))

/* --- @mp_drop@ --- *
 *
 * Arguments:	@mp *m@ = pointer to a multiprecision integer
 *
 * Returns:	---
 *
 * Use:		Drops a reference to an integer which isn't wanted any more.
 *		If there are no more references, the integer is destroyed.
 */

extern void mp_drop(mp */*m*/);

#define MP_DROP(m) do {							\
  mp *_mm = (m);							\
  if (_mm->ref > 1)							\
    _mm->ref--;								\
  else if (!(_mm->f & MP_CONST))					\
    mp_destroy(_mm);							\
} while (0)
		   
/* --- @mp_split@ --- *
 *
 * Arguments:	@mp *m@ = pointer to a multiprecision integer
 *
 * Returns:	A reference to the same integer, possibly with a different
 *		address.
 *
 * Use:		Splits off a modifiable version of the integer referred to.
 */

extern mp *mp_split(mp */*m*/);

#define MP_SPLIT(m) do {						\
  mp *_mm = (m);							\
  if ((_mm->f & MP_CONST) || _mm->ref != 1) {				\
    mp *_dd = mp_create(_mm->sz);					\
    _dd->vl = _dd->v + MP_LEN(_mm);					\
    _dd->f = _mm->f & (MP_NEG | MP_BURN);				\
    memcpy(_dd->v, _mm->v, MPWS(MP_LEN(_mm)));				\
    _dd->ref = 1;							\
    _mm->ref--;								\
    (m) = _dd;								\
  }									\
} while (0)

/* --- @mp_resize@ --- *
 *
 * Arguments:	@mp *m@ = pointer to a multiprecision integer
 *		@size_t sz@ = new size
 *
 * Returns:	---
 *
 * Use:		Resizes the vector containing the integer's digits.  The new
 *		size must be at least as large as the current integer's
 *		length.  The integer's length is increased and new digits are
 *		filled with zeroes.  This isn't really intended for client
 *		use.
 */

extern void mp_resize(mp */*m*/, size_t /*sz*/);

#define MP_RESIZE(m, ssz) do {						\
  mp *_m = (m);								\
  size_t _sz = (ssz);							\
  size_t _len = MP_LEN(_m);						\
  mpw *_v = MP_ALLOC(_sz);						\
  if (!(_m->f & MP_UNDEF))						\
    memcpy(_v, _m->v, MPWS(_len));					\
  if (_m->f & MP_BURN)							\
    memset(_m->v, 0, MPWS(_m->sz));					\
  MP_FREE(_m->v);							\
  _m->v = _v;								\
  _m->vl = _v + _len;							\
  _m->sz = _sz;								\
} while (0)

/* --- @mp_ensure@ --- *
 *
 * Arguments:	@mp *m@ = pointer to a multiprecision integer
 *		@size_t sz@ = required size
 *
 * Returns:	---
 *
 * Use:		Ensures that the integer has enough space for @sz@ digits.
 *		The value is not changed.
 */

extern void mp_ensure(mp */*m*/, size_t /*sz*/);

#define MP_ENSURE(m, ssz) do {						\
  mp *_mm = (m);							\
  size_t _ssz = (ssz);							\
  size_t _len = MP_LEN(_mm);						\
  if (_ssz > _mm->sz)							\
    MP_RESIZE(_mm, _ssz);						\
  if (!(_mm->f & MP_UNDEF) && _ssz > _len)				\
    memset(_mm->vl, 0, MPWS(_ssz - _len));				\
  _mm->vl = _mm->v + _ssz;						\
} while (0)

/* --- @mp_modify@ --- *
 *
 * Arguments:	@mp *m@ = pointer to a multiprecision integer
 *		@size_t sz@ = size required
 *
 * Returns:	Pointer to the integer (possibly different).
 *
 * Use:		Prepares an integer to be overwritten.  It's split off from
 *		other references to the same integer, and sufficient space is
 *		allocated.
 */

extern mp *mp_modify(mp */*m*/, size_t /*sz*/);

#define MP_MODIFY(m, sz) do {						\
  size_t _rq = (sz);							\
  mp *_m = (m);								\
  if (_m == MP_NEW || m->ref > 1 || (_m->f & MP_CONST))	{		\
    if (_m)								\
      MP_DROP(_m);							\
    _m = mp_create(_rq);						\
  } else								\
    MP_ENSURE(_m, _rq);							\
  (m) = _m;								\
} while (0)

/*----- Size manipulation -------------------------------------------------*/

/* --- @mp_shrink@ --- *
 *
 * Arguments:	@mp *m@ = pointer to a multiprecision integer
 *
 * Returns:	---
 *
 * Use:		Reduces the recorded length of an integer.  This doesn't
 *		reduce the amount of memory used, although it can improve
 *		performance a bit.  To reduce memory, use @mp_minimize@
 *		instead.  This can't change the value of an integer, and is
 *		therefore safe to use even when there are multiple
 *		references.
 */

extern void mp_shrink(mp */*m*/);

#define MP_SHRINK(m) do {						\
  mp *_mm = (m);							\
  MPX_SHRINK(_mm->v, _mm->vl);						\
  if (!MP_LEN(_mm))							\
    _mm->f &= ~MP_NEG;							\
} while (0)

/* --- @mp_minimize@ --- *
 *
 * Arguments:	@mp *m@ = pointer to a multiprecision integer
 *
 * Returns:	---
 *
 * Use:		Reduces the amount of memory an integer uses.  It's best to
 *		do this to numbers which aren't going to change in the
 *		future.
 */

extern void mp_minimize(mp */*m*/);

/*----- Bit scanning ------------------------------------------------------*/

#ifndef CATACOMB_MPSCAN_H
#  include "mpscan.h"
#endif

/* --- @mp_scan@ --- *
 *
 * Arguments:	@mpscan *sc@ = pointer to bitscanner block
 *		@const mp *m@ = pointer to a multiprecision integer
 *
 * Returns:	---
 *
 * Use:		Initializes a bitscanner on a multiprecision integer.
 */

extern void mp_scan(mpscan */*sc*/, const mp */*m*/);

#define MP_SCAN(sc, m) do {						\
  const mp *_mm = (m);							\
  mpscan *_sc = (sc);							\
  MPSCAN_INITX(_sc, _mm->v, _mm->vl);					\
} while (0)

/* --- Other bitscanning aliases --- */

#define mp_step mpscan_step
#define mp_bit mpscan_bit

#define MP_STEP MPSCAN_STEP
#define MP_BIT MPSCAN_BIT

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

/* --- @mp_octets@ --- *
 *
 * Arguments:	@const mp *m@ = a multiprecision integer
 *
 * Returns:	The number of octets required to represent @m@.
 *
 * Use:		Calculates the external storage required for a multiprecision
 *		integer.
 */

extern size_t mp_octets(const mp */*m*/);

/* --- @mp_bits@ --- *
 *
 * Arguments:	@const mp *m@ = a multiprecision integer
 *
 * Returns:	The number of bits required to represent @m@.
 *
 * Use:		Calculates the external storage required for a multiprecision
 *		integer.
 */

extern unsigned long mp_bits(const mp */*m*/);

/* --- @mp_loadl@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@const void *pv@ = pointer to source data
 *		@size_t sz@ = size of the source data
 *
 * Returns:	Resulting multiprecision number.
 *
 * Use:		Loads a multiprecision number from an array of octets.  The
 *		first byte in the array is the least significant.  More
 *		formally, if the bytes are %$b_0, b_1, \ldots, b_{n-1}$%
 *		then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%.
 */

extern mp *mp_loadl(mp */*d*/, const void */*pv*/, size_t /*sz*/);

/* --- @mp_storel@ --- *
 *
 * Arguments:	@const mp *m@ = source
 *		@void *pv@ = pointer to output array
 *		@size_t sz@ = size of the output array
 *
 * Returns:	---
 *
 * Use:		Stores a multiprecision number in an array of octets.  The
 *		first byte in the array is the least significant.  If the
 *		array is too small to represent the number, high-order bits
 *		are truncated; if the array is too large, high order bytes
 *		are filled with zeros.  More formally, if the number is
 *		%$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%,
 *		then the array is %$b_0, b_1, \ldots, b_{n-1}$%.
 */

extern void mp_storel(const mp */*m*/, void */*pv*/, size_t /*sz*/);

/* --- @mp_loadb@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@const void *pv@ = pointer to source data
 *		@size_t sz@ = size of the source data
 *
 * Returns:	Resulting multiprecision number.
 *
 * Use:		Loads a multiprecision number from an array of octets.  The
 *		last byte in the array is the least significant.  More
 *		formally, if the bytes are %$b_{n-1}, b_{n-2}, \ldots, b_0$%
 *		then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%.
 */

extern mp *mp_loadb(mp */*d*/, const void */*pv*/, size_t /*sz*/);

/* --- @mp_storeb@ --- *
 *
 * Arguments:	@const mp *m@ = source
 *		@void *pv@ = pointer to output array
 *		@size_t sz@ = size of the output array
 *
 * Returns:	---
 *
 * Use:		Stores a multiprecision number in an array of octets.  The
 *		last byte in the array is the least significant.  If the
 *		array is too small to represent the number, high-order bits
 *		are truncated; if the array is too large, high order bytes
 *		are filled with zeros.  More formally, if the number is
 *		%$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%,
 *		then the array is %$b_{n-1}, b_{n-2}, \ldots, b_0$%.
 */

extern void mp_storeb(const mp */*m*/, void */*pv*/, size_t /*sz*/);

/*----- Simple arithmetic -------------------------------------------------*/

/* --- @mp_2c@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a@ = source
 *
 * Returns:	Result, @a@ converted to two's complement notation.
 */

extern mp *mp_2c(mp */*d*/, mp */*a*/);

/* --- @mp_sm@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a@ = source
 *
 * Returns:	Result, @a@ converted to the native signed-magnitude
 *		notation.
 */

extern mp *mp_sm(mp */*d*/, mp */*a*/);

/* --- @mp_lsl@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a@ = source
 *		@size_t n@ = number of bits to move
 *
 * Returns:	Result, @a@ shifted left by @n@.
 */

extern mp *mp_lsl(mp */*d*/, mp */*a*/, size_t /*n*/);

/* --- @mp_lsr@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a@ = source
 *		@size_t n@ = number of bits to move
 *
 * Returns:	Result, @a@ shifted left by @n@.
 */

extern mp *mp_lsr(mp */*d*/, mp */*a*/, size_t /*n*/);

/* --- @mp_cmp@ --- *
 *
 * Arguments:	@const mp *a, *b@ = two numbers
 *
 * Returns:	Less than, equal to or greater than zero, according to
 *		whether @a@ is less than, equal to or greater than @b@.
 */

extern int mp_cmp(const mp */*a*/, const mp */*b*/);

#define MP_CMP(a, op, b) (mp_cmp((a), (b)) op 0)

/* --- @mp_add@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a, *b@ = sources
 *
 * Returns:	Result, @a@ added to @b@.
 */

extern mp *mp_add(mp */*d*/, mp */*a*/, mp */*b*/);

/* --- @mp_sub@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a, *b@ = sources
 *
 * Returns:	Result, @b@ subtracted from @a@.
 */

extern mp *mp_sub(mp */*d*/, mp */*a*/, mp */*b*/);

/* --- @mp_mul@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a, *b@ = sources
 *
 * Returns:	Result, @a@ multiplied by @b@.
 */

extern mp *mp_mul(mp */*d*/, mp */*a*/, mp */*b*/);

/* --- @mp_sqr@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a@ = source
 *
 * Returns:	Result, @a@ squared.
 */

extern mp *mp_sqr(mp */*d*/, mp */*a*/);

/* --- @mp_div@ --- *
 *
 * Arguments:	@mp **qq, **rr@ = destination, quotient and remainder
 *		@mp *a, *b@ = sources
 *
 * Use:		Calculates the quotient and remainder when @a@ is divided by
 *		@b@.
 */

extern void mp_div(mp **/*qq*/, mp **/*rr*/, mp */*a*/, mp */*b*/);

/*----- More advanced algorithms ------------------------------------------*/

/* --- @mp_gcd@ --- *
 *
 * Arguments:	@mp **gcd, **xx, **yy@ = where to write the results
 *		@mp *a, *b@ = sources (must be nonzero)
 *
 * Returns:	---
 *
 * Use:		Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that
 *		@ax + by = gcd(a, b)@.  This is useful for computing modular
 *		inverses.  Neither @a@ nor @b@ may be zero.
 */

extern void mp_gcd(mp **/*gcd*/, mp **/*xx*/, mp **/*yy*/,
		   mp */*a*/, mp */*b*/);

/* --- @mp_jacobi@ --- *
 *
 * Arguments:	@mp *a@ = an integer less than @n@
 *		@mp *n@ = an odd integer
 *
 * Returns:	@-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
 *
 * Use:		Computes the Jacobi symbol.  If @n@ is prime, this is the
 *		Legendre symbol and is equal to 1 if and only if @a@ is a
 *		quadratic residue mod @n@.  The result is zero if and only if
 *		@a@ and @n@ have a common factor greater than one.
 */

int mp_jacobi(mp */*a*/, mp */*n*/);

/*----- Test harness support ----------------------------------------------*/

#include <mLib/testrig.h>

#ifndef CATACOMB_MPTEXT_H
#  include "mptext.h"
#endif

extern const test_type type_mp;

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

#ifdef __cplusplus
  }
#endif

#endif