1 # User Documentation for the IMath Library
3 Author: [M. J. Fromberger](https://github.com/creachadair)
7 1. Edit Makefile to select compiler and options. The default is to use gcc.
8 You may want to change CC to `clang` instead of `gcc` (and on macOS that
9 what you will get anyway), but you should be able to use the default GCC
12 By default, the Makefile assumes you can use 64-bit integer types, even
13 though they were not standard in ANSI C90. If you cannot, add
14 `-DUSE_32BIT_WORDS` to the compiler options.
16 2. Type `make` or `make test` to build the test driver and run the unit tests.
17 None of these should fail. If they do, see below for how you can report
20 To build with debugging enabled (and optimization disabled), run `make
21 DEBUG=Y`. This sets the preprocessor macro `DEBUG` to 1, and several other
22 things (see Makefile for details).
24 To use the library in your code, include "imath.h" wherever you intend to use
25 the library's routines. The integer library is just a single source file, so
26 you can compile it into your project in whatever way makes sense. If you wish
27 to use rational arithmetic, you will also need to include "imrat.h".
31 The basic types defined by the imath library are `mpz_t`, an arbitrary
32 precision signed integer, and `mpq_t`, an arbitrary precision signed rational
33 number. The type `mp_int` is a pointer to an `mpz_t`, and `mp_rat` is a
34 pointer to an `mpq_t`.
36 Most of the functions in the imath library return a value of type `mp_result`.
37 This is a signed integer type which can be used to convey status information
38 and also return small values. Any negative value is considered to be a status
39 message. The following constants are defined for processing these:
41 | Status | Description |
42 | ----------- | -------------------------------------------- |
43 | `MP_OK` | operation successful, all is well (= 0) |
44 | `MP_FALSE` | boolean false (= `MP_OK`) |
45 | `MP_TRUE` | boolean true |
46 | `MP_MEMORY` | out of memory |
47 | `MP_RANGE` | parameter out of range |
48 | `MP_UNDEF` | result is undefined (e.g., division by zero) |
49 | `MP_TRUNC` | output value was truncated |
50 | `MP_BADARG` | an invalid parameter was passed |
52 If you obtain a zero or negative value of an `mp_result`, you can use the
53 `mp_error_string()` routine to obtain a pointer to a brief human-readable
54 string describing the error. These strings are statically allocated, so they
55 need not be freed by the caller; the same strings are re-used from call to
58 Unless otherwise noted, it is legal to use the same parameter for both inputs
59 and output with most of the functions in this library. For example, you can
60 add a number to itself and replace the original by writing:
62 mp_int_add(a, a, a); /* a = a + a */
64 Any cases in which this is not legal will be noted in the function summaries
65 below (if you discover that this is not so, please report it as a bug; I will
66 fix either the function or the documentation :)
70 Each of the API functions is documented here. The general format of the
75 > return_type function_name(parameters ...)
77 > - English description.
79 Unless otherwise noted, any API function that returns `mp_result` may be
80 expected to return `MP_OK`, `MP_BADARG`, or `MP_MEMORY`. Other return values
81 should be documented in the description. Please let me know if you discover
84 The following macros are defined in "imath.h", to define the sizes of the
85 various data types used in the library:
87 | Constant | Description
88 | --------------- | ----------------------------------------
89 | `MP_DIGIT_BIT` | the number of bits in a single `mpz_t` digit.
90 | `MP_WORD_BIT` | the number of bits in a `mpz_t` word.
91 | `MP_SMALL_MIN` | the minimum value representable by an `mp_small`.
92 | `MP_SMALL_MAX` | the maximum value representable by an `mp_small`.
93 | `MP_USMALL_MAX` | the maximum value representable by an `mp_usmall`.
94 | `MP_MIN_RADIX` | the minimum radix accepted for base conversion.
95 | `MP_MAX_RADIX` | the maximum radix accepted for base conversion.
99 An `mp_int` must be initialized before use. By default, an `mp_int` is
100 initialized with a certain minimum amount of storage for digits, and the
101 storage is expanded automatically as needed. To initialize an `mp_int`, use
102 the following functions:
105 mp_int_init mp_int_alloc mp_int_init_size
112 When you are finished with an `mp_int`, you must free the memory it uses:
114 {{insert "imath.h" mp_int_clear mp_int_free}}
118 To set an `mp_int` which has already been initialized to a small integer value,
121 {{insert "imath.h" mp_int_set_value mp_int_set_uvalue}}
123 To copy one initialized `mp_int` to another, use:
125 {{insert "imath.h" mp_int_copy}}
127 ### Arithmetic Functions
130 mp_int_is_odd mp_int_is_even
134 mp_int_add mp_int_add_value
135 mp_int_sub mp_int_sub_value
136 mp_int_mul mp_int_mul_value mp_int_mul_pow2
138 mp_int_root mp_int_sqrt
139 mp_int_div mp_int_div_value mp_int_div_pow2
140 mp_int_mod mp_int_mod_value
141 mp_int_expt mp_int_expt_value mp_int_expt_full
144 ### Comparison Functions
146 Unless otherwise specified, comparison between values `x` and `y` returns a
147 **comparator**, an integer value < 0 if `x` is less than `y`, 0 if `x` is equal
148 to `y`, and > 0 if `x` is greater than `y`.
151 mp_int_compare mp_int_compare_unsigned mp_int_compare_zero
152 mp_int_compare_value mp_int_compare_uvalue
153 mp_int_divisible_value mp_int_is_pow2
156 ### Modular Operations
159 mp_int_exptmod mp_int_exptmod_evalue mp_int_exptmod_bvalue
160 mp_int_exptmod_known mp_int_redux_const
162 mp_int_gcd mp_int_egcd mp_int_lcm
165 ### Conversion of Values
168 mp_int_to_int mp_int_to_uint
169 mp_int_to_string mp_int_string_len
170 mp_int_read_string mp_int_read_cstring
172 mp_int_to_binary mp_int_read_binary mp_int_binary_len
173 mp_int_to_unsigned mp_int_read_unsigned mp_int_unsigned_len
178 Ordinarily, integer multiplication and squaring are done using the simple
179 quadratic "schoolbook" algorithm. However, for sufficiently large values,
180 there is a more efficient algorithm usually attributed to Karatsuba and Ofman
181 that is usually faster. See Knuth Vol. 2 for more details about how this
184 The breakpoint between the "normal" and the recursive algorithm is controlled
185 by a static digit threshold defined in `imath.c`. Values with fewer significant
186 digits use the standard algorithm. This value can be modified by calling
187 `mp_int_multiply_threshold(n)`. The `imtimer` program and the
188 `findthreshold.py` script (Python) can help you find a suitable value for for
189 your particular platform.
191 {{insert "imath.h" mp_error_string}}
193 ## Rational Arithmetic
197 ## Representation Details
199 > NOTE: You do not need to read this section to use IMath. This is provided
200 > for the benefit of developers wishing to extend or modify the internals of
203 IMath uses a signed magnitude representation for arbitrary precision integers.
204 The magnitude is represented as an array of radix-R digits in increasing order
205 of significance; the value of R is chosen to be half the size of the largest
206 available unsigned integer type, so typically 16 or 32 bits. Digits are
207 represented as mp_digit, which must be an unsigned integral type.
209 Digit arrays are allocated using `malloc(3)` and `realloc(3)`. Because this
210 can be an expensive operation, the library takes pains to avoid allocation as
211 much as possible. For this reason, the `mpz_t` structure distinguishes between
212 how many digits are allocated and how many digits are actually consumed by the
213 representation. The fields of an `mpz_t` are:
215 mp_digit single; /* single-digit value (see note) */
216 mp_digit *digits; /* array of digits */
217 mp_size alloc; /* how many digits are allocated */
218 mp_size used; /* how many digits are in use */
219 mp_sign sign; /* the sign of the value */
221 The elements of `digits` at indices less than `used` are the significant
222 figures of the value; the elements at indices greater than or equal to `used`
223 are undefined (and may contain garbage). At all times, `used` must be at least
224 1 and at most `alloc`.
226 To avoid interaction with the memory allocator, single-digit values are stored
227 directly in the `mpz_t` structure, in the `single` field. The semantics of
228 access are the same as the more general case.
230 The number of digits allocated for an `mpz_t` is referred to in the library
231 documentation as its "precision". Operations that affect an `mpz_t` cause
232 precision to increase as needed. In any case, all allocations are measured in
233 digits, and rounded up to the nearest `mp_word` boundary. There is a default
234 minimum precision stored as a static constant default_precision (`imath.c`).
235 This value can be set using `mp_int_default_precision(n)`.
237 Note that the allocated size of an `mpz_t` can only grow; the library never
238 reallocates in order to decrease the size. A simple way to do so explicitly is
239 to use `mp_int_init_copy()`, as in:
245 mp_int_init_copy(&new, &big);
246 mp_int_swap(&new, &big);
250 The value of `sign` is 0 for positive values and zero, 1 for negative values.
251 Constants `MP_ZPOS` and `MP_NEG` are defined for these; no other sign values
254 If you are adding to this library, you should be careful to preserve the
255 convention that inputs and outputs can overlap, as described above. So, for
256 example, `mp_int_add(a, a, a)` is legal. Often, this means you must maintain
257 one or more temporary mpz_t structures for intermediate values. The private
258 macros `DECLARE_TEMP(N)`, `CLEANUP_TEMP()`, and `TEMP(K)` can be used to
259 maintain a conventional structure like this:
263 /* Declare how many temp values you need.
264 Use TEMP(i) to access the ith value (0-indexed). */
268 /* Perform actions that must return MP_OK or fail. */
269 REQUIRE(mp_int_copy(x, TEMP(1)));
271 REQUIRE(mp_int_expt(TEMP(1), TEMP(2), TEMP(3)));
274 /* You can also use REQUIRE directly for more complex cases. */
275 if (some_difficult_question(TEMP(3)) != answer(x)) {
276 REQUIRE(MP_RANGE); /* falls through to cleanup (below) */
279 /* Ensure temporary values are cleaned up at exit.
281 If control reaches here via a REQUIRE failure, the code below
282 the cleanup will not be executed.
289 Under the covers, these macros are just maintaining an array of `mpz_t` values,
290 and a jump label to handle cleanup. You may only have one `DECLARE_TEMP` and
291 its corresponding `CLEANUP_TEMP` per function body.
293 "Small" integer values are represented by the types `mp_small` and `mp_usmall`,
294 which are mapped to appropriately-sized types on the host system. The default
295 for `mp_small` is "long" and the default for `mp_usmall` is "unsigned long".
296 You may change these, provided you insure that `mp_small` is signed and
297 `mp_usmall` is unsigned. You will also need to adjust the size macros:
299 MP_SMALL_MIN, MP_SMALL_MAX
300 MP_USMALL_MIN, MP_USMALL_MAX
302 ... which are defined in `<imath.h>`, if you change these.
304 Rational numbers are represented using a pair of arbitrary precision integers,
305 with the convention that the sign of the numerator is the sign of the rational
306 value, and that the result of any rational operation is always represented in
307 lowest terms. The canonical representation for rational zero is 0/1. See
310 ## Testing and Reporting of Bugs
312 Test vectors are included in the `tests/` subdirectory of the imath
313 distribution. When you run `make test`, it builds the `imtest` program and
314 runs all available test vectors. If any tests fail, you will get a line like
319 Here, _x_ is the line number of the test which failed, _y_ is index of the test
320 within the file, and _v_ is the value(s) actually computed. The name of the
321 file is printed at the beginning of each test, so you can find out what test
322 vector failed by executing the following (with x, y, and v replaced by the
323 above values, and where "foo.t" is the name of the test file that was being
324 processed at the time):
326 % tail +x tests/foo.t | head -1
328 None of the tests should fail (but see [Note 2](#note2)); if any do, it
329 probably indicates a bug in the library (or at the very least, some assumption
330 I made which I shouldn't have). Please [file an
331 issue](https://github.com/creachadair/imath/issues/new), including the `FAILED`
332 test line(s), as well as the output of the above `tail` command (so I know what
333 inputs caused the failure).
335 If you build with the preprocessor symbol `DEBUG` defined as a positive
336 integer, the digit allocators (`s_alloc`, `s_realloc`) fill all new buffers
337 with the value `0xdeadbeefabad1dea`, or as much of it as will fit in a digit,
338 so that you can more easily catch uninitialized reads in the debugger.
342 1. <a name="note1"></a>You can generally use the same variables for both input
343 and output. One exception is that you may not use the same variable for
344 both the quotient and the remainder of `mp_int_div()`.
346 2. <a name="note2"></a>Many of the tests for this library were written under
347 the assumption that the `mp_small` type is 32 bits or more. If you compile
348 with a smaller type, you may see `MP_RANGE` errors in some of the tests that
349 otherwise pass (due to conversion failures). Also, the pi generator (pi.c)
350 will not work correctly if `mp_small` is too short, as its algorithm for arc
351 tangent is fairly simple-minded.
355 The IMath library was written by Michael J. Fromberger.
357 If you discover any bugs or testing failures, please [open an
358 issue](https://github.com/creachadair/imath/issues/new). Please be sure to
359 include a complete description of what went wrong, and if possible, a test
360 vector for `imtest` and/or a minimal test program that will demonstrate the bug
361 on your system. Please also let me know what hardware, operating system, and
362 compiler you're using.
366 The algorithms used in this library came from Vol. 2 of Donald Knuth's "The Art
367 of Computer Programming" (Seminumerical Algorithms). Thanks to Nelson Bolyard,
368 Bryan Olson, Tom St. Denis, Tushar Udeshi, and Eric Silva for excellent
369 feedback on earlier versions of this code. Special thanks to Jonathan Shapiro
370 for some very helpful design advice, as well as feedback and some clever ideas
371 for improving performance in some common use cases.
373 ## License and Disclaimers
375 IMath is Copyright 2002-2009 Michael J. Fromberger
376 You may use it subject to the following Licensing Terms:
378 Permission is hereby granted, free of charge, to any person obtaining a copy of
379 this software and associated documentation files (the "Software"), to deal in
380 the Software without restriction, including without limitation the rights to
381 use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
382 of the Software, and to permit persons to whom the Software is furnished to do
383 so, subject to the following conditions:
385 The above copyright notice and this permission notice shall be included in all
386 copies or substantial portions of the Software.
388 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
389 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
390 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
391 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
392 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
393 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE