2 This file was generated from "doc.md.in" by mkdoc.py
6 # User Documentation for the IMath Library
8 Author: [M. J. Fromberger](https://github.com/creachadair)
12 1. Edit Makefile to select compiler and options. The default is to use gcc.
13 You may want to change CC to `clang` instead of `gcc` (and on macOS that
14 what you will get anyway), but you should be able to use the default GCC
17 By default, the Makefile assumes you can use 64-bit integer types, even
18 though they were not standard in ANSI C90. If you cannot, add
19 `-DUSE_32BIT_WORDS` to the compiler options.
21 2. Type `make` or `make test` to build the test driver and run the unit tests.
22 None of these should fail. If they do, see below for how you can report
25 To build with debugging enabled (and optimization disabled), run `make
26 DEBUG=Y`. This sets the preprocessor macro `DEBUG` to 1, and several other
27 things (see Makefile for details).
29 To use the library in your code, include "imath.h" wherever you intend to use
30 the library's routines. The integer library is just a single source file, so
31 you can compile it into your project in whatever way makes sense. If you wish
32 to use rational arithmetic, you will also need to include "imrat.h".
36 The basic types defined by the imath library are `mpz_t`, an arbitrary
37 precision signed integer, and `mpq_t`, an arbitrary precision signed rational
38 number. The type `mp_int` is a pointer to an `mpz_t`, and `mp_rat` is a
39 pointer to an `mpq_t`.
41 Most of the functions in the imath library return a value of type `mp_result`.
42 This is a signed integer type which can be used to convey status information
43 and also return small values. Any negative value is considered to be a status
44 message. The following constants are defined for processing these:
46 | Status | Description |
47 | ----------- | -------------------------------------------- |
48 | `MP_OK` | operation successful, all is well (= 0) |
49 | `MP_FALSE` | boolean false (= `MP_OK`) |
50 | `MP_TRUE` | boolean true |
51 | `MP_MEMORY` | out of memory |
52 | `MP_RANGE` | parameter out of range |
53 | `MP_UNDEF` | result is undefined (e.g., division by zero) |
54 | `MP_TRUNC` | output value was truncated |
55 | `MP_BADARG` | an invalid parameter was passed |
57 If you obtain a zero or negative value of an `mp_result`, you can use the
58 `mp_error_string()` routine to obtain a pointer to a brief human-readable
59 string describing the error. These strings are statically allocated, so they
60 need not be freed by the caller; the same strings are re-used from call to
63 Unless otherwise noted, it is legal to use the same parameter for both inputs
64 and output with most of the functions in this library. For example, you can
65 add a number to itself and replace the original by writing:
67 mp_int_add(a, a, a); /* a = a + a */
69 Any cases in which this is not legal will be noted in the function summaries
70 below (if you discover that this is not so, please report it as a bug; I will
71 fix either the function or the documentation :)
75 Each of the API functions is documented here. The general format of the
80 > return_type function_name(parameters ...)
82 > - English description.
84 Unless otherwise noted, any API function that returns `mp_result` may be
85 expected to return `MP_OK`, `MP_BADARG`, or `MP_MEMORY`. Other return values
86 should be documented in the description. Please let me know if you discover
89 The following macros are defined in "imath.h", to define the sizes of the
90 various data types used in the library:
92 | Constant | Description
93 | --------------- | ----------------------------------------
94 | `MP_DIGIT_BIT` | the number of bits in a single `mpz_t` digit.
95 | `MP_WORD_BIT` | the number of bits in a `mpz_t` word.
96 | `MP_SMALL_MIN` | the minimum value representable by an `mp_small`.
97 | `MP_SMALL_MAX` | the maximum value representable by an `mp_small`.
98 | `MP_USMALL_MAX` | the maximum value representable by an `mp_usmall`.
99 | `MP_MIN_RADIX` | the minimum radix accepted for base conversion.
100 | `MP_MAX_RADIX` | the maximum radix accepted for base conversion.
104 An `mp_int` must be initialized before use. By default, an `mp_int` is
105 initialized with a certain minimum amount of storage for digits, and the
106 storage is expanded automatically as needed. To initialize an `mp_int`, use
107 the following functions:
110 <a id="mp_int_init"></a><pre>
111 mp_result <a href="imath.h#L115">mp_int_init</a>(mp_int z);
113 - Initializes `z` with 1-digit precision and sets it to zero. This function
114 cannot fail unless `z == NULL`.
117 <a id="mp_int_alloc"></a><pre>
118 mp_int <a href="imath.h#L119">mp_int_alloc</a>(void);
120 - Allocates a fresh zero-valued `mpz_t` on the heap, returning NULL in case
121 of error. The only possible error is out-of-memory.
124 <a id="mp_int_init_size"></a><pre>
125 mp_result <a href="imath.h#L124">mp_int_init_size</a>(mp_int z, mp_size prec);
127 - Initializes `z` with at least `prec` digits of storage, and sets it to
128 zero. If `prec` is zero, the default precision is used. In either case the
129 size is rounded up to the nearest multiple of the word size.
132 <a id="mp_int_init_copy"></a><pre>
133 mp_result <a href="imath.h#L128">mp_int_init_copy</a>(mp_int z, mp_int old);
135 - Initializes `z` to be a copy of an already-initialized value in `old`. The
136 new copy does not share storage with the original.
139 <a id="mp_int_init_value"></a><pre>
140 mp_result <a href="imath.h#L131">mp_int_init_value</a>(mp_int z, mp_small value);
142 - Initializes `z` to the specified signed `value` at default precision.
148 When you are finished with an `mp_int`, you must free the memory it uses:
151 <a id="mp_int_clear"></a><pre>
152 void <a href="imath.h#L143">mp_int_clear</a>(mp_int z);
154 - Releases the storage used by `z`.
157 <a id="mp_int_free"></a><pre>
158 void <a href="imath.h#L147">mp_int_free</a>(mp_int z);
160 - Releases the storage used by `z` and also `z` itself.
161 This should only be used for `z` allocated by `mp_int_alloc()`.
167 To set an `mp_int` which has already been initialized to a small integer value,
171 <a id="mp_int_set_value"></a><pre>
172 mp_result <a href="imath.h#L137">mp_int_set_value</a>(mp_int z, mp_small value);
174 - Sets `z` to the value of the specified signed `value`.
177 <a id="mp_int_set_uvalue"></a><pre>
178 mp_result <a href="imath.h#L140">mp_int_set_uvalue</a>(mp_int z, mp_usmall uvalue);
180 - Sets `z` to the value of the specified unsigned `value`.
184 To copy one initialized `mp_int` to another, use:
187 <a id="mp_int_copy"></a><pre>
188 mp_result <a href="imath.h#L151">mp_int_copy</a>(mp_int a, mp_int c);
190 - Replaces the value of `c` with a copy of the value of `a`. No new memory is
191 allocated unless `a` has more significant digits than `c` has allocated.
195 ### Arithmetic Functions
198 <a id="mp_int_is_odd"></a><pre>
199 static inline bool <a href="imath.h#L108">mp_int_is_odd</a>(mp_int z);
201 - Reports whether `z` is odd, having remainder 1 when divided by 2.
204 <a id="mp_int_is_even"></a><pre>
205 static inline bool <a href="imath.h#L111">mp_int_is_even</a>(mp_int z);
207 - Reports whether `z` is even, having remainder 0 when divided by 2.
210 <a id="mp_int_zero"></a><pre>
211 void <a href="imath.h#L157">mp_int_zero</a>(mp_int z);
213 - Sets `z` to zero. The allocated storage of `z` is not changed.
216 <a id="mp_int_abs"></a><pre>
217 mp_result <a href="imath.h#L160">mp_int_abs</a>(mp_int a, mp_int c);
219 - Sets `c` to the absolute value of `a`.
222 <a id="mp_int_neg"></a><pre>
223 mp_result <a href="imath.h#L163">mp_int_neg</a>(mp_int a, mp_int c);
225 - Sets `c` to the additive inverse (negation) of `a`.
228 <a id="mp_int_add"></a><pre>
229 mp_result <a href="imath.h#L166">mp_int_add</a>(mp_int a, mp_int b, mp_int c);
231 - Sets `c` to the sum of `a` and `b`.
234 <a id="mp_int_add_value"></a><pre>
235 mp_result <a href="imath.h#L169">mp_int_add_value</a>(mp_int a, mp_small value, mp_int c);
237 - Sets `c` to the sum of `a` and `value`.
240 <a id="mp_int_sub"></a><pre>
241 mp_result <a href="imath.h#L172">mp_int_sub</a>(mp_int a, mp_int b, mp_int c);
243 - Sets `c` to the difference of `a` less `b`.
246 <a id="mp_int_sub_value"></a><pre>
247 mp_result <a href="imath.h#L175">mp_int_sub_value</a>(mp_int a, mp_small value, mp_int c);
249 - Sets `c` to the difference of `a` less `value`.
252 <a id="mp_int_mul"></a><pre>
253 mp_result <a href="imath.h#L178">mp_int_mul</a>(mp_int a, mp_int b, mp_int c);
255 - Sets `c` to the product of `a` and `b`.
258 <a id="mp_int_mul_value"></a><pre>
259 mp_result <a href="imath.h#L181">mp_int_mul_value</a>(mp_int a, mp_small value, mp_int c);
261 - Sets `c` to the product of `a` and `value`.
264 <a id="mp_int_mul_pow2"></a><pre>
265 mp_result <a href="imath.h#L184">mp_int_mul_pow2</a>(mp_int a, mp_small p2, mp_int c);
267 - Sets `c` to the product of `a` and `2^p2`. Requires `p2 >= 0`.
270 <a id="mp_int_sqr"></a><pre>
271 mp_result <a href="imath.h#L187">mp_int_sqr</a>(mp_int a, mp_int c);
273 - Sets `c` to the square of `a`.
276 <a id="mp_int_root"></a><pre>
277 mp_result <a href="imath.h#L306">mp_int_root</a>(mp_int a, mp_small b, mp_int c);
279 - Sets `c` to the greatest integer not less than the `b`th root of `a`,
280 using Newton's root-finding algorithm.
281 It returns `MP_UNDEF` if `a < 0` and `b` is even.
284 <a id="mp_int_sqrt"></a><pre>
285 static inline mp_result <a href="imath.h#L310">mp_int_sqrt</a>(mp_int a, mp_int c);
287 - Sets `c` to the greatest integer not less than the square root of `a`.
288 This is a special case of `mp_int_root()`.
291 <a id="mp_int_div"></a><pre>
292 mp_result <a href="imath.h#L195">mp_int_div</a>(mp_int a, mp_int b, mp_int q, mp_int r);
294 - Sets `q` and `r` to the quotent and remainder of `a / b`. Division by
295 powers of 2 is detected and handled efficiently. The remainder is pinned
298 Either of `q` or `r` may be NULL, but not both, and `q` and `r` may not
299 point to the same value.
302 <a id="mp_int_div_value"></a><pre>
303 mp_result <a href="imath.h#L200">mp_int_div_value</a>(mp_int a, mp_small value, mp_int q, mp_small *r);
305 - Sets `q` and `*r` to the quotent and remainder of `a / value`. Division by
306 powers of 2 is detected and handled efficiently. The remainder is pinned to
307 `0 <= *r < b`. Either of `q` or `r` may be NULL.
310 <a id="mp_int_div_pow2"></a><pre>
311 mp_result <a href="imath.h#L206">mp_int_div_pow2</a>(mp_int a, mp_small p2, mp_int q, mp_int r);
313 - Sets `q` and `r` to the quotient and remainder of `a / 2^p2`. This is a
314 special case for division by powers of two that is more efficient than
315 using ordinary division. Note that `mp_int_div()` will automatically handle
316 this case, this function is for cases where you have only the exponent.
319 <a id="mp_int_mod"></a><pre>
320 mp_result <a href="imath.h#L210">mp_int_mod</a>(mp_int a, mp_int m, mp_int c);
322 - Sets `c` to the remainder of `a / m`.
323 The remainder is pinned to `0 <= c < m`.
326 <a id="mp_int_mod_value"></a><pre>
327 static inline mp_result <a href="imath.h#L226">mp_int_mod_value</a>(mp_int a, mp_small value, mp_small* r);
329 - Sets `*r` to the remainder of `a / value`.
330 The remainder is pinned to `0 <= r < value`.
333 <a id="mp_int_expt"></a><pre>
334 mp_result <a href="imath.h#L214">mp_int_expt</a>(mp_int a, mp_small b, mp_int c);
336 - Sets `c` to the value of `a` raised to the `b` power.
337 It returns `MP_RANGE` if `b < 0`.
340 <a id="mp_int_expt_value"></a><pre>
341 mp_result <a href="imath.h#L218">mp_int_expt_value</a>(mp_small a, mp_small b, mp_int c);
343 - Sets `c` to the value of `a` raised to the `b` power.
344 It returns `MP_RANGE` if `b < 0`.
347 <a id="mp_int_expt_full"></a><pre>
348 mp_result <a href="imath.h#L222">mp_int_expt_full</a>(mp_int a, mp_int b, mp_int c);
350 - Sets `c` to the value of `a` raised to the `b` power.
351 It returns `MP_RANGE`) if `b < 0`.
355 ### Comparison Functions
357 Unless otherwise specified, comparison between values `x` and `y` returns a
358 **comparator**, an integer value < 0 if `x` is less than `y`, 0 if `x` is equal
359 to `y`, and > 0 if `x` is greater than `y`.
362 <a id="mp_int_compare"></a><pre>
363 int <a href="imath.h#L232">mp_int_compare</a>(mp_int a, mp_int b);
365 - Returns the comparator of `a` and `b`.
368 <a id="mp_int_compare_unsigned"></a><pre>
369 int <a href="imath.h#L236">mp_int_compare_unsigned</a>(mp_int a, mp_int b);
371 - Returns the comparator of the magnitudes of `a` and `b`, disregarding their
372 signs. Neither `a` nor `b` is modified by the comparison.
375 <a id="mp_int_compare_zero"></a><pre>
376 int <a href="imath.h#L239">mp_int_compare_zero</a>(mp_int z);
378 - Returns the comparator of `z` and zero.
381 <a id="mp_int_compare_value"></a><pre>
382 int <a href="imath.h#L242">mp_int_compare_value</a>(mp_int z, mp_small v);
384 - Returns the comparator of `z` and the signed value `v`.
387 <a id="mp_int_compare_uvalue"></a><pre>
388 int <a href="imath.h#L245">mp_int_compare_uvalue</a>(mp_int z, mp_usmall uv);
390 - Returns the comparator of `z` and the unsigned value `uv`.
393 <a id="mp_int_divisible_value"></a><pre>
394 bool <a href="imath.h#L248">mp_int_divisible_value</a>(mp_int a, mp_small v);
396 - Reports whether `a` is divisible by `v`.
399 <a id="mp_int_is_pow2"></a><pre>
400 int <a href="imath.h#L252">mp_int_is_pow2</a>(mp_int z);
402 - Returns `k >= 0` such that `z` is `2^k`, if such a `k` exists. If no such
403 `k` exists, the function returns -1.
407 ### Modular Operations
410 <a id="mp_int_exptmod"></a><pre>
411 mp_result <a href="imath.h#L256">mp_int_exptmod</a>(mp_int a, mp_int b, mp_int m, mp_int c);
413 - Sets `c` to the value of `a` raised to the `b` power, reduced modulo `m`.
414 It returns `MP_RANGE` if `b < 0` or `MP_UNDEF` if `m == 0`.
417 <a id="mp_int_exptmod_evalue"></a><pre>
418 mp_result <a href="imath.h#L260">mp_int_exptmod_evalue</a>(mp_int a, mp_small value, mp_int m, mp_int c);
420 - Sets `c` to the value of `a` raised to the `value` power, modulo `m`.
421 It returns `MP_RANGE` if `value < 0` or `MP_UNDEF` if `m == 0`.
424 <a id="mp_int_exptmod_bvalue"></a><pre>
425 mp_result <a href="imath.h#L264">mp_int_exptmod_bvalue</a>(mp_small value, mp_int b, mp_int m, mp_int c);
427 - Sets `c` to the value of `value` raised to the `b` power, modulo `m`.
428 It returns `MP_RANGE` if `b < 0` or `MP_UNDEF` if `m == 0`.
431 <a id="mp_int_exptmod_known"></a><pre>
432 mp_result <a href="imath.h#L271">mp_int_exptmod_known</a>(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c);
434 - Sets `c` to the value of `a` raised to the `b` power, reduced modulo `m`,
435 given a precomputed reduction constant `mu` defined for Barrett's modular
438 It returns `MP_RANGE` if `b < 0` or `MP_UNDEF` if `m == 0`.
441 <a id="mp_int_redux_const"></a><pre>
442 mp_result <a href="imath.h#L275">mp_int_redux_const</a>(mp_int m, mp_int c);
444 - Sets `c` to the reduction constant for Barrett reduction by modulus `m`.
445 Requires that `c` and `m` point to distinct locations.
448 <a id="mp_int_invmod"></a><pre>
449 mp_result <a href="imath.h#L282">mp_int_invmod</a>(mp_int a, mp_int m, mp_int c);
451 - Sets `c` to the multiplicative inverse of `a` modulo `m`, if it exists.
452 The least non-negative representative of the congruence class is computed.
454 It returns `MP_UNDEF` if the inverse does not exist, or `MP_RANGE` if `a ==
458 <a id="mp_int_gcd"></a><pre>
459 mp_result <a href="imath.h#L288">mp_int_gcd</a>(mp_int a, mp_int b, mp_int c);
461 - Sets `c` to the greatest common divisor of `a` and `b`.
463 It returns `MP_UNDEF` if the GCD is undefined, such as for example if `a`
464 and `b` are both zero.
467 <a id="mp_int_egcd"></a><pre>
468 mp_result <a href="imath.h#L295">mp_int_egcd</a>(mp_int a, mp_int b, mp_int c, mp_int x, mp_int y);
470 - Sets `c` to the greatest common divisor of `a` and `b`, and sets `x` and
471 `y` to values satisfying Bezout's identity `gcd(a, b) = ax + by`.
473 It returns `MP_UNDEF` if the GCD is undefined, such as for example if `a`
474 and `b` are both zero.
477 <a id="mp_int_lcm"></a><pre>
478 mp_result <a href="imath.h#L301">mp_int_lcm</a>(mp_int a, mp_int b, mp_int c);
480 - Sets `c` to the least common multiple of `a` and `b`.
482 It returns `MP_UNDEF` if the LCM is undefined, such as for example if `a`
483 and `b` are both zero.
487 ### Conversion of Values
490 <a id="mp_int_to_int"></a><pre>
491 mp_result <a href="imath.h#L315">mp_int_to_int</a>(mp_int z, mp_small *out);
493 - Returns `MP_OK` if `z` is representable as `mp_small`, else `MP_RANGE`.
494 If `out` is not NULL, `*out` is set to the value of `z` when `MP_OK`.
497 <a id="mp_int_to_uint"></a><pre>
498 mp_result <a href="imath.h#L319">mp_int_to_uint</a>(mp_int z, mp_usmall *out);
500 - Returns `MP_OK` if `z` is representable as `mp_usmall`, or `MP_RANGE`.
501 If `out` is not NULL, `*out` is set to the value of `z` when `MP_OK`.
504 <a id="mp_int_to_string"></a><pre>
505 mp_result <a href="imath.h#L327">mp_int_to_string</a>(mp_int z, mp_size radix, char *str, int limit);
507 - Converts `z` to a zero-terminated string of characters in the specified
508 `radix`, writing at most `limit` characters to `str` including the
509 terminating NUL value. A leading `-` is used to indicate a negative value.
511 Returns `MP_TRUNC` if `limit` was to small to write all of `z`.
512 Requires `MP_MIN_RADIX <= radix <= MP_MAX_RADIX`.
515 <a id="mp_int_string_len"></a><pre>
516 mp_result <a href="imath.h#L332">mp_int_string_len</a>(mp_int z, mp_size radix);
518 - Reports the minimum number of characters required to represent `z` as a
519 zero-terminated string in the given `radix`.
520 Requires `MP_MIN_RADIX <= radix <= MP_MAX_RADIX`.
523 <a id="mp_int_read_string"></a><pre>
524 mp_result <a href="imath.h#L347">mp_int_read_string</a>(mp_int z, mp_size radix, const char *str);
526 - Reads a string of ASCII digits in the specified `radix` from the zero
527 terminated `str` provided into `z`. For values of `radix > 10`, the letters
528 `A`..`Z` or `a`..`z` are accepted. Letters are interpreted without respect
531 Leading whitespace is ignored, and a leading `+` or `-` is interpreted as a
532 sign flag. Processing stops when a NUL or any other character out of range
533 for a digit in the given radix is encountered.
535 If the whole string was consumed, `MP_OK` is returned; otherwise
536 `MP_TRUNC`. is returned.
538 Requires `MP_MIN_RADIX <= radix <= MP_MAX_RADIX`.
541 <a id="mp_int_read_cstring"></a><pre>
542 mp_result <a href="imath.h#L365">mp_int_read_cstring</a>(mp_int z, mp_size radix, const char *str, char **end);
544 - Reads a string of ASCII digits in the specified `radix` from the zero
545 terminated `str` provided into `z`. For values of `radix > 10`, the letters
546 `A`..`Z` or `a`..`z` are accepted. Letters are interpreted without respect
549 Leading whitespace is ignored, and a leading `+` or `-` is interpreted as a
550 sign flag. Processing stops when a NUL or any other character out of range
551 for a digit in the given radix is encountered.
553 If the whole string was consumed, `MP_OK` is returned; otherwise
554 `MP_TRUNC`. is returned. If `end` is not NULL, `*end` is set to point to
555 the first unconsumed byte of the input string (the NUL byte if the whole
556 string was consumed). This emulates the behavior of the standard C
559 Requires `MP_MIN_RADIX <= radix <= MP_MAX_RADIX`.
562 <a id="mp_int_count_bits"></a><pre>
563 mp_result <a href="imath.h#L368">mp_int_count_bits</a>(mp_int z);
565 - Returns the number of significant bits in `z`.
568 <a id="mp_int_to_binary"></a><pre>
569 mp_result <a href="imath.h#L383">mp_int_to_binary</a>(mp_int z, unsigned char *buf, int limit);
571 - Converts `z` to 2's complement binary, writing at most `limit` bytes into
572 the given `buf`. Returns `MP_TRUNC` if the buffer limit was too small to
573 contain the whole value. If this occurs, the contents of buf will be
574 effectively garbage, as the function uses the buffer as scratch space.
576 The binary representation of `z` is in base-256 with digits ordered from
577 most significant to least significant (network byte ordering). The
578 high-order bit of the first byte is set for negative values, clear for
581 As a result, non-negative values will be padded with a leading zero byte if
582 the high-order byte of the base-256 magnitude is set. This extra byte is
583 accounted for by the `mp_int_binary_len()` function.
586 <a id="mp_int_read_binary"></a><pre>
587 mp_result <a href="imath.h#L388">mp_int_read_binary</a>(mp_int z, unsigned char *buf, int len);
589 - Reads a 2's complement binary value from `buf` into `z`, where `len` is the
590 length of the buffer. The contents of `buf` may be overwritten during
591 processing, although they will be restored when the function returns.
594 <a id="mp_int_binary_len"></a><pre>
595 mp_result <a href="imath.h#L391">mp_int_binary_len</a>(mp_int z);
597 - Returns the number of bytes to represent `z` in 2's complement binary.
600 <a id="mp_int_to_unsigned"></a><pre>
601 mp_result <a href="imath.h#L402">mp_int_to_unsigned</a>(mp_int z, unsigned char *buf, int limit);
603 - Converts the magnitude of `z` to unsigned binary, writing at most `limit`
604 bytes into the given `buf`. The sign of `z` is ignored, but `z` is not
605 modified. Returns `MP_TRUNC` if the buffer limit was too small to contain
606 the whole value. If this occurs, the contents of `buf` will be effectively
607 garbage, as the function uses the buffer as scratch space during
610 The binary representation of `z` is in base-256 with digits ordered from
611 most significant to least significant (network byte ordering).
614 <a id="mp_int_read_unsigned"></a><pre>
615 mp_result <a href="imath.h#L407">mp_int_read_unsigned</a>(mp_int z, unsigned char *buf, int len);
617 - Reads an unsigned binary value from `buf` into `z`, where `len` is the
618 length of the buffer. The contents of `buf` are not modified during
622 <a id="mp_int_unsigned_len"></a><pre>
623 mp_result <a href="imath.h#L411">mp_int_unsigned_len</a>(mp_int z);
625 - Returns the number of bytes required to represent `z` as an unsigned binary
632 Ordinarily, integer multiplication and squaring are done using the simple
633 quadratic "schoolbook" algorithm. However, for sufficiently large values,
634 there is a more efficient algorithm usually attributed to Karatsuba and Ofman
635 that is usually faster. See Knuth Vol. 2 for more details about how this
638 The breakpoint between the "normal" and the recursive algorithm is controlled
639 by a static digit threshold defined in `imath.c`. Values with fewer significant
640 digits use the standard algorithm. This value can be modified by calling
641 `mp_int_multiply_threshold(n)`. The `imtimer` program and the
642 `findthreshold.py` script (Python) can help you find a suitable value for for
643 your particular platform.
646 <a id="mp_error_string"></a><pre>
647 const char *<a href="imath.h#L416">mp_error_string</a>(mp_result res);
649 - Returns a pointer to a brief, human-readable, zero-terminated string
650 describing `res`. The returned string is statically allocated and must not
651 be freed by the caller.
655 ## Rational Arithmetic
658 <a id="mp_rat_init"></a><pre>
659 mp_result <a href="imrat.h#L59">mp_rat_init</a>(mp_rat r);
661 - Initializes `r` with 1-digit precision and sets it to zero. This function
662 cannot fail unless `r` is NULL.
665 <a id="mp_rat_alloc"></a><pre>
666 mp_rat <a href="imrat.h#L63">mp_rat_alloc</a>(void);
668 - Allocates a fresh zero-valued `mpq_t` on the heap, returning NULL in case
669 of error. The only possible error is out-of-memory.
672 <a id="mp_rat_reduce"></a><pre>
673 mp_result <a href="imrat.h#L69">mp_rat_reduce</a>(mp_rat r);
675 - Reduces `r` in-place to lowest terms and canonical form.
677 Zero is represented as 0/1, one as 1/1, and signs are adjusted so that the
678 sign of the value is carried by the numerator.
681 <a id="mp_rat_init_size"></a><pre>
682 mp_result <a href="imrat.h#L76">mp_rat_init_size</a>(mp_rat r, mp_size n_prec, mp_size d_prec);
684 - Initializes `r` with at least `n_prec` digits of storage for the numerator
685 and `d_prec` digits of storage for the denominator, and value zero.
687 If either precision is zero, the default precision is used, rounded up to
688 the nearest word size.
691 <a id="mp_rat_init_copy"></a><pre>
692 mp_result <a href="imrat.h#L80">mp_rat_init_copy</a>(mp_rat r, mp_rat old);
694 - Initializes `r` to be a copy of an already-initialized value in `old`. The
695 new copy does not share storage with the original.
698 <a id="mp_rat_set_value"></a><pre>
699 mp_result <a href="imrat.h#L84">mp_rat_set_value</a>(mp_rat r, mp_small numer, mp_small denom);
701 - Sets the value of `r` to the ratio of signed `numer` to signed `denom`. It
702 returns `MP_UNDEF` if `denom` is zero.
705 <a id="mp_rat_set_uvalue"></a><pre>
706 mp_result <a href="imrat.h#L88">mp_rat_set_uvalue</a>(mp_rat r, mp_usmall numer, mp_usmall denom);
708 - Sets the value of `r` to the ratio of unsigned `numer` to unsigned
709 `denom`. It returns `MP_UNDEF` if `denom` is zero.
712 <a id="mp_rat_clear"></a><pre>
713 void <a href="imrat.h#L91">mp_rat_clear</a>(mp_rat r);
715 - Releases the storage used by `r`.
718 <a id="mp_rat_free"></a><pre>
719 void <a href="imrat.h#L95">mp_rat_free</a>(mp_rat r);
721 - Releases the storage used by `r` and also `r` itself.
722 This should only be used for `r` allocated by `mp_rat_alloc()`.
725 <a id="mp_rat_numer"></a><pre>
726 mp_result <a href="imrat.h#L98">mp_rat_numer</a>(mp_rat r, mp_int z);
728 - Sets `z` to a copy of the numerator of `r`.
731 <a id="mp_rat_numer_ref"></a><pre>
732 mp_int <a href="imrat.h#L101">mp_rat_numer_ref</a>(mp_rat r);
734 - Returns a pointer to the numerator of `r`.
737 <a id="mp_rat_denom"></a><pre>
738 mp_result <a href="imrat.h#L104">mp_rat_denom</a>(mp_rat r, mp_int z);
740 - Sets `z` to a copy of the denominator of `r`.
743 <a id="mp_rat_denom_ref"></a><pre>
744 mp_int <a href="imrat.h#L107">mp_rat_denom_ref</a>(mp_rat r);
746 - Returns a pointer to the denominator of `r`.
749 <a id="mp_rat_sign"></a><pre>
750 mp_sign <a href="imrat.h#L110">mp_rat_sign</a>(mp_rat r);
752 - Reports the sign of `r`.
755 <a id="mp_rat_copy"></a><pre>
756 mp_result <a href="imrat.h#L115">mp_rat_copy</a>(mp_rat a, mp_rat c);
758 - Sets `c` to a copy of the value of `a`. No new memory is allocated unless a
759 term of `a` has more significant digits than the corresponding term of `c`
763 <a id="mp_rat_zero"></a><pre>
764 void <a href="imrat.h#L118">mp_rat_zero</a>(mp_rat r);
766 - Sets `r` to zero. The allocated storage of `r` is not changed.
769 <a id="mp_rat_abs"></a><pre>
770 mp_result <a href="imrat.h#L121">mp_rat_abs</a>(mp_rat a, mp_rat c);
772 - Sets `c` to the absolute value of `a`.
775 <a id="mp_rat_neg"></a><pre>
776 mp_result <a href="imrat.h#L124">mp_rat_neg</a>(mp_rat a, mp_rat c);
778 - Sets `c` to the absolute value of `a`.
781 <a id="mp_rat_recip"></a><pre>
782 mp_result <a href="imrat.h#L128">mp_rat_recip</a>(mp_rat a, mp_rat c);
784 - Sets `c` to the reciprocal of `a` if the reciprocal is defined.
785 It returns `MP_UNDEF` if `a` is zero.
788 <a id="mp_rat_add"></a><pre>
789 mp_result <a href="imrat.h#L131">mp_rat_add</a>(mp_rat a, mp_rat b, mp_rat c);
791 - Sets `c` to the sum of `a` and `b`.
794 <a id="mp_rat_sub"></a><pre>
795 mp_result <a href="imrat.h#L134">mp_rat_sub</a>(mp_rat a, mp_rat b, mp_rat c);
797 - Sets `c` to the difference of `a` less `b`.
800 <a id="mp_rat_mul"></a><pre>
801 mp_result <a href="imrat.h#L137">mp_rat_mul</a>(mp_rat a, mp_rat b, mp_rat c);
803 - Sets `c` to the product of `a` and `b`.
806 <a id="mp_rat_div"></a><pre>
807 mp_result <a href="imrat.h#L141">mp_rat_div</a>(mp_rat a, mp_rat b, mp_rat c);
809 - Sets `c` to the ratio `a / b` if that ratio is defined.
810 It returns `MP_UNDEF` if `b` is zero.
813 <a id="mp_rat_add_int"></a><pre>
814 mp_result <a href="imrat.h#L144">mp_rat_add_int</a>(mp_rat a, mp_int b, mp_rat c);
816 - Sets `c` to the sum of `a` and integer `b`.
819 <a id="mp_rat_sub_int"></a><pre>
820 mp_result <a href="imrat.h#L147">mp_rat_sub_int</a>(mp_rat a, mp_int b, mp_rat c);
822 - Sets `c` to the difference of `a` less integer `b`.
825 <a id="mp_rat_mul_int"></a><pre>
826 mp_result <a href="imrat.h#L150">mp_rat_mul_int</a>(mp_rat a, mp_int b, mp_rat c);
828 - Sets `c` to the product of `a` and integer `b`.
831 <a id="mp_rat_div_int"></a><pre>
832 mp_result <a href="imrat.h#L154">mp_rat_div_int</a>(mp_rat a, mp_int b, mp_rat c);
834 - Sets `c` to the ratio `a / b` if that ratio is defined.
835 It returns `MP_UNDEF` if `b` is zero.
838 <a id="mp_rat_expt"></a><pre>
839 mp_result <a href="imrat.h#L158">mp_rat_expt</a>(mp_rat a, mp_small b, mp_rat c);
841 - Sets `c` to the value of `a` raised to the `b` power.
842 It returns `MP_RANGE` if `b < 0`.
845 <a id="mp_rat_compare"></a><pre>
846 int <a href="imrat.h#L161">mp_rat_compare</a>(mp_rat a, mp_rat b);
848 - Returns the comparator of `a` and `b`.
851 <a id="mp_rat_compare_unsigned"></a><pre>
852 int <a href="imrat.h#L165">mp_rat_compare_unsigned</a>(mp_rat a, mp_rat b);
854 - Returns the comparator of the magnitudes of `a` and `b`, disregarding their
855 signs. Neither `a` nor `b` is modified by the comparison.
858 <a id="mp_rat_compare_zero"></a><pre>
859 int <a href="imrat.h#L168">mp_rat_compare_zero</a>(mp_rat r);
861 - Returns the comparator of `r` and zero.
864 <a id="mp_rat_compare_value"></a><pre>
865 int <a href="imrat.h#L172">mp_rat_compare_value</a>(mp_rat r, mp_small n, mp_small d);
867 - Returns the comparator of `r` and the signed ratio `n / d`.
868 It returns `MP_UNDEF` if `d` is zero.
871 <a id="mp_rat_is_integer"></a><pre>
872 bool <a href="imrat.h#L175">mp_rat_is_integer</a>(mp_rat r);
874 - Reports whether `r` is an integer, having canonical denominator 1.
877 <a id="mp_rat_to_ints"></a><pre>
878 mp_result <a href="imrat.h#L180">mp_rat_to_ints</a>(mp_rat r, mp_small *num, mp_small *den);
880 - Reports whether the numerator and denominator of `r` can be represented as
881 small signed integers, and if so stores the corresponding values to `num`
882 and `den`. It returns `MP_RANGE` if either cannot be so represented.
885 <a id="mp_rat_to_string"></a><pre>
886 mp_result <a href="imrat.h#L186">mp_rat_to_string</a>(mp_rat r, mp_size radix, char *str, int limit);
888 - Converts `r` to a zero-terminated string of the format `"n/d"` with `n` and
889 `d` in the specified radix and writing no more than `limit` bytes to the
890 given output buffer `str`. The output of the numerator includes a sign flag
891 if `r` is negative. Requires `MP_MIN_RADIX <= radix <= MP_MAX_RADIX`.
894 <a id="mp_rat_to_decimal"></a><pre>
895 mp_result <a href="imrat.h#L215">mp_rat_to_decimal</a>(mp_rat r, mp_size radix, mp_size prec, mp_round_mode round, char *str, int limit);
897 - Converts the value of `r` to a string in decimal-point notation with the
898 specified radix, writing no more than `limit` bytes of data to the given
899 output buffer. It generates `prec` digits of precision, and requires
900 `MP_MIN_RADIX <= radix <= MP_MAX_RADIX`.
902 Ratios usually must be rounded when they are being converted for output as
903 a decimal value. There are four rounding modes currently supported:
907 Truncates the value toward zero.
908 Example: 12.009 to 2dp becomes 12.00
913 Rounds the value away from zero:
914 Example: 12.001 to 2dp becomes 12.01, but
915 12.000 to 2dp remains 12.00
920 Rounds the value to nearest digit, half goes toward zero.
921 Example: 12.005 to 2dp becomes 12.00, but
922 12.006 to 2dp becomes 12.01
927 Rounds the value to nearest digit, half rounds upward.
928 Example: 12.005 to 2dp becomes 12.01, but
929 12.004 to 2dp becomes 12.00
933 <a id="mp_rat_string_len"></a><pre>
934 mp_result <a href="imrat.h#L221">mp_rat_string_len</a>(mp_rat r, mp_size radix);
936 - Reports the minimum number of characters required to represent `r` as a
937 zero-terminated string in the given `radix`.
938 Requires `MP_MIN_RADIX <= radix <= MP_MAX_RADIX`.
941 <a id="mp_rat_decimal_len"></a><pre>
942 mp_result <a href="imrat.h#L226">mp_rat_decimal_len</a>(mp_rat r, mp_size radix, mp_size prec);
944 - Reports the length in bytes of the buffer needed to convert `r` using the
945 `mp_rat_to_decimal()` function with the specified `radix` and `prec`. The
946 buffer size estimate may slightly exceed the actual required capacity.
949 <a id="mp_rat_read_string"></a><pre>
950 mp_result <a href="imrat.h#L231">mp_rat_read_string</a>(mp_rat r, mp_size radix, const char *str);
952 - Sets `r` to the value represented by a zero-terminated string `str` in the
953 format `"n/d"` including a sign flag. It returns `MP_UNDEF` if the encoded
954 denominator has value zero.
957 <a id="mp_rat_read_cstring"></a><pre>
958 mp_result <a href="imrat.h#L238">mp_rat_read_cstring</a>(mp_rat r, mp_size radix, const char *str, char **end);
960 - Sets `r` to the value represented by a zero-terminated string `str` in the
961 format `"n/d"` including a sign flag. It returns `MP_UNDEF` if the encoded
962 denominator has value zero. If `end` is not NULL then `*end` is set to
963 point to the first unconsumed character in the string, after parsing.
966 <a id="mp_rat_read_ustring"></a><pre>
967 mp_result <a href="imrat.h#L252">mp_rat_read_ustring</a>(mp_rat r, mp_size radix, const char *str, char **end);
969 - Sets `r` to the value represented by a zero-terminated string `str` having
970 one of the following formats, each with an optional leading sign flag:
973 n : integer format, e.g. "123"
974 n/d : ratio format, e.g., "-12/5"
975 z.ffff : decimal format, e.g., "1.627"
978 It returns `MP_UNDEF` if the effective denominator is zero. If `end` is not
979 NULL then `*end` is set to point to the first unconsumed character in the
980 string, after parsing.
983 <a id="mp_rat_read_decimal"></a><pre>
984 mp_result <a href="imrat.h#L258">mp_rat_read_decimal</a>(mp_rat r, mp_size radix, const char *str);
986 - Sets `r` to the value represented by a zero-terminated string `str` in the
987 format `"z.ffff"` including a sign flag. It returns `MP_UNDEF` if the
988 effective denominator.
991 <a id="mp_rat_read_cdecimal"></a><pre>
992 mp_result <a href="imrat.h#L264">mp_rat_read_cdecimal</a>(mp_rat r, mp_size radix, const char *str, char **end);
994 - Sets `r` to the value represented by a zero-terminated string `str` in the
995 format `"z.ffff"` including a sign flag. It returns `MP_UNDEF` if the
996 effective denominator. If `end` is not NULL then `*end` is set to point to
997 the first unconsumed character in the string, after parsing.
1001 ## Representation Details
1003 > NOTE: You do not need to read this section to use IMath. This is provided
1004 > for the benefit of developers wishing to extend or modify the internals of
1007 IMath uses a signed magnitude representation for arbitrary precision integers.
1008 The magnitude is represented as an array of radix-R digits in increasing order
1009 of significance; the value of R is chosen to be half the size of the largest
1010 available unsigned integer type, so typically 16 or 32 bits. Digits are
1011 represented as mp_digit, which must be an unsigned integral type.
1013 Digit arrays are allocated using `malloc(3)` and `realloc(3)`. Because this
1014 can be an expensive operation, the library takes pains to avoid allocation as
1015 much as possible. For this reason, the `mpz_t` structure distinguishes between
1016 how many digits are allocated and how many digits are actually consumed by the
1017 representation. The fields of an `mpz_t` are:
1019 mp_digit single; /* single-digit value (see note) */
1020 mp_digit *digits; /* array of digits */
1021 mp_size alloc; /* how many digits are allocated */
1022 mp_size used; /* how many digits are in use */
1023 mp_sign sign; /* the sign of the value */
1025 The elements of `digits` at indices less than `used` are the significant
1026 figures of the value; the elements at indices greater than or equal to `used`
1027 are undefined (and may contain garbage). At all times, `used` must be at least
1028 1 and at most `alloc`.
1030 To avoid interaction with the memory allocator, single-digit values are stored
1031 directly in the `mpz_t` structure, in the `single` field. The semantics of
1032 access are the same as the more general case.
1034 The number of digits allocated for an `mpz_t` is referred to in the library
1035 documentation as its "precision". Operations that affect an `mpz_t` cause
1036 precision to increase as needed. In any case, all allocations are measured in
1037 digits, and rounded up to the nearest `mp_word` boundary. There is a default
1038 minimum precision stored as a static constant default_precision (`imath.c`).
1039 This value can be set using `mp_int_default_precision(n)`.
1041 Note that the allocated size of an `mpz_t` can only grow; the library never
1042 reallocates in order to decrease the size. A simple way to do so explicitly is
1043 to use `mp_int_init_copy()`, as in:
1049 mp_int_init_copy(&new, &big);
1050 mp_int_swap(&new, &big);
1054 The value of `sign` is 0 for positive values and zero, 1 for negative values.
1055 Constants `MP_ZPOS` and `MP_NEG` are defined for these; no other sign values
1058 If you are adding to this library, you should be careful to preserve the
1059 convention that inputs and outputs can overlap, as described above. So, for
1060 example, `mp_int_add(a, a, a)` is legal. Often, this means you must maintain
1061 one or more temporary mpz_t structures for intermediate values. The private
1062 macros `DECLARE_TEMP(N)`, `CLEANUP_TEMP()`, and `TEMP(K)` can be used to
1063 maintain a conventional structure like this:
1067 /* Declare how many temp values you need.
1068 Use TEMP(i) to access the ith value (0-indexed). */
1072 /* Perform actions that must return MP_OK or fail. */
1073 REQUIRE(mp_int_copy(x, TEMP(1)));
1075 REQUIRE(mp_int_expt(TEMP(1), TEMP(2), TEMP(3)));
1078 /* You can also use REQUIRE directly for more complex cases. */
1079 if (some_difficult_question(TEMP(3)) != answer(x)) {
1080 REQUIRE(MP_RANGE); /* falls through to cleanup (below) */
1083 /* Ensure temporary values are cleaned up at exit.
1085 If control reaches here via a REQUIRE failure, the code below
1086 the cleanup will not be executed.
1093 Under the covers, these macros are just maintaining an array of `mpz_t` values,
1094 and a jump label to handle cleanup. You may only have one `DECLARE_TEMP` and
1095 its corresponding `CLEANUP_TEMP` per function body.
1097 "Small" integer values are represented by the types `mp_small` and `mp_usmall`,
1098 which are mapped to appropriately-sized types on the host system. The default
1099 for `mp_small` is "long" and the default for `mp_usmall` is "unsigned long".
1100 You may change these, provided you insure that `mp_small` is signed and
1101 `mp_usmall` is unsigned. You will also need to adjust the size macros:
1103 MP_SMALL_MIN, MP_SMALL_MAX
1104 MP_USMALL_MIN, MP_USMALL_MAX
1106 ... which are defined in `<imath.h>`, if you change these.
1108 Rational numbers are represented using a pair of arbitrary precision integers,
1109 with the convention that the sign of the numerator is the sign of the rational
1110 value, and that the result of any rational operation is always represented in
1111 lowest terms. The canonical representation for rational zero is 0/1. See
1114 ## Testing and Reporting of Bugs
1116 Test vectors are included in the `tests/` subdirectory of the imath
1117 distribution. When you run `make test`, it builds the `imtest` program and
1118 runs all available test vectors. If any tests fail, you will get a line like
1123 Here, _x_ is the line number of the test which failed, _y_ is index of the test
1124 within the file, and _v_ is the value(s) actually computed. The name of the
1125 file is printed at the beginning of each test, so you can find out what test
1126 vector failed by executing the following (with x, y, and v replaced by the
1127 above values, and where "foo.t" is the name of the test file that was being
1128 processed at the time):
1130 % tail +x tests/foo.t | head -1
1132 None of the tests should fail (but see [Note 2](#note2)); if any do, it
1133 probably indicates a bug in the library (or at the very least, some assumption
1134 I made which I shouldn't have). Please [file an
1135 issue](https://github.com/creachadair/imath/issues/new), including the `FAILED`
1136 test line(s), as well as the output of the above `tail` command (so I know what
1137 inputs caused the failure).
1139 If you build with the preprocessor symbol `DEBUG` defined as a positive
1140 integer, the digit allocators (`s_alloc`, `s_realloc`) fill all new buffers
1141 with the value `0xdeadbeefabad1dea`, or as much of it as will fit in a digit,
1142 so that you can more easily catch uninitialized reads in the debugger.
1146 1. <a name="note1"></a>You can generally use the same variables for both input
1147 and output. One exception is that you may not use the same variable for
1148 both the quotient and the remainder of `mp_int_div()`.
1150 2. <a name="note2"></a>Many of the tests for this library were written under
1151 the assumption that the `mp_small` type is 32 bits or more. If you compile
1152 with a smaller type, you may see `MP_RANGE` errors in some of the tests that
1153 otherwise pass (due to conversion failures). Also, the pi generator (pi.c)
1154 will not work correctly if `mp_small` is too short, as its algorithm for arc
1155 tangent is fairly simple-minded.
1159 The IMath library was written by Michael J. Fromberger.
1161 If you discover any bugs or testing failures, please [open an
1162 issue](https://github.com/creachadair/imath/issues/new). Please be sure to
1163 include a complete description of what went wrong, and if possible, a test
1164 vector for `imtest` and/or a minimal test program that will demonstrate the bug
1165 on your system. Please also let me know what hardware, operating system, and
1166 compiler you're using.
1170 The algorithms used in this library came from Vol. 2 of Donald Knuth's "The Art
1171 of Computer Programming" (Seminumerical Algorithms). Thanks to Nelson Bolyard,
1172 Bryan Olson, Tom St. Denis, Tushar Udeshi, and Eric Silva for excellent
1173 feedback on earlier versions of this code. Special thanks to Jonathan Shapiro
1174 for some very helpful design advice, as well as feedback and some clever ideas
1175 for improving performance in some common use cases.
1177 ## License and Disclaimers
1179 IMath is Copyright 2002-2009 Michael J. Fromberger
1180 You may use it subject to the following Licensing Terms:
1182 Permission is hereby granted, free of charge, to any person obtaining a copy of
1183 this software and associated documentation files (the "Software"), to deal in
1184 the Software without restriction, including without limitation the rights to
1185 use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
1186 of the Software, and to permit persons to whom the Software is furnished to do
1187 so, subject to the following conditions:
1189 The above copyright notice and this permission notice shall be included in all
1190 copies or substantial portions of the Software.
1192 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
1193 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1194 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
1195 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
1196 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
1197 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE