45 KiB

Raw Blame History

User Documentation for the IMath Library

Installation

Edit Makefile to select compiler and options. The default is to use gcc. You may want to change CC to clang instead of gcc (and on macOS that what you will get anyway), but you should be able to use the default GCC settings for either.

By default, the Makefile assumes you can use 64-bit integer types, even though they were not standard in ANSI C90. If you cannot, add -DUSE_32BIT_WORDS to the compiler options.
Type make or make test to build the test driver and run the unit tests. None of these should fail. If they do, see below for how you can report bugs.

To build with debugging enabled (and optimization disabled), run make DEBUG=Y. This sets the preprocessor macro DEBUG to 1, and several other things (see Makefile for details).

To use the library in your code, include "imath.h" wherever you intend to use the library's routines. The integer library is just a single source file, so you can compile it into your project in whatever way makes sense. If you wish to use rational arithmetic, you will also need to include "imrat.h".

Background

The basic types defined by the imath library are mpz_t, an arbitrary precision signed integer, and mpq_t, an arbitrary precision signed rational number. The type mp_int is a pointer to an mpz_t, and mp_rat is a pointer to an mpq_t.

Most of the functions in the imath library return a value of type mp_result. This is a signed integer type which can be used to convey status information and also return small values. Any negative value is considered to be a status message. The following constants are defined for processing these:

Status	Description
`MP_OK`	operation successful, all is well (= 0)
`MP_FALSE`	boolean false (= `MP_OK`)
`MP_TRUE`	boolean true
`MP_MEMORY`	out of memory
`MP_RANGE`	parameter out of range
`MP_UNDEF`	result is undefined (e.g., division by zero)
`MP_TRUNC`	output value was truncated
`MP_BADARG`	an invalid parameter was passed

If you obtain a zero or negative value of an mp_result, you can use the mp_error_string() routine to obtain a pointer to a brief human-readable string describing the error. These strings are statically allocated, so they need not be freed by the caller; the same strings are re-used from call to call.

Unless otherwise noted, it is legal to use the same parameter for both inputs and output with most of the functions in this library. For example, you can add a number to itself and replace the original by writing:

mp_int_add(a, a, a);  /* a = a + a */

Any cases in which this is not legal will be noted in the function summaries below (if you discover that this is not so, please report it as a bug; I will fix either the function or the documentation :)

The IMath API

Each of the API functions is documented here. The general format of the entries is:

return_type function_name(parameters ...)
English description.

Unless otherwise noted, any API function that returns mp_result may be expected to return MP_OK, MP_BADARG, or MP_MEMORY. Other return values should be documented in the description. Please let me know if you discover this is not the case.

The following macros are defined in "imath.h", to define the sizes of the various data types used in the library:

Constant	Description
`MP_DIGIT_BIT`	the number of bits in a single `mpz_t` digit.
`MP_WORD_BIT`	the number of bits in a `mpz_t` word.
`MP_SMALL_MIN`	the minimum value representable by an `mp_small`.
`MP_SMALL_MAX`	the maximum value representable by an `mp_small`.
`MP_USMALL_MAX`	the maximum value representable by an `mp_usmall`.
`MP_MIN_RADIX`	the minimum radix accepted for base conversion.
`MP_MAX_RADIX`	the maximum radix accepted for base conversion.

Initialization

An mp_int must be initialized before use. By default, an mp_int is initialized with a certain minimum amount of storage for digits, and the storage is expanded automatically as needed. To initialize an mp_int, use the following functions:

mp_result mp_int_init(mp_int z);

Initializes z with 1-digit precision and sets it to zero. This function cannot fail unless z == NULL.

mp_int mp_int_alloc(void);

Allocates a fresh zero-valued mpz_t on the heap, returning NULL in case of error. The only possible error is out-of-memory.

mp_result mp_int_init_size(mp_int z, mp_size prec);

Initializes z with at least prec digits of storage, and sets it to zero. If prec is zero, the default precision is used. In either case the size is rounded up to the nearest multiple of the word size.

mp_result mp_int_init_copy(mp_int z, mp_int old);

Initializes z to be a copy of an already-initialized value in old. The new copy does not share storage with the original.

mp_result mp_int_init_value(mp_int z, mp_small value);

Initializes z to the specified signed value at default precision.

Cleanup

When you are finished with an mp_int, you must free the memory it uses:

void mp_int_clear(mp_int z);

Releases the storage used by z.

void mp_int_free(mp_int z);

Releases the storage used by z and also z itself. This should only be used for z allocated by mp_int_alloc().

Setting Values

To set an mp_int which has already been initialized to a small integer value, use:

mp_result mp_int_set_value(mp_int z, mp_small value);

Sets z to the value of the specified signed value.

mp_result mp_int_set_uvalue(mp_int z, mp_usmall uvalue);

Sets z to the value of the specified unsigned value.

To copy one initialized mp_int to another, use:

mp_result mp_int_copy(mp_int a, mp_int c);

Replaces the value of c with a copy of the value of a. No new memory is allocated unless a has more significant digits than c has allocated.

Arithmetic Functions

static inline bool mp_int_is_odd(mp_int z);

Reports whether z is odd, having remainder 1 when divided by 2.

static inline bool mp_int_is_even(mp_int z);

Reports whether z is even, having remainder 0 when divided by 2.

void mp_int_zero(mp_int z);

Sets z to zero. The allocated storage of z is not changed.

mp_result mp_int_abs(mp_int a, mp_int c);

Sets c to the absolute value of a.

mp_result mp_int_neg(mp_int a, mp_int c);

Sets c to the additive inverse (negation) of a.

mp_result mp_int_add(mp_int a, mp_int b, mp_int c);

Sets c to the sum of a and b.

mp_result mp_int_add_value(mp_int a, mp_small value, mp_int c);

Sets c to the sum of a and value.

mp_result mp_int_sub(mp_int a, mp_int b, mp_int c);

Sets c to the difference of a less b.

mp_result mp_int_sub_value(mp_int a, mp_small value, mp_int c);

Sets c to the difference of a less value.

mp_result mp_int_mul(mp_int a, mp_int b, mp_int c);

Sets c to the product of a and b.

mp_result mp_int_mul_value(mp_int a, mp_small value, mp_int c);

Sets c to the product of a and value.

mp_result mp_int_mul_pow2(mp_int a, mp_small p2, mp_int c);

Sets c to the product of a and 2^p2. Requires p2 >= 0.

mp_result mp_int_sqr(mp_int a, mp_int c);

Sets c to the square of a.

mp_result mp_int_root(mp_int a, mp_small b, mp_int c);

Sets c to the greatest integer not less than the bth root of a, using Newton's root-finding algorithm. It returns MP_UNDEF if a < 0 and b is even.

static inline mp_result mp_int_sqrt(mp_int a, mp_int c);

Sets c to the greatest integer not less than the square root of a. This is a special case of mp_int_root().

mp_result mp_int_div(mp_int a, mp_int b, mp_int q, mp_int r);

Sets q and r to the quotent and remainder of a / b. Division by powers of 2 is detected and handled efficiently. The remainder is pinned to 0 <= r < b.

Either of q or r may be NULL, but not both, and q and r may not point to the same value.

mp_result mp_int_div_value(mp_int a, mp_small value, mp_int q, mp_small *r);

Sets q and *r to the quotent and remainder of a / value. Division by powers of 2 is detected and handled efficiently. The remainder is pinned to 0 <= *r < b. Either of q or r may be NULL.

mp_result mp_int_div_pow2(mp_int a, mp_small p2, mp_int q, mp_int r);

Sets q and r to the quotient and remainder of a / 2^p2. This is a special case for division by powers of two that is more efficient than using ordinary division. Note that mp_int_div() will automatically handle this case, this function is for cases where you have only the exponent.

mp_result mp_int_mod(mp_int a, mp_int m, mp_int c);

Sets c to the remainder of a / m. The remainder is pinned to 0 <= c < m.

