progs/perftest.c: Use from Glibc syscall numbers.

[catacomb] / math / mp.h

/* -*-c-*-
 *
 * 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.
 */

#ifndef CATACOMB_MP_H
#define CATACOMB_MP_H

#ifdef __cplusplus
  extern "C" {
#endif

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

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

#include <mLib/macros.h>
#include <mLib/sub.h>

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

#ifndef CATACOMB_ARENA_H
#  include "arena.h"
#endif

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

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

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

/* --- A multiprecision integer --- */

typedef struct mp {
  mpw *v, *vl;				/* Vector of digits, current limit */
  size_t sz;				/* Size of digit buffer in words */
  mparena *a;				/* Arena for buffer allocation */
  unsigned f;				/* Flags (see below) */
  unsigned ref;				/* Reference counter */
} mp;

#define MP_NEG 1u			/* Negative (signed magnitude) */
#define MP_BURN 2u			/* Secret (viral flag) */
#define MP_CONST 4u			/* Uses strange memory allocation */
#define MP_UNDEF 8u			/* Contains nothing interesting */
#define MP_DESTROYED 16u		/* Has been destroyed */

/* --- A factor for simultaneous exponentation --- *
 *
 * Used by the Montgomery and Barrett exponentiators.
 */

typedef struct mp_expfactor {
  mp *base;
  mp *exp;
} mp_expfactor;

/*----- 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_256	 (&mp_const[7])
#define MP_MONE	 (&mp_const[8])

#define MP_NEW ((mp *)0)
#define MP_NEWSEC (&mp_const[9])

/*----- 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_new@ --- *
 *
 * Arguments:	@size_t sz@ = size of vector required
 *		@unsigned f@ = flags to set
 *
 * Returns:	Pointer to a new MP structure.
 *
 * Use:		Allocates a new multiprecision integer.  The data space is
 *		allocated from either the standard global or secret arena,
 *		depending on the initial flags requested.
 */

extern mp *mp_new(size_t /*sz*/, unsigned /*f*/);

/* --- @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_createsecure@ --- *
 *
 * 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.  The integer's
 *		data space is allocated from the secure arena.  Its burn flag
 *		is set.
 */

extern mp *mp_createsecure(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) MUFFLE_WARNINGS_EXPR(GCC_WARNING("-Wunused-value"),	\
					((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);							\
  _mm->ref--;								\
  if (_mm->ref == 0 && !(_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 *_m = (m);								\
  if ((_m->f & MP_CONST) || _m->ref > 1) {				\
    size_t _len = MP_LEN(_m);						\
    mp *_mm = mp_new(_len, _m->f);					\
    if (!(_m->f & MP_UNDEF))						\
      memcpy(_mm->v, _m->v, MPWS(_len));				\
    _m->ref--;								\
    _m = _mm;								\
  }									\
  (m) = _m;								\
} 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.  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);							\
  mparena *_a = (_m->f & MP_BURN) ? MPARENA_SECURE : MPARENA_GLOBAL;	\
  mpw *_v;								\
  size_t _len = MP_LEN(_m);						\
  assert(((void)"can't make size less than length", _sz >= _len));	\
  _v = mpalloc(_a, _sz);						\
  if (!(_m->f & MP_UNDEF))						\
    memcpy(_v, _m->v, MPWS(_len));					\
  if (_m->f & MP_BURN)							\
    memset(_m->v, 0, MPWS(_m->sz));					\
  mpfree(_m->a, _m->v);							\
  _m->a = _a;								\
  _m->v = _v;								\
  _m->vl = _v + _len;							\
} 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 *_m = (m);								\
  size_t _ssz = (ssz);							\
  size_t _len = MP_LEN(_m);						\
  if (_ssz >= _len) {							\
    if (_ssz > _m->sz)							\
      mp_resize(_m, _ssz);						\
    if (!(_m->f & MP_UNDEF) && _ssz > _len)				\
      memset(_m->vl, 0, MPWS(_ssz - _len));				\
    _m->vl = _m->v + _ssz;						\
  }									\
} while (0)

/* --- @mp_dest@ --- *
 *
 * Arguments:	@mp *m@ = a suggested destination integer
 *		@size_t sz@ = size required for result, in digits
 *		@unsigned f@ = various flags
 *
 * Returns:	A pointer to an appropriate destination.
 *
 * Use:		Converts a suggested destination into a real destination with
 *		the required properties.  If the real destination is @d@,
 *		then the following properties will hold:
 *
 *		  * @d@ will have exactly one reference.
 *
 *		  * If @m@ is not @MP_NEW@, then the contents of @m@ will not
 *		    change, unless @f@ has the @MP_UNDEF@ flag set.
 *
 *		  * If @m@ is not @MP_NEW@, then he reference count of @m@ on
 *		    entry is equal to the sum of the counts of @d@ and @m@ on
 *		    exit.
 *
 *		  * The size of @d@ will be at least @sz@.
 *
 *		  * If @f@ has the @MP_BURN@ flag set, then @d@ will be
 *		    allocated from @MPARENA_SECURE@.
 *
 *		Understanding this function is crucial to using Catacomb's
 *		multiprecision integer library effectively.
 */

extern mp *mp_dest(mp */*m*/, size_t /*sz*/, unsigned /*f*/);

#define MP_DEST(m, ssz, f) do {						\
  mp *_m = (m);								\
  size_t _ssz = (ssz);							\
  unsigned _f = (f);							\
  _m = mp_dest(_m, _ssz, _f);						\
  (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_ZEROP(_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)

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

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

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

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

#define mp_step mpscan_step
#define mp_bit mpscan_bit
#define mp_rstep mpscan_rstep
#define mp_rbit mpscan_rbit

#define MP_STEP MPSCAN_STEP
#define MP_BIT MPSCAN_BIT
#define MP_RSTEP MPSCAN_RSTEP
#define MP_RBIT MPSCAN_RBIT

/*----- 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_octets2c@ --- *
 *
 * 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 represented as two's complement.
 */

extern size_t mp_octets2c(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*/);

/* --- @mp_loadl2c@ --- *
 *
 * 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 as
 *		two's complement.  The first byte in the array is the least
 *		significant.
 */

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

/* --- @mp_storel2c@ --- *
 *
 * 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 as two's
 *		complement.  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 sign-extended.
 */

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

/* --- @mp_loadb2c@ --- *
 *
 * 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 as
 *		two's complement.  The last byte in the array is the least
 *		significant.
 */

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

/* --- @mp_storeb2c@ --- *
 *
 * 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, as
 *		two's complement.  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 sign-extended.
 */

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

/*----- Bit operations ----------------------------------------------------*/

/* --- @mp_not@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a@ = source
 *
 * Returns:	The bitwise complement of the source.
 */

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

/* --- @mp_bitop@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a, *b@ = sources
 *
 * Returns:	The result of the given bitwise operation.  These functions
 *		don't handle negative numbers at all sensibly.  For that, use
 *		the @...2c@ variants.  The functions are named after the
 *		truth tables they generate:
 *
 *			a:	0011
 *			b:	0101
 *			@mpx_bitXXXX@
 */

#define MP_BITDECL(string)						\
  extern mp *mp_bit##string(mp */*d*/, mp */*a*/, mp */*b*/);
MPX_DOBIN(MP_BITDECL)

/* --- @mp_[n]and@, @mp_[n]or@, @mp_[n]xor@, @mp_not@ --- *
 *
 * Synonyms for the commonly-used functions.
 */

#define mp_and	mp_bit0001
#define mp_or	mp_bit0111
#define mp_nand mp_bit1110
#define mp_nor	mp_bit1000
#define mp_xor	mp_bit0110

/* --- @mp_testbit@ --- *
 *
 * Arguments:	@mp *x@ = a large integer
 *		@unsigned long n@ = which bit to test
 *
 * Returns:	Nonzero if the bit is set, zero if not.
 */

extern int mp_testbit(mp */*x*/, unsigned long /*n*/);

/* --- @mp_setbit@, @mp_clearbit@ --- *
 *
 * Arguments:	@mp *d@ = a destination
 *		@mp *x@ = a large integer
 *		@unsigned long n@ = which bit to modify
 *
 * Returns:	The argument @x@, with the appropriate bit set or cleared.
 */

extern mp *mp_setbit(mp */*d*/, mp */*x*/, unsigned long /*n*/);
extern mp *mp_clearbit(mp */*d*/, mp */*x*/, unsigned long /*n*/);

/* --- @mp_lsl@, @mp_lslc@, @mp_lsr@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a@ = source
 *		@size_t n@ = number of bits to move
 *
 * Returns:	Result, @a@ shifted left or right by @n@.
 *
 * Use:		Bitwise shift operators.  @mp_lslc@ fills the bits introduced
 *		on the right with ones instead of zeroes: it's used
 *		internally by @mp_lsl2c@, though it may be useful on its
 *		own.
 */

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

/* --- @mp_not2c@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a@ = source
 *
 * Returns:	The sign-extended complement of the argument.
 */

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

/* --- @mp_bitop2c@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a, *b@ = sources
 *
 * Returns:	The result of the given bitwise operation.  Negative numbers
 *		are treated as two's complement, sign-extended infinitely to
 *		the left.  The functions are named after the truth tables
 *		they generate:
 *
 *			a:	0011
 *			b:	0101
 *			@mpx_bitXXXX@
 */

#define MP_BIT2CDECL(string)						\
  extern mp *mp_bit##string##2c(mp */*d*/, mp */*a*/, mp */*b*/);
MPX_DOBIN(MP_BIT2CDECL)

/* --- @mp_[n]and@, @mp_[n]or@, @mp_[n]xor@, @mp_not@ --- *
 *
 * Synonyms for the commonly-used functions.
 */

#define mp_and2c  mp_bit00012c
#define mp_or2c	  mp_bit01112c
#define mp_nand2c mp_bit11102c
#define mp_nor2c  mp_bit10002c
#define mp_xor2c  mp_bit01102c

/* --- @mp_lsl2c@, @mp_lsr2c@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a@ = source
 *		@size_t n@ = number of bits to move
 *
 * Returns:	Result, @a@ shifted left or right by @n@.  Handles the
 *		pretence of sign-extension for negative numbers.
 */

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

/* --- @mp_testbit2c@ --- *
 *
 * Arguments:	@mp *x@ = a large integer
 *		@unsigned long n@ = which bit to test
 *
 * Returns:	Nonzero if the bit is set, zero if not.  Fakes up two's
 *		complement representation.
 */

extern int mp_testbit2c(mp */*x*/, unsigned long /*n*/);

/* --- @mp_setbit2c@, @mp_clearbit2c@ --- *
 *
 * Arguments:	@mp *d@ = a destination
 *		@mp *x@ = a large integer
 *		@unsigned long n@ = which bit to modify
 *
 * Returns:	The argument @x@, with the appropriate bit set or cleared.
 *		Fakes up two's complement representation.
 */

extern mp *mp_setbit2c(mp */*d*/, mp */*x*/, unsigned long /*n*/);
extern mp *mp_clearbit2c(mp */*d*/, mp */*x*/, unsigned long /*n*/);

/*----- Comparisons -------------------------------------------------------*/

/* --- @mp_eq@ --- *
 *
 * Arguments:	@const mp *a, *b@ = two numbers
 *
 * Returns:	Nonzero if the numbers are equal.
 */

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

#define MP_EQ(a, b)							\
  ((((a)->f ^ (b)->f) & MP_NEG) == 0 &&					\
   mpx_ueq((a)->v, (a)->vl, (b)->v, (b)->vl))

/* --- @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)

/* --- Other handy macros --- */

#define MP_NEGP(x) ((x)->f & MP_NEG)
#define MP_ZEROP(x) (!MP_LEN(x))
#define MP_POSP(x) (!MP_NEGP(x) && !MP_ZEROP(x))
#define MP_ODDP(x) (!MP_ZEROP(x) && ((x)->v[0] & 1u))
#define MP_EVENP(x) (!MP_ODDP(x))

/*----- Arithmetic operations ---------------------------------------------*/

/* --- @mp_neg@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *a@ = argument
 *
 * Returns:	The negation of the argument.
 *
 * Use:		Negates its argument.
 */

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

/* --- @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*/);

/* --- @mp_exp@ --- *
 *
 * Arguments:	@mp *d@ = fake destination
 *		@mp *a@ = base
 *		@mp *e@ = exponent
 *
 * Returns:	Result, %$a^e$%.
 */

extern mp *mp_exp(mp */*d*/, mp */*a*/, mp */*e*/);

/* --- @mp_odd@ --- *
 *
 * Arguments:	@mp *d@ = pointer to destination integer
 *		@mp *m@ = pointer to source integer
 *		@size_t *s@ = where to store the power of 2
 *
 * Returns:	An odd integer integer %$t$% such that %$m = 2^s t$%.
 *
 * Use:		Computes a power of two and an odd integer which, when
 *		multiplied, give a specified result.  This sort of thing is
 *		useful in number theory quite often.
 */

extern mp *mp_odd(mp */*d*/, mp */*m*/, size_t */*s*/);

/* --- @mp_leastcongruent@ --- *
 *
 * Arguments:	@mp *d@ = pointer to destination
 *		@mp *b@ = lower bound
 *		@mp *r@ = representative
 *		@mp *m@ = modulus
 *
 * Returns:	The smallest integer %$x \equiv r \pmod{m}$% such that
 *		%$x \ge b$%.
 */

extern mp *mp_leastcongruent(mp */*d*/, mp */*b*/, mp */*r*/, mp */*m*/);

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

/* --- @mp_sqrt@ --- *
 *
 * Arguments:	@mp *d@ = pointer to destination integer
 *		@mp *a@ = (nonnegative) integer to take square root of
 *
 * Returns:	The largest integer %$x$% such that %$x^2 \le a$%.
 *
 * Use:		Computes integer square roots.
 *
 *		The current implementation isn't very good: it uses the
 *		Newton-Raphson method to find an approximation to %$a$%.  If
 *		there's any demand for a better version, I'll write one.
 */

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

/* --- @mp_nthrt@ --- *
 *
 * Arguments:	@mp *d@ = fake destination
 *		@mp *a@ = an integer
 *		@mp *n@ = a strictly positive integer
 *		@int *exectp_out@ = set nonzero if an exact solution is found
 *
 * Returns:	The integer %$\bigl\lfloor \sqrt[n]{a} \bigr\rfloor$%.
 *
 * Use:		Return an approximation to the %$n$%th root of %$a$%.
 *		Specifically, it returns the largest integer %$x$% such that
 *		%$x^n \le a$%.  If %$x^n = a$% then @*exactp_out@ is set
 *		nonzero; otherwise it is set zero.  (If @exactp_out@ is null
 *		then this information is discarded.)
 *
 *		The exponent %$n$% must be strictly positive: it's not clear
 *		to me what the right answer is for %$n \le 0$%.  If %$a$% is
 *		negative then %$n$% must be odd; otherwise there is no real
 *		solution.
 */

extern mp *mp_nthrt(mp */*d*/, mp */*a*/, mp */*n*/, int */*exactp_out*/);

/* --- @mp_perfect_power_p@ --- *
 *
 * Arguments:	@mp **x@ = where to write the base
 *		@mp **n@ = where to write the exponent
 *		@mp *a@ = an integer
 *
 * Returns:	Nonzero if %$a$% is a perfect power.
 *
 * Use:		Returns whether an integer %$a$% is a perfect power, i.e.,
 *		whether it can be written in the form %$a = x^n$% where
 *		%$|x| > 1$% and %$n > 1$% are integers.  If this is possible,
 *		then (a) store %$x$% and the largest such %$n$% in @*x@ and
 *		@*n@, and return nonzero; otherwise, store %$x = a$% and
 *		%$n = 1$% and return zero.  (Either @x@ or @n@, or both, may
 *		be null to discard these outputs.)
 *
 *		Note that %$-1$%, %$0$% and %$1$% are not considered perfect
 *		powers by this definition.  (The exponent is not well-defined
 *		in these cases, but it seemed better to implement a function
 *		which worked for all integers.)  Note also that %$-4$% is not
 *		a perfect power since it has no real square root.
 */

extern int mp_perfect_power_p(mp **/*x*/, mp **/*n*/, mp */*a*/);

/* --- @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_modinv@ --- *
 *
 * Arguments:	@mp *d@ = destination
 *		@mp *x@ = argument
 *		@mp *p@ = modulus
 *
 * Returns:	The inverse %$x^{-1} \bmod p$%.
 *
 * Use:		Computes a modular inverse.    An assertion fails if %$p$%
 *		has no inverse.
 */

extern mp *mp_modinv(mp */*d*/, mp */*x*/, mp */*p*/);

/* --- @mp_jacobi@ --- *
 *
 * Arguments:	@mp *a@ = an integer
 *		@mp *n@ = another integer
 *
 * Returns:	@-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
 *
 * Use:		Computes the Kronecker symbol %$\jacobi{a}{n}$%.  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.
 *
 *		If @n@ is composite, then this computes the Kronecker symbol
 *
 *		  %$\jacobi{a}{n}=\jacobi{a}{u}\prod_i\jacobi{a}{p_i}^{e_i}$%
 *
 *		where %$n = u p_0^{e_0} \ldots p_{n-1}^{e_{n-1}}$% is the
 *		prime factorization of %$n$%.  The missing bits are:
 *
 *		  * %$\jacobi{a}{1} = 1$%;
 *		  * %$\jacobi{a}{-1} = 1$% if @a@ is negative, or 1 if
 *		    positive;
 *		  * %$\jacobi{a}{0} = 0$%;
 *		  * %$\jacobi{a}{2}$ is 0 if @a@ is even, 1 if @a@ is
 *		    congruent to 1 or 7 (mod 8), or %$-1$% otherwise.
 *
 *		If %$n$% is positive and odd, then this is the Jacobi
 *		symbol.  (The Kronecker symbol is a consistent domain
 *		extension; the Jacobi symbol was implemented first, and the
 *		name stuck.)
 */

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

/* --- @mp_modsqrt@ --- *
 *
 * Arguments:	@mp *d@ = destination integer
 *		@mp *a@ = source integer
 *		@mp *p@ = modulus (must be prime)
 *
 * Returns:	If %$a$% is a quadratic residue, a square root of %$a$%; else
 *		a null pointer.
 *
 * Use:		Returns an integer %$x$% such that %$x^2 \equiv a \pmod{p}$%,
 *		if one exists; else a null pointer.  This function will not
 *		work if %$p$% is composite: you must factor the modulus, take
 *		a square root mod each factor, and recombine the results
 *		using the Chinese Remainder Theorem.
 *
 *		We guarantee that the square root returned is the smallest
 *		one (i.e., the `positive' square root).
 */

extern mp *mp_modsqrt(mp */*d*/, mp */*a*/, mp */*p*/);

/* --- @mp_modexp@ --- *
 *
 * Arguments:	@mp *d@ = fake destination
 *		@mp *x@ = base of exponentiation
 *		@mp *e@ = exponent
 *		@mp *n@ = modulus (must be positive)
 *
 * Returns:	The value %$x^e \bmod n$%.
 */

extern mp *mp_modexp(mp */*d*/, mp */*x*/, mp */*e*/, 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