bigint.c - tomo

(586 lines)
   1 // Integer type infos and methods
   2 
   3 #include <stdio.h> // Must be before gmp.h
   4 
   5 #include <ctype.h>
   6 #include <gc.h>
   7 #include <gmp.h>
   8 #include <stdbool.h>
   9 #include <stdint.h>
  10 #include <stdio.h>
  11 #include <stdlib.h>
  12 
  13 #include "datatypes.h"
  14 #include "integers.h"
  15 #include "optionals.h"
  16 #include "print.h"
  17 #include "siphash.h"
  18 #include "text.h"
  19 #include "types.h"
  20 
  21 public
  22 int Int$print(FILE *f, Int_t i) {
  23     if (likely(i.small & 1L)) {
  24         return _print_int(f, (int64_t)((i.small) >> 2L));
  25     } else {
  26         return gmp_fprintf(f, "%Zd", i.big);
  27     }
  28 }
  29 
  30 static Text_t _int64_to_text(int64_t n) {
  31     if (n == INT64_MIN) return Text("-9223372036854775808");
  32 
  33     char buf[21] = {[20] = 0}; // Big enough for INT64_MIN + '\0'
  34     char *p = &buf[19];
  35     bool negative = n < 0;
  36     if (negative) n = -n; // Safe to do because we checked for INT64_MIN earlier
  37 
  38     do {
  39         *(p--) = '0' + (n % 10);
  40         n /= 10;
  41     } while (n > 0);
  42 
  43     if (negative) *(p--) = '-';
  44 
  45     return Text$from_strn(p + 1, (size_t)(&buf[19] - p));
  46 }
  47 
  48 public
  49 Text_t Int$value_as_text(Int_t i) {
  50     if (likely(i.small & 1L)) {
  51         return _int64_to_text(i.small >> 2L);
  52     } else {
  53         char *str = mpz_get_str(NULL, 10, i.big);
  54         return Text$from_str(str);
  55     }
  56 }
  57 
  58 public
  59 Text_t Int$as_text(const void *i, bool colorize, const TypeInfo_t *info) {
  60     (void)info;
  61     if (!i) return Text("Int");
  62     Text_t text = Int$value_as_text(*(Int_t *)i);
  63     if (colorize) text = Text$concat(Text("\x1b[35m"), text, Text("\x1b[m"));
  64     return text;
  65 }
  66 
  67 static bool Int$is_none(const void *i, const TypeInfo_t *info) {
  68     (void)info;
  69     return ((Int_t *)i)->small == 0L;
  70 }
  71 
  72 public
  73 PUREFUNC int32_t Int$compare_value(const Int_t x, const Int_t y) {
  74     if (likely(x.small & y.small & 1L)) return (x.small > y.small) - (x.small < y.small);
  75     else if (x.small & 1) return -mpz_cmp_si(y.big, (x.small >> 2));
  76     else if (y.small & 1) return mpz_cmp_si(x.big, (y.small >> 2));
  77     else return x.big == y.big ? 0 : mpz_cmp(x.big, y.big);
  78 }
  79 
  80 public
  81 PUREFUNC int32_t Int$compare(const void *x, const void *y, const TypeInfo_t *info) {
  82     (void)info;
  83     return Int$compare_value(*(Int_t *)x, *(Int_t *)y);
  84 }
  85 
  86 public
  87 PUREFUNC bool Int$equal_value(const Int_t x, const Int_t y) {
  88     if (likely((x.small | y.small) & 1L)) return x.small == y.small;
  89     else return x.big == y.big ? 0 : (mpz_cmp(x.big, y.big) == 0);
  90 }
  91 
  92 public
  93 PUREFUNC bool Int$equal(const void *x, const void *y, const TypeInfo_t *info) {
  94     (void)info;
  95     return Int$equal_value(*(Int_t *)x, *(Int_t *)y);
  96 }
  97 
  98 public
  99 CONSTFUNC Int_t Int$clamped(Int_t x, Int_t low, Int_t high) {
 100     return (Int$compare(&x, &low, &Int$info) <= 0) ? low : (Int$compare(&x, &high, &Int$info) >= 0 ? high : x);
 101 }
 102 
 103 public
 104 CONSTFUNC bool Int$is_between(const Int_t x, const Int_t low, const Int_t high) {
 105     int32_t low_cmp = Int$compare_value(x, low);
 106     int32_t high_cmp = Int$compare_value(x, high);
 107     return (low_cmp >= 0 && high_cmp <= 0) || (low_cmp <= 0 && high_cmp >= 0);
 108 }
 109 
 110 public
 111 PUREFUNC uint64_t Int$hash(const void *vx, const TypeInfo_t *info) {
 112     (void)info;
 113     Int_t *x = (Int_t *)vx;
 114     if (likely(x->small & 1L)) {
 115         return siphash24((void *)x, sizeof(Int_t));
 116     } else {
 117         char *str = mpz_get_str(NULL, 16, x->big);
 118         return siphash24((void *)str, strlen(str));
 119     }
 120 }
 121 
 122 public
 123 Text_t Int$hex(Int_t i, Int_t digits_int, bool uppercase, bool prefix) {
 124     if (Int$is_negative(i)) return Text$concat(Text("-"), Int$hex(Int$negative(i), digits_int, uppercase, prefix));
 125 
 126     if (likely(i.small & 1L)) {
 127         uint64_t u64 = (uint64_t)(i.small >> 2);
 128         return Text$from_str(String(
 129             hex(u64, .no_prefix = !prefix, .digits = Int32$from_int(digits_int, false), .uppercase = uppercase)));
 130     } else {
 131         char *str = mpz_get_str(NULL, 16, i.big);
 132         if (uppercase) {
 133             for (char *c = str; *c; c++)
 134                 *c = (char)toupper(*c);
 135         }
 136         int64_t digits = Int64$from_int(digits_int, false);
 137         int64_t needed_zeroes = digits - (int64_t)strlen(str);
 138         if (needed_zeroes <= 0) return prefix ? Text$concat(Text("0x"), Text$from_str(str)) : Text$from_str(str);
 139 
 140         char *zeroes = GC_MALLOC_ATOMIC((size_t)(needed_zeroes));
 141         memset(zeroes, '0', (size_t)(needed_zeroes));
 142         if (prefix) return Text$concat(Text("0x"), Text$from_str(zeroes), Text$from_str(str));
 143         else return Text$concat(Text$from_str(zeroes), Text$from_str(str));
 144     }
 145 }
 146 
 147 public
 148 Text_t Int$octal(Int_t i, Int_t digits_int, bool prefix) {
 149     if (Int$is_negative(i)) return Text$concat(Text("-"), Int$octal(Int$negative(i), digits_int, prefix));
 150 
 151     if (likely(i.small & 1L)) {
 152         uint64_t u64 = (uint64_t)(i.small >> 2);
 153         return Text$from_str(String(oct(u64, .no_prefix = !prefix, .digits = Int32$from_int(digits_int, false))));
 154     } else {
 155         int64_t digits = Int64$from_int(digits_int, false);
 156         char *str = mpz_get_str(NULL, 8, i.big);
 157         int64_t needed_zeroes = digits - (int64_t)strlen(str);
 158         if (needed_zeroes <= 0) return prefix ? Text$concat(Text("0o"), Text$from_str(str)) : Text$from_str(str);
 159 
 160         char *zeroes = GC_MALLOC_ATOMIC((size_t)(needed_zeroes));
 161         memset(zeroes, '0', (size_t)(needed_zeroes));
 162         if (prefix) return Text$concat(Text("0o"), Text$from_str(zeroes), Text$from_str(str));
 163         else return Text$concat(Text$from_str(zeroes), Text$from_str(str));
 164     }
 165 }
 166 
 167 public
 168 Int_t Int$slow_plus(Int_t x, Int_t y) {
 169     mpz_t result;
 170     mpz_init_set_int(result, x);
 171     if (y.small & 1L) {
 172         if (y.small < 0L) mpz_sub_ui(result, result, (uint64_t)(-(y.small >> 2L)));
 173         else mpz_add_ui(result, result, (uint64_t)(y.small >> 2L));
 174     } else {
 175         mpz_add(result, result, y.big);
 176     }
 177     return Int$from_mpz(result);
 178 }
 179 
 180 public
 181 Int_t Int$slow_minus(Int_t x, Int_t y) {
 182     mpz_t result;
 183     mpz_init_set_int(result, x);
 184     if (y.small & 1L) {
 185         if (y.small < 0L) mpz_add_ui(result, result, (uint64_t)(-(y.small >> 2L)));
 186         else mpz_sub_ui(result, result, (uint64_t)(y.small >> 2L));
 187     } else {
 188         mpz_sub(result, result, y.big);
 189     }
 190     return Int$from_mpz(result);
 191 }
 192 
 193 public
 194 Int_t Int$slow_times(Int_t x, Int_t y) {
 195     mpz_t result;
 196     mpz_init_set_int(result, x);
 197     if (y.small & 1L) mpz_mul_si(result, result, y.small >> 2L);
 198     else mpz_mul(result, result, y.big);
 199     return Int$from_mpz(result);
 200 }
 201 
 202 public
 203 Int_t Int$slow_divided_by(Int_t dividend, Int_t divisor) {
 204     // Euclidean division, see:
 205     // https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/divmodnote-letter.pdf
 206     mpz_t quotient, remainder;
 207     mpz_init_set_int(quotient, dividend);
 208     mpz_init_set_int(remainder, divisor);
 209     mpz_tdiv_qr(quotient, remainder, quotient, remainder);
 210     if (mpz_sgn(remainder) < 0) {
 211         bool d_positive = likely(divisor.small & 1L) ? divisor.small > 0x1L : mpz_sgn(divisor.big) > 0;
 212         if (d_positive) mpz_sub_ui(quotient, quotient, 1);
 213         else mpz_add_ui(quotient, quotient, 1);
 214     }
 215     return Int$from_mpz(quotient);
 216 }
 217 
 218 public
 219 Int_t Int$slow_modulo(Int_t x, Int_t modulus) {
 220     mpz_t result;
 221     mpz_init_set_int(result, x);
 222     mpz_t divisor;
 223     mpz_init_set_int(divisor, modulus);
 224     mpz_mod(result, result, divisor);
 225     return Int$from_mpz(result);
 226 }
 227 
 228 public
 229 Int_t Int$slow_modulo1(Int_t x, Int_t modulus) {
 230     mpz_t result;
 231     mpz_init_set_int(result, x);
 232     mpz_sub_ui(result, result, 1);
 233     mpz_t divisor;
 234     mpz_init_set_int(divisor, modulus);
 235     mpz_mod(result, result, divisor);
 236     mpz_add_ui(result, result, 1);
 237     return Int$from_mpz(result);
 238 }
 239 
 240 public
 241 Int_t Int$slow_left_shifted(Int_t x, Int_t y) {
 242     mp_bitcnt_t bits = (mp_bitcnt_t)Int64$from_int(y, false);
 243     mpz_t result;
 244     mpz_init_set_int(result, x);
 245     mpz_mul_2exp(result, result, bits);
 246     return Int$from_mpz(result);
 247 }
 248 
 249 public
 250 Int_t Int$slow_right_shifted(Int_t x, Int_t y) {
 251     mp_bitcnt_t bits = (mp_bitcnt_t)Int64$from_int(y, false);
 252     mpz_t result;
 253     mpz_init_set_int(result, x);
 254     mpz_tdiv_q_2exp(result, result, bits);
 255     return Int$from_mpz(result);
 256 }
 257 
 258 public
 259 Int_t Int$slow_bit_and(Int_t x, Int_t y) {
 260     mpz_t result;
 261     mpz_init_set_int(result, x);
 262     mpz_t y_mpz;
 263     mpz_init_set_int(y_mpz, y);
 264     mpz_and(result, result, y_mpz);
 265     return Int$from_mpz(result);
 266 }
 267 
 268 public
 269 Int_t Int$slow_bit_or(Int_t x, Int_t y) {
 270     mpz_t result;
 271     mpz_init_set_int(result, x);
 272     mpz_t y_mpz;
 273     mpz_init_set_int(y_mpz, y);
 274     mpz_ior(result, result, y_mpz);
 275     return Int$from_mpz(result);
 276 }
 277 
 278 public
 279 Int_t Int$slow_bit_xor(Int_t x, Int_t y) {
 280     mpz_t result;
 281     mpz_init_set_int(result, x);
 282     mpz_t y_mpz;
 283     mpz_init_set_int(y_mpz, y);
 284     mpz_xor(result, result, y_mpz);
 285     return Int$from_mpz(result);
 286 }
 287 
 288 public
 289 Int_t Int$slow_negated(Int_t x) {
 290     mpz_t result;
 291     mpz_init_set_int(result, x);
 292     mpz_neg(result, result);
 293     mpz_sub_ui(result, result, 1);
 294     return Int$from_mpz(result);
 295 }
 296 
 297 public
 298 Int_t Int$slow_negative(Int_t x) {
 299     if (likely(x.small & 1L)) return (Int_t){.small = 4L * -((x.small) >> 2L) + 1L};
 300 
 301     mpz_t result;
 302     mpz_init_set_int(result, x);
 303     mpz_neg(result, result);
 304     return Int$from_mpz(result);
 305 }
 306 
 307 public
 308 Int_t Int$abs(Int_t x) {
 309     if (likely(x.small & 1L)) return (Int_t){.small = 4L * labs((x.small) >> 2L) + 1L};
 310 
 311     mpz_t result;
 312     mpz_init_set_int(result, x);
 313     mpz_abs(result, result);
 314     return Int$from_mpz(result);
 315 }
 316 
 317 public
 318 Int_t Int$power(Int_t base, Int_t exponent) {
 319     int64_t exp = Int64$from_int(exponent, false);
 320     if (unlikely(exp < 0)) fail("Cannot take a negative power of an integer!");
 321     mpz_t result;
 322     mpz_init_set_int(result, base);
 323     mpz_pow_ui(result, result, (uint64_t)exp);
 324     return Int$from_mpz(result);
 325 }
 326 
 327 public
 328 Int_t Int$gcd(Int_t x, Int_t y) {
 329     if (likely(x.small & y.small & 0x1L)) return I_small(Int32$gcd(x.small >> 2L, y.small >> 2L));
 330 
 331     mpz_t result;
 332     mpz_init(result);
 333     if (x.small & 0x1L) mpz_gcd_ui(result, y.big, (uint64_t)labs(x.small >> 2L));
 334     else if (y.small & 0x1L) mpz_gcd_ui(result, x.big, (uint64_t)labs(y.small >> 2L));
 335     else mpz_gcd(result, x.big, y.big);
 336     return Int$from_mpz(result);
 337 }
 338 
 339 public
 340 OptionalInt_t Int$sqrt(Int_t i) {
 341     if (Int$compare_value(i, I(0)) < 0) return NONE_INT;
 342     mpz_t result;
 343     mpz_init_set_int(result, i);
 344     mpz_sqrt(result, result);
 345     return Int$from_mpz(result);
 346 }
 347 
 348 public
 349 bool Int$get_bit(Int_t x, Int_t bit_index) {
 350     mpz_t i;
 351     mpz_init_set_int(i, x);
 352     if (Int$compare_value(bit_index, I(1)) < 0) fail("Invalid bit index (expected 1 or higher): ", bit_index);
 353     if (Int$compare_value(bit_index, Int$from_int64(INT64_MAX)) > 0) fail("Bit index is too large! ", bit_index);
 354 
 355     int is_bit_set = mpz_tstbit(i, (mp_bitcnt_t)(Int64$from_int(bit_index, true) - 1));
 356     return (bool)is_bit_set;
 357 }
 358 
 359 typedef struct {
 360     OptionalInt_t current, last;
 361     Int_t step;
 362 } IntRange_t;
 363 
 364 static OptionalInt_t _next_int(IntRange_t *info) {
 365     OptionalInt_t i = info->current;
 366     if (!Int$is_none(&i, &Int$info)) {
 367         Int_t next = Int$plus(i, info->step);
 368         if (!Int$is_none(&info->last, &Int$info)
 369             && Int$compare_value(next, info->last) == Int$compare_value(info->step, I(0)))
 370             next = NONE_INT;
 371         info->current = next;
 372     }
 373     return i;
 374 }
 375 
 376 public
 377 PUREFUNC Closure_t Int$to(Int_t first, Int_t last, OptionalInt_t step) {
 378     IntRange_t *range = GC_MALLOC(sizeof(IntRange_t));
 379     range->current = first;
 380     range->last = last;
 381     range->step = Int$is_none(&step, &Int$info) ? Int$compare_value(last, first) >= 0
 382                                                       ? (Int_t){.small = (1L << 2L) | 1L}
 383                                                       : (Int_t){.small = (-1L >> 2L) | 1L}
 384                                                 : step;
 385     return (Closure_t){.fn = _next_int, .userdata = range};
 386 }
 387 
 388 public
 389 PUREFUNC Closure_t Int$onward(Int_t first, Int_t step) {
 390     IntRange_t *range = GC_MALLOC(sizeof(IntRange_t));
 391     range->current = first;
 392     range->last = NONE_INT;
 393     range->step = step;
 394     return (Closure_t){.fn = _next_int, .userdata = range};
 395 }
 396 
 397 public
 398 Int_t Int$from_str(const char *str) {
 399     return Int$parse(Text$from_str(str), NONE_INT, NULL);
 400 }
 401 
 402 public
 403 OptionalInt_t Int$parse(Text_t text, OptionalInt_t base, Text_t *remainder) {
 404     const char *str = Text$as_c_string(text);
 405     bool negative = (*str == '-');
 406     if (negative || *str == '+') str += 1;
 407     const char *end = str;
 408     int32_t base32;
 409     if (base.small != 0) {
 410         base32 = Int32$from_int(base, false);
 411         switch (base32) {
 412         case 16:
 413             if (strncmp(str, "0x", 2) == 0) {
 414             base16_prefix:
 415                 str += 2;
 416             }
 417             end = str + strspn(str, "0123456789abcdefABCDEF");
 418             break;
 419         case 10:
 420         base10:
 421             end = str + strspn(str, "0123456789");
 422             break;
 423         case 8:
 424             if (strncmp(str, "0o", 2) == 0) {
 425             base8_prefix:
 426                 str += 2;
 427             }
 428             end = str + strspn(str, "01234567");
 429             break;
 430         case 2:
 431             if (strncmp(str, "0b", 2) == 0) {
 432             base2_prefix:
 433                 str += 2;
 434             }
 435             end = str + strspn(str, "01");
 436             break;
 437         case 1: {
 438             str += strspn(str, "0");
 439             size_t n = strspn(str, "1");
 440             end = str + n;
 441             if (remainder) *remainder = Text$from_str(end);
 442             else if (*end != '\0') return NONE_INT;
 443             return Int$from_int64((int64_t)n);
 444         }
 445         default: {
 446             if (base32 < 1 || base32 > 36) {
 447                 if (remainder) *remainder = text;
 448                 return NONE_INT;
 449             }
 450             for (; *end; end++) {
 451                 char c = *end;
 452                 int32_t digit;
 453                 if ('0' <= c && c <= '9') {
 454                     digit = (c - (int)'0');
 455                 } else if ('a' <= c && c <= 'z') {
 456                     digit = (c - (int)'a');
 457                 } else if ('A' <= c && c <= 'Z') {
 458                     digit = (c - (int)'A');
 459                 } else {
 460                     break;
 461                 }
 462                 if (digit >= base32) break;
 463             }
 464         }
 465         }
 466     } else if (strncmp(str, "0x", 2) == 0) {
 467         base32 = 16;
 468         goto base16_prefix;
 469     } else if (strncmp(str, "0o", 2) == 0) {
 470         base32 = 8;
 471         goto base8_prefix;
 472     } else if (strncmp(str, "0b", 2) == 0) {
 473         base32 = 2;
 474         goto base2_prefix;
 475     } else {
 476         base32 = 10;
 477         goto base10;
 478     }
 479 
 480     if (remainder) *remainder = Text$from_str(end);
 481     else if (*end != '\0') return NONE_INT;
 482 
 483     mpz_t i;
 484     int result = mpz_init_set_str(i, String(string_slice(str, (size_t)(end - str))), base32);
 485     if (result != 0) {
 486         if (remainder) *remainder = text;
 487         return NONE_INT;
 488     }
 489     if (negative) {
 490         mpz_neg(i, i);
 491     }
 492     return Int$from_mpz(i);
 493 }
 494 
 495 public
 496 bool Int$is_prime(Int_t x, Int_t reps) {
 497     mpz_t p;
 498     mpz_init_set_int(p, x);
 499     if (unlikely(Int$compare_value(reps, I(9999)) > 0))
 500         fail("Number of prime-test repetitions should not be above 9999");
 501     int reps_int = Int32$from_int(reps, false);
 502     return (mpz_probab_prime_p(p, reps_int) != 0);
 503 }
 504 
 505 public
 506 Int_t Int$next_prime(Int_t x) {
 507     mpz_t p;
 508     mpz_init_set_int(p, x);
 509     mpz_nextprime(p, p);
 510     return Int$from_mpz(p);
 511 }
 512 
 513 public
 514 Int_t Int$choose(Int_t n, Int_t k) {
 515     if unlikely (Int$compare_value(n, I_small(0)) < 0) fail("Negative inputs are not supported for choose()");
 516 
 517     mpz_t ret;
 518     mpz_init(ret);
 519 
 520     int64_t k_i64 = Int64$from_int(k, false);
 521     if unlikely (k_i64 < 0) fail("Negative inputs are not supported for choose()");
 522 
 523     if likely (n.small & 1L) {
 524         mpz_bin_uiui(ret, (unsigned long)(n.small >> 2L), (unsigned long)k_i64);
 525     } else {
 526         mpz_t n_mpz;
 527         mpz_init_set_int(n_mpz, n);
 528         mpz_bin_ui(ret, n_mpz, (unsigned long)k_i64);
 529     }
 530     return Int$from_mpz(ret);
 531 }
 532 
 533 public
 534 Int_t Int$factorial(Int_t n) {
 535     mpz_t ret;
 536     mpz_init(ret);
 537     int64_t n_i64 = Int64$from_int(n, false);
 538     if unlikely (n_i64 < 0) fail("Factorials are not defined for negative numbers");
 539     mpz_fac_ui(ret, (unsigned long)n_i64);
 540     return Int$from_mpz(ret);
 541 }
 542 
 543 static void Int$serialize(const void *obj, FILE *out, Table_t *pointers, const TypeInfo_t *info) {
 544     (void)info;
 545     Int_t i = *(Int_t *)obj;
 546     if (likely(i.small & 1L)) {
 547         fputc(0, out);
 548         int64_t i64 = i.small >> 2L;
 549         Int64$serialize(&i64, out, pointers, &Int64$info);
 550     } else {
 551         fputc(1, out);
 552         mpz_t n;
 553         mpz_init_set_int(n, *(Int_t *)obj);
 554         mpz_out_raw(out, n);
 555     }
 556 }
 557 
 558 static void Int$deserialize(FILE *in, void *obj, List_t *pointers, const TypeInfo_t *info) {
 559     (void)info;
 560     if (fgetc(in) == 0) {
 561         int64_t i = 0;
 562         Int64$deserialize(in, &i, pointers, &Int64$info);
 563         *(Int_t *)obj = (Int_t){.small = (i << 2L) | 1L};
 564     } else {
 565         mpz_t n;
 566         mpz_init(n);
 567         mpz_inp_raw(n, in);
 568         *(Int_t *)obj = Int$from_mpz(n);
 569     }
 570 }
 571 
 572 public
 573 const TypeInfo_t Int$info = {
 574     .size = sizeof(Int_t),
 575     .align = __alignof__(Int_t),
 576     .metamethods =
 577         {
 578             .compare = Int$compare,
 579             .equal = Int$equal,
 580             .hash = Int$hash,
 581             .as_text = Int$as_text,
 582             .is_none = Int$is_none,
 583             .serialize = Int$serialize,
 584             .deserialize = Int$deserialize,
 585         },
 586 };