static inline mp_result mp_int_mod_value(mp_int a, mp_small value, mp_small* r);

Sets *r to the remainder of a / value. The remainder is pinned to 0 <= r < value.

mp_result mp_int_expt(mp_int a, mp_small b, mp_int c);

Sets c to the value of a raised to the b power. It returns MP_RANGE if b < 0.

mp_result mp_int_expt_value(mp_small a, mp_small b, mp_int c);

Sets c to the value of a raised to the b power. It returns MP_RANGE if b < 0.

mp_result mp_int_expt_full(mp_int a, mp_int b, mp_int c);

Sets c to the value of a raised to the b power. It returns MP_RANGE) if b < 0.

Comparison Functions

Unless otherwise specified, comparison between values x and y returns a comparator, an integer value < 0 if x is less than y, 0 if x is equal to y, and > 0 if x is greater than y.

int mp_int_compare(mp_int a, mp_int b);

Returns the comparator of a and b.

int mp_int_compare_unsigned(mp_int a, mp_int b);

Returns the comparator of the magnitudes of a and b, disregarding their signs. Neither a nor b is modified by the comparison.

int mp_int_compare_zero(mp_int z);

Returns the comparator of z and zero.

int mp_int_compare_value(mp_int z, mp_small v);

Returns the comparator of z and the signed value v.

int mp_int_compare_uvalue(mp_int z, mp_usmall uv);

Returns the comparator of z and the unsigned value uv.

bool mp_int_divisible_value(mp_int a, mp_small v);

Reports whether a is divisible by v.

int mp_int_is_pow2(mp_int z);

Returns k >= 0 such that z is 2^k, if such a k exists. If no such k exists, the function returns -1.

Modular Operations

mp_result mp_int_exptmod(mp_int a, mp_int b, mp_int m, mp_int c);

Sets c to the value of a raised to the b power, reduced modulo m. It returns MP_RANGE if b < 0 or MP_UNDEF if m == 0.

mp_result mp_int_exptmod_evalue(mp_int a, mp_small value, mp_int m, mp_int c);

Sets c to the value of a raised to the value power, modulo m. It returns MP_RANGE if value < 0 or MP_UNDEF if m == 0.

mp_result mp_int_exptmod_bvalue(mp_small value, mp_int b, mp_int m, mp_int c);

Sets c to the value of value raised to the b power, modulo m. It returns MP_RANGE if b < 0 or MP_UNDEF if m == 0.

mp_result mp_int_exptmod_known(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c);

Sets c to the value of a raised to the b power, reduced modulo m, given a precomputed reduction constant mu defined for Barrett's modular reduction algorithm.

It returns MP_RANGE if b < 0 or MP_UNDEF if m == 0.

mp_result mp_int_redux_const(mp_int m, mp_int c);

Sets c to the reduction constant for Barrett reduction by modulus m. Requires that c and m point to distinct locations.

mp_result mp_int_invmod(mp_int a, mp_int m, mp_int c);

Sets c to the multiplicative inverse of a modulo m, if it exists. The least non-negative representative of the congruence class is computed.

It returns MP_UNDEF if the inverse does not exist, or MP_RANGE if a == 0 or m <= 0.

mp_result mp_int_gcd(mp_int a, mp_int b, mp_int c);

Sets c to the greatest common divisor of a and b.

It returns MP_UNDEF if the GCD is undefined, such as for example if a and b are both zero.

mp_result mp_int_egcd(mp_int a, mp_int b, mp_int c, mp_int x, mp_int y);

Sets c to the greatest common divisor of a and b, and sets x and y to values satisfying Bezout's identity gcd(a, b) = ax + by.

It returns MP_UNDEF if the GCD is undefined, such as for example if a and b are both zero.

mp_result mp_int_lcm(mp_int a, mp_int b, mp_int c);

Sets c to the least common multiple of a and b.

It returns MP_UNDEF if the LCM is undefined, such as for example if a and b are both zero.

Conversion of Values

mp_result mp_int_to_int(mp_int z, mp_small *out);

