#include "libtommath/tommath.h" #if 0 void *MP_MALLOC(size_t size) { return malloc(size); } void *MP_REALLOC(void *mem, size_t oldsize, size_t newsize) { return realloc(mem, newsize); } void *MP_CALLOC(size_t nmemb, size_t size) { return calloc(nmemb, size); } void MP_FREE(void *mem, size_t size) { free(mem); } #else void *MP_MALLOC(size_t size) { return gb_alloc(permanent_allocator(), cast(isize)size); } void *MP_REALLOC(void *mem, size_t oldsize, size_t newsize) { return gb_resize(permanent_allocator(), mem, cast(isize)oldsize, cast(isize)newsize); } void *MP_CALLOC(size_t nmemb, size_t size) { size_t total = nmemb*size; return gb_alloc(permanent_allocator(), cast(isize)total); } void MP_FREE(void *mem, size_t size) { // DO NOTHING } #endif #ifndef MAX_BIG_INT_SHIFT #define MAX_BIG_INT_SHIFT 1024 #endif typedef mp_int BigInt; gb_internal void big_int_from_u64(BigInt *dst, u64 x); gb_internal void big_int_from_i64(BigInt *dst, i64 x); gb_internal void big_int_init (BigInt *dst, BigInt const *src); gb_internal void big_int_from_string(BigInt *dst, String const &s, bool *success); gb_internal void big_int_dealloc(BigInt *dst) { mp_clear(dst); } gb_internal BigInt big_int_make(BigInt const *b, bool abs=false); gb_internal BigInt big_int_make_abs(BigInt const *b); gb_internal BigInt big_int_make_u64(u64 x); gb_internal BigInt big_int_make_i64(i64 x); gb_internal u64 big_int_to_u64 (BigInt const *x); gb_internal i64 big_int_to_i64 (BigInt const *x); gb_internal f64 big_int_to_f64 (BigInt const *x); gb_internal String big_int_to_string(gbAllocator allocator, BigInt const *x, u64 base = 10); gb_internal void big_int_add (BigInt *dst, BigInt const *x, BigInt const *y); gb_internal void big_int_sub (BigInt *dst, BigInt const *x, BigInt const *y); gb_internal void big_int_shl (BigInt *dst, BigInt const *x, BigInt const *y); gb_internal void big_int_shr (BigInt *dst, BigInt const *x, BigInt const *y); gb_internal void big_int_mul (BigInt *dst, BigInt const *x, BigInt const *y); gb_internal void big_int_mul_u64(BigInt *dst, BigInt const *x, u64 y); gb_internal void big_int_exp_u64(BigInt *dst, BigInt const *x, u64 y, bool *success); gb_internal void big_int_quo_rem(BigInt const *x, BigInt const *y, BigInt *q, BigInt *r); gb_internal void big_int_quo (BigInt *z, BigInt const *x, BigInt const *y); gb_internal void big_int_rem (BigInt *z, BigInt const *x, BigInt const *y); gb_internal void big_int_and (BigInt *dst, BigInt const *x, BigInt const *y); gb_internal void big_int_and_not(BigInt *dst, BigInt const *x, BigInt const *y); gb_internal void big_int_xor (BigInt *dst, BigInt const *x, BigInt const *y); gb_internal void big_int_or (BigInt *dst, BigInt const *x, BigInt const *y); gb_internal void big_int_not (BigInt *dst, BigInt const *x, i32 bit_count, bool is_signed); gb_internal void big_int_add_eq(BigInt *dst, BigInt const *x); gb_internal void big_int_sub_eq(BigInt *dst, BigInt const *x); gb_internal void big_int_shl_eq(BigInt *dst, BigInt const *x); gb_internal void big_int_shr_eq(BigInt *dst, BigInt const *x); gb_internal void big_int_mul_eq(BigInt *dst, BigInt const *x); gb_internal void big_int_quo_eq(BigInt *dst, BigInt const *x); gb_internal void big_int_rem_eq(BigInt *dst, BigInt const *x); gb_internal bool big_int_is_neg(BigInt const *x); gb_internal void big_int_neg(BigInt *dst, BigInt const *x); gb_internal void big_int_add_eq(BigInt *dst, BigInt const *x) { BigInt res = {}; big_int_init(&res, dst); big_int_add(dst, &res, x); } gb_internal void big_int_sub_eq(BigInt *dst, BigInt const *x) { BigInt res = {}; big_int_init(&res, dst); big_int_sub(dst, &res, x); } gb_internal void big_int_shl_eq(BigInt *dst, BigInt const *x) { BigInt res = {}; big_int_init(&res, dst); big_int_shl(dst, &res, x); } gb_internal void big_int_shr_eq(BigInt *dst, BigInt const *x) { BigInt res = {}; big_int_init(&res, dst); big_int_shr(dst, &res, x); } gb_internal void big_int_mul_eq(BigInt *dst, BigInt const *x) { BigInt res = {}; big_int_init(&res, dst); big_int_mul(dst, &res, x); } gb_internal void big_int_quo_eq(BigInt *dst, BigInt const *x) { BigInt res = {}; big_int_init(&res, dst); big_int_quo(dst, &res, x); } gb_internal void big_int_rem_eq(BigInt *dst, BigInt const *x) { BigInt res = {}; big_int_init(&res, dst); big_int_rem(dst, &res, x); } gb_internal i64 big_int_sign(BigInt const *x) { if (mp_iszero(x)) { return 0; } return x->sign == MP_ZPOS ? +1 : -1; } gb_internal void big_int_from_u64(BigInt *dst, u64 x) { mp_init_u64(dst, x); } gb_internal void big_int_from_i64(BigInt *dst, i64 x) { mp_init_i64(dst, x); } gb_internal void big_int_init(BigInt *dst, BigInt const *src) { if (dst == src) { return; } mp_init_copy(dst, src); } gb_internal BigInt big_int_make(BigInt const *b, bool abs) { BigInt i = {}; big_int_init(&i, b); if (abs) mp_abs(&i, &i); return i; } gb_internal BigInt big_int_make_abs(BigInt const *b) { return big_int_make(b, true); } gb_internal BigInt big_int_make_u64(u64 x) { BigInt i = {}; big_int_from_u64(&i, x); return i; } gb_internal BigInt big_int_make_i64(i64 x) { BigInt i = {}; big_int_from_i64(&i, x); return i; } gb_internal void big_int_from_string(BigInt *dst, String const &s, bool *success) { *success = true; bool is_negative = false; u64 base = 10; bool has_prefix = false; if (s.len > 2 && s[0] == '0') { switch (s[1]) { case 'b': base = 2; has_prefix = true; break; case 'o': base = 8; has_prefix = true; break; case 'd': base = 10; has_prefix = true; break; case 'z': base = 12; has_prefix = true; break; case 'x': base = 16; has_prefix = true; break; case 'h': base = 16; has_prefix = true; break; } } u8 *text = s.text; isize len = s.len; if (has_prefix) { text += 2; len -= 2; } BigInt b = {}; big_int_from_u64(&b, base); mp_zero(dst); isize i = 0; for (; i < len; i++) { Rune r = cast(Rune)text[i]; if (r == '-') { if (is_negative) { // NOTE(Jeroen): Can't have a doubly negative number. *success = false; return; } is_negative = true; continue; } if (r == '_') { continue; } u64 v = u64_digit_value(r); if (v >= base) { // NOTE(Jeroen): Can still be a valid integer if the next character is an `e` or `E`. if (r != 'e' && r != 'E') { *success = false; } break; } BigInt val = big_int_make_u64(v); big_int_mul_eq(dst, &b); big_int_add_eq(dst, &val); } if (i < len && (text[i] == 'e' || text[i] == 'E')) { i += 1; GB_ASSERT(base == 10); GB_ASSERT(text[i] != '-'); if (text[i] == '+') { i += 1; } u64 exp = 0; for (; i < len; i++) { char r = cast(char)text[i]; if (r == '_') { continue; } u64 v = 0; if (gb_char_is_digit(r)) { v = u64_digit_value(r); } else { *success = false; break; } exp *= 10; exp += v; } BigInt tmp = {}; mp_init(&tmp); big_int_exp_u64(&tmp, &b, exp, success); big_int_mul_eq(dst, &tmp); } if (is_negative) { big_int_neg(dst, dst); } } gb_internal u64 big_int_to_u64(BigInt const *x) { GB_ASSERT(x->sign == 0); return mp_get_u64(x); } gb_internal i64 big_int_to_i64(BigInt const *x) { return mp_get_i64(x); } gb_internal f64 big_int_to_f64(BigInt const *x) { return mp_get_double(x); } gb_internal void big_int_neg(BigInt *dst, BigInt const *x) { mp_neg(x, dst); } gb_internal int big_int_cmp(BigInt const *x, BigInt const *y) { return mp_cmp(x, y); } gb_internal int big_int_cmp_zero(BigInt const *x) { if (mp_iszero(x)) { return 0; } return x->sign ? -1 : +1; } gb_internal bool big_int_is_zero(BigInt const *x) { return mp_iszero(x); } gb_internal void big_int_add(BigInt *dst, BigInt const *x, BigInt const *y) { mp_add(x, y, dst); } gb_internal void big_int_sub(BigInt *dst, BigInt const *x, BigInt const *y) { mp_sub(x, y, dst); } gb_internal void big_int_shl(BigInt *dst, BigInt const *x, BigInt const *y) { u32 yy = mp_get_u32(y); mp_mul_2d(x, yy, dst); } gb_internal void big_int_shr(BigInt *dst, BigInt const *x, BigInt const *y) { u32 yy = mp_get_u32(y); BigInt d = {}; mp_div_2d(x, yy, dst, &d); big_int_dealloc(&d); } gb_internal void big_int_mul_u64(BigInt *dst, BigInt const *x, u64 y) { BigInt d = {}; big_int_from_u64(&d, y); mp_mul(x, &d, dst); big_int_dealloc(&d); } gb_internal void big_int_exp_u64(BigInt *dst, BigInt const *x, u64 y, bool *success) { if (y > INT_MAX) { *success = false; return; } // Note: The cutoff for square-multiply being faster than the naive // for loop is when exp > 4, but it probably isn't worth adding // a fast path. mp_err err = mp_expt_n(x, int(y), dst); *success = err == MP_OKAY; } gb_internal void big_int_mul(BigInt *dst, BigInt const *x, BigInt const *y) { mp_mul(x, y, dst); } gb_internal u64 leading_zeros_u64(u64 x) { #if defined(GB_COMPILER_MSVC) #if defined(GB_ARCH_64_BIT) return __lzcnt64(x); #else u64 y, n; n = 0; y = x; L: if (x < 0) { return n; } if (y == 0) { return 64-n; } n++; x <<= 1; y >>= 1; goto L; #endif #else return cast(u64)__builtin_clzll(cast(unsigned long long)x); #endif } // `big_int_quo_rem` sets z to the quotient x/y and r to the remainder x%y // and returns the pair (z, r) for y != 0. // if y == 0, a division-by-zero run-time panic occurs. // // q = x/y with the result truncated to zero // r = x - y*q gb_internal void big_int_quo_rem(BigInt const *x, BigInt const *y, BigInt *q_, BigInt *r_) { mp_div(x, y, q_, r_); } gb_internal void big_int_quo(BigInt *z, BigInt const *x, BigInt const *y) { BigInt r = {}; big_int_quo_rem(x, y, z, &r); big_int_dealloc(&r); } gb_internal void big_int_rem(BigInt *z, BigInt const *x, BigInt const *y) { BigInt q = {}; big_int_quo_rem(x, y, &q, z); big_int_dealloc(&q); } gb_internal void big_int_euclidean_mod(BigInt *z, BigInt const *x, BigInt const *y) { BigInt y0 = {}; big_int_init(&y0, y); BigInt q = {}; big_int_quo_rem(x, y, &q, z); if (z->sign) { if (y0.sign) { big_int_sub(z, z, &y0); } else { big_int_add(z, z, &y0); } } } gb_internal void big_int_and(BigInt *dst, BigInt const *x, BigInt const *y) { mp_and(x, y, dst); } gb_internal void big_int_and_not(BigInt *dst, BigInt const *x, BigInt const *y) { if (mp_iszero(x)) { big_int_init(dst, y); return; } if (mp_iszero(y)) { big_int_init(dst, x); return; } if (x->sign == y->sign) { if (x->sign) { // (-x) &~ (-y) == ~(x-1) &~ ~(y-1) == ~(x-1) & (y-1) == (y-1) &~ (x-1) BigInt x1 = big_int_make_abs(x); BigInt y1 = big_int_make_abs(y); mp_decr(&x1); mp_decr(&y1); BigInt ny1 = {}; mp_complement(&y1, &ny1); mp_and(&x1, &ny1, dst); big_int_dealloc(&x1); big_int_dealloc(&y1); big_int_dealloc(&ny1); return; } BigInt ny = {}; mp_complement(y, &ny); mp_and(x, &ny, dst); big_int_dealloc(&ny); return; } if (x->sign) { // (-x) &~ y == ~(x-1) &~ y == ~(x-1) & ~y == ~((x-1) | y) == -(((x-1) | y) + 1) BigInt x1 = big_int_make_abs(x); BigInt y1 = big_int_make_abs(y); mp_decr(&x1); BigInt z1 = {}; big_int_or(&z1, &x1, &y1); mp_add_d(&z1, 1, dst); big_int_dealloc(&x1); big_int_dealloc(&y1); big_int_dealloc(&z1); return; } // x &~ (-y) == x &~ ~(y-1) == x & (y-1) BigInt x1 = big_int_make_abs(x); BigInt y1 = big_int_make_abs(y); mp_decr(&y1); big_int_and(dst, &x1, &y1); big_int_dealloc(&x1); big_int_dealloc(&y1); return; } gb_internal void big_int_xor(BigInt *dst, BigInt const *x, BigInt const *y) { mp_xor(x, y, dst); } gb_internal void big_int_or(BigInt *dst, BigInt const *x, BigInt const *y) { mp_or(x, y, dst); } gb_internal void debug_print_big_int(BigInt const *x) { TEMPORARY_ALLOCATOR_GUARD(); String s = big_int_to_string(temporary_allocator(), x, 10); gb_printf_err("[DEBUG] %.*s\n", LIT(s)); } gb_internal void big_int_not(BigInt *dst, BigInt const *x, i32 bit_count, bool is_signed) { GB_ASSERT(bit_count >= 0); if (bit_count == 0) { big_int_from_u64(dst, 0); return; } if (big_int_is_neg(x)) { // ~x == -x - 1 big_int_neg(dst, x); mp_decr(dst); mp_mod_2d(dst, bit_count, dst); return; } BigInt pow2b = {}; mp_2expt(&pow2b, bit_count); BigInt mask = {}; mp_2expt(&mask, bit_count); mp_decr(&mask); BigInt v = {}; mp_init_copy(&v, x); mp_mod_2d(&v, bit_count, &v); mp_xor(&v, &mask, dst); if (is_signed) { BigInt pmask = {}; BigInt pmask_minus_one = {}; mp_2expt(&pmask, bit_count-1); mp_sub_d(&pmask, 1, &pmask_minus_one); BigInt a = {}; BigInt b = {}; big_int_and(&a, dst, &pmask_minus_one); big_int_and(&b, dst, &pmask); big_int_sub(dst, &a, &b); big_int_dealloc(&a); big_int_dealloc(&b); } big_int_dealloc(&pow2b); big_int_dealloc(&mask); big_int_dealloc(&v); } gb_internal bool big_int_is_neg(BigInt const *x) { if (x == nullptr) { return false; } return x->sign != MP_ZPOS; } gb_internal char digit_to_char(u8 digit) { GB_ASSERT(digit < 16); if (digit <= 9) { return digit + '0'; } else if (digit <= 15) { return digit + 'a'; } return '0'; } gb_internal String big_int_to_string(gbAllocator allocator, BigInt const *x, u64 base) { GB_ASSERT(base <= 16); if (mp_iszero(x)) { u8 *buf = gb_alloc_array(allocator, u8, 1); buf[0] = '0'; return make_string(buf, 1); } Array buf = {}; array_init(&buf, allocator, 0, 32); if (x->used >= 498) { // 2^498 ~ 10^150 mp_int val = {}; mp_abs(x, &val); int exp = 0; mp_err err = mp_log_n(&val, 10, &exp); GB_ASSERT(err == MP_OKAY); GB_ASSERT(exp >= 100); mp_int thousand_below = {}; mp_int thousand_above = {}; mp_init_i32(&thousand_below, 10); mp_expt_n(&thousand_below, exp-3, &thousand_below); mp_div(&val, &thousand_below, &thousand_above, nullptr); double mant = 1.0e-3 * mp_get_double(&thousand_above); char val_buf[256] = {}; isize n = gb_snprintf(val_buf, gb_size_of(val_buf)-1, "~ %s%.fe%u", (x->sign ? "-" : ""), mant, exp); array_add_elems(&buf, val_buf, n-1); } else { BigInt v = {}; mp_init_copy(&v, x); if (v.sign) { array_add(&buf, '-'); mp_abs(&v, &v); } isize first_word_idx = buf.count; BigInt r = {}; BigInt b = {}; big_int_from_u64(&b, base); u8 digit = 0; while (big_int_cmp(&v, &b) >= 0) { big_int_quo_rem(&v, &b, &v, &r); digit = cast(u8)big_int_to_u64(&r); array_add(&buf, digit_to_char(digit)); } big_int_rem(&r, &v, &b); digit = cast(u8)big_int_to_u64(&r); array_add(&buf, digit_to_char(digit)); big_int_dealloc(&r); big_int_dealloc(&b); for (isize i = first_word_idx; i < buf.count/2; i++) { isize j = buf.count + first_word_idx - i - 1; char tmp = buf[i]; buf[i] = buf[j]; buf[j] = tmp; } } return make_string(cast(u8 *)buf.data, buf.count); } gb_internal int big_int_log2(BigInt const *x) { return mp_count_bits(x) - 1; }