1# User Documentation for the IMath Library 2 3Author: [M. J. Fromberger](https://github.com/creachadair) 4 5## Installation 6 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 that 9 what you will get anyway), but you should be able to use the default GCC 10 settings for either. 11 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. 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 report 18 bugs. 19 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). 23 24To use the library in your code, include "imath.h" wherever you intend to use 25the library's routines. The integer library is just a single source file, so 26you can compile it into your project in whatever way makes sense. If you wish 27to use rational arithmetic, you will also need to include "imrat.h". 28 29## Background 30 31The basic types defined by the imath library are `mpz_t`, an arbitrary 32precision signed integer, and `mpq_t`, an arbitrary precision signed rational 33number. The type `mp_int` is a pointer to an `mpz_t`, and `mp_rat` is a 34pointer 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 information 38and also return small values. Any negative value is considered to be a status 39message. 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 the 53`mp_error_string()` routine to obtain a pointer to a brief human-readable 54string describing the error. These strings are statically allocated, so they 55need not be freed by the caller; the same strings are re-used from call to 56call. 57 58Unless otherwise noted, it is legal to use the same parameter for both inputs 59and output with most of the functions in this library. For example, you can 60add 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 summaries 65below (if you discover that this is not so, please report it as a bug; I will 66fix either the function or the documentation :) 67 68## The IMath API 69 70Each of the API functions is documented here. The general format of the 71entries 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 be 80expected to return `MP_OK`, `MP_BADARG`, or `MP_MEMORY`. Other return values 81should be documented in the description. Please let me know if you discover 82this is not the case. 83 84The following macros are defined in "imath.h", to define the sizes of the 85various data types used in the library: 86 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. 96 97#### Initialization 98 99An `mp_int` must be initialized before use. By default, an `mp_int` is 100initialized with a certain minimum amount of storage for digits, and the 101storage is expanded automatically as needed. To initialize an `mp_int`, use 102the following functions: 103 104{{insert "imath.h" 105 mp_int_init mp_int_alloc mp_int_init_size 106 mp_int_init_copy 107 mp_int_init_value 108}} 109 110#### Cleanup 111 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 Values 117 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 Functions 128 129{{insert "imath.h" 130 mp_int_is_odd mp_int_is_even 131 mp_int_zero 132 mp_int_abs 133 mp_int_neg 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 137 mp_int_sqr 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 142}} 143 144### Comparison Functions 145 146Unless 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 148to `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_zero 152 mp_int_compare_value mp_int_compare_uvalue 153 mp_int_divisible_value mp_int_is_pow2 154}} 155 156### Modular Operations 157 158{{insert "imath.h" 159 mp_int_exptmod mp_int_exptmod_evalue mp_int_exptmod_bvalue 160 mp_int_exptmod_known mp_int_redux_const 161 mp_int_invmod 162 mp_int_gcd mp_int_egcd mp_int_lcm 163}} 164 165### Conversion of Values 166 167{{insert "imath.h" 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 171 mp_int_count_bits 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 174}} 175 176### Other Functions 177 178Ordinarily, integer multiplication and squaring are done using the simple 179quadratic "schoolbook" algorithm. However, for sufficiently large values, 180there is a more efficient algorithm usually attributed to Karatsuba and Ofman 181that is usually faster. See Knuth Vol. 2 for more details about how this 182algorithm works. 183 184The breakpoint between the "normal" and the recursive algorithm is controlled 185by a static digit threshold defined in `imath.c`. Values with fewer significant 186digits 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 189your particular platform. 190 191{{insert "imath.h" mp_error_string}} 192 193## Rational Arithmetic 194 195{{insert "imrat.h"}} 196 197## Representation Details 198 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 201> 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 order 205of significance; the value of R is chosen to be half the size of the largest 206available unsigned integer type, so typically 16 or 32 bits. Digits are 207represented as mp_digit, which must be an unsigned integral type. 208 209Digit arrays are allocated using `malloc(3)` and `realloc(3)`. Because this 210can be an expensive operation, the library takes pains to avoid allocation as 211much as possible. For this reason, the `mpz_t` structure distinguishes between 212how many digits are allocated and how many digits are actually consumed by the 213representation. 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 significant 222figures 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 least 2241 and at most `alloc`. 225 226To avoid interaction with the memory allocator, single-digit values are stored 227directly in the `mpz_t` structure, in the `single` field. The semantics of 228access are the same as the more general case. 229 230The number of digits allocated for an `mpz_t` is referred to in the library 231documentation as its "precision". Operations that affect an `mpz_t` cause 232precision to increase as needed. In any case, all allocations are measured in 233digits, and rounded up to the nearest `mp_word` boundary. There is a default 234minimum 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 never 238reallocates in order to decrease the size. A simple way to do so explicitly is 239to 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 values 252are used. 253 254If you are adding to this library, you should be careful to preserve the 255convention that inputs and outputs can overlap, as described above. So, for 256example, `mp_int_add(a, a, a)` is legal. Often, this means you must maintain 257one or more temporary mpz_t structures for intermediate values. The private 258macros `DECLARE_TEMP(N)`, `CLEANUP_TEMP()`, and `TEMP(K)` can be used to 259maintain a conventional structure like this: 260 261```c 262{ 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 below 282 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` and 291its 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 default 295for `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 and 297`mp_usmall` is unsigned. You will also need to adjust the size macros: 298 299 MP_SMALL_MIN, MP_SMALL_MAX 300 MP_USMALL_MIN, MP_USMALL_MAX 301 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 rational 306value, and that the result of any rational operation is always represented in 307lowest terms. The canonical representation for rational zero is 0/1. See 308"imrat.h". 309 310## Testing and Reporting of Bugs 311 312Test vectors are included in the `tests/` subdirectory of the imath 313distribution. When you run `make test`, it builds the `imtest` program and 314runs all available test vectors. If any tests fail, you will get a line like 315this: 316 317 x y FAILED v 318 319Here, _x_ is the line number of the test which failed, _y_ is index of the test 320within the file, and _v_ is the value(s) actually computed. The name of the 321file is printed at the beginning of each test, so you can find out what test 322vector failed by executing the following (with x, y, and v replaced by the 323above values, and where "foo.t" is the name of the test file that was being 324processed at the time): 325 326 % tail +x tests/foo.t | head -1 327 328None of the tests should fail (but see [Note 2](#note2)); if any do, it 329probably indicates a bug in the library (or at the very least, some assumption 330I made which I shouldn't have). Please [file an 331issue](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 what 333inputs caused the failure). 334 335If you build with the preprocessor symbol `DEBUG` defined as a positive 336integer, the digit allocators (`s_alloc`, `s_realloc`) fill all new buffers 337with 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## Notes 341 3421. <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()`. 345 3462. <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. 352 353## Contacts 354 355The IMath library was written by Michael J. Fromberger. 356 357If you discover any bugs or testing failures, please [open an 358issue](https://github.com/creachadair/imath/issues/new). Please be sure to 359include a complete description of what went wrong, and if possible, a test 360vector for `imtest` and/or a minimal test program that will demonstrate the bug 361on your system. Please also let me know what hardware, operating system, and 362compiler you're using. 363 364## Acknowledgements 365 366The algorithms used in this library came from Vol. 2 of Donald Knuth's "The Art 367of Computer Programming" (Seminumerical Algorithms). Thanks to Nelson Bolyard, 368Bryan Olson, Tom St. Denis, Tushar Udeshi, and Eric Silva for excellent 369feedback on earlier versions of this code. Special thanks to Jonathan Shapiro 370for some very helpful design advice, as well as feedback and some clever ideas 371for improving performance in some common use cases. 372 373## License and Disclaimers 374 375IMath is Copyright 2002-2009 Michael J. Fromberger 376You may use it subject to the following Licensing Terms: 377 378Permission is hereby granted, free of charge, to any person obtaining a copy of 379this software and associated documentation files (the "Software"), to deal in 380the Software without restriction, including without limitation the rights to 381use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies 382of the Software, and to permit persons to whom the Software is furnished to do 383so, subject to the following conditions: 384 385The above copyright notice and this permission notice shall be included in all 386copies or substantial portions of the Software. 387 388THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 389IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 390FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 391AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 392LIABILITY, 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 THE 394SOFTWARE. 395