From bac188ce07b957807d4c649cb5d4e5e253360278 Mon Sep 17 00:00:00 2001 From: Bruce Hill Date: Fri, 16 Aug 2024 14:24:20 -0400 Subject: Change division and modulus to use euclidean division, plus fix up a few integer bugs --- builtins/integers.c | 27 +++++++++++++----------- builtins/integers.h | 59 ++++++++++++++++++++++++++++++++++++++++++++--------- 2 files changed, 64 insertions(+), 22 deletions(-) (limited to 'builtins') diff --git a/builtins/integers.c b/builtins/integers.c index bb82fab6..0bf7dc22 100644 --- a/builtins/integers.c +++ b/builtins/integers.c @@ -152,17 +152,20 @@ public Int_t Int$slow_times(Int_t x, Int_t y) { return Int$from_mpz(result); } -public Int_t Int$slow_divided_by(Int_t x, Int_t y) { - mpz_t result; - mpz_init_set_int(result, x); - if (y.small & 1) { - mpz_t y_mpz; - mpz_init_set_si(y_mpz, y.small >> 2); - mpz_cdiv_q(result, result, y_mpz); - } else { - mpz_cdiv_q(result, result, *y.big); +public Int_t Int$slow_divided_by(Int_t dividend, Int_t divisor) { + // Euclidean division, see: https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/divmodnote-letter.pdf + mpz_t quotient, remainder; + mpz_init_set_int(quotient, dividend); + mpz_init_set_int(remainder, divisor); + mpz_tdiv_qr(quotient, remainder, quotient, remainder); + if (mpz_sgn(remainder) < 0) { + bool d_positive = __builtin_expect(divisor.small & 1, 1) ? divisor.small > 0x1 : mpz_sgn(*divisor.big) > 0; + if (d_positive) + mpz_sub_ui(quotient, quotient, 1); + else + mpz_add_ui(quotient, quotient, 1); } - return Int$from_mpz(result); + return Int$from_mpz(quotient); } public Int_t Int$slow_modulo(Int_t x, Int_t modulus) @@ -359,7 +362,7 @@ public const TypeInfo $Int = { } \ public CORD KindOfInt ## $format(c_type i, Int_t digits_int) { \ int64_t digits = Int_to_Int64(digits_int, false); \ - return CORD_asprintf("%0*" fmt, (int)digits, i); \ + return CORD_asprintf("%0*ld", (int)digits, (int64_t)i); \ } \ public CORD KindOfInt ## $hex(c_type i, Int_t digits_int, bool uppercase, bool prefix) { \ int64_t digits = Int_to_Int64(digits_int, false); \ @@ -427,7 +430,7 @@ public const TypeInfo $Int = { .CustomInfo={.compare=(void*)KindOfInt##$compare, .as_text=(void*)KindOfInt##$as_text}, \ }; -DEFINE_INT_TYPE(int64_t, Int64, "ld", INT64_MIN, INT64_MAX); +DEFINE_INT_TYPE(int64_t, Int64, "ld_i64", INT64_MIN, INT64_MAX); DEFINE_INT_TYPE(int32_t, Int32, "d_i32", INT32_MIN, INT32_MAX); DEFINE_INT_TYPE(int16_t, Int16, "d_i16", INT16_MIN, INT16_MAX); DEFINE_INT_TYPE(int8_t, Int8, "d_i8", INT8_MIN, INT8_MAX); diff --git a/builtins/integers.h b/builtins/integers.h index 898469e2..e6b5b1fb 100644 --- a/builtins/integers.h +++ b/builtins/integers.h @@ -35,7 +35,26 @@ Range_t type_name ## $to(c_type from, c_type to); \ c_type type_name ## $from_text(CORD text, CORD *the_rest); \ extern const c_type type_name ## $min, type_name##$max; \ - extern const TypeInfo $ ## type_name; + extern const TypeInfo $ ## type_name; \ + static inline c_type type_name ## $divided_by(c_type D, c_type d) { \ + c_type q = D/d, r = D%d; \ + if (r < 0) { \ + if (d > 0) q = q-1; \ + else q = q+1; \ + } \ + return q; \ + } \ + static inline c_type type_name ## $modulo(c_type D, c_type d) { \ + c_type r = D%d; \ + if (r < 0) { \ + if (d > 0) r = r + d; \ + else r = r - d; \ + } \ + return r; \ + } \ + static inline c_type type_name ## $modulo1(c_type D, c_type d) { \ + return type_name ## $modulo(D-1, d) + 1; \ + } DEFINE_INT_TYPE(int64_t, Int64); DEFINE_INT_TYPE(int32_t, Int32); @@ -128,27 +147,47 @@ static inline Int_t Int$times(Int_t x, Int_t y) { static inline Int_t Int$divided_by(Int_t x, Int_t y) { if (__builtin_expect(((x.small & y.small) & 1) != 0, 1)) { - const int64_t z = ((x.small>>1) / (y.small>>1)) << 2; - if (__builtin_expect(z == (int32_t)z, 1)) - return (Int_t){.small=z|1}; + // Euclidean division, see: https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/divmodnote-letter.pdf + const int64_t D = (x.small>>2); + const int64_t d = (y.small>>2); + int64_t q = D/d; + int64_t r = D%d; + if (r < 0) { + if (d > 0) q = q-1; + else q = q+1; + } + if (__builtin_expect(q == (int32_t)q, 1)) + return (Int_t){.small=(q<<2)|1}; } return Int$slow_divided_by(x, y); } static inline Int_t Int$modulo(Int_t x, Int_t y) { if (__builtin_expect(((x.small & y.small) & 1) != 0, 1)) { - int64_t mod = (x.small>>2) % (y.small>>2); - if (mod < 0) mod += (y.small>>2); - return (Int_t){.small=(mod<<2)+1}; + // Euclidean modulus, see: https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/divmodnote-letter.pdf + const int64_t D = (x.small>>2); + const int64_t d = (y.small>>2); + int64_t r = D%d; + if (r < 0) { + if (d > 0) r = r + d; + else r = r - d; + } + return (Int_t){.small=(r<<2)|1}; } return Int$slow_modulo(x, y); } static inline Int_t Int$modulo1(Int_t x, Int_t y) { if (__builtin_expect(((x.small & y.small) & 1) != 0, 1)) { - int64_t mod = ((x.small>>2)-1) % (y.small>>2); - if (mod < 0) mod += (y.small>>2); - return (Int_t){.small=((mod+1)<<2)+1}; + // Euclidean modulus, see: https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/divmodnote-letter.pdf + const int64_t D = (x.small>>2)-1; + const int64_t d = (y.small>>2); + int64_t r = D%d; + if (r < 0) { + if (d > 0) r = r + d; + else r = r - d; + } + return (Int_t){.small=((r+1)<<2)|1}; } return Int$slow_modulo1(x, y); } -- cgit v1.2.3