395 lines · plain
1# User Documentation for the IMath Library2 3Author: [M. J. Fromberger](https://github.com/creachadair)4 5## Installation6 71. 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 that9 what you will get anyway), but you should be able to use the default GCC10 settings for either.11 12 By default, the Makefile assumes you can use 64-bit integer types, even13 though they were not standard in ANSI C90. If you cannot, add14 `-DUSE_32BIT_WORDS` to the compiler options.15 162. 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 report18 bugs.19 20 To build with debugging enabled (and optimization disabled), run `make21 DEBUG=Y`. This sets the preprocessor macro `DEBUG` to 1, and several other22 things (see Makefile for details).23 24To use the library in your code, include "imath.h" wherever you intend to use25the library's routines. The integer library is just a single source file, so26you can compile it into your project in whatever way makes sense. If you wish27to use rational arithmetic, you will also need to include "imrat.h".28 29## Background30 31The basic types defined by the imath library are `mpz_t`, an arbitrary32precision signed integer, and `mpq_t`, an arbitrary precision signed rational33number. The type `mp_int` is a pointer to an `mpz_t`, and `mp_rat` is a34pointer to an `mpq_t`.35 36Most of the functions in the imath library return a value of type `mp_result`.37This is a signed integer type which can be used to convey status information38and also return small values. Any negative value is considered to be a status39message. The following constants are defined for processing these:40 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 |51 52If you obtain a zero or negative value of an `mp_result`, you can use the53`mp_error_string()` routine to obtain a pointer to a brief human-readable54string describing the error. These strings are statically allocated, so they55need not be freed by the caller; the same strings are re-used from call to56call.57 58Unless otherwise noted, it is legal to use the same parameter for both inputs59and output with most of the functions in this library. For example, you can60add a number to itself and replace the original by writing:61 62 mp_int_add(a, a, a); /* a = a + a */63 64Any cases in which this is not legal will be noted in the function summaries65below (if you discover that this is not so, please report it as a bug; I will66fix either the function or the documentation :)67 68## The IMath API69 70Each of the API functions is documented here. The general format of the71entries is:72 73> ------------74> <pre>75> return_type function_name(parameters ...)76> </pre>77> - English description.78 79Unless otherwise noted, any API function that returns `mp_result` may be80expected to return `MP_OK`, `MP_BADARG`, or `MP_MEMORY`. Other return values81should be documented in the description. Please let me know if you discover82this is not the case.83 84The following macros are defined in "imath.h", to define the sizes of the85various data types used in the library:86 87| Constant | Description88| --------------- | ----------------------------------------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.96 97#### Initialization98 99An `mp_int` must be initialized before use. By default, an `mp_int` is100initialized with a certain minimum amount of storage for digits, and the101storage is expanded automatically as needed. To initialize an `mp_int`, use102the following functions:103 104{{insert "imath.h"105 mp_int_init mp_int_alloc mp_int_init_size106 mp_int_init_copy107 mp_int_init_value108}}109 110#### Cleanup111 112When you are finished with an `mp_int`, you must free the memory it uses:113 114{{insert "imath.h" mp_int_clear mp_int_free}}115 116#### Setting Values117 118To set an `mp_int` which has already been initialized to a small integer value,119use:120 121{{insert "imath.h" mp_int_set_value mp_int_set_uvalue}}122 123To copy one initialized `mp_int` to another, use:124 125{{insert "imath.h" mp_int_copy}}126 127### Arithmetic Functions128 129{{insert "imath.h"130 mp_int_is_odd mp_int_is_even131 mp_int_zero132 mp_int_abs133 mp_int_neg134 mp_int_add mp_int_add_value135 mp_int_sub mp_int_sub_value136 mp_int_mul mp_int_mul_value mp_int_mul_pow2137 mp_int_sqr138 mp_int_root mp_int_sqrt139 mp_int_div mp_int_div_value mp_int_div_pow2140 mp_int_mod mp_int_mod_value141 mp_int_expt mp_int_expt_value mp_int_expt_full142}}143 144### Comparison Functions145 146Unless otherwise specified, comparison between values `x` and `y` returns a147**comparator**, an integer value < 0 if `x` is less than `y`, 0 if `x` is equal148to `y`, and > 0 if `x` is greater than `y`.149 150{{insert "imath.h"151 mp_int_compare mp_int_compare_unsigned mp_int_compare_zero152 mp_int_compare_value mp_int_compare_uvalue153 mp_int_divisible_value mp_int_is_pow2154}}155 156### Modular Operations157 158{{insert "imath.h"159 mp_int_exptmod mp_int_exptmod_evalue mp_int_exptmod_bvalue160 mp_int_exptmod_known mp_int_redux_const161 mp_int_invmod162 mp_int_gcd mp_int_egcd mp_int_lcm163}}164 165### Conversion of Values166 167{{insert "imath.h"168 mp_int_to_int mp_int_to_uint169 mp_int_to_string mp_int_string_len170 mp_int_read_string mp_int_read_cstring171 mp_int_count_bits172 mp_int_to_binary mp_int_read_binary mp_int_binary_len173 mp_int_to_unsigned mp_int_read_unsigned mp_int_unsigned_len174}}175 176### Other Functions177 178Ordinarily, integer multiplication and squaring are done using the simple179quadratic "schoolbook" algorithm. However, for sufficiently large values,180there is a more efficient algorithm usually attributed to Karatsuba and Ofman181that is usually faster. See Knuth Vol. 2 for more details about how this182algorithm works.183 184The breakpoint between the "normal" and the recursive algorithm is controlled185by a static digit threshold defined in `imath.c`. Values with fewer significant186digits use the standard algorithm. This value can be modified by calling187`mp_int_multiply_threshold(n)`. The `imtimer` program and the188`findthreshold.py` script (Python) can help you find a suitable value for for189your particular platform.190 191{{insert "imath.h" mp_error_string}}192 193## Rational Arithmetic194 195{{insert "imrat.h"}}196 197## Representation Details198 199> NOTE: You do not need to read this section to use IMath. This is provided200> for the benefit of developers wishing to extend or modify the internals of201> the library.202 203IMath uses a signed magnitude representation for arbitrary precision integers.204The magnitude is represented as an array of radix-R digits in increasing order205of significance; the value of R is chosen to be half the size of the largest206available unsigned integer type, so typically 16 or 32 bits. Digits are207represented as mp_digit, which must be an unsigned integral type.208 209Digit arrays are allocated using `malloc(3)` and `realloc(3)`. Because this210can be an expensive operation, the library takes pains to avoid allocation as211much as possible. For this reason, the `mpz_t` structure distinguishes between212how many digits are allocated and how many digits are actually consumed by the213representation. The fields of an `mpz_t` are:214 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 */220 221The elements of `digits` at indices less than `used` are the significant222figures of the value; the elements at indices greater than or equal to `used`223are undefined (and may contain garbage). At all times, `used` must be at least2241 and at most `alloc`.225 226To avoid interaction with the memory allocator, single-digit values are stored227directly in the `mpz_t` structure, in the `single` field. The semantics of228access are the same as the more general case.229 230The number of digits allocated for an `mpz_t` is referred to in the library231documentation as its "precision". Operations that affect an `mpz_t` cause232precision to increase as needed. In any case, all allocations are measured in233digits, and rounded up to the nearest `mp_word` boundary. There is a default234minimum precision stored as a static constant default_precision (`imath.c`).235This value can be set using `mp_int_default_precision(n)`.236 237Note that the allocated size of an `mpz_t` can only grow; the library never238reallocates in order to decrease the size. A simple way to do so explicitly is239to use `mp_int_init_copy()`, as in:240 241```242mpz_t big, new;243 244/* ... */245mp_int_init_copy(&new, &big);246mp_int_swap(&new, &big);247mp_int_clear(&new);248```249 250The value of `sign` is 0 for positive values and zero, 1 for negative values.251Constants `MP_ZPOS` and `MP_NEG` are defined for these; no other sign values252are used.253 254If you are adding to this library, you should be careful to preserve the255convention that inputs and outputs can overlap, as described above. So, for256example, `mp_int_add(a, a, a)` is legal. Often, this means you must maintain257one or more temporary mpz_t structures for intermediate values. The private258macros `DECLARE_TEMP(N)`, `CLEANUP_TEMP()`, and `TEMP(K)` can be used to259maintain a conventional structure like this:260 261```c262{263 /* Declare how many temp values you need.264 Use TEMP(i) to access the ith value (0-indexed). */265 DECLARE_TEMP(8);266 ...267 268 /* Perform actions that must return MP_OK or fail. */269 REQUIRE(mp_int_copy(x, TEMP(1)));270 ...271 REQUIRE(mp_int_expt(TEMP(1), TEMP(2), TEMP(3)));272 ...273 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) */277 }278 279 /* Ensure temporary values are cleaned up at exit.280 281 If control reaches here via a REQUIRE failure, the code below282 the cleanup will not be executed.283 */284 CLEANUP_TEMP();285 return MP_OK;286}287```288 289Under the covers, these macros are just maintaining an array of `mpz_t` values,290and a jump label to handle cleanup. You may only have one `DECLARE_TEMP` and291its corresponding `CLEANUP_TEMP` per function body.292 293"Small" integer values are represented by the types `mp_small` and `mp_usmall`,294which are mapped to appropriately-sized types on the host system. The default295for `mp_small` is "long" and the default for `mp_usmall` is "unsigned long".296You may change these, provided you insure that `mp_small` is signed and297`mp_usmall` is unsigned. You will also need to adjust the size macros:298 299 MP_SMALL_MIN, MP_SMALL_MAX300 MP_USMALL_MIN, MP_USMALL_MAX301 302... which are defined in `<imath.h>`, if you change these.303 304Rational numbers are represented using a pair of arbitrary precision integers,305with the convention that the sign of the numerator is the sign of the rational306value, and that the result of any rational operation is always represented in307lowest terms. The canonical representation for rational zero is 0/1. See308"imrat.h".309 310## Testing and Reporting of Bugs311 312Test vectors are included in the `tests/` subdirectory of the imath313distribution. When you run `make test`, it builds the `imtest` program and314runs all available test vectors. If any tests fail, you will get a line like315this:316 317 x y FAILED v318 319Here, _x_ is the line number of the test which failed, _y_ is index of the test320within the file, and _v_ is the value(s) actually computed. The name of the321file is printed at the beginning of each test, so you can find out what test322vector failed by executing the following (with x, y, and v replaced by the323above values, and where "foo.t" is the name of the test file that was being324processed at the time):325 326 % tail +x tests/foo.t | head -1327 328None of the tests should fail (but see [Note 2](#note2)); if any do, it329probably indicates a bug in the library (or at the very least, some assumption330I made which I shouldn't have). Please [file an331issue](https://github.com/creachadair/imath/issues/new), including the `FAILED`332test line(s), as well as the output of the above `tail` command (so I know what333inputs caused the failure).334 335If you build with the preprocessor symbol `DEBUG` defined as a positive336integer, the digit allocators (`s_alloc`, `s_realloc`) fill all new buffers337with the value `0xdeadbeefabad1dea`, or as much of it as will fit in a digit,338so that you can more easily catch uninitialized reads in the debugger.339 340## Notes341 3421. <a name="note1"></a>You can generally use the same variables for both input343 and output. One exception is that you may not use the same variable for344 both the quotient and the remainder of `mp_int_div()`.345 3462. <a name="note2"></a>Many of the tests for this library were written under347 the assumption that the `mp_small` type is 32 bits or more. If you compile348 with a smaller type, you may see `MP_RANGE` errors in some of the tests that349 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 arc351 tangent is fairly simple-minded.352 353## Contacts354 355The IMath library was written by Michael J. Fromberger.356 357If you discover any bugs or testing failures, please [open an358issue](https://github.com/creachadair/imath/issues/new). Please be sure to359include a complete description of what went wrong, and if possible, a test360vector for `imtest` and/or a minimal test program that will demonstrate the bug361on your system. Please also let me know what hardware, operating system, and362compiler you're using.363 364## Acknowledgements365 366The algorithms used in this library came from Vol. 2 of Donald Knuth's "The Art367of Computer Programming" (Seminumerical Algorithms). Thanks to Nelson Bolyard,368Bryan Olson, Tom St. Denis, Tushar Udeshi, and Eric Silva for excellent369feedback on earlier versions of this code. Special thanks to Jonathan Shapiro370for some very helpful design advice, as well as feedback and some clever ideas371for improving performance in some common use cases.372 373## License and Disclaimers374 375IMath is Copyright 2002-2009 Michael J. Fromberger376You may use it subject to the following Licensing Terms:377 378Permission is hereby granted, free of charge, to any person obtaining a copy of379this software and associated documentation files (the "Software"), to deal in380the Software without restriction, including without limitation the rights to381use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies382of the Software, and to permit persons to whom the Software is furnished to do383so, subject to the following conditions:384 385The above copyright notice and this permission notice shall be included in all386copies or substantial portions of the Software.387 388THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR389IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,390FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE391AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER392LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,393OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE394SOFTWARE.395