polly/lib/External/isl/imath/gmp_compat.c

   1 /*
   2   Name:     gmp_compat.c
   3   Purpose:  Provide GMP compatiable routines for imath library
   4   Author:   David Peixotto
   5
   6   Copyright (c) 2012 Qualcomm Innovation Center, Inc. All rights reserved.
   7
   8   Permission is hereby granted, free of charge, to any person obtaining a copy
   9   of this software and associated documentation files (the "Software"), to deal
  10   in the Software without restriction, including without limitation the rights
  11   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  12   copies of the Software, and to permit persons to whom the Software is
  13   furnished to do so, subject to the following conditions:
  14
  15   The above copyright notice and this permission notice shall be included in
  16   all copies or substantial portions of the Software.
  17
  18   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  19   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  20   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
  21   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  22   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  23   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  24   SOFTWARE.
  25  */
  26 #include "gmp_compat.h"
  27 #include <assert.h>
  28 #include <ctype.h>
  29 #include <stdio.h>
  30 #include <stdlib.h>
  31 #include <string.h>
  32
  33 #if defined(_MSC_VER)
  34 #include <BaseTsd.h>
  35 typedef SSIZE_T ssize_t;
  36 #else
  37 #include <sys/types.h>
  38 #endif
  39
  40 #ifdef NDEBUG
  41 #define CHECK(res) (res)
  42 #else
  43 #define CHECK(res) assert(((res) == MP_OK) && "expected MP_OK")
  44 #endif
  45
  46 /* *(signed char *)&endian_test will thus either be:
  47  *     0b00000001 =  1 on big-endian
  48  *     0b11111111 = -1 on little-endian */
  49 static const uint16_t endian_test = 0x1FF;
  50 #define HOST_ENDIAN (*(signed char *)&endian_test)
  51
  52 /*************************************************************************
  53  *
  54  * Functions with direct translations
  55  *
  56  *************************************************************************/
  57 /* gmp: mpq_clear */
  58 void GMPQAPI(clear)(mp_rat x) { mp_rat_clear(x); }
  59
  60 /* gmp: mpq_cmp */
  61 int GMPQAPI(cmp)(mp_rat op1, mp_rat op2) { return mp_rat_compare(op1, op2); }
  62
  63 /* gmp: mpq_init */
  64 void GMPQAPI(init)(mp_rat x) { CHECK(mp_rat_init(x)); }
  65
  66 /* gmp: mpq_mul */
  67 void GMPQAPI(mul)(mp_rat product, mp_rat multiplier, mp_rat multiplicand) {
  68   CHECK(mp_rat_mul(multiplier, multiplicand, product));
  69 }
  70
  71 /* gmp: mpq_set */
  72 void GMPQAPI(set)(mp_rat rop, mp_rat op) { CHECK(mp_rat_copy(op, rop)); }
  73
  74 /* gmp: mpz_abs */
  75 void GMPZAPI(abs)(mp_int rop, mp_int op) { CHECK(mp_int_abs(op, rop)); }
  76
  77 /* gmp: mpz_add */
  78 void GMPZAPI(add)(mp_int rop, mp_int op1, mp_int op2) {
  79   CHECK(mp_int_add(op1, op2, rop));
  80 }
  81
  82 /* gmp: mpz_clear */
  83 void GMPZAPI(clear)(mp_int x) { mp_int_clear(x); }
  84
  85 /* gmp: mpz_cmp_si */
  86 int GMPZAPI(cmp_si)(mp_int op1, long op2) {
  87   return mp_int_compare_value(op1, op2);
  88 }
  89
  90 /* gmp: mpz_cmpabs */
  91 int GMPZAPI(cmpabs)(mp_int op1, mp_int op2) {
  92   return mp_int_compare_unsigned(op1, op2);
  93 }
  94
  95 /* gmp: mpz_cmp */
  96 int GMPZAPI(cmp)(mp_int op1, mp_int op2) { return mp_int_compare(op1, op2); }
  97
  98 /* gmp: mpz_init */
  99 void GMPZAPI(init)(mp_int x) { CHECK(mp_int_init(x)); }
 100
 101 /* gmp: mpz_mul */
 102 void GMPZAPI(mul)(mp_int rop, mp_int op1, mp_int op2) {
 103   CHECK(mp_int_mul(op1, op2, rop));
 104 }
 105
 106 /* gmp: mpz_neg */
 107 void GMPZAPI(neg)(mp_int rop, mp_int op) { CHECK(mp_int_neg(op, rop)); }
 108
 109 /* gmp: mpz_set_si */
 110 void GMPZAPI(set_si)(mp_int rop, long op) { CHECK(mp_int_set_value(rop, op)); }
 111
 112 /* gmp: mpz_set */
 113 void GMPZAPI(set)(mp_int rop, mp_int op) { CHECK(mp_int_copy(op, rop)); }
 114
 115 /* gmp: mpz_sub */
 116 void GMPZAPI(sub)(mp_int rop, mp_int op1, mp_int op2) {
 117   CHECK(mp_int_sub(op1, op2, rop));
 118 }
 119
 120 /* gmp: mpz_swap */
 121 void GMPZAPI(swap)(mp_int rop1, mp_int rop2) { mp_int_swap(rop1, rop2); }
 122
 123 /* gmp: mpq_sgn */
 124 int GMPQAPI(sgn)(mp_rat op) { return mp_rat_compare_zero(op); }
 125
 126 /* gmp: mpz_sgn */
 127 int GMPZAPI(sgn)(mp_int op) { return mp_int_compare_zero(op); }
 128
 129 /* gmp: mpq_set_ui */
 130 void GMPQAPI(set_ui)(mp_rat rop, unsigned long op1, unsigned long op2) {
 131   CHECK(mp_rat_set_uvalue(rop, op1, op2));
 132 }
 133
 134 /* gmp: mpz_set_ui */
 135 void GMPZAPI(set_ui)(mp_int rop, unsigned long op) {
 136   CHECK(mp_int_set_uvalue(rop, op));
 137 }
 138
 139 /* gmp: mpq_den_ref */
 140 mp_int GMPQAPI(denref)(mp_rat op) { return mp_rat_denom_ref(op); }
 141
 142 /* gmp: mpq_num_ref */
 143 mp_int GMPQAPI(numref)(mp_rat op) { return mp_rat_numer_ref(op); }
 144
 145 /* gmp: mpq_canonicalize */
 146 void GMPQAPI(canonicalize)(mp_rat op) { CHECK(mp_rat_reduce(op)); }
 147
 148 /*
 149  * Functions that can be implemented as a combination of imath functions
 150  */
 151
 152 /* gmp: mpz_addmul */
 153 /* gmp: rop = rop + (op1 * op2) */
 154 void GMPZAPI(addmul)(mp_int rop, mp_int op1, mp_int op2) {
 155   mpz_t tempz;
 156   mp_int temp = &tempz;
 157   mp_int_init(temp);
 158
 159   CHECK(mp_int_mul(op1, op2, temp));
 160   CHECK(mp_int_add(rop, temp, rop));
 161   mp_int_clear(temp);
 162 }
 163
 164 /* gmp: mpz_divexact */
 165 /* gmp: only produces correct results when d divides n */
 166 void GMPZAPI(divexact)(mp_int q, mp_int n, mp_int d) {
 167   CHECK(mp_int_div(n, d, q, NULL));
 168 }
 169
 170 /* gmp: mpz_divisible_p */
 171 /* gmp: return 1 if d divides n, 0 otherwise */
 172 /* gmp: 0 is considered to divide only 0 */
 173 int GMPZAPI(divisible_p)(mp_int n, mp_int d) {
 174   /* variables to hold remainder */
 175   mpz_t rz;
 176   mp_int r = &rz;
 177   int r_is_zero;
 178
 179   /* check for d = 0 */
 180   int n_is_zero = mp_int_compare_zero(n) == 0;
 181   int d_is_zero = mp_int_compare_zero(d) == 0;
 182   if (d_is_zero) return n_is_zero;
 183
 184   /* return true if remainder is 0 */
 185   CHECK(mp_int_init(r));
 186   CHECK(mp_int_div(n, d, NULL, r));
 187   r_is_zero = mp_int_compare_zero(r) == 0;
 188   mp_int_clear(r);
 189
 190   return r_is_zero;
 191 }
 192
 193 /* gmp: mpz_submul */
 194 /* gmp: rop = rop - (op1 * op2) */
 195 void GMPZAPI(submul)(mp_int rop, mp_int op1, mp_int op2) {
 196   mpz_t tempz;
 197   mp_int temp = &tempz;
 198   mp_int_init(temp);
 199
 200   CHECK(mp_int_mul(op1, op2, temp));
 201   CHECK(mp_int_sub(rop, temp, rop));
 202
 203   mp_int_clear(temp);
 204 }
 205
 206 /* gmp: mpz_add_ui */
 207 void GMPZAPI(add_ui)(mp_int rop, mp_int op1, unsigned long op2) {
 208   mpz_t tempz;
 209   mp_int temp = &tempz;
 210   CHECK(mp_int_init_uvalue(temp, op2));
 211
 212   CHECK(mp_int_add(op1, temp, rop));
 213
 214   mp_int_clear(temp);
 215 }
 216
 217 /* gmp: mpz_divexact_ui */
 218 /* gmp: only produces correct results when d divides n */
 219 void GMPZAPI(divexact_ui)(mp_int q, mp_int n, unsigned long d) {
 220   mpz_t tempz;
 221   mp_int temp = &tempz;
 222   CHECK(mp_int_init_uvalue(temp, d));
 223
 224   CHECK(mp_int_div(n, temp, q, NULL));
 225
 226   mp_int_clear(temp);
 227 }
 228
 229 /* gmp: mpz_mul_ui */
 230 void GMPZAPI(mul_ui)(mp_int rop, mp_int op1, unsigned long op2) {
 231   mpz_t tempz;
 232   mp_int temp = &tempz;
 233   CHECK(mp_int_init_uvalue(temp, op2));
 234
 235   CHECK(mp_int_mul(op1, temp, rop));
 236
 237   mp_int_clear(temp);
 238 }
 239
 240 /* gmp: mpz_pow_ui */
 241 /* gmp: 0^0 = 1 */
 242 void GMPZAPI(pow_ui)(mp_int rop, mp_int base, unsigned long exp) {
 243   mpz_t tempz;
 244   mp_int temp = &tempz;
 245
 246   /* check for 0^0 */
 247   if (exp == 0 && mp_int_compare_zero(base) == 0) {
 248     CHECK(mp_int_set_value(rop, 1));
 249     return;
 250   }
 251
 252   /* rop = base^exp */
 253   CHECK(mp_int_init_uvalue(temp, exp));
 254   CHECK(mp_int_expt_full(base, temp, rop));
 255   mp_int_clear(temp);
 256 }
 257
 258 /* gmp: mpz_sub_ui */
 259 void GMPZAPI(sub_ui)(mp_int rop, mp_int op1, unsigned long op2) {
 260   mpz_t tempz;
 261   mp_int temp = &tempz;
 262   CHECK(mp_int_init_uvalue(temp, op2));
 263
 264   CHECK(mp_int_sub(op1, temp, rop));
 265
 266   mp_int_clear(temp);
 267 }
 268
 269 /*************************************************************************
 270  *
 271  * Functions with different behavior in corner cases
 272  *
 273  *************************************************************************/
 274
 275 /* gmp: mpz_gcd */
 276 void GMPZAPI(gcd)(mp_int rop, mp_int op1, mp_int op2) {
 277   int op1_is_zero = mp_int_compare_zero(op1) == 0;
 278   int op2_is_zero = mp_int_compare_zero(op2) == 0;
 279
 280   if (op1_is_zero && op2_is_zero) {
 281     mp_int_zero(rop);
 282     return;
 283   }
 284
 285   CHECK(mp_int_gcd(op1, op2, rop));
 286 }
 287
 288 /* gmp: mpz_get_str */
 289 char *GMPZAPI(get_str)(char *str, int radix, mp_int op) {
 290   int i, r, len;
 291
 292   /* Support negative radix like gmp */
 293   r = radix;
 294   if (r < 0) r = -r;
 295
 296   /* Compute the length of the string needed to hold the int */
 297   len = mp_int_string_len(op, r);
 298   if (str == NULL) {
 299     str = malloc(len);
 300   }
 301
 302   /* Convert to string using imath function */
 303   CHECK(mp_int_to_string(op, r, str, len));
 304
 305   /* Change case to match gmp */
 306   for (i = 0; i < len - 1; i++) {
 307     if (radix < 0) {
 308       str[i] = toupper(str[i]);
 309     } else {
 310       str[i] = tolower(str[i]);
 311     }
 312   }
 313   return str;
 314 }
 315
 316 /* gmp: mpq_get_str */
 317 char *GMPQAPI(get_str)(char *str, int radix, mp_rat op) {
 318   int i, r, len;
 319
 320   /* Only print numerator if it is a whole number */
 321   if (mp_int_compare_value(mp_rat_denom_ref(op), 1) == 0)
 322     return GMPZAPI(get_str)(str, radix, mp_rat_numer_ref(op));
 323
 324   /* Support negative radix like gmp */
 325   r = radix;
 326   if (r < 0) r = -r;
 327
 328   /* Compute the length of the string needed to hold the int */
 329   len = mp_rat_string_len(op, r);
 330   if (str == NULL) {
 331     str = malloc(len);
 332   }
 333
 334   /* Convert to string using imath function */
 335   CHECK(mp_rat_to_string(op, r, str, len));
 336
 337   /* Change case to match gmp */
 338   for (i = 0; i < len; i++) {
 339     if (radix < 0) {
 340       str[i] = toupper(str[i]);
 341     } else {
 342       str[i] = tolower(str[i]);
 343     }
 344   }
 345
 346   return str;
 347 }
 348
 349 /* gmp: mpz_set_str */
 350 int GMPZAPI(set_str)(mp_int rop, char *str, int base) {
 351   mp_result res = mp_int_read_string(rop, base, str);
 352   return ((res == MP_OK) ? 0 : -1);
 353 }
 354
 355 /* gmp: mpq_set_str */
 356 int GMPQAPI(set_str)(mp_rat rop, char *s, int base) {
 357   char *slash;
 358   char *str;
 359   mp_result resN;
 360   mp_result resD;
 361   int res = 0;
 362
 363   /* Copy string to temporary storage so we can modify it below */
 364   str = malloc(strlen(s) + 1);
 365   strcpy(str, s);
 366
 367   /* Properly format the string as an int by terminating at the / */
 368   slash = strchr(str, '/');
 369   if (slash) *slash = '\0';
 370
 371   /* Parse numerator */
 372   resN = mp_int_read_string(mp_rat_numer_ref(rop), base, str);
 373
 374   /* Parse denominator if given or set to 1 if not */
 375   if (slash) {
 376     resD = mp_int_read_string(mp_rat_denom_ref(rop), base, slash + 1);
 377   } else {
 378     resD = mp_int_set_uvalue(mp_rat_denom_ref(rop), 1);
 379   }
 380
 381   /* Return failure if either parse failed */
 382   if (resN != MP_OK || resD != MP_OK) {
 383     res = -1;
 384   }
 385
 386   free(str);
 387   return res;
 388 }
 389
 390 static unsigned long get_long_bits(mp_int op) {
 391   /* Deal with integer that does not fit into unsigned long. We want to grab
 392    * the least significant digits that will fit into the long.  Read the digits
 393    * into the long starting at the most significant digit that fits into a
 394    * long. The long is shifted over by MP_DIGIT_BIT before each digit is added.
 395    *
 396    * The shift is decomposed into two steps (following the pattern used in the
 397    * rest of the imath library) to accommodate architectures that don't deal
 398    * well with 32-bit shifts.
 399    */
 400   mp_size digits_to_copy =
 401       (sizeof(unsigned long) + sizeof(mp_digit) - 1) / sizeof(mp_digit);
 402   if (digits_to_copy > MP_USED(op)) {
 403     digits_to_copy = MP_USED(op);
 404   }
 405
 406   mp_digit *digits = MP_DIGITS(op);
 407   unsigned long out = 0;
 408
 409   for (int i = digits_to_copy - 1; i >= 0; i--) {
 410     out <<= (MP_DIGIT_BIT / 2);
 411     out <<= (MP_DIGIT_BIT / 2);
 412     out |= digits[i];
 413   }
 414
 415   return out;
 416 }
 417
 418 /* gmp: mpz_get_ui */
 419 unsigned long GMPZAPI(get_ui)(mp_int op) {
 420   unsigned long out;
 421
 422   /* Try a standard conversion that fits into an unsigned long */
 423   mp_result res = mp_int_to_uint(op, &out);
 424   if (res == MP_OK) return out;
 425
 426   /* Abort the try if we don't have a range error in the conversion.
 427    * The range error indicates that the value cannot fit into a long. */
 428   CHECK(res == MP_RANGE ? MP_OK : MP_RANGE);
 429   if (res != MP_RANGE) return 0;
 430
 431   return get_long_bits(op);
 432 }
 433
 434 /* gmp: mpz_get_si */
 435 long GMPZAPI(get_si)(mp_int op) {
 436   long out;
 437   unsigned long uout;
 438   int long_msb;
 439
 440   /* Try a standard conversion that fits into a long */
 441   mp_result res = mp_int_to_int(op, &out);
 442   if (res == MP_OK) return out;
 443
 444   /* Abort the try if we don't have a range error in the conversion.
 445    * The range error indicates that the value cannot fit into a long. */
 446   CHECK(res == MP_RANGE ? MP_OK : MP_RANGE);
 447   if (res != MP_RANGE) return 0;
 448
 449   /* get least significant bits into an unsigned long */
 450   uout = get_long_bits(op);
 451
 452   /* clear the top bit */
 453   long_msb = (sizeof(unsigned long) * 8) - 1;
 454   uout &= (~(1UL << long_msb));
 455
 456   /* convert to negative if needed based on sign of op */
 457   if (MP_SIGN(op) == MP_NEG) {
 458     uout = 0 - uout;
 459   }
 460
 461   out = (long)uout;
 462   return out;
 463 }
 464
 465 /* gmp: mpz_lcm */
 466 void GMPZAPI(lcm)(mp_int rop, mp_int op1, mp_int op2) {
 467   int op1_is_zero = mp_int_compare_zero(op1) == 0;
 468   int op2_is_zero = mp_int_compare_zero(op2) == 0;
 469
 470   if (op1_is_zero || op2_is_zero) {
 471     mp_int_zero(rop);
 472     return;
 473   }
 474
 475   CHECK(mp_int_lcm(op1, op2, rop));
 476   CHECK(mp_int_abs(rop, rop));
 477 }
 478
 479 /* gmp: mpz_mul_2exp */
 480 /* gmp: allow big values for op2 when op1 == 0 */
 481 void GMPZAPI(mul_2exp)(mp_int rop, mp_int op1, unsigned long op2) {
 482   if (mp_int_compare_zero(op1) == 0)
 483     mp_int_zero(rop);
 484   else
 485     CHECK(mp_int_mul_pow2(op1, op2, rop));
 486 }
 487
 488 /*
 489  * Functions needing expanded functionality
 490  */
 491 /* [Note]Overview of division implementation
 492
 493     All division operations (N / D) compute q and r such that
 494
 495       N = q * D + r, with 0 <= abs(r) < abs(d)
 496
 497     The q and r values are not uniquely specified by N and D. To specify which q
 498     and r values should be used, GMP implements three different rounding modes
 499     for integer division:
 500
 501       ceiling  - round q twords +infinity, r has opposite sign as d
 502       floor    - round q twords -infinity, r has same sign as d
 503       truncate - round q twords zero,      r has same sign as n
 504
 505     The imath library only supports truncate as a rounding mode. We need to
 506     implement the other rounding modes in terms of truncating division. We first
 507     perform the division in trucate mode and then adjust q accordingly. Once we
 508     know q, we can easily compute the correct r according the the formula above
 509     by computing:
 510
 511       r = N - q * D
 512
 513     The main task is to compute q. We can compute the correct q from a trucated
 514     version as follows.
 515
 516     For ceiling rounding mode, if q is less than 0 then the truncated rounding
 517     mode is the same as the ceiling rounding mode.  If q is greater than zero
 518     then we need to round q up by one because the truncated version was rounded
 519     down to zero. If q equals zero then check to see if the result of the
 520     divison is positive. A positive result needs to increment q to one.
 521
 522     For floor rounding mode, if q is greater than 0 then the trucated rounding
 523     mode is the same as the floor rounding mode. If q is less than zero then we
 524     need to round q down by one because the trucated mode rounded q up by one
 525     twords zero. If q is zero then we need to check to see if the result of the
 526     division is negative. A negative result needs to decrement q to negative
 527     one.
 528  */
 529
 530 /* gmp: mpz_cdiv_q */
 531 void GMPZAPI(cdiv_q)(mp_int q, mp_int n, mp_int d) {
 532   mpz_t rz;
 533   mp_int r = &rz;
 534   int qsign, rsign, nsign, dsign;
 535   CHECK(mp_int_init(r));
 536
 537   /* save signs before division because q can alias with n or d */
 538   nsign = mp_int_compare_zero(n);
 539   dsign = mp_int_compare_zero(d);
 540
 541   /* truncating division */
 542   CHECK(mp_int_div(n, d, q, r));
 543
 544   /* see: [Note]Overview of division implementation */
 545   qsign = mp_int_compare_zero(q);
 546   rsign = mp_int_compare_zero(r);
 547   if (qsign > 0) {    /* q > 0 */
 548     if (rsign != 0) { /* r != 0 */
 549       CHECK(mp_int_add_value(q, 1, q));
 550     }
 551   } else if (qsign == 0) { /* q == 0 */
 552     if (rsign != 0) {      /* r != 0 */
 553       if ((nsign > 0 && dsign > 0) || (nsign < 0 && dsign < 0)) {
 554         CHECK(mp_int_set_value(q, 1));
 555       }
 556     }
 557   }
 558   mp_int_clear(r);
 559 }
 560
 561 /* gmp: mpz_fdiv_q */
 562 void GMPZAPI(fdiv_q)(mp_int q, mp_int n, mp_int d) {
 563   mpz_t rz;
 564   mp_int r = &rz;
 565   int qsign, rsign, nsign, dsign;
 566   CHECK(mp_int_init(r));
 567
 568   /* save signs before division because q can alias with n or d */
 569   nsign = mp_int_compare_zero(n);
 570   dsign = mp_int_compare_zero(d);
 571
 572   /* truncating division */
 573   CHECK(mp_int_div(n, d, q, r));
 574
 575   /* see: [Note]Overview of division implementation */
 576   qsign = mp_int_compare_zero(q);
 577   rsign = mp_int_compare_zero(r);
 578   if (qsign < 0) {    /* q  < 0 */
 579     if (rsign != 0) { /* r != 0 */
 580       CHECK(mp_int_sub_value(q, 1, q));
 581     }
 582   } else if (qsign == 0) { /* q == 0 */
 583     if (rsign != 0) {      /* r != 0 */
 584       if ((nsign < 0 && dsign > 0) || (nsign > 0 && dsign < 0)) {
 585         CHECK(mp_int_set_value(q, -1));
 586       }
 587     }
 588   }
 589   mp_int_clear(r);
 590 }
 591
 592 /* gmp: mpz_fdiv_r */
 593 void GMPZAPI(fdiv_r)(mp_int r, mp_int n, mp_int d) {
 594   mpz_t qz;
 595   mpz_t tempz;
 596   mpz_t orig_dz;
 597   mpz_t orig_nz;
 598   mp_int q = &qz;
 599   mp_int temp = &tempz;
 600   mp_int orig_d = &orig_dz;
 601   mp_int orig_n = &orig_nz;
 602   CHECK(mp_int_init(q));
 603   CHECK(mp_int_init(temp));
 604   /* Make a copy of n in case n and d in case they overlap with q */
 605   CHECK(mp_int_init_copy(orig_d, d));
 606   CHECK(mp_int_init_copy(orig_n, n));
 607
 608   /* floor division */
 609   GMPZAPI(fdiv_q)(q, n, d);
 610
 611   /* see: [Note]Overview of division implementation */
 612   /* n = q * d + r  ==>  r = n - q * d */
 613   mp_int_mul(q, orig_d, temp);
 614   mp_int_sub(orig_n, temp, r);
 615
 616   mp_int_clear(q);
 617   mp_int_clear(temp);
 618   mp_int_clear(orig_d);
 619   mp_int_clear(orig_n);
 620 }
 621
 622 /* gmp: mpz_tdiv_q */
 623 void GMPZAPI(tdiv_q)(mp_int q, mp_int n, mp_int d) {
 624   /* truncating division*/
 625   CHECK(mp_int_div(n, d, q, NULL));
 626 }
 627
 628 /* gmp: mpz_fdiv_q_ui */
 629 unsigned long GMPZAPI(fdiv_q_ui)(mp_int q, mp_int n, unsigned long d) {
 630   mpz_t tempz;
 631   mp_int temp = &tempz;
 632   mpz_t rz;
 633   mp_int r = &rz;
 634   mpz_t orig_nz;
 635   mp_int orig_n = &orig_nz;
 636   unsigned long rl;
 637   CHECK(mp_int_init_uvalue(temp, d));
 638   CHECK(mp_int_init(r));
 639   /* Make a copy of n in case n and q overlap */
 640   CHECK(mp_int_init_copy(orig_n, n));
 641
 642   /* use floor division mode to compute q and r */
 643   GMPZAPI(fdiv_q)(q, n, temp);
 644   GMPZAPI(fdiv_r)(r, orig_n, temp);
 645   CHECK(mp_int_to_uint(r, &rl));
 646
 647   mp_int_clear(temp);
 648   mp_int_clear(r);
 649   mp_int_clear(orig_n);
 650
 651   return rl;
 652 }
 653
 654 /* gmp: mpz_export */
 655 void *GMPZAPI(export)(void *rop, size_t *countp, int order, size_t size,
 656                       int endian, size_t nails, mp_int op) {
 657   size_t i, j;
 658   size_t num_used_bytes;
 659   size_t num_words, num_missing_bytes;
 660   ssize_t word_offset;
 661   unsigned char *dst;
 662   mp_digit *src;
 663   int src_bits;
 664
 665   /* We do not have a complete implementation. Assert to ensure our
 666    * restrictions are in place.
 667    */
 668   assert(nails == 0 && "Do not support non-full words");
 669   assert(endian == 1 || endian == 0 || endian == -1);
 670   assert(order == 1 || order == -1);
 671
 672   /* Test for zero */
 673   if (mp_int_compare_zero(op) == 0) {
 674     if (countp) *countp = 0;
 675     return rop;
 676   }
 677
 678   /* Calculate how many words we need */
 679   num_used_bytes = mp_int_unsigned_len(op);
 680   num_words = (num_used_bytes + (size - 1)) / size; /* ceil division */
 681   assert(num_used_bytes > 0);
 682
 683   /* Check to see if we will have missing bytes in the last word.
 684
 685      Missing bytes can only occur when the size of words we output is
 686      greater than the size of words used internally by imath. The number of
 687      missing bytes is the number of bytes needed to fill out the last word. If
 688      this number is greater than the size of a single mp_digit, then we need to
 689      pad the word with extra zeros. Otherwise, the missing bytes can be filled
 690      directly from the zeros in the last digit in the number.
 691    */
 692   num_missing_bytes = (size * num_words) - num_used_bytes;
 693   assert(num_missing_bytes < size);
 694
 695   /* Allocate space for the result if needed */
 696   if (rop == NULL) {
 697     rop = malloc(num_words * size);
 698   }
 699
 700   if (endian == 0) {
 701     endian = HOST_ENDIAN;
 702   }
 703
 704   /* Initialize dst and src pointers */
 705   dst = (unsigned char *)rop + (order >= 0 ? (num_words - 1) * size : 0) +
 706         (endian >= 0 ? size - 1 : 0);
 707   src = MP_DIGITS(op);
 708   src_bits = MP_DIGIT_BIT;
 709
 710   word_offset = (endian >= 0 ? size : -size) + (order < 0 ? size : -size);
 711
 712   for (i = 0; i < num_words; i++) {
 713     for (j = 0; j < size && i * size + j < num_used_bytes; j++) {
 714       if (src_bits == 0) {
 715         ++src;
 716         src_bits = MP_DIGIT_BIT;
 717       }
 718       *dst = (*src >> (MP_DIGIT_BIT - src_bits)) & 0xFF;
 719       src_bits -= 8;
 720       dst -= endian;
 721     }
 722     for (; j < size; j++) {
 723       *dst = 0;
 724       dst -= endian;
 725     }
 726     dst += word_offset;
 727   }
 728
 729   if (countp) *countp = num_words;
 730   return rop;
 731 }
 732
 733 /* gmp: mpz_import */
 734 void GMPZAPI(import)(mp_int rop, size_t count, int order, size_t size,
 735                      int endian, size_t nails, const void *op) {
 736   mpz_t tmpz;
 737   mp_int tmp = &tmpz;
 738   size_t total_size;
 739   size_t num_digits;
 740   ssize_t word_offset;
 741   const unsigned char *src;
 742   mp_digit *dst;
 743   int dst_bits;
 744   size_t i, j;
 745   if (count == 0 || op == NULL) return;
 746
 747   /* We do not have a complete implementation. Assert to ensure our
 748    * restrictions are in place. */
 749   assert(nails == 0 && "Do not support non-full words");
 750   assert(endian == 1 || endian == 0 || endian == -1);
 751   assert(order == 1 || order == -1);
 752
 753   if (endian == 0) {
 754     endian = HOST_ENDIAN;
 755   }
 756
 757   /* Compute number of needed digits by ceil division */
 758   total_size = count * size;
 759   num_digits = (total_size + sizeof(mp_digit) - 1) / sizeof(mp_digit);
 760
 761   /* Init temporary */
 762   mp_int_init_size(tmp, num_digits);
 763   for (i = 0; i < num_digits; i++) tmp->digits[i] = 0;
 764
 765   /* Copy bytes */
 766   src = (const unsigned char *)op + (order >= 0 ? (count - 1) * size : 0) +
 767         (endian >= 0 ? size - 1 : 0);
 768   dst = MP_DIGITS(tmp);
 769   dst_bits = 0;
 770
 771   word_offset = (endian >= 0 ? size : -size) + (order < 0 ? size : -size);
 772
 773   for (i = 0; i < count; i++) {
 774     for (j = 0; j < size; j++) {
 775       if (dst_bits == MP_DIGIT_BIT) {
 776         ++dst;
 777         dst_bits = 0;
 778       }
 779       *dst |= ((mp_digit)*src) << dst_bits;
 780       dst_bits += 8;
 781       src -= endian;
 782     }
 783     src += word_offset;
 784   }
 785
 786   tmp->used = num_digits;
 787
 788   /* Remove leading zeros from number */
 789   {
 790     mp_size uz_ = tmp->used;
 791     mp_digit *dz_ = MP_DIGITS(tmp) + uz_ - 1;
 792     while (uz_ > 1 && (*dz_-- == 0)) --uz_;
 793     tmp->used = uz_;
 794   }
 795
 796   /* Copy to destination */
 797   mp_int_copy(tmp, rop);
 798   mp_int_clear(tmp);
 799 }
 800
 801 /* gmp: mpz_sizeinbase */
 802 size_t GMPZAPI(sizeinbase)(mp_int op, int base) {
 803   mp_result res;
 804   size_t size;
 805
 806   /* If op == 0, return 1 */
 807   if (mp_int_compare_zero(op) == 0) return 1;
 808
 809   /* Compute string length in base */
 810   res = mp_int_string_len(op, base);
 811   CHECK((res > 0) == MP_OK);
 812
 813   /* Now adjust the final size by getting rid of string artifacts */
 814   size = res;
 815
 816   /* subtract one for the null terminator */
 817   size -= 1;
 818
 819   /* subtract one for the negative sign */
 820   if (mp_int_compare_zero(op) < 0) size -= 1;
 821
 822   return size;
 823 }