Returns MP_OK if z is representable as mp_small, else MP_RANGE. If out is not NULL, *out is set to the value of z when MP_OK.

mp_result mp_int_to_uint(mp_int z, mp_usmall *out);

Returns MP_OK if z is representable as mp_usmall, or MP_RANGE. If out is not NULL, *out is set to the value of z when MP_OK.

mp_result mp_int_to_string(mp_int z, mp_size radix, char *str, int limit);

Converts z to a zero-terminated string of characters in the specified radix, writing at most limit characters to str including the terminating NUL value. A leading - is used to indicate a negative value.

Returns MP_TRUNC if limit was to small to write all of z. Requires MP_MIN_RADIX <= radix <= MP_MAX_RADIX.

mp_result mp_int_string_len(mp_int z, mp_size radix);

Reports the minimum number of characters required to represent z as a zero-terminated string in the given radix. Requires MP_MIN_RADIX <= radix <= MP_MAX_RADIX.

mp_result mp_int_read_string(mp_int z, mp_size radix, const char *str);

Reads a string of ASCII digits in the specified radix from the zero terminated str provided into z. For values of radix > 10, the letters A..Z or a..z are accepted. Letters are interpreted without respect to case.

Leading whitespace is ignored, and a leading + or - is interpreted as a sign flag. Processing stops when a NUL or any other character out of range for a digit in the given radix is encountered.

If the whole string was consumed, MP_OK is returned; otherwise MP_TRUNC. is returned.

Requires MP_MIN_RADIX <= radix <= MP_MAX_RADIX.

mp_result mp_int_read_cstring(mp_int z, mp_size radix, const char *str, char **end);

Reads a string of ASCII digits in the specified radix from the zero terminated str provided into z. For values of radix > 10, the letters A..Z or a..z are accepted. Letters are interpreted without respect to case.

Leading whitespace is ignored, and a leading + or - is interpreted as a sign flag. Processing stops when a NUL or any other character out of range for a digit in the given radix is encountered.

If the whole string was consumed, MP_OK is returned; otherwise MP_TRUNC. is returned. If end is not NULL, *end is set to point to the first unconsumed byte of the input string (the NUL byte if the whole string was consumed). This emulates the behavior of the standard C strtol() function.

Requires MP_MIN_RADIX <= radix <= MP_MAX_RADIX.

mp_result mp_int_count_bits(mp_int z);

Returns the number of significant bits in z.

mp_result mp_int_to_binary(mp_int z, unsigned char *buf, int limit);

Converts z to 2's complement binary, writing at most limit bytes into the given buf. Returns MP_TRUNC if the buffer limit was too small to contain the whole value. If this occurs, the contents of buf will be effectively garbage, as the function uses the buffer as scratch space.

The binary representation of z is in base-256 with digits ordered from most significant to least significant (network byte ordering). The high-order bit of the first byte is set for negative values, clear for non-negative values.

As a result, non-negative values will be padded with a leading zero byte if the high-order byte of the base-256 magnitude is set. This extra byte is accounted for by the mp_int_binary_len() function.

mp_result mp_int_read_binary(mp_int z, unsigned char *buf, int len);

Reads a 2's complement binary value from buf into z, where len is the length of the buffer. The contents of buf may be overwritten during processing, although they will be restored when the function returns.

mp_result mp_int_binary_len(mp_int z);

Returns the number of bytes to represent z in 2's complement binary.

mp_result mp_int_to_unsigned(mp_int z, unsigned char *buf, int limit);

Converts the magnitude of z to unsigned binary, writing at most limit bytes into the given buf. The sign of z is ignored, but z is not modified. Returns MP_TRUNC if the buffer limit was too small to contain the whole value. If this occurs, the contents of buf will be effectively garbage, as the function uses the buffer as scratch space during conversion.

The binary representation of z is in base-256 with digits ordered from most significant to least significant (network byte ordering).

mp_result mp_int_read_unsigned(mp_int z, unsigned char *buf, int len);

Reads an unsigned binary value from buf into z, where len is the length of the buffer. The contents of buf are not modified during processing.

