# User Documentation for the IMath Library
Author: [M. J. Fromberger](https://github.com/creachadair)
## Installation
1. 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.
2. 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 `b`th 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:
```c
{
/* 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 ``, 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](#note2)); 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](https://github.com/creachadair/imath/issues/new), 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
1. 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()`.
2. 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](https://github.com/creachadair/imath/issues/new). 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
IMath is Copyright 2002-2009 Michael J. Fromberger
You may use it subject to the following Licensing Terms:
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.