1255 lines · c
1/*2 * Copyright 2015 INRIA Paris-Rocquencourt3 *4 * Use of this software is governed by the MIT license5 *6 * Written by Michael Kruse, INRIA Paris-Rocquencourt,7 * Domaine de Voluceau, Rocquenqourt, B.P. 105,8 * 78153 Le Chesnay Cedex France9 */10#ifndef ISL_INT_SIOIMATH_H11#define ISL_INT_SIOIMATH_H12 13#include <inttypes.h>14#include <limits.h>15#include <stdint.h>16#include <stdlib.h>17 18#include <isl_imath.h>19#include <isl/hash.h>20 21#define ARRAY_SIZE(array) (sizeof(array)/sizeof(*array))22 23/* Visual Studio before VS2015 does not support the inline keyword when24 * compiling in C mode because it was introduced in C99 which it does not25 * officially support. Instead, it has a proprietary extension using __inline.26 */27#if defined(_MSC_VER) && (_MSC_VER < 1900)28#define inline __inline29#endif30 31/* The type to represent integers optimized for small values. It is either a32 * pointer to an mp_int ( = mpz_t*; big representation) or an int32_t (small33 * represenation) with a discriminator at the least significant bit. In big34 * representation it will be always zero because of heap alignment. It is set35 * to 1 for small representation and use the 32 most significant bits for the36 * int32_t.37 *38 * Structure on 64 bit machines, with 8-byte aligment (3 bits):39 *40 * Big representation:41 * MSB LSB42 * |------------------------------------------------------------00043 * | mpz_t* |44 * | != NULL |45 *46 * Small representation:47 * MSB 32 LSB48 * |------------------------------|0000000000000000000000000000000149 * | int32_t |50 * | 2147483647 ... -2147483647 |51 * ^52 * |53 * discriminator bit54 *55 * On 32 bit machines isl_sioimath type is blown up to 8 bytes, i.e.56 * isl_sioimath is guaranteed to be at least 8 bytes. This is to ensure the57 * int32_t can be hidden in that type without data loss. In the future we might58 * optimize this to use 31 hidden bits in a 32 bit pointer. We may also use 6359 * bits on 64 bit machines, but this comes with the cost of additional overflow60 * checks because there is no standardized 128 bit integer we could expand to.61 *62 * We use native integer types and avoid union structures to avoid assumptions63 * on the machine's endianness.64 *65 * This implementation makes the following assumptions:66 * - long can represent any int32_t67 * - mp_small is signed long68 * - mp_usmall is unsigned long69 * - adresses returned by malloc are aligned to 2-byte boundaries (leastmost70 * bit is zero)71 */72#if UINT64_MAX > UINTPTR_MAX73typedef uint64_t isl_sioimath;74#else75typedef uintptr_t isl_sioimath;76#endif77 78/* The negation of the smallest possible number in int32_t, INT32_MIN79 * (0x80000000u, -2147483648), cannot be represented in an int32_t, therefore80 * every operation that may produce this value needs to special-case it.81 * The operations are:82 * abs(INT32_MIN)83 * -INT32_MIN (negation)84 * -1 * INT32_MIN (multiplication)85 * INT32_MIN/-1 (any division: divexact, fdiv, cdiv, tdiv)86 * To avoid checking these cases, we exclude INT32_MIN from small87 * representation.88 */89#define ISL_SIOIMATH_SMALL_MIN (-INT32_MAX)90 91/* Largest possible number in small representation */92#define ISL_SIOIMATH_SMALL_MAX INT32_MAX93 94/* Used for function parameters the function modifies. */95typedef isl_sioimath *isl_sioimath_ptr;96 97/* Used for function parameters that are read-only. */98typedef isl_sioimath isl_sioimath_src;99 100/* Return whether the argument is stored in small representation.101 */102inline int isl_sioimath_is_small(isl_sioimath val)103{104 return val & 0x00000001;105}106 107/* Return whether the argument is stored in big representation.108 */109inline int isl_sioimath_is_big(isl_sioimath val)110{111 return !isl_sioimath_is_small(val);112}113 114/* Get the number of an isl_int in small representation. Result is undefined if115 * val is not stored in that format.116 */117inline int32_t isl_sioimath_get_small(isl_sioimath val)118{119 return val >> 32;120}121 122/* Get the number of an in isl_int in big representation. Result is undefined if123 * val is not stored in that format.124 */125inline mp_int isl_sioimath_get_big(isl_sioimath val)126{127 return (mp_int)(uintptr_t) val;128}129 130/* Return 1 if val is stored in small representation and store its value to131 * small. We rely on the compiler to optimize the isl_sioimath_get_small such132 * that the shift is moved into the branch that executes in case of small133 * representation. If there is no such branch, then a single shift is still134 * cheaper than introducing branching code.135 */136inline int isl_sioimath_decode_small(isl_sioimath val, int32_t *small)137{138 *small = isl_sioimath_get_small(val);139 return isl_sioimath_is_small(val);140}141 142/* Return 1 if val is stored in big representation and store its value to big.143 */144inline int isl_sioimath_decode_big(isl_sioimath val, mp_int *big)145{146 *big = isl_sioimath_get_big(val);147 return isl_sioimath_is_big(val);148}149 150/* Encode a small representation into an isl_int.151 */152inline isl_sioimath isl_sioimath_encode_small(int32_t val)153{154 return ((isl_sioimath) val) << 32 | 0x00000001;155}156 157/* Encode a big representation.158 */159inline isl_sioimath isl_sioimath_encode_big(mp_int val)160{161 return (isl_sioimath)(uintptr_t) val;162}163 164/* A common situation is to call an IMath function with at least one argument165 * that is currently in small representation or an integer parameter, i.e. a big166 * representation of the same number is required. Promoting the original167 * argument comes with multiple problems, such as modifying a read-only168 * argument, the responsibility of deallocation and the execution cost. Instead,169 * we make a copy by 'faking' the IMath internal structure.170 *171 * We reserve the maximum number of required digits on the stack to avoid heap172 * allocations.173 *174 * mp_digit can be uint32_t or uint16_t. This code must work for little and big175 * endian digits. The structure for an uint64_t argument and 32-bit mp_digits is176 * sketched below.177 *178 * |----------------------------|179 * uint64_t180 *181 * |-------------||-------------|182 * mp_digit mp_digit183 * digits[1] digits[0]184 * Most sig digit Least sig digit185 */186typedef struct {187 mpz_t big;188 mp_digit digits[(sizeof(uintmax_t) + sizeof(mp_digit) - 1) /189 sizeof(mp_digit)];190} isl_sioimath_scratchspace_t;191 192/* Convert a native integer to IMath's digit representation. A native integer193 * might be big- or little endian, but IMath always stores the least significant194 * digit in the lowest array indices. memcpy therefore is not possible.195 *196 * We also have to consider that long and mp_digit can be of different sizes,197 * depending on the compiler (LP64, LLP64) and IMath's USE_64BIT_WORDS. This198 * macro should work for all of them.199 *200 * "used" is set to the number of written digits. It must be minimal (IMath201 * checks zeroness using the used field), but always at least one. Also note202 * that the result of num>>(sizeof(num)*CHAR_BIT) is undefined.203 */204#define ISL_SIOIMATH_TO_DIGITS(num, digits, used) \205 do { \206 int i = 0; \207 do { \208 (digits)[i] = \209 ((num) >> (sizeof(mp_digit) * CHAR_BIT * i)); \210 i += 1; \211 if (i >= (sizeof(num) + sizeof(mp_digit) - 1) / \212 sizeof(mp_digit)) \213 break; \214 if (((num) >> (sizeof(mp_digit) * CHAR_BIT * i)) == 0) \215 break; \216 } while (1); \217 (used) = i; \218 } while (0)219 220inline void isl_siomath_uint32_to_digits(uint32_t num, mp_digit *digits,221 mp_size *used)222{223 ISL_SIOIMATH_TO_DIGITS(num, digits, *used);224}225 226inline void isl_siomath_ulong_to_digits(unsigned long num, mp_digit *digits,227 mp_size *used)228{229 ISL_SIOIMATH_TO_DIGITS(num, digits, *used);230}231 232inline void isl_siomath_uint64_to_digits(uint64_t num, mp_digit *digits,233 mp_size *used)234{235 ISL_SIOIMATH_TO_DIGITS(num, digits, *used);236}237 238/* Get the IMath representation of an isl_int without modifying it.239 * For the case it is not in big representation yet, pass some scratch space we240 * can use to store the big representation in.241 * In order to avoid requiring init and free on the scratch space, we directly242 * modify the internal representation.243 *244 * The name derives from its indented use: getting the big representation of an245 * input (src) argument.246 */247inline mp_int isl_sioimath_bigarg_src(isl_sioimath arg,248 isl_sioimath_scratchspace_t *scratch)249{250 mp_int big;251 int32_t small;252 uint32_t num;253 254 if (isl_sioimath_decode_big(arg, &big))255 return big;256 257 small = isl_sioimath_get_small(arg);258 scratch->big.digits = scratch->digits;259 scratch->big.alloc = ARRAY_SIZE(scratch->digits);260 if (small >= 0) {261 scratch->big.sign = MP_ZPOS;262 num = small;263 } else {264 scratch->big.sign = MP_NEG;265 num = -small;266 }267 268 isl_siomath_uint32_to_digits(num, scratch->digits, &scratch->big.used);269 return &scratch->big;270}271 272/* Create a temporary IMath mp_int for a signed long.273 */274inline mp_int isl_sioimath_siarg_src(signed long arg,275 isl_sioimath_scratchspace_t *scratch)276{277 unsigned long num;278 279 scratch->big.digits = scratch->digits;280 scratch->big.alloc = ARRAY_SIZE(scratch->digits);281 if (arg >= 0) {282 scratch->big.sign = MP_ZPOS;283 num = arg;284 } else {285 scratch->big.sign = MP_NEG;286 num = (arg == LONG_MIN) ? ((unsigned long) LONG_MAX) + 1 : -arg;287 }288 289 isl_siomath_ulong_to_digits(num, scratch->digits, &scratch->big.used);290 return &scratch->big;291}292 293/* Create a temporary IMath mp_int for an int64_t.294 */295inline mp_int isl_sioimath_si64arg_src(int64_t arg,296 isl_sioimath_scratchspace_t *scratch)297{298 uint64_t num;299 300 scratch->big.digits = scratch->digits;301 scratch->big.alloc = ARRAY_SIZE(scratch->digits);302 if (arg >= 0) {303 scratch->big.sign = MP_ZPOS;304 num = arg;305 } else {306 scratch->big.sign = MP_NEG;307 num = (arg == INT64_MIN) ? ((uint64_t) INT64_MAX) + 1 : -arg;308 }309 310 isl_siomath_uint64_to_digits(num, scratch->digits, &scratch->big.used);311 return &scratch->big;312}313 314/* Create a temporary IMath mp_int for an unsigned long.315 */316inline mp_int isl_sioimath_uiarg_src(unsigned long arg,317 isl_sioimath_scratchspace_t *scratch)318{319 scratch->big.digits = scratch->digits;320 scratch->big.alloc = ARRAY_SIZE(scratch->digits);321 scratch->big.sign = MP_ZPOS;322 323 isl_siomath_ulong_to_digits(arg, scratch->digits, &scratch->big.used);324 return &scratch->big;325}326 327/* Ensure big representation. Does not preserve the current number.328 * Callers may use the fact that the value _is_ preserved if the presentation329 * was big before.330 */331inline mp_int isl_sioimath_reinit_big(isl_sioimath_ptr ptr)332{333 if (isl_sioimath_is_small(*ptr))334 *ptr = isl_sioimath_encode_big(mp_int_alloc());335 return isl_sioimath_get_big(*ptr);336}337 338/* Set ptr to a number in small representation.339 */340inline void isl_sioimath_set_small(isl_sioimath_ptr ptr, int32_t val)341{342 if (isl_sioimath_is_big(*ptr))343 mp_int_free(isl_sioimath_get_big(*ptr));344 *ptr = isl_sioimath_encode_small(val);345}346 347/* Set ptr to val, choosing small representation if possible.348 */349inline void isl_sioimath_set_int32(isl_sioimath_ptr ptr, int32_t val)350{351 if (ISL_SIOIMATH_SMALL_MIN <= val && val <= ISL_SIOIMATH_SMALL_MAX) {352 isl_sioimath_set_small(ptr, val);353 return;354 }355 356 mp_int_init_value(isl_sioimath_reinit_big(ptr), val);357}358 359/* Assign an int64_t number using small representation if possible.360 */361inline void isl_sioimath_set_int64(isl_sioimath_ptr ptr, int64_t val)362{363 if (ISL_SIOIMATH_SMALL_MIN <= val && val <= ISL_SIOIMATH_SMALL_MAX) {364 isl_sioimath_set_small(ptr, val);365 return;366 }367 368 isl_sioimath_scratchspace_t scratch;369 mp_int_copy(isl_sioimath_si64arg_src(val, &scratch),370 isl_sioimath_reinit_big(ptr));371}372 373/* Convert to big representation while preserving the current number.374 */375inline void isl_sioimath_promote(isl_sioimath_ptr dst)376{377 int32_t small;378 379 if (isl_sioimath_is_big(*dst))380 return;381 382 small = isl_sioimath_get_small(*dst);383 mp_int_set_value(isl_sioimath_reinit_big(dst), small);384}385 386/* Convert to small representation while preserving the current number. Does387 * nothing if dst doesn't fit small representation.388 */389inline void isl_sioimath_try_demote(isl_sioimath_ptr dst)390{391 mp_small small;392 393 if (isl_sioimath_is_small(*dst))394 return;395 396 if (mp_int_to_int(isl_sioimath_get_big(*dst), &small) != MP_OK)397 return;398 399 if (ISL_SIOIMATH_SMALL_MIN <= small && small <= ISL_SIOIMATH_SMALL_MAX)400 isl_sioimath_set_small(dst, small);401}402 403/* Initialize an isl_int. The implicit value is 0 in small representation.404 */405inline void isl_sioimath_init(isl_sioimath_ptr dst)406{407 *dst = isl_sioimath_encode_small(0);408}409 410/* Free the resources taken by an isl_int.411 */412inline void isl_sioimath_clear(isl_sioimath_ptr dst)413{414 if (isl_sioimath_is_small(*dst))415 return;416 417 mp_int_free(isl_sioimath_get_big(*dst));418}419 420/* Copy the value of one isl_int to another.421 */422inline void isl_sioimath_set(isl_sioimath_ptr dst, isl_sioimath_src val)423{424 if (isl_sioimath_is_small(val)) {425 isl_sioimath_set_small(dst, isl_sioimath_get_small(val));426 return;427 }428 429 mp_int_copy(isl_sioimath_get_big(val), isl_sioimath_reinit_big(dst));430}431 432/* Store a signed long into an isl_int.433 */434inline void isl_sioimath_set_si(isl_sioimath_ptr dst, long val)435{436 if (ISL_SIOIMATH_SMALL_MIN <= val && val <= ISL_SIOIMATH_SMALL_MAX) {437 isl_sioimath_set_small(dst, val);438 return;439 }440 441 mp_int_set_value(isl_sioimath_reinit_big(dst), val);442}443 444/* Store an unsigned long into an isl_int.445 */446inline void isl_sioimath_set_ui(isl_sioimath_ptr dst, unsigned long val)447{448 if (val <= ISL_SIOIMATH_SMALL_MAX) {449 isl_sioimath_set_small(dst, val);450 return;451 }452 453 mp_int_set_uvalue(isl_sioimath_reinit_big(dst), val);454}455 456/* Return whether a number can be represented by a signed long.457 */458inline int isl_sioimath_fits_slong(isl_sioimath_src val)459{460 mp_small dummy;461 462 if (isl_sioimath_is_small(val))463 return 1;464 465 return mp_int_to_int(isl_sioimath_get_big(val), &dummy) == MP_OK;466}467 468/* Return a number as signed long. Result is undefined if the number cannot be469 * represented as long.470 */471inline long isl_sioimath_get_si(isl_sioimath_src val)472{473 mp_small result;474 475 if (isl_sioimath_is_small(val))476 return isl_sioimath_get_small(val);477 478 mp_int_to_int(isl_sioimath_get_big(val), &result);479 return result;480}481 482/* Return whether a number can be represented as unsigned long.483 */484inline int isl_sioimath_fits_ulong(isl_sioimath_src val)485{486 mp_usmall dummy;487 488 if (isl_sioimath_is_small(val))489 return isl_sioimath_get_small(val) >= 0;490 491 return mp_int_to_uint(isl_sioimath_get_big(val), &dummy) == MP_OK;492}493 494/* Return a number as unsigned long. Result is undefined if the number cannot be495 * represented as unsigned long.496 */497inline unsigned long isl_sioimath_get_ui(isl_sioimath_src val)498{499 mp_usmall result;500 501 if (isl_sioimath_is_small(val))502 return isl_sioimath_get_small(val);503 504 mp_int_to_uint(isl_sioimath_get_big(val), &result);505 return result;506}507 508/* Return a number as floating point value.509 */510inline double isl_sioimath_get_d(isl_sioimath_src val)511{512 mp_int big;513 double result = 0;514 int i;515 516 if (isl_sioimath_is_small(val))517 return isl_sioimath_get_small(val);518 519 big = isl_sioimath_get_big(val);520 for (i = 0; i < big->used; ++i)521 result = result * (double) ((uintmax_t) MP_DIGIT_MAX + 1) +522 (double) big->digits[i];523 524 if (big->sign == MP_NEG)525 result = -result;526 527 return result;528}529 530/* Format a number as decimal string.531 *532 * The largest possible string from small representation is 12 characters533 * ("-2147483647").534 */535inline char *isl_sioimath_get_str(isl_sioimath_src val)536{537 char *result;538 539 if (isl_sioimath_is_small(val)) {540 result = malloc(12);541 snprintf(result, 12, "%" PRIi32, isl_sioimath_get_small(val));542 return result;543 }544 545 return impz_get_str(NULL, 10, isl_sioimath_get_big(val));546}547 548/* Return the absolute value.549 */550inline void isl_sioimath_abs(isl_sioimath_ptr dst, isl_sioimath_src arg)551{552 if (isl_sioimath_is_small(arg)) {553 isl_sioimath_set_small(dst, labs(isl_sioimath_get_small(arg)));554 return;555 }556 557 mp_int_abs(isl_sioimath_get_big(arg), isl_sioimath_reinit_big(dst));558}559 560/* Return the negation of a number.561 */562inline void isl_sioimath_neg(isl_sioimath_ptr dst, isl_sioimath_src arg)563{564 if (isl_sioimath_is_small(arg)) {565 isl_sioimath_set_small(dst, -isl_sioimath_get_small(arg));566 return;567 }568 569 mp_int_neg(isl_sioimath_get_big(arg), isl_sioimath_reinit_big(dst));570}571 572/* Swap two isl_ints.573 *574 * isl_sioimath can be copied bytewise; nothing depends on its address. It can575 * also be stored in a CPU register.576 */577inline void isl_sioimath_swap(isl_sioimath_ptr lhs, isl_sioimath_ptr rhs)578{579 isl_sioimath tmp = *lhs;580 *lhs = *rhs;581 *rhs = tmp;582}583 584/* Add an unsigned long to the number.585 *586 * On LP64 unsigned long exceeds the range of an int64_t, therefore we check in587 * advance whether small representation possibly overflows.588 */589inline void isl_sioimath_add_ui(isl_sioimath_ptr dst, isl_sioimath lhs,590 unsigned long rhs)591{592 int32_t smalllhs;593 isl_sioimath_scratchspace_t lhsscratch;594 595 if (isl_sioimath_decode_small(lhs, &smalllhs) &&596 (rhs <= (uint64_t) INT64_MAX - (uint64_t) ISL_SIOIMATH_SMALL_MAX)) {597 isl_sioimath_set_int64(dst, (int64_t) smalllhs + rhs);598 return;599 }600 601 impz_add_ui(isl_sioimath_reinit_big(dst),602 isl_sioimath_bigarg_src(lhs, &lhsscratch), rhs);603 isl_sioimath_try_demote(dst);604}605 606/* Subtract an unsigned long.607 *608 * On LP64 unsigned long exceeds the range of an int64_t. If609 * ISL_SIOIMATH_SMALL_MIN-rhs>=INT64_MIN we can do the calculation using int64_t610 * without risking an overflow.611 */612inline void isl_sioimath_sub_ui(isl_sioimath_ptr dst, isl_sioimath lhs,613 unsigned long rhs)614{615 int32_t smalllhs;616 isl_sioimath_scratchspace_t lhsscratch;617 618 if (isl_sioimath_decode_small(lhs, &smalllhs) &&619 (rhs < (uint64_t) INT64_MIN - (uint64_t) ISL_SIOIMATH_SMALL_MIN)) {620 isl_sioimath_set_int64(dst, (int64_t) smalllhs - rhs);621 return;622 }623 624 impz_sub_ui(isl_sioimath_reinit_big(dst),625 isl_sioimath_bigarg_src(lhs, &lhsscratch), rhs);626 isl_sioimath_try_demote(dst);627}628 629/* Sum of two isl_ints.630 */631inline void isl_sioimath_add(isl_sioimath_ptr dst, isl_sioimath_src lhs,632 isl_sioimath_src rhs)633{634 isl_sioimath_scratchspace_t scratchlhs, scratchrhs;635 int32_t smalllhs, smallrhs;636 637 if (isl_sioimath_decode_small(lhs, &smalllhs) &&638 isl_sioimath_decode_small(rhs, &smallrhs)) {639 isl_sioimath_set_int64(640 dst, (int64_t) smalllhs + (int64_t) smallrhs);641 return;642 }643 644 mp_int_add(isl_sioimath_bigarg_src(lhs, &scratchlhs),645 isl_sioimath_bigarg_src(rhs, &scratchrhs),646 isl_sioimath_reinit_big(dst));647 isl_sioimath_try_demote(dst);648}649 650/* Subtract two isl_ints.651 */652inline void isl_sioimath_sub(isl_sioimath_ptr dst, isl_sioimath_src lhs,653 isl_sioimath_src rhs)654{655 isl_sioimath_scratchspace_t scratchlhs, scratchrhs;656 int32_t smalllhs, smallrhs;657 658 if (isl_sioimath_decode_small(lhs, &smalllhs) &&659 isl_sioimath_decode_small(rhs, &smallrhs)) {660 isl_sioimath_set_int64(661 dst, (int64_t) smalllhs - (int64_t) smallrhs);662 return;663 }664 665 mp_int_sub(isl_sioimath_bigarg_src(lhs, &scratchlhs),666 isl_sioimath_bigarg_src(rhs, &scratchrhs),667 isl_sioimath_reinit_big(dst));668 isl_sioimath_try_demote(dst);669}670 671/* Multiply two isl_ints.672 */673inline void isl_sioimath_mul(isl_sioimath_ptr dst, isl_sioimath_src lhs,674 isl_sioimath_src rhs)675{676 isl_sioimath_scratchspace_t scratchlhs, scratchrhs;677 int32_t smalllhs, smallrhs;678 679 if (isl_sioimath_decode_small(lhs, &smalllhs) &&680 isl_sioimath_decode_small(rhs, &smallrhs)) {681 isl_sioimath_set_int64(682 dst, (int64_t) smalllhs * (int64_t) smallrhs);683 return;684 }685 686 mp_int_mul(isl_sioimath_bigarg_src(lhs, &scratchlhs),687 isl_sioimath_bigarg_src(rhs, &scratchrhs),688 isl_sioimath_reinit_big(dst));689 isl_sioimath_try_demote(dst);690}691 692/* Shift lhs by rhs bits to the left and store the result in dst. Effectively,693 * this operation computes 'lhs * 2^rhs'.694 */695inline void isl_sioimath_mul_2exp(isl_sioimath_ptr dst, isl_sioimath lhs,696 unsigned long rhs)697{698 isl_sioimath_scratchspace_t scratchlhs;699 int32_t smalllhs;700 701 if (isl_sioimath_decode_small(lhs, &smalllhs) && (rhs <= 32ul)) {702 isl_sioimath_set_int64(dst, ((int64_t) smalllhs) << rhs);703 return;704 }705 706 mp_int_mul_pow2(isl_sioimath_bigarg_src(lhs, &scratchlhs), rhs,707 isl_sioimath_reinit_big(dst));708}709 710/* Multiply an isl_int and a signed long.711 */712inline void isl_sioimath_mul_si(isl_sioimath_ptr dst, isl_sioimath lhs,713 signed long rhs)714{715 isl_sioimath_scratchspace_t scratchlhs, scratchrhs;716 int32_t smalllhs;717 718 if (isl_sioimath_decode_small(lhs, &smalllhs) && (rhs > LONG_MIN) &&719 (labs(rhs) <= UINT32_MAX)) {720 isl_sioimath_set_int64(dst, (int64_t) smalllhs * (int64_t) rhs);721 return;722 }723 724 mp_int_mul(isl_sioimath_bigarg_src(lhs, &scratchlhs),725 isl_sioimath_siarg_src(rhs, &scratchrhs),726 isl_sioimath_reinit_big(dst));727 isl_sioimath_try_demote(dst);728}729 730/* Multiply an isl_int and an unsigned long.731 */732inline void isl_sioimath_mul_ui(isl_sioimath_ptr dst, isl_sioimath lhs,733 unsigned long rhs)734{735 isl_sioimath_scratchspace_t scratchlhs, scratchrhs;736 int32_t smalllhs;737 738 if (isl_sioimath_decode_small(lhs, &smalllhs) && (rhs <= UINT32_MAX)) {739 isl_sioimath_set_int64(dst, (int64_t) smalllhs * (int64_t) rhs);740 return;741 }742 743 mp_int_mul(isl_sioimath_bigarg_src(lhs, &scratchlhs),744 isl_sioimath_uiarg_src(rhs, &scratchrhs),745 isl_sioimath_reinit_big(dst));746 isl_sioimath_try_demote(dst);747}748 749/* Compute the power of an isl_int to an unsigned long.750 * Always let IMath do it; the result is unlikely to be small except in some751 * special cases.752 * Note: 0^0 == 1753 */754inline void isl_sioimath_pow_ui(isl_sioimath_ptr dst, isl_sioimath_src lhs,755 unsigned long rhs)756{757 isl_sioimath_scratchspace_t scratchlhs, scratchrhs;758 int32_t smalllhs;759 760 switch (rhs) {761 case 0:762 isl_sioimath_set_small(dst, 1);763 return;764 case 1:765 isl_sioimath_set(dst, lhs);766 return;767 case 2:768 isl_sioimath_mul(dst, lhs, lhs);769 return;770 }771 772 if (isl_sioimath_decode_small(lhs, &smalllhs)) {773 switch (smalllhs) {774 case 0:775 isl_sioimath_set_small(dst, 0);776 return;777 case 1:778 isl_sioimath_set_small(dst, 1);779 return;780 case 2:781 isl_sioimath_set_small(dst, 1);782 isl_sioimath_mul_2exp(dst, *dst, rhs);783 return;784 default:785 if ((MP_SMALL_MIN <= rhs) && (rhs <= MP_SMALL_MAX)) {786 mp_int_expt_value(smalllhs, rhs,787 isl_sioimath_reinit_big(dst));788 isl_sioimath_try_demote(dst);789 return;790 }791 }792 }793 794 mp_int_expt_full(isl_sioimath_bigarg_src(lhs, &scratchlhs),795 isl_sioimath_uiarg_src(rhs, &scratchrhs),796 isl_sioimath_reinit_big(dst));797 isl_sioimath_try_demote(dst);798}799 800/* Fused multiply-add.801 */802inline void isl_sioimath_addmul(isl_sioimath_ptr dst, isl_sioimath_src lhs,803 isl_sioimath_src rhs)804{805 isl_sioimath tmp;806 isl_sioimath_init(&tmp);807 isl_sioimath_mul(&tmp, lhs, rhs);808 isl_sioimath_add(dst, *dst, tmp);809 isl_sioimath_clear(&tmp);810}811 812/* Fused multiply-add with an unsigned long.813 */814inline void isl_sioimath_addmul_ui(isl_sioimath_ptr dst, isl_sioimath_src lhs,815 unsigned long rhs)816{817 isl_sioimath tmp;818 isl_sioimath_init(&tmp);819 isl_sioimath_mul_ui(&tmp, lhs, rhs);820 isl_sioimath_add(dst, *dst, tmp);821 isl_sioimath_clear(&tmp);822}823 824/* Fused multiply-subtract.825 */826inline void isl_sioimath_submul(isl_sioimath_ptr dst, isl_sioimath_src lhs,827 isl_sioimath_src rhs)828{829 isl_sioimath tmp;830 isl_sioimath_init(&tmp);831 isl_sioimath_mul(&tmp, lhs, rhs);832 isl_sioimath_sub(dst, *dst, tmp);833 isl_sioimath_clear(&tmp);834}835 836/* Fused multiply-add with an unsigned long.837 */838inline void isl_sioimath_submul_ui(isl_sioimath_ptr dst, isl_sioimath_src lhs,839 unsigned long rhs)840{841 isl_sioimath tmp;842 isl_sioimath_init(&tmp);843 isl_sioimath_mul_ui(&tmp, lhs, rhs);844 isl_sioimath_sub(dst, *dst, tmp);845 isl_sioimath_clear(&tmp);846}847 848void isl_sioimath_gcd(isl_sioimath_ptr dst, isl_sioimath_src lhs,849 isl_sioimath_src rhs);850void isl_sioimath_lcm(isl_sioimath_ptr dst, isl_sioimath_src lhs,851 isl_sioimath_src rhs);852 853/* Divide lhs by rhs, rounding to zero (Truncate).854 */855inline void isl_sioimath_tdiv_q(isl_sioimath_ptr dst, isl_sioimath_src lhs,856 isl_sioimath_src rhs)857{858 isl_sioimath_scratchspace_t lhsscratch, rhsscratch;859 int32_t lhssmall, rhssmall;860 861 if (isl_sioimath_decode_small(lhs, &lhssmall) &&862 isl_sioimath_decode_small(rhs, &rhssmall)) {863 isl_sioimath_set_small(dst, lhssmall / rhssmall);864 return;865 }866 867 mp_int_div(isl_sioimath_bigarg_src(lhs, &lhsscratch),868 isl_sioimath_bigarg_src(rhs, &rhsscratch),869 isl_sioimath_reinit_big(dst), NULL);870 isl_sioimath_try_demote(dst);871 return;872}873 874/* Divide lhs by an unsigned long rhs, rounding to zero (Truncate).875 */876inline void isl_sioimath_tdiv_q_ui(isl_sioimath_ptr dst, isl_sioimath_src lhs,877 unsigned long rhs)878{879 isl_sioimath_scratchspace_t lhsscratch, rhsscratch;880 int32_t lhssmall;881 882 if (isl_sioimath_is_small(lhs) && (rhs <= (unsigned long) INT32_MAX)) {883 lhssmall = isl_sioimath_get_small(lhs);884 isl_sioimath_set_small(dst, lhssmall / (int32_t) rhs);885 return;886 }887 888 if (rhs <= MP_SMALL_MAX) {889 mp_int_div_value(isl_sioimath_bigarg_src(lhs, &lhsscratch), rhs,890 isl_sioimath_reinit_big(dst), NULL);891 isl_sioimath_try_demote(dst);892 return;893 }894 895 mp_int_div(isl_sioimath_bigarg_src(lhs, &lhsscratch),896 isl_sioimath_uiarg_src(rhs, &rhsscratch),897 isl_sioimath_reinit_big(dst), NULL);898 isl_sioimath_try_demote(dst);899}900 901/* Divide lhs by rhs, rounding to positive infinity (Ceil).902 */903inline void isl_sioimath_cdiv_q(isl_sioimath_ptr dst, isl_sioimath_src lhs,904 isl_sioimath_src rhs)905{906 int32_t lhssmall, rhssmall;907 isl_sioimath_scratchspace_t lhsscratch, rhsscratch;908 int32_t q;909 910 if (isl_sioimath_decode_small(lhs, &lhssmall) &&911 isl_sioimath_decode_small(rhs, &rhssmall)) {912 if ((lhssmall >= 0) && (rhssmall >= 0))913 q = ((int64_t) lhssmall + (int64_t) rhssmall - 1) /914 rhssmall;915 else if ((lhssmall < 0) && (rhssmall < 0))916 q = ((int64_t) lhssmall + (int64_t) rhssmall + 1) /917 rhssmall;918 else919 q = lhssmall / rhssmall;920 isl_sioimath_set_small(dst, q);921 return;922 }923 924 impz_cdiv_q(isl_sioimath_reinit_big(dst),925 isl_sioimath_bigarg_src(lhs, &lhsscratch),926 isl_sioimath_bigarg_src(rhs, &rhsscratch));927 isl_sioimath_try_demote(dst);928}929 930/* Compute the division of lhs by a rhs of type unsigned long, rounding towards931 * positive infinity (Ceil).932 */933inline void isl_sioimath_cdiv_q_ui(isl_sioimath_ptr dst, isl_sioimath_src lhs,934 unsigned long rhs)935{936 isl_sioimath_scratchspace_t lhsscratch, rhsscratch;937 int32_t lhssmall, q;938 939 if (isl_sioimath_decode_small(lhs, &lhssmall) && (rhs <= INT32_MAX)) {940 if (lhssmall >= 0)941 q = ((int64_t) lhssmall + ((int64_t) rhs - 1)) /942 (int64_t) rhs;943 else944 q = lhssmall / (int32_t) rhs;945 isl_sioimath_set_small(dst, q);946 return;947 }948 949 impz_cdiv_q(isl_sioimath_reinit_big(dst),950 isl_sioimath_bigarg_src(lhs, &lhsscratch),951 isl_sioimath_uiarg_src(rhs, &rhsscratch));952 isl_sioimath_try_demote(dst);953}954 955/* Divide lhs by rhs, rounding to negative infinity (Floor).956 */957inline void isl_sioimath_fdiv_q(isl_sioimath_ptr dst, isl_sioimath_src lhs,958 isl_sioimath_src rhs)959{960 isl_sioimath_scratchspace_t lhsscratch, rhsscratch;961 int32_t lhssmall, rhssmall;962 int32_t q;963 964 if (isl_sioimath_decode_small(lhs, &lhssmall) &&965 isl_sioimath_decode_small(rhs, &rhssmall)) {966 if ((lhssmall < 0) && (rhssmall >= 0))967 q = ((int64_t) lhssmall - ((int64_t) rhssmall - 1)) /968 rhssmall;969 else if ((lhssmall >= 0) && (rhssmall < 0))970 q = ((int64_t) lhssmall - ((int64_t) rhssmall + 1)) /971 rhssmall;972 else973 q = lhssmall / rhssmall;974 isl_sioimath_set_small(dst, q);975 return;976 }977 978 impz_fdiv_q(isl_sioimath_reinit_big(dst),979 isl_sioimath_bigarg_src(lhs, &lhsscratch),980 isl_sioimath_bigarg_src(rhs, &rhsscratch));981 isl_sioimath_try_demote(dst);982}983 984/* Compute the division of lhs by a rhs of type unsigned long, rounding towards985 * negative infinity (Floor).986 */987inline void isl_sioimath_fdiv_q_ui(isl_sioimath_ptr dst, isl_sioimath_src lhs,988 unsigned long rhs)989{990 isl_sioimath_scratchspace_t lhsscratch, rhsscratch;991 int32_t lhssmall, q;992 993 if (isl_sioimath_decode_small(lhs, &lhssmall) && (rhs <= INT32_MAX)) {994 if (lhssmall >= 0)995 q = (uint32_t) lhssmall / rhs;996 else997 q = ((int64_t) lhssmall - ((int64_t) rhs - 1)) /998 (int64_t) rhs;999 isl_sioimath_set_small(dst, q);1000 return;1001 }1002 1003 impz_fdiv_q(isl_sioimath_reinit_big(dst),1004 isl_sioimath_bigarg_src(lhs, &lhsscratch),1005 isl_sioimath_uiarg_src(rhs, &rhsscratch));1006 isl_sioimath_try_demote(dst);1007}1008 1009/* Get the remainder of: lhs divided by rhs rounded towards negative infinite1010 * (Floor).1011 */1012inline void isl_sioimath_fdiv_r(isl_sioimath_ptr dst, isl_sioimath_src lhs,1013 isl_sioimath_src rhs)1014{1015 isl_sioimath_scratchspace_t lhsscratch, rhsscratch;1016 int64_t lhssmall, rhssmall;1017 int32_t r;1018 1019 if (isl_sioimath_is_small(lhs) && isl_sioimath_is_small(rhs)) {1020 lhssmall = isl_sioimath_get_small(lhs);1021 rhssmall = isl_sioimath_get_small(rhs);1022 r = (rhssmall + lhssmall % rhssmall) % rhssmall;1023 isl_sioimath_set_small(dst, r);1024 return;1025 }1026 1027 impz_fdiv_r(isl_sioimath_reinit_big(dst),1028 isl_sioimath_bigarg_src(lhs, &lhsscratch),1029 isl_sioimath_bigarg_src(rhs, &rhsscratch));1030 isl_sioimath_try_demote(dst);1031}1032 1033void isl_sioimath_read(isl_sioimath_ptr dst, const char *str);1034 1035/* Return:1036 * +1 for a positive number1037 * -1 for a negative number1038 * 0 if the number is zero1039 */1040inline int isl_sioimath_sgn(isl_sioimath_src arg)1041{1042 int32_t small;1043 1044 if (isl_sioimath_decode_small(arg, &small))1045 return (small > 0) - (small < 0);1046 1047 return mp_int_compare_zero(isl_sioimath_get_big(arg));1048}1049 1050/* Return:1051 * +1 if lhs > rhs1052 * -1 if lhs < rhs1053 * 0 if lhs = rhs1054 */1055inline int isl_sioimath_cmp(isl_sioimath_src lhs, isl_sioimath_src rhs)1056{1057 isl_sioimath_scratchspace_t lhsscratch, rhsscratch;1058 int32_t lhssmall, rhssmall;1059 1060 if (isl_sioimath_decode_small(lhs, &lhssmall) &&1061 isl_sioimath_decode_small(rhs, &rhssmall))1062 return (lhssmall > rhssmall) - (lhssmall < rhssmall);1063 1064 if (isl_sioimath_decode_small(rhs, &rhssmall))1065 return mp_int_compare_value(1066 isl_sioimath_bigarg_src(lhs, &lhsscratch), rhssmall);1067 1068 if (isl_sioimath_decode_small(lhs, &lhssmall))1069 return -mp_int_compare_value(1070 isl_sioimath_bigarg_src(rhs, &rhsscratch), lhssmall);1071 1072 return mp_int_compare(1073 isl_sioimath_get_big(lhs), isl_sioimath_get_big(rhs));1074}1075 1076/* As isl_sioimath_cmp, but with signed long rhs.1077 */1078inline int isl_sioimath_cmp_si(isl_sioimath_src lhs, signed long rhs)1079{1080 int32_t lhssmall;1081 1082 if (isl_sioimath_decode_small(lhs, &lhssmall))1083 return (lhssmall > rhs) - (lhssmall < rhs);1084 1085 return mp_int_compare_value(isl_sioimath_get_big(lhs), rhs);1086}1087 1088/* Return:1089 * +1 if |lhs| > |rhs|1090 * -1 if |lhs| < |rhs|1091 * 0 if |lhs| = |rhs|1092 */1093inline int isl_sioimath_abs_cmp(isl_sioimath_src lhs, isl_sioimath_src rhs)1094{1095 isl_sioimath_scratchspace_t lhsscratch, rhsscratch;1096 int32_t lhssmall, rhssmall;1097 1098 if (isl_sioimath_decode_small(lhs, &lhssmall) &&1099 isl_sioimath_decode_small(rhs, &rhssmall)) {1100 lhssmall = labs(lhssmall);1101 rhssmall = labs(rhssmall);1102 return (lhssmall > rhssmall) - (lhssmall < rhssmall);1103 }1104 1105 return mp_int_compare_unsigned(1106 isl_sioimath_bigarg_src(lhs, &lhsscratch),1107 isl_sioimath_bigarg_src(rhs, &rhsscratch));1108}1109 1110/* Return whether lhs is divisible by rhs.1111 * In particular, can rhs be multiplied by some integer to result in lhs?1112 * If rhs is zero, then this means lhs has to be zero too.1113 */1114inline int isl_sioimath_is_divisible_by(isl_sioimath_src lhs,1115 isl_sioimath_src rhs)1116{1117 isl_sioimath_scratchspace_t lhsscratch, rhsscratch;1118 int32_t lhssmall, rhssmall;1119 mpz_t rem;1120 int cmp;1121 1122 if (isl_sioimath_sgn(rhs) == 0)1123 return isl_sioimath_sgn(lhs) == 0;1124 1125 if (isl_sioimath_decode_small(lhs, &lhssmall) &&1126 isl_sioimath_decode_small(rhs, &rhssmall))1127 return lhssmall % rhssmall == 0;1128 1129 if (isl_sioimath_decode_small(rhs, &rhssmall))1130 return mp_int_divisible_value(1131 isl_sioimath_bigarg_src(lhs, &lhsscratch), rhssmall);1132 1133 mp_int_init(&rem);1134 mp_int_div(isl_sioimath_bigarg_src(lhs, &lhsscratch),1135 isl_sioimath_bigarg_src(rhs, &rhsscratch), NULL, &rem);1136 cmp = mp_int_compare_zero(&rem);1137 mp_int_clear(&rem);1138 return cmp == 0;1139}1140 1141/* Return a hash code of an isl_sioimath.1142 * The hash code for a number in small and big representation must be identical1143 * on the same machine because small representation if not obligatory if fits.1144 */1145inline uint32_t isl_sioimath_hash(isl_sioimath_src arg, uint32_t hash)1146{1147 int32_t small;1148 int i;1149 uint32_t num;1150 mp_digit digits[(sizeof(uint32_t) + sizeof(mp_digit) - 1) /1151 sizeof(mp_digit)];1152 mp_size used;1153 const unsigned char *digitdata = (const unsigned char *) &digits;1154 1155 if (isl_sioimath_decode_small(arg, &small)) {1156 if (small < 0)1157 isl_hash_byte(hash, 0xFF);1158 num = labs(small);1159 1160 isl_siomath_uint32_to_digits(num, digits, &used);1161 for (i = 0; i < used * sizeof(mp_digit); i += 1)1162 isl_hash_byte(hash, digitdata[i]);1163 return hash;1164 }1165 1166 return isl_imath_hash(isl_sioimath_get_big(arg), hash);1167}1168 1169/* Return the number of digits in a number of the given base or more, i.e. the1170 * string length without sign and null terminator.1171 *1172 * Current implementation for small representation returns the maximal number1173 * of binary digits in that representation, which can be much larger than the1174 * smallest possible solution.1175 */1176inline size_t isl_sioimath_sizeinbase(isl_sioimath_src arg, int base)1177{1178 int32_t small;1179 1180 if (isl_sioimath_decode_small(arg, &small))1181 return sizeof(int32_t) * CHAR_BIT - 1;1182 1183 return impz_sizeinbase(isl_sioimath_get_big(arg), base);1184}1185 1186void isl_sioimath_print(FILE *out, isl_sioimath_src i, int width);1187void isl_sioimath_dump(isl_sioimath_src arg);1188 1189typedef isl_sioimath isl_int[1];1190#define isl_int_init(i) isl_sioimath_init((i))1191#define isl_int_clear(i) isl_sioimath_clear((i))1192 1193#define isl_int_set(r, i) isl_sioimath_set((r), *(i))1194#define isl_int_set_si(r, i) isl_sioimath_set_si((r), i)1195#define isl_int_set_ui(r, i) isl_sioimath_set_ui((r), i)1196#define isl_int_fits_slong(r) isl_sioimath_fits_slong(*(r))1197#define isl_int_get_si(r) isl_sioimath_get_si(*(r))1198#define isl_int_fits_ulong(r) isl_sioimath_fits_ulong(*(r))1199#define isl_int_get_ui(r) isl_sioimath_get_ui(*(r))1200#define isl_int_get_d(r) isl_sioimath_get_d(*(r))1201#define isl_int_get_str(r) isl_sioimath_get_str(*(r))1202#define isl_int_abs(r, i) isl_sioimath_abs((r), *(i))1203#define isl_int_neg(r, i) isl_sioimath_neg((r), *(i))1204#define isl_int_swap(i, j) isl_sioimath_swap((i), (j))1205#define isl_int_swap_or_set(i, j) isl_sioimath_swap((i), (j))1206#define isl_int_add_ui(r, i, j) isl_sioimath_add_ui((r), *(i), j)1207#define isl_int_sub_ui(r, i, j) isl_sioimath_sub_ui((r), *(i), j)1208 1209#define isl_int_add(r, i, j) isl_sioimath_add((r), *(i), *(j))1210#define isl_int_sub(r, i, j) isl_sioimath_sub((r), *(i), *(j))1211#define isl_int_mul(r, i, j) isl_sioimath_mul((r), *(i), *(j))1212#define isl_int_mul_2exp(r, i, j) isl_sioimath_mul_2exp((r), *(i), j)1213#define isl_int_mul_si(r, i, j) isl_sioimath_mul_si((r), *(i), j)1214#define isl_int_mul_ui(r, i, j) isl_sioimath_mul_ui((r), *(i), j)1215#define isl_int_pow_ui(r, i, j) isl_sioimath_pow_ui((r), *(i), j)1216#define isl_int_addmul(r, i, j) isl_sioimath_addmul((r), *(i), *(j))1217#define isl_int_addmul_ui(r, i, j) isl_sioimath_addmul_ui((r), *(i), j)1218#define isl_int_submul(r, i, j) isl_sioimath_submul((r), *(i), *(j))1219#define isl_int_submul_ui(r, i, j) isl_sioimath_submul_ui((r), *(i), j)1220 1221#define isl_int_gcd(r, i, j) isl_sioimath_gcd((r), *(i), *(j))1222#define isl_int_lcm(r, i, j) isl_sioimath_lcm((r), *(i), *(j))1223#define isl_int_divexact(r, i, j) isl_sioimath_tdiv_q((r), *(i), *(j))1224#define isl_int_divexact_ui(r, i, j) isl_sioimath_tdiv_q_ui((r), *(i), j)1225#define isl_int_tdiv_q(r, i, j) isl_sioimath_tdiv_q((r), *(i), *(j))1226#define isl_int_cdiv_q(r, i, j) isl_sioimath_cdiv_q((r), *(i), *(j))1227#define isl_int_cdiv_q_ui(r, i, j) isl_sioimath_cdiv_q_ui((r), *(i), j)1228#define isl_int_fdiv_q(r, i, j) isl_sioimath_fdiv_q((r), *(i), *(j))1229#define isl_int_fdiv_r(r, i, j) isl_sioimath_fdiv_r((r), *(i), *(j))1230#define isl_int_fdiv_q_ui(r, i, j) isl_sioimath_fdiv_q_ui((r), *(i), j)1231 1232#define isl_int_read(r, s) isl_sioimath_read((r), s)1233#define isl_int_sgn(i) isl_sioimath_sgn(*(i))1234#define isl_int_cmp(i, j) isl_sioimath_cmp(*(i), *(j))1235#define isl_int_cmp_si(i, si) isl_sioimath_cmp_si(*(i), si)1236#define isl_int_eq(i, j) (isl_sioimath_cmp(*(i), *(j)) == 0)1237#define isl_int_ne(i, j) (isl_sioimath_cmp(*(i), *(j)) != 0)1238#define isl_int_lt(i, j) (isl_sioimath_cmp(*(i), *(j)) < 0)1239#define isl_int_le(i, j) (isl_sioimath_cmp(*(i), *(j)) <= 0)1240#define isl_int_gt(i, j) (isl_sioimath_cmp(*(i), *(j)) > 0)1241#define isl_int_ge(i, j) (isl_sioimath_cmp(*(i), *(j)) >= 0)1242#define isl_int_abs_cmp(i, j) isl_sioimath_abs_cmp(*(i), *(j))1243#define isl_int_abs_eq(i, j) (isl_sioimath_abs_cmp(*(i), *(j)) == 0)1244#define isl_int_abs_ne(i, j) (isl_sioimath_abs_cmp(*(i), *(j)) != 0)1245#define isl_int_abs_lt(i, j) (isl_sioimath_abs_cmp(*(i), *(j)) < 0)1246#define isl_int_abs_gt(i, j) (isl_sioimath_abs_cmp(*(i), *(j)) > 0)1247#define isl_int_abs_ge(i, j) (isl_sioimath_abs_cmp(*(i), *(j)) >= 0)1248#define isl_int_is_divisible_by(i, j) isl_sioimath_is_divisible_by(*(i), *(j))1249 1250#define isl_int_hash(v, h) isl_sioimath_hash(*(v), h)1251#define isl_int_free_str(s) free(s)1252#define isl_int_print(out, i, width) isl_sioimath_print(out, *(i), width)1253 1254#endif /* ISL_INT_SIOIMATH_H */1255