mp_result mp_int_unsigned_len(mp_int z);

Returns the number of bytes required to represent z as an unsigned binary value in base 256.

Other Functions

Ordinarily, integer multiplication and squaring are done using the simple quadratic "schoolbook" algorithm. However, for sufficiently large values, there is a more efficient algorithm usually attributed to Karatsuba and Ofman that is usually faster. See Knuth Vol. 2 for more details about how this algorithm works.

The breakpoint between the "normal" and the recursive algorithm is controlled by a static digit threshold defined in imath.c. Values with fewer significant digits use the standard algorithm. This value can be modified by calling mp_int_multiply_threshold(n). The imtimer program and the findthreshold.py script (Python) can help you find a suitable value for for your particular platform.

const char *mp_error_string(mp_result res);

Returns a pointer to a brief, human-readable, zero-terminated string describing res. The returned string is statically allocated and must not be freed by the caller.

Rational Arithmetic

mp_result mp_rat_init(mp_rat r);

Initializes r with 1-digit precision and sets it to zero. This function cannot fail unless r is NULL.

mp_rat mp_rat_alloc(void);

Allocates a fresh zero-valued mpq_t on the heap, returning NULL in case of error. The only possible error is out-of-memory.

mp_result mp_rat_reduce(mp_rat r);

Reduces r in-place to lowest terms and canonical form.

Zero is represented as 0/1, one as 1/1, and signs are adjusted so that the sign of the value is carried by the numerator.

mp_result mp_rat_init_size(mp_rat r, mp_size n_prec, mp_size d_prec);

Initializes r with at least n_prec digits of storage for the numerator and d_prec digits of storage for the denominator, and value zero.

If either precision is zero, the default precision is used, rounded up to the nearest word size.

mp_result mp_rat_init_copy(mp_rat r, mp_rat old);

Initializes r to be a copy of an already-initialized value in old. The new copy does not share storage with the original.

mp_result mp_rat_set_value(mp_rat r, mp_small numer, mp_small denom);

Sets the value of r to the ratio of signed numer to signed denom. It returns MP_UNDEF if denom is zero.

mp_result mp_rat_set_uvalue(mp_rat r, mp_usmall numer, mp_usmall denom);

Sets the value of r to the ratio of unsigned numer to unsigned denom. It returns MP_UNDEF if denom is zero.

void mp_rat_clear(mp_rat r);

Releases the storage used by r.

void mp_rat_free(mp_rat r);

Releases the storage used by r and also r itself. This should only be used for r allocated by mp_rat_alloc().

mp_result mp_rat_numer(mp_rat r, mp_int z);

Sets z to a copy of the numerator of r.

mp_int mp_rat_numer_ref(mp_rat r);

Returns a pointer to the numerator of r.

mp_result mp_rat_denom(mp_rat r, mp_int z);

Sets z to a copy of the denominator of r.

mp_int mp_rat_denom_ref(mp_rat r);

Returns a pointer to the denominator of r.

mp_sign mp_rat_sign(mp_rat r);

Reports the sign of r.

mp_result mp_rat_copy(mp_rat a, mp_rat c);

Sets c to a copy of the value of a. No new memory is allocated unless a term of a has more significant digits than the corresponding term of c has allocated.

void mp_rat_zero(mp_rat r);

Sets r to zero. The allocated storage of r is not changed.

mp_result mp_rat_abs(mp_rat a, mp_rat c);

Sets c to the absolute value of a.

mp_result mp_rat_neg(mp_rat a, mp_rat c);

Sets c to the absolute value of a.

mp_result mp_rat_recip(mp_rat a, mp_rat c);

Sets c to the reciprocal of a if the reciprocal is defined. It returns MP_UNDEF if a is zero.

mp_result mp_rat_add(mp_rat a, mp_rat b, mp_rat c);

Sets c to the sum of a and b.

mp_result mp_rat_sub(mp_rat a, mp_rat b, mp_rat c);

Sets c to the difference of a less b.

mp_result mp_rat_mul(mp_rat a, mp_rat b, mp_rat c);

