From a51d48201b5245dc2cc2bfb00e0ac8e7b52203d9 Mon Sep 17 00:00:00 2001 From: Bruce Hill Date: Fri, 11 Jul 2025 15:38:42 -0400 Subject: Use _Decimal64 instead of mpdecimal --- src/stdlib/datatypes.h | 5 +- src/stdlib/decimals.c | 245 +++++++++++++++++++++++++++---------------------- src/stdlib/decimals.h | 41 ++++----- src/stdlib/print.c | 93 ++++++++++++++++++- src/stdlib/print.h | 4 + src/stdlib/stdlib.c | 12 --- 6 files changed, 252 insertions(+), 148 deletions(-) (limited to 'src/stdlib') diff --git a/src/stdlib/datatypes.h b/src/stdlib/datatypes.h index 00a9411f..8ca2876d 100644 --- a/src/stdlib/datatypes.h +++ b/src/stdlib/datatypes.h @@ -3,7 +3,6 @@ // Common datastructures (lists, tables, closures) #include -#include #include #include #include @@ -23,8 +22,8 @@ #define Num_t double #define Num32_t float -#define Dec_t mpd_t* -#define OptionalDec_t mpd_t* +#define Dec_t _Decimal64 +#define OptionalDec_t _Decimal64 #define Int64_t int64_t #define Int32_t int32_t diff --git a/src/stdlib/decimals.c b/src/stdlib/decimals.c index 8bdc2907..4a3f0c08 100644 --- a/src/stdlib/decimals.c +++ b/src/stdlib/decimals.c @@ -3,7 +3,7 @@ #include #include -#include +#include #include #include #include @@ -17,45 +17,21 @@ #include "nums.h" #include "optionals.h" #include "print.h" -#include "siphash.h" #include "text.h" #include "types.h" -static mpd_context_t ctx = { - .prec=30, - .emax=25, - .emin=-25, - .traps=MPD_Division_by_zero | MPD_Overflow | MPD_Subnormal | MPD_Underflow, - .status=0, - .newtrap=0, - .round=MPD_ROUND_HALF_EVEN, - .clamp=1, - .allcr=1, -}; - -static inline char *Dec$as_str(Dec_t d) { - char *str = mpd_format(d, "f", &ctx); - char *point = strchr(str, '.'); - if (point == NULL) return str; - char *p; - for (p = point + strlen(point)-1; p > point && *p == '0'; p--) - *p = '\0'; - if (*p == '.') *p = '\0'; - return str; -} - public int Dec$print(FILE *f, Dec_t d) { - return fputs(Dec$as_str(d), f); + return fprint(f, d); } public Text_t Dec$value_as_text(Dec_t d) { - return Text$from_str(Dec$as_str(d)); + return Text$from_str(String(d)); } public Text_t Dec$as_text(const void *d, bool colorize, const TypeInfo_t *info) { (void)info; if (!d) return Text("Dec"); - Text_t text = Text$from_str(Dec$as_str(*(Dec_t*)d)); + Text_t text = Text$from_str(String(*(Dec_t*)d)); if (colorize) text = Text$concat(Text("\x1b[35m"), text, Text("\x1b[m")); return text; } @@ -63,84 +39,81 @@ public Text_t Dec$as_text(const void *d, bool colorize, const TypeInfo_t *info) static bool Dec$is_none(const void *d, const TypeInfo_t *info) { (void)info; - return *(Dec_t*)d == NULL; + return *(int64_t*)d == -1; } -public PUREFUNC int32_t Dec$compare_value(const Dec_t x, const Dec_t y) { - return mpd_compare(D(0), x, y, &ctx); +public CONSTFUNC int32_t Dec$compare_value(const Dec_t x, const Dec_t y) { + return (x > y) - (x < y); } -public PUREFUNC int32_t Dec$compare(const void *x, const void *y, const TypeInfo_t *info) { +public CONSTFUNC int32_t Dec$compare(const void *x, const void *y, const TypeInfo_t *info) { (void)info; - return mpd_compare(D(0), *(Dec_t*)x, *(Dec_t*)y, &ctx); + return Dec$compare_value(*(Dec_t*)x, *(Dec_t*)y); } -public PUREFUNC bool Dec$equal_value(const Dec_t x, const Dec_t y) { - return Dec$compare_value(x, y) == 0; -} - -public PUREFUNC bool Dec$equal(const void *x, const void *y, const TypeInfo_t *info) { - (void)info; - return Dec$compare(x, y, info) == 0; +public CONSTFUNC bool Dec$equal_value(const Dec_t x, const Dec_t y) { + return x == y; } -public PUREFUNC uint64_t Dec$hash(const void *vx, const TypeInfo_t *info) { +public CONSTFUNC bool Dec$equal(const void *x, const void *y, const TypeInfo_t *info) { (void)info; - Dec_t d = *(Dec_t*)vx; - char *str = Dec$as_str(d); - return siphash24((void*)str, strlen(str)); + return *(_Decimal64*)x == *(_Decimal64*)y; } -public Dec_t Dec$plus(Dec_t x, Dec_t y) { - Dec_t result = mpd_new(&ctx); - mpd_add(result, x, y, &ctx); - return result; +public CONSTFUNC Dec_t Dec$plus(Dec_t x, Dec_t y) { + return x + y; } -public Dec_t Dec$negative(Dec_t x) { - Dec_t result = mpd_new(&ctx); - mpd_minus(result, x, &ctx); - return result; +public CONSTFUNC Dec_t Dec$negative(Dec_t x) { + return -x; } -public Dec_t Dec$minus(Dec_t x, Dec_t y) { - Dec_t result = mpd_new(&ctx); - mpd_sub(result, x, y, &ctx); - return result; +public CONSTFUNC Dec_t Dec$minus(Dec_t x, Dec_t y) { + return x - y; } -public Dec_t Dec$times(Dec_t x, Dec_t y) { - Dec_t result = mpd_new(&ctx); - mpd_mul(result, x, y, &ctx); - return result; +public CONSTFUNC Dec_t Dec$times(Dec_t x, Dec_t y) { + return x * y; } -public Dec_t Dec$divided_by(Dec_t x, Dec_t y) { - Dec_t result = mpd_new(&ctx); - mpd_div(result, x, y, &ctx); - return result; +public CONSTFUNC Dec_t Dec$divided_by(Dec_t x, Dec_t y) { + return x / y; } -public Dec_t Dec$modulo(Dec_t x, Dec_t modulus) { - Dec_t result = mpd_new(&ctx); - mpd_rem(result, x, modulus, &ctx); - return result; +public CONSTFUNC Dec_t Dec$modulo(Dec_t x, Dec_t modulus) { + // TODO: improve the accuracy of this approach: + return (Dec_t)Num$mod((double)x, (double)modulus); } -public Dec_t Dec$modulo1(Dec_t x, Dec_t modulus) { - return Dec$plus(Dec$modulo(Dec$minus(x, D(1)), modulus), D(1)); +public CONSTFUNC Dec_t Dec$modulo1(Dec_t x, Dec_t modulus) { + // TODO: improve the accuracy of this approach: + return (Dec_t)Num$mod1((double)x, (double)modulus); } -public Dec_t Dec$from_str(const char *str) { - Dec_t result = mpd_new(&ctx); - mpd_set_string(result, str, &ctx); - return result; +public PUREFUNC OptionalDec_t Dec$from_str(const char *str) { + _Decimal64 n = 0.0DD; + const char *p = str; + bool negative = (*p == '-'); + if (negative) + p += 1; + for (; *p; p++) { + if (*p == '.') break; + if (*p == '_') continue; + if (!isdigit(*p)) return NONE_DEC; + n = 10.0DD * n + (_Decimal64)(*p - '0'); + } + _Decimal64 denominator = 1.0DD; + for (; *p; p++) { + if (*p == '_') continue; + if (!isdigit(*p)) return NONE_DEC; + n = 10.0DD * n + (_Decimal64)(*p - '0'); + denominator *= 0.1DD; + } + return n * denominator; } -public Dec_t Dec$from_int64(int64_t i) { - Dec_t result = mpd_new(&ctx); - mpd_set_i64(result, i, &ctx); - return result; +public CONSTFUNC Dec_t Dec$from_int64(int64_t i) { + return (_Decimal64)i; } public Dec_t Dec$from_int(Int_t i) { @@ -149,21 +122,15 @@ public Dec_t Dec$from_int(Int_t i) { } Text_t text = Int$value_as_text(i); const char *str = Text$as_c_string(text); - Dec_t result = mpd_new(&ctx); - mpd_set_string(result, str, &ctx); - return result; + return Dec$from_str(str); } -public Dec_t Dec$from_num(double n) { - Text_t text = Num$as_text(&n, false, &Num$info); - const char *str = Text$as_c_string(text); - Dec_t result = mpd_new(&ctx); - mpd_set_string(result, str, &ctx); - return result; +CONSTFUNC public Dec_t Dec$from_num(double n) { + return (_Decimal64)n; } public Int_t Dec$as_int(Dec_t d, bool truncate) { - char *str = Dec$as_str(d); + char *str = String(d); char *decimal = strchr(str, '.'); if (!truncate && decimal) fail("Could not convert to an integer without truncation: ", str); @@ -191,43 +158,106 @@ public Byte_t Dec$as_byte(Dec_t d, bool truncate) { return Byte$from_int(Dec$as_int(d, truncate), truncate); } -public bool Dec$as_bool(Dec_t d) { - return !mpd_iszero(d); +CONSTFUNC public bool Dec$as_bool(Dec_t d) { + return d != 0.0DD; } public double Dec$as_num(Dec_t d) { - const char *str = Dec$as_str(d); + const char *str = String(d); return strtod(str, NULL); } -public Dec_t Dec$power(Dec_t base, Dec_t exponent) { - Dec_t result = mpd_new(&ctx); - mpd_pow(result, base, exponent, &ctx); - return result; +#define NAN_MASK 0x7C00000000000000UL +#define INF_MASK 0x7800000000000000UL +#define DEC_BITS(n) ((union { uint64_t bits; _Decimal64 d; }){.d=n}).bits + +static bool Dec$isfinite(Dec_t d) { + uint64_t bits = DEC_BITS(d); + return (((bits & NAN_MASK) != NAN_MASK) && + ((bits & INF_MASK) != INF_MASK)); +} + +static bool Dec$isnan(Dec_t d) { + uint64_t bits = DEC_BITS(d); + return ((bits & NAN_MASK) == NAN_MASK); +} + +CONSTFUNC static Dec_t Dec$int_power(Dec_t x, int64_t exponent) +{ + if (exponent == 0) { + return 1.DD; + } else if (exponent == 1) { + return x; + } else if (exponent % 2 == 0) { + Dec_t y = Dec$int_power(x, exponent/2); + return y*y; + } else { + return x * Dec$int_power(x, exponent - 1); + } +} + +public Dec_t Dec$power(Dec_t x, Dec_t y) { + if (x == 0.DD && y < 0.DD) + fail("The following math operation is not supported: ", x, "^", y); + + /* For any y, including a NaN. */ + if (x == 1.DD) + return x; + + if (Dec$isnan(x) || Dec$isnan(y)) + return NONE_DEC; + + if (y == 0.DD) + return 1.DD; + + if (x < 0.DD && y < 0.DD) { + return NONE_DEC; + } else if (x == 0.DD) { + return y < 0.DD ? NONE_DEC : 0.DD; + } else if (!Dec$isfinite(x)) { + return y < 0.DD ? 0.DD : x; + } + + int64_t int_y = (int64_t)y; + if ((Dec_t)int_y == y) + return Dec$int_power(x, int_y); + + // TODO: improve the accuracy of this approach: + return (Dec_t)powl((long double)x, (long double)y); } public Dec_t Dec$round(Dec_t d, Int_t digits) { - Dec_t result = mpd_new(&ctx); - if (digits.small != 1L) - d = Dec$times(d, Dec$power(D(10), Dec$from_int(digits))); - mpd_round_to_int(result, d, &ctx); - if (digits.small != 1L) - result = Dec$divided_by(result, Dec$power(D(10), Dec$from_int(digits))); - return result; + int64_t digits64 = Int64$from_int(digits, false); + if (digits.small != 1L) { + for (int64_t i = digits64; i > 0; i--) + d *= 10.0DD; + for (int64_t i = digits64; i < 0; i++) + d *= 0.1DD; + } + _Decimal64 truncated = (_Decimal64)(int64_t)d; + _Decimal64 difference = (d - truncated); + _Decimal64 rounded; + if (difference < 0.0DD) { + rounded = (difference < -0.5DD) ? truncated - 1.0DD : truncated; + } else { + rounded = (difference >= 0.5DD) ? truncated + 1.0DD : truncated; + } + for (int64_t i = digits64; i > 0; i--) + rounded *= 0.1DD; + for (int64_t i = digits64; i < 0; i++) + rounded *= 10.0DD; + return rounded; } public OptionalDec_t Dec$parse(Text_t text) { - Dec_t result = mpd_new(&ctx); - uint32_t status = 0; - mpd_qset_string(result, Text$as_c_string(text), &ctx, &status); - return status == 0 ? result : NONE_DEC; + return Dec$from_str(Text$as_c_string(text)); } static void Dec$serialize(const void *obj, FILE *out, Table_t *pointers, const TypeInfo_t *info) { (void)info; Dec_t d = *(Dec_t*)obj; - char *str = Dec$as_str(d); + char *str = String(d); int64_t len = (int64_t)strlen(str); Int64$serialize(&len, out, pointers, &Int64$info); if (fwrite(str, sizeof(char), (size_t)len, out) != (size_t)len) @@ -253,7 +283,6 @@ public const TypeInfo_t Dec$info = { .metamethods={ .compare=Dec$compare, .equal=Dec$equal, - .hash=Dec$hash, .as_text=Dec$as_text, .is_none=Dec$is_none, .serialize=Dec$serialize, diff --git a/src/stdlib/decimals.h b/src/stdlib/decimals.h index 83e69fce..d6eeb332 100644 --- a/src/stdlib/decimals.h +++ b/src/stdlib/decimals.h @@ -5,7 +5,6 @@ #include #include #include -#include #include "print.h" #include "datatypes.h" @@ -13,34 +12,30 @@ #include "types.h" #include "util.h" -#define NONE_DEC ((mpd_t*)NULL) - -// #define D(i) Dec$from_int64(i) -#define D(i) ((mpd_t[1]){MPD_STATIC|MPD_CONST_DATA, 0, 1, 1, 1, (mpd_uint_t[]){i}}) +#define NONE_DEC (((union { int64_t i; _Decimal64 d; }){.i=-1}).d) int Dec$print(FILE *f, Dec_t d); Text_t Dec$value_as_text(Dec_t d); Text_t Dec$as_text(const void *d, bool colorize, const TypeInfo_t *info); -PUREFUNC int32_t Dec$compare_value(const Dec_t x, const Dec_t y); -PUREFUNC int32_t Dec$compare(const void *x, const void *y, const TypeInfo_t *info); -PUREFUNC bool Dec$equal_value(const Dec_t x, const Dec_t y); -PUREFUNC bool Dec$equal(const void *x, const void *y, const TypeInfo_t *info); -PUREFUNC uint64_t Dec$hash(const void *vx, const TypeInfo_t *info); -Dec_t Dec$round(Dec_t d, Int_t digits); -Dec_t Dec$power(Dec_t base, Dec_t exponent); -Dec_t Dec$plus(Dec_t x, Dec_t y); -Dec_t Dec$negative(Dec_t x); -Dec_t Dec$minus(Dec_t x, Dec_t y); -Dec_t Dec$times(Dec_t x, Dec_t y); -Dec_t Dec$divided_by(Dec_t x, Dec_t y); -Dec_t Dec$modulo(Dec_t x, Dec_t modulus); -Dec_t Dec$modulo1(Dec_t x, Dec_t modulus); -Dec_t Dec$from_str(const char *str); +CONSTFUNC int32_t Dec$compare_value(const Dec_t x, const Dec_t y); +CONSTFUNC int32_t Dec$compare(const void *x, const void *y, const TypeInfo_t *info); +CONSTFUNC bool Dec$equal_value(const Dec_t x, const Dec_t y); +CONSTFUNC bool Dec$equal(const void *x, const void *y, const TypeInfo_t *info); +CONSTFUNC Dec_t Dec$round(Dec_t d, Int_t digits); +CONSTFUNC Dec_t Dec$power(Dec_t base, Dec_t exponent); +CONSTFUNC Dec_t Dec$plus(Dec_t x, Dec_t y); +CONSTFUNC Dec_t Dec$negative(Dec_t x); +CONSTFUNC Dec_t Dec$minus(Dec_t x, Dec_t y); +CONSTFUNC Dec_t Dec$times(Dec_t x, Dec_t y); +CONSTFUNC Dec_t Dec$divided_by(Dec_t x, Dec_t y); +CONSTFUNC Dec_t Dec$modulo(Dec_t x, Dec_t modulus); +CONSTFUNC Dec_t Dec$modulo1(Dec_t x, Dec_t modulus); +PUREFUNC Dec_t Dec$from_str(const char *str); OptionalDec_t Dec$parse(Text_t text); -Dec_t Dec$from_int64(int64_t i); +CONSTFUNC Dec_t Dec$from_int64(int64_t i); Dec_t Dec$from_int(Int_t i); -Dec_t Dec$from_num(double n); +CONSTFUNC Dec_t Dec$from_num(double n); #define Dec$from_num32(n) Dec$from_num((double)n) #define Dec$from_int32(i) Dec$from_int64((int64_t)i) #define Dec$from_int16(i) Dec$from_int64((int64_t)i) @@ -54,7 +49,7 @@ int32_t Dec$as_int32(Dec_t d, bool truncate); int16_t Dec$as_int16(Dec_t d, bool truncate); int8_t Dec$as_int8(Dec_t d, bool truncate); Byte_t Dec$as_byte(Dec_t d, bool truncate); -bool Dec$as_bool(Dec_t d); +CONSTFUNC bool Dec$as_bool(Dec_t d); double Dec$as_num(Dec_t d); #define Dec$as_num32(d) ((float)Dec$as_num(d)) diff --git a/src/stdlib/print.c b/src/stdlib/print.c index 3270c765..a2ef7d4e 100644 --- a/src/stdlib/print.c +++ b/src/stdlib/print.c @@ -22,7 +22,7 @@ public int _print_int(FILE *f, int64_t n) if (negative) *(p--) = '-'; - return fwrite(p + 1, sizeof(char), (size_t)(&buf[19] - p), f); + return (int)fwrite(p + 1, sizeof(char), (size_t)(&buf[19] - p), f); } public int _print_uint(FILE *f, uint64_t n) @@ -35,7 +35,7 @@ public int _print_uint(FILE *f, uint64_t n) n /= 10; } while (n > 0); - return fwrite(p + 1, sizeof(char), (size_t)(&buf[19] - p), f); + return (int)fwrite(p + 1, sizeof(char), (size_t)(&buf[19] - p), f); } public int _print_hex(FILE *f, hex_format_t hex) @@ -178,6 +178,95 @@ public int _print_char(FILE *f, char c) #undef ESC } +public int _print_decimal64(FILE *f , _Decimal64 x) +{ + union { + _Decimal64 decimal; + struct { + uint64_t mantissa:52; + uint64_t exponent:11; + bool negative:1; + }; + uint64_t bits; + } info = {.decimal = x}; + + if ((info.bits >> 58 & 0x1F) == 0x1E) + return fputs(info.negative ? "-INF" : "INF", f); + else if ((info.bits >> 58 & 0x1F) == 0x1F) + return fputs("NAN", f); + + // determine exponent e, and mantissa m + // where e and m are depend on the bits in m2 + uint64_t e; + uint64_t m; + uint64_t m2 = info.bits >> 61 & 0x3; + if (m2 == 0x3) { + e = info.bits >> 51 & 0x3FF; + m = 0x20000000000000 | (info.bits & 0x7FFFFFFFFFFFF); + } else { + e = info.bits >> 53 & 0x3FF; + m = info.bits & 0x1FFFFFFFFFFFFF; + } + + if (m == 0) return fputs("0", f); + + char buf[64] = {[63]=0}; + char *p = &buf[62]; + int64_t exponent = (int64_t)e - 398; + + uint64_t n = m; + do { + *(p--) = '0' + (n % 10); + n /= 10; + } while (n > 0); + + const char *digit_str = p + 1; + + int printed = 0; + int64_t digits = (int64_t)(&buf[63] - digit_str); + + if (info.negative) + printed += fputc('-', f); + + while (exponent < 0 && digits > 1 && digit_str[digits-1] == '0') { + digits -= 1; + exponent += 1; + } + + if (exponent >= 0) { + printed += (int)fwrite(digit_str, sizeof(char), (size_t)digits, f); + for (int64_t i = 0; i < exponent; i++) + printed += (int)fwrite("0", sizeof(char), 1, f); + } else { + int64_t digits_above_zero = MAX(digits + exponent, 0); + if (digits_above_zero > 0) { + printed += (int)fwrite(digit_str, sizeof(char), (size_t)digits_above_zero, f); + for (int64_t i = -digits; i > exponent; i--) + printed += fputc('0', f); + } else { + printed += fputc('0', f); + } + + int64_t digits_below_zero = digits - digits_above_zero; + if (digits_below_zero > 0) { + const char *rest = digit_str + digits_above_zero; + if (*rest) { + printed += fputc('.', f); + for (int64_t i = digits_below_zero; i < -exponent; i++) + printed += fputc('0', f); + printed += fputs(digit_str+digits_above_zero, f); + } + } + } + + return printed; +} + +public int _print_decimal32(FILE *f, _Decimal32 x) +{ + return _print_decimal64(f, (_Decimal64)x); +} + public int _print_quoted(FILE *f, quoted_t quoted) { #define ESC(e) "\\" e diff --git a/src/stdlib/print.h b/src/stdlib/print.h index 5ef5b6ed..5a4a56f2 100644 --- a/src/stdlib/print.h +++ b/src/stdlib/print.h @@ -83,6 +83,8 @@ int _print_double(FILE *f, double x); int _print_hex(FILE *f, hex_format_t hex); int _print_hex_double(FILE *f, hex_double_t hex); int _print_oct(FILE *f, oct_format_t oct); +int _print_decimal32(FILE *f, _Decimal32 d); +int _print_decimal64(FILE *f, _Decimal64 d); PRINT_FN _print_float(FILE *f, float x) { return _print_double(f, (double)x); } PRINT_FN _print_pointer(FILE *f, void *p) { return _print_hex(f, hex((uint64_t)p)); } PRINT_FN _print_bool(FILE *f, bool b) { return fputs(b ? "yes" : "no", f); } @@ -116,6 +118,8 @@ extern int Int$print(FILE *f, Int_t i); uint8_t: _print_uint, \ float: _print_float, \ double: _print_double, \ + _Decimal32: _print_decimal32, \ + _Decimal64: _print_decimal64, \ hex_format_t: _print_hex, \ hex_double_t: _print_hex_double, \ oct_format_t: _print_oct, \ diff --git a/src/stdlib/stdlib.c b/src/stdlib/stdlib.c index 472be526..82470c41 100644 --- a/src/stdlib/stdlib.c +++ b/src/stdlib/stdlib.c @@ -5,7 +5,6 @@ #include #include #include -#include #include #include #include @@ -71,12 +70,6 @@ static _Noreturn void fpe_handler(int sig, siginfo_t *info, void *userdata) _exit(1); } - -static void *GC_calloc(size_t n, size_t size) -{ - return GC_malloc(n*size); -} - public void tomo_init(void) { GC_INIT(); @@ -87,11 +80,6 @@ public void tomo_init(void) setlocale(LC_ALL, ""); assert(getrandom(TOMO_HASH_KEY, sizeof(TOMO_HASH_KEY), 0) == sizeof(TOMO_HASH_KEY)); - mpd_mallocfunc = GC_malloc; - mpd_callocfunc = GC_calloc; - mpd_reallocfunc = GC_realloc; - mpd_free = GC_free; - struct sigaction ill_sigaction; ill_sigaction.sa_sigaction = signal_handler; sigemptyset(&ill_sigaction.sa_mask); -- cgit v1.2.3