Sets c to the product of a and b.

mp_result mp_rat_div(mp_rat a, mp_rat b, mp_rat c);

Sets c to the ratio a / b if that ratio is defined. It returns MP_UNDEF if b is zero.

mp_result mp_rat_add_int(mp_rat a, mp_int b, mp_rat c);

Sets c to the sum of a and integer b.

mp_result mp_rat_sub_int(mp_rat a, mp_int b, mp_rat c);

Sets c to the difference of a less integer b.

mp_result mp_rat_mul_int(mp_rat a, mp_int b, mp_rat c);

Sets c to the product of a and integer b.

mp_result mp_rat_div_int(mp_rat a, mp_int b, mp_rat c);

Sets c to the ratio a / b if that ratio is defined. It returns MP_UNDEF if b is zero.

mp_result mp_rat_expt(mp_rat a, mp_small b, mp_rat c);

Sets c to the value of a raised to the b power. It returns MP_RANGE if b < 0.

int mp_rat_compare(mp_rat a, mp_rat b);

Returns the comparator of a and b.

int mp_rat_compare_unsigned(mp_rat a, mp_rat b);

Returns the comparator of the magnitudes of a and b, disregarding their signs. Neither a nor b is modified by the comparison.

int mp_rat_compare_zero(mp_rat r);

Returns the comparator of r and zero.

int mp_rat_compare_value(mp_rat r, mp_small n, mp_small d);

Returns the comparator of r and the signed ratio n / d. It returns MP_UNDEF if d is zero.

bool mp_rat_is_integer(mp_rat r);

Reports whether r is an integer, having canonical denominator 1.

mp_result mp_rat_to_ints(mp_rat r, mp_small *num, mp_small *den);

Reports whether the numerator and denominator of r can be represented as small signed integers, and if so stores the corresponding values to num and den. It returns MP_RANGE if either cannot be so represented.

mp_result mp_rat_to_string(mp_rat r, mp_size radix, char *str, int limit);

Converts r to a zero-terminated string of the format "n/d" with n and d in the specified radix and writing no more than limit bytes to the given output buffer str. The output of the numerator includes a sign flag if r is negative. Requires MP_MIN_RADIX <= radix <= MP_MAX_RADIX.

mp_result mp_rat_to_decimal(mp_rat r, mp_size radix, mp_size prec, mp_round_mode round, char *str, int limit);

Converts the value of r to a string in decimal-point notation with the specified radix, writing no more than limit bytes of data to the given output buffer. It generates prec digits of precision, and requires MP_MIN_RADIX <= radix <= MP_MAX_RADIX.

Ratios usually must be rounded when they are being converted for output as a decimal value. There are four rounding modes currently supported:

  MP_ROUND_DOWN
    Truncates the value toward zero.
    Example:  12.009 to 2dp becomes 12.00

  MP_ROUND_UP
    Rounds the value away from zero:
    Example:  12.001 to 2dp becomes 12.01, but
              12.000 to 2dp remains 12.00

  MP_ROUND_HALF_DOWN
     Rounds the value to nearest digit, half goes toward zero.
     Example:  12.005 to 2dp becomes 12.00, but
               12.006 to 2dp becomes 12.01

  MP_ROUND_HALF_UP
     Rounds the value to nearest digit, half rounds upward.
     Example:  12.005 to 2dp becomes 12.01, but
               12.004 to 2dp becomes 12.00

mp_result mp_rat_string_len(mp_rat r, mp_size radix);

Reports the minimum number of characters required to represent r as a zero-terminated string in the given radix. Requires MP_MIN_RADIX <= radix <= MP_MAX_RADIX.

mp_result mp_rat_decimal_len(mp_rat r, mp_size radix, mp_size prec);

Reports the length in bytes of the buffer needed to convert r using the mp_rat_to_decimal() function with the specified radix and prec. The buffer size estimate may slightly exceed the actual required capacity.

mp_result mp_rat_read_string(mp_rat r, mp_size radix, const char *str);

Sets r to the value represented by a zero-terminated string str in the format "n/d" including a sign flag. It returns MP_UNDEF if the encoded denominator has value zero.

mp_result mp_rat_read_cstring(mp_rat r, mp_size radix, const char *str, char **end);

Sets r to the value represented by a zero-terminated string str in the format "n/d" including a sign flag. It returns MP_UNDEF if the encoded denominator has value zero. If end is not NULL then *end is set to point to the first unconsumed character in the string, after parsing.

mp_result mp_rat_read_ustring(mp_rat r, mp_size radix, const char *str, char **end);

Sets r to the value represented by a zero-terminated string str having one of the following formats, each with an optional leading sign flag:
```
   n         : integer format, e.g. "123"
   n/d       : ratio format, e.g., "-12/5"
   z.ffff    : decimal format, e.g., "1.627"
```
It returns MP_UNDEF if the effective denominator is zero. If end is not NULL then *end is set to point to the first unconsumed character in the string, after parsing.

mp_result mp_rat_read_decimal(mp_rat r, mp_size radix, const char *str);

Sets r to the value represented by a zero-terminated string str in the format "z.ffff" including a sign flag. It returns MP_UNDEF if the effective denominator.

mp_result mp_rat_read_cdecimal(mp_rat r, mp_size radix, const char *str, char **end);

Sets r to the value represented by a zero-terminated string str in the format "z.ffff" including a sign flag. It returns MP_UNDEF if the effective denominator. If end is not NULL then *end is set to point to the first unconsumed character in the string, after parsing.

Representation Details

NOTE: You do not need to read this section to use IMath. This is provided for the benefit of developers wishing to extend or modify the internals of the library.

IMath uses a signed magnitude representation for arbitrary precision integers. The magnitude is represented as an array of radix-R digits in increasing order of significance; the value of R is chosen to be half the size of the largest available unsigned integer type, so typically 16 or 32 bits. Digits are represented as mp_digit, which must be an unsigned integral type.

Digit arrays are allocated using malloc(3) and realloc(3). Because this can be an expensive operation, the library takes pains to avoid allocation as much as possible. For this reason, the mpz_t structure distinguishes between how many digits are allocated and how many digits are actually consumed by the representation. The fields of an mpz_t are:

mp_digit    single;  /* single-digit value (see note) */
mp_digit   *digits;  /* array of digits               */
mp_size     alloc;   /* how many digits are allocated */
mp_size     used;    /* how many digits are in use    */
mp_sign     sign;    /* the sign of the value         */

The elements of digits at indices less than used are the significant figures of the value; the elements at indices greater than or equal to used are undefined (and may contain garbage). At all times, used must be at least 1 and at most alloc.

To avoid interaction with the memory allocator, single-digit values are stored directly in the mpz_t structure, in the single field. The semantics of access are the same as the more general case.

The number of digits allocated for an mpz_t is referred to in the library documentation as its "precision". Operations that affect an mpz_t cause precision to increase as needed. In any case, all allocations are measured in digits, and rounded up to the nearest mp_word boundary. There is a default minimum precision stored as a static constant default_precision (imath.c). This value can be set using mp_int_default_precision(n).

Note that the allocated size of an mpz_t can only grow; the library never reallocates in order to decrease the size. A simple way to do so explicitly is to use mp_int_init_copy(), as in:

mpz_t big, new;

/* ... */
mp_int_init_copy(&new, &big);
mp_int_swap(&new, &big);
mp_int_clear(&new);

The value of sign is 0 for positive values and zero, 1 for negative values. Constants MP_ZPOS and MP_NEG are defined for these; no other sign values are used.

If you are adding to this library, you should be careful to preserve the convention that inputs and outputs can overlap, as described above. So, for example, mp_int_add(a, a, a) is legal. Often, this means you must maintain one or more temporary mpz_t structures for intermediate values. The private macros DECLARE_TEMP(N), CLEANUP_TEMP(), and TEMP(K) can be used to maintain a conventional structure like this:

{
  /* Declare how many temp values you need.
	 Use TEMP(i) to access the ith value (0-indexed). */
  DECLARE_TEMP(8);
  ...

  /* Perform actions that must return MP_OK or fail. */
  REQUIRE(mp_int_copy(x, TEMP(1)));
  ...
  REQUIRE(mp_int_expt(TEMP(1), TEMP(2), TEMP(3)));
  ...

  /* You can also use REQUIRE directly for more complex cases. */
  if (some_difficult_question(TEMP(3)) != answer(x)) {
	REQUIRE(MP_RANGE);  /* falls through to cleanup (below) */
  }

  /* Ensure temporary values are cleaned up at exit.

     If control reaches here via a REQUIRE failure, the code below
	 the cleanup will not be executed.
   */
  CLEANUP_TEMP();
  return MP_OK;
}

Under the covers, these macros are just maintaining an array of mpz_t values, and a jump label to handle cleanup. You may only have one DECLARE_TEMP and its corresponding CLEANUP_TEMP per function body.

"Small" integer values are represented by the types mp_small and mp_usmall, which are mapped to appropriately-sized types on the host system. The default for mp_small is "long" and the default for mp_usmall is "unsigned long". You may change these, provided you insure that mp_small is signed and mp_usmall is unsigned. You will also need to adjust the size macros:

MP_SMALL_MIN, MP_SMALL_MAX
MP_USMALL_MIN, MP_USMALL_MAX

... which are defined in <imath.h>, if you change these.

Rational numbers are represented using a pair of arbitrary precision integers, with the convention that the sign of the numerator is the sign of the rational value, and that the result of any rational operation is always represented in lowest terms. The canonical representation for rational zero is 0/1. See "imrat.h".

Testing and Reporting of Bugs

Test vectors are included in the tests/ subdirectory of the imath distribution. When you run make test, it builds the imtest program and runs all available test vectors. If any tests fail, you will get a line like this:

x    y    FAILED      v

Here, x is the line number of the test which failed, y is index of the test within the file, and v is the value(s) actually computed. The name of the file is printed at the beginning of each test, so you can find out what test vector failed by executing the following (with x, y, and v replaced by the above values, and where "foo.t" is the name of the test file that was being processed at the time):

% tail +x tests/foo.t | head -1

None of the tests should fail (but see Note 2); if any do, it probably indicates a bug in the library (or at the very least, some assumption I made which I shouldn't have). Please file an issue, including the FAILED test line(s), as well as the output of the above tail command (so I know what inputs caused the failure).

If you build with the preprocessor symbol DEBUG defined as a positive integer, the digit allocators (s_alloc, s_realloc) fill all new buffers with the value 0xdeadbeefabad1dea, or as much of it as will fit in a digit, so that you can more easily catch uninitialized reads in the debugger.

Notes

You can generally use the same variables for both input and output. One exception is that you may not use the same variable for both the quotient and the remainder of mp_int_div().
Many of the tests for this library were written under the assumption that the mp_small type is 32 bits or more. If you compile with a smaller type, you may see MP_RANGE errors in some of the tests that otherwise pass (due to conversion failures). Also, the pi generator (pi.c) will not work correctly if mp_small is too short, as its algorithm for arc tangent is fairly simple-minded.

Contacts

The IMath library was written by Michael J. Fromberger.

If you discover any bugs or testing failures, please open an issue. Please be sure to include a complete description of what went wrong, and if possible, a test vector for imtest and/or a minimal test program that will demonstrate the bug on your system. Please also let me know what hardware, operating system, and compiler you're using.

Acknowledgements

The algorithms used in this library came from Vol. 2 of Donald Knuth's "The Art of Computer Programming" (Seminumerical Algorithms). Thanks to Nelson Bolyard, Bryan Olson, Tom St. Denis, Tushar Udeshi, and Eric Silva for excellent feedback on earlier versions of this code. Special thanks to Jonathan Shapiro for some very helpful design advice, as well as feedback and some clever ideas for improving performance in some common use cases.

License and Disclaimers

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

45 KiB Raw Blame History