From 2f342c3be266b46c7b345ee12d69322f31294bf7 Mon Sep 17 00:00:00 2001 From: "russell@unturf.com" Date: Sat, 6 Jun 2026 20:23:37 -0400 Subject: [PATCH] =?UTF-8?q?c-tier=20bignum=20=E2=80=94=20arbitrary-precisi?= =?UTF-8?q?on=20integers=20unblock=20secp256k1=20widths?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds tagged bignum support alongside the existing 48-bit fixnum on the C tier. Tag 6 = bignum, heap struct sign-magnitude with u64 little-endian limbs. Reader emits bignums for any literal past the fixnum range; +, -, *, quotient, remainder, modulo, expt, =, <, >, abs, odd?, even?, integer?, exact?, number->string, string->number all promote fixnum → bignum on overflow & demote back when results fit. Boehm GC owns every allocation. Schoolbook O(n²) mul + shift-subtract divmod is sufficient at our 4-limb / 256-bit scale. Before: (expt 2 48) = 0, (expt 2 256) = 0, secp256k1-p = -4294968273. After: all three return their exact arbitrary-precision values, matching Python tier byte-for-byte. Validated: - c/test.c — 85/85 pass (+2 new bignum unit tests). - tests/functional.lsp — 205/205 pass on both C & Python tiers. - tests/bignum-cross-tier.lsp — 33/33 pass byte-identical on both tiers (diff produces no output). - ecdsa/runs/lumbda-sweep-003/c-tier-bignum-probe.lsp — all four assertions now match the Python oracle. - ecdsa Phase B byte-identity sweep inside QEMU guest: n+1=9 p=251 sha256 c668bbe3... — matches Python oracle. n+1=18 p=131071 sha256 8a031f96... — matches Python oracle. n+1=33 p=2³²-5 sha256 0bc56905... — matches Python oracle. Previously the n+1=33 C tier emitted sha256 b024d6d9... (26,078 fewer Toffolis due to silent fixnum wrap). Bignums close that gate. secp256k1 production-width emit (n+1=257) is now structurally unblocked on C tier; downstream agent (#55) drives that next-step on the ecdsa side. Asm tier inherits in a follow-up port. --- c/Makefile | 2 +- c/bignum.c | 496 ++++++++++++++++++++++++++++++++++++ c/builtins.c | 65 ++++- c/eval.c | 6 +- c/lumbda.h | 77 +++++- c/printer.c | 14 + c/reader.c | 23 +- c/test.c | 47 ++++ c/types.c | 69 ++++- c/vm.c | 2 +- tests/bignum-cross-tier.lsp | 146 +++++++++++ 11 files changed, 918 insertions(+), 29 deletions(-) create mode 100644 c/bignum.c create mode 100644 tests/bignum-cross-tier.lsp diff --git a/c/Makefile b/c/Makefile index fa21583..15a9c1e 100644 --- a/c/Makefile +++ b/c/Makefile @@ -17,7 +17,7 @@ CFLAGS += -DUSE_BOEHM_GC LDFLAGS += -lgc endif -SRCS = types.c reader.c printer.c eval.c builtins.c vm.c jit.c portal.c +SRCS = types.c bignum.c reader.c printer.c eval.c builtins.c vm.c jit.c portal.c OBJS = $(SRCS:.c=.o) .PHONY: all clean test bench diff --git a/c/bignum.c b/c/bignum.c new file mode 100644 index 0000000..808b827 --- /dev/null +++ b/c/bignum.c @@ -0,0 +1,496 @@ +/* + * bignum.c — Arbitrary-precision integer arithmetic for lumbda C tier. + * + * Design (full notes in lumbda.h after the Rational block): + * - Sign-magnitude representation: int32_t sign (-1, 0, +1), + * uint32_t n_limbs, uint64_t *limbs (little-endian magnitude). + * - bignum_normalize() trims leading-zero limbs & demotes to fixnum + * when |value| <= FIXNUM_MAX so callers never branch on tag. + * - All arithmetic uses 64-bit limbs with 128-bit products (__int128) + * for schoolbook multiply. Division uses long-division by single + * limb when divisor fits one limb, otherwise schoolbook with + * normalization (sufficient for our scale: secp256k1 needs at most + * 4 limbs / 256 bits, far below where Karatsuba helps). + * - Boehm GC owns the Bignum struct & its limbs array. + * + * Cross-tier oracle: Python tier (native arbitrary int) is the byte- + * identity standard. Every operation here must agree with Python on + * matching inputs — see tests/bignum-cross-tier.lsp. + */ +#include "lumbda.h" + +/* ═══════════════════════════════════════════════════════════════════════════ + * Allocation helpers + * ═══════════════════════════════════════════════════════════════════════════ */ + +static Bignum *alloc_bignum(uint32_t n_limbs) { + Bignum *b = (Bignum *)ul_malloc(sizeof(Bignum)); + b->hdr.type = OBJ_BIGNUM; + b->sign = 0; + b->n_limbs = n_limbs; + if (n_limbs == 0) { + b->limbs = NULL; + } else { + b->limbs = (uint64_t *)ul_malloc(sizeof(uint64_t) * n_limbs); + memset(b->limbs, 0, sizeof(uint64_t) * n_limbs); + } + return b; +} + +static Value wrap_bignum(Bignum *b) { + return NANBOX(TAG_BIGNUM, (uintptr_t)b); +} + +/* Trim leading zero limbs in place. Zero magnitude => sign 0. */ +static void trim(Bignum *b) { + while (b->n_limbs > 0 && b->limbs[b->n_limbs - 1] == 0) + b->n_limbs--; + if (b->n_limbs == 0) b->sign = 0; +} + +/* ═══════════════════════════════════════════════════════════════════════════ + * Constructors + * ═══════════════════════════════════════════════════════════════════════════ */ + +Value make_bignum_from_i64(int64_t v) { + if (FITS_FIXNUM(v)) return VAL_INT(v); + Bignum *b = alloc_bignum(1); + if (v < 0) { + b->sign = -1; + /* Avoid -INT64_MIN UB: cast then negate via unsigned. */ + b->limbs[0] = (uint64_t)(-(v + 1)) + 1ULL; + } else if (v > 0) { + b->sign = 1; + b->limbs[0] = (uint64_t)v; + } else { + b->sign = 0; + b->n_limbs = 0; + } + return wrap_bignum(b); +} + +Value bignum_from_limbs(int32_t sign, uint64_t *limbs, uint32_t n) { + Bignum *b = alloc_bignum(n); + if (n > 0) memcpy(b->limbs, limbs, sizeof(uint64_t) * n); + b->sign = sign; + trim(b); + return wrap_bignum(b); +} + +/* Demote bignum-shape to fixnum when |value| fits 48 signed bits. + * On demote we return VAL_INT(...); otherwise allocate a fresh Bignum. */ +Value bignum_normalize(int32_t sign, uint64_t *limbs, uint32_t n) { + /* Strip top zeros (compute effective limb count). */ + while (n > 0 && limbs[n - 1] == 0) n--; + if (n == 0) return VAL_INT(0); + if (n == 1) { + uint64_t lo = limbs[0]; + /* Signed 48-bit range: -2^47 .. 2^47 - 1. */ + if (sign >= 0 && lo <= (uint64_t)FIXNUM_MAX) + return VAL_INT((int64_t)lo); + if (sign < 0 && lo <= (uint64_t)(1ULL << 47)) + return VAL_INT(-(int64_t)lo); + } + Bignum *b = alloc_bignum(n); + memcpy(b->limbs, limbs, sizeof(uint64_t) * n); + b->sign = sign; + return wrap_bignum(b); +} + +/* Promote a fixnum to a Bignum struct (caller-owned). Bignum input passes + * through (the same struct is returned — Bignum structs are immutable from + * the perspective of arithmetic so sharing is safe). */ +Bignum *as_bignum(Value v) { + if (IS_BIGNUM(v)) return AS_BIGNUM(v); + if (IS_INT(v)) { + int64_t x = as_int(v); + Bignum *b = alloc_bignum(x == 0 ? 0 : 1); + if (x > 0) { b->sign = 1; b->limbs[0] = (uint64_t)x; } + else if (x < 0) { + b->sign = -1; + b->limbs[0] = (uint64_t)(-(x + 1)) + 1ULL; + } + return b; + } + if (IS_DOUBLE(v)) { + double d = as_double(v); + return AS_BIGNUM(make_bignum_from_i64((int64_t)d)); + } + lisp_error("not an integer"); + return NULL; +} + +Value to_bignum(Value v) { + if (IS_BIGNUM(v)) return v; + return wrap_bignum(as_bignum(v)); +} + +/* ═══════════════════════════════════════════════════════════════════════════ + * Parsing (string → bignum) + * ═══════════════════════════════════════════════════════════════════════════ */ + +/* Multiply magnitude by single limb; add another limb. Magnitude grows. */ +static void mag_mul_add(uint64_t *limbs, uint32_t *n, uint64_t mul, uint64_t add) { + __uint128_t carry = add; + for (uint32_t i = 0; i < *n; i++) { + __uint128_t prod = (__uint128_t)limbs[i] * mul + carry; + limbs[i] = (uint64_t)prod; + carry = prod >> 64; + } + while (carry) { + limbs[(*n)++] = (uint64_t)carry; + carry >>= 64; + } +} + +Value make_bignum_from_str(const char *s, int base) { + int32_t sign = 1; + if (*s == '-') { sign = -1; s++; } + else if (*s == '+') s++; + /* Conservative limb count: log2(10) ≈ 3.33 bits per digit → roughly + * len/19 limbs for base 10; len/16 limbs for base 16. Round up + 1. */ + size_t len = strlen(s); + uint32_t cap = (uint32_t)((len * 4) / 19) + 2; + uint64_t *limbs = (uint64_t *)ul_malloc(sizeof(uint64_t) * cap); + memset(limbs, 0, sizeof(uint64_t) * cap); + uint32_t n = 0; + for (; *s; s++) { + int d; + if (*s >= '0' && *s <= '9') d = *s - '0'; + else if (*s >= 'a' && *s <= 'f') d = *s - 'a' + 10; + else if (*s >= 'A' && *s <= 'F') d = *s - 'A' + 10; + else { ul_free(limbs); lisp_error("bad digit '%c' in number", *s); } + if (d >= base) { ul_free(limbs); lisp_error("digit out of base"); } + mag_mul_add(limbs, &n, (uint64_t)base, (uint64_t)d); + if (n >= cap - 1) { + cap *= 2; + limbs = (uint64_t *)ul_realloc(limbs, sizeof(uint64_t) * cap); + memset(limbs + n, 0, sizeof(uint64_t) * (cap - n)); + } + } + Value v = bignum_normalize(sign, limbs, n); + ul_free(limbs); + return v; +} + +/* ═══════════════════════════════════════════════════════════════════════════ + * Magnitude-level primitives + * ═══════════════════════════════════════════════════════════════════════════ */ + +/* Compare magnitudes. Returns -1, 0, +1. */ +static int mag_cmp(const uint64_t *a, uint32_t na, const uint64_t *b, uint32_t nb) { + if (na != nb) return na < nb ? -1 : 1; + for (int32_t i = (int32_t)na - 1; i >= 0; i--) { + if (a[i] != b[i]) return a[i] < b[i] ? -1 : 1; + } + return 0; +} + +/* out = a + b magnitudes. out_cap must be at least max(na, nb) + 1. */ +static uint32_t mag_add(const uint64_t *a, uint32_t na, + const uint64_t *b, uint32_t nb, + uint64_t *out) { + uint32_t n = na > nb ? na : nb; + __uint128_t carry = 0; + for (uint32_t i = 0; i < n; i++) { + __uint128_t s = carry; + if (i < na) s += a[i]; + if (i < nb) s += b[i]; + out[i] = (uint64_t)s; + carry = s >> 64; + } + if (carry) out[n++] = (uint64_t)carry; + return n; +} + +/* out = a - b magnitudes, assuming a >= b. out_cap >= na. */ +static uint32_t mag_sub(const uint64_t *a, uint32_t na, + const uint64_t *b, uint32_t nb, + uint64_t *out) { + int64_t borrow = 0; + for (uint32_t i = 0; i < na; i++) { + __int128 s = (__int128)a[i] - borrow; + if (i < nb) s -= b[i]; + if (s < 0) { s += ((__int128)1 << 64); borrow = 1; } + else borrow = 0; + out[i] = (uint64_t)s; + } + while (na > 0 && out[na - 1] == 0) na--; + return na; +} + +/* out = a * b magnitudes. out_cap >= na + nb. */ +static uint32_t mag_mul(const uint64_t *a, uint32_t na, + const uint64_t *b, uint32_t nb, + uint64_t *out) { + memset(out, 0, sizeof(uint64_t) * (na + nb)); + for (uint32_t i = 0; i < na; i++) { + __uint128_t carry = 0; + for (uint32_t j = 0; j < nb; j++) { + __uint128_t prod = (__uint128_t)a[i] * b[j] + out[i + j] + carry; + out[i + j] = (uint64_t)prod; + carry = prod >> 64; + } + out[i + nb] += (uint64_t)carry; + } + uint32_t n = na + nb; + while (n > 0 && out[n - 1] == 0) n--; + return n; +} + +/* Divide magnitude a by single 64-bit limb d. Quotient into q (cap >= na), + * remainder returned. */ +static uint64_t mag_divmod_u64(const uint64_t *a, uint32_t na, + uint64_t d, uint64_t *q, uint32_t *nq) { + __uint128_t r = 0; + for (int32_t i = (int32_t)na - 1; i >= 0; i--) { + __uint128_t cur = (r << 64) | a[i]; + q[i] = (uint64_t)(cur / d); + r = cur % d; + } + *nq = na; + while (*nq > 0 && q[*nq - 1] == 0) (*nq)--; + return (uint64_t)r; +} + +/* Schoolbook division: a / b → (q, r) where a = q*b + r, 0 <= r < b. + * For our scale (operands <= 4 limbs commonly, up to a few dozen for + * intermediates in expt), a basic Knuth Algorithm D is overkill — we use + * shift-subtract on bits. O(na * nb * 64). */ +static void mag_divmod(const uint64_t *a, uint32_t na, + const uint64_t *b, uint32_t nb, + uint64_t *q, uint32_t *nq, + uint64_t *r, uint32_t *nr) { + if (nb == 0) lisp_error("division by zero"); + if (mag_cmp(a, na, b, nb) < 0) { + /* a < b: q=0, r=a */ + *nq = 0; + memcpy(r, a, sizeof(uint64_t) * na); + *nr = na; + return; + } + if (nb == 1) { + uint64_t rem = mag_divmod_u64(a, na, b[0], q, nq); + if (rem) { r[0] = rem; *nr = 1; } else { *nr = 0; } + return; + } + /* Shift-subtract long division: walk bits MSB → LSB of dividend. */ + memset(q, 0, sizeof(uint64_t) * na); + memset(r, 0, sizeof(uint64_t) * (nb + 1)); + uint32_t rn = 0; + for (int64_t bit = (int64_t)na * 64 - 1; bit >= 0; bit--) { + /* r <<= 1 */ + uint64_t carry = 0; + for (uint32_t i = 0; i <= nb; i++) { + uint64_t new_carry = r[i] >> 63; + r[i] = (r[i] << 1) | carry; + carry = new_carry; + } + if (carry || r[nb]) rn = nb + 1; + else { while (rn > 0 && r[rn - 1] == 0) rn--; } + /* r |= bit of a */ + uint64_t bit_val = (a[bit / 64] >> (bit % 64)) & 1ULL; + r[0] |= bit_val; + if (bit_val && rn == 0) rn = 1; + /* if r >= b: r -= b; set bit in q */ + if (mag_cmp(r, rn, b, nb) >= 0) { + uint32_t new_rn = mag_sub(r, rn, b, nb, r); + rn = new_rn; + q[bit / 64] |= 1ULL << (bit % 64); + } + } + *nq = na; + while (*nq > 0 && q[*nq - 1] == 0) (*nq)--; + *nr = rn; +} + +/* ═══════════════════════════════════════════════════════════════════════════ + * Signed arithmetic (sign-magnitude wrappers) + * ═══════════════════════════════════════════════════════════════════════════ */ + +Value big_add(Value a, Value b) { + Bignum *ba = as_bignum(a), *bb = as_bignum(b); + if (ba->sign == 0) return bignum_normalize(bb->sign, bb->limbs, bb->n_limbs); + if (bb->sign == 0) return bignum_normalize(ba->sign, ba->limbs, ba->n_limbs); + uint32_t cap = (ba->n_limbs > bb->n_limbs ? ba->n_limbs : bb->n_limbs) + 1; + uint64_t *out = (uint64_t *)ul_malloc(sizeof(uint64_t) * cap); + Value res; + if (ba->sign == bb->sign) { + uint32_t n = mag_add(ba->limbs, ba->n_limbs, bb->limbs, bb->n_limbs, out); + res = bignum_normalize(ba->sign, out, n); + } else { + int cmp = mag_cmp(ba->limbs, ba->n_limbs, bb->limbs, bb->n_limbs); + if (cmp == 0) { res = VAL_INT(0); } + else if (cmp > 0) { + uint32_t n = mag_sub(ba->limbs, ba->n_limbs, bb->limbs, bb->n_limbs, out); + res = bignum_normalize(ba->sign, out, n); + } else { + uint32_t n = mag_sub(bb->limbs, bb->n_limbs, ba->limbs, ba->n_limbs, out); + res = bignum_normalize(bb->sign, out, n); + } + } + ul_free(out); + return res; +} + +Value big_neg(Value a) { + if (IS_INT(a)) { + int64_t x = as_int(a); + return make_bignum_from_i64(-x); + } + Bignum *ba = as_bignum(a); + if (ba->sign == 0) return VAL_INT(0); + return bignum_normalize(-ba->sign, ba->limbs, ba->n_limbs); +} + +Value big_sub(Value a, Value b) { + return big_add(a, big_neg(b)); +} + +Value big_mul(Value a, Value b) { + Bignum *ba = as_bignum(a), *bb = as_bignum(b); + if (ba->sign == 0 || bb->sign == 0) return VAL_INT(0); + uint32_t cap = ba->n_limbs + bb->n_limbs; + uint64_t *out = (uint64_t *)ul_malloc(sizeof(uint64_t) * cap); + uint32_t n = mag_mul(ba->limbs, ba->n_limbs, bb->limbs, bb->n_limbs, out); + Value res = bignum_normalize(ba->sign * bb->sign, out, n); + ul_free(out); + return res; +} + +bool big_is_zero(Value v) { + if (IS_INT(v)) return as_int(v) == 0; + if (IS_BIGNUM(v)) return AS_BIGNUM(v)->sign == 0; + return false; +} + +bool big_is_odd(Value v) { + if (IS_INT(v)) return (as_int(v) & 1) != 0; + if (IS_BIGNUM(v)) { + Bignum *b = AS_BIGNUM(v); + return b->n_limbs > 0 && (b->limbs[0] & 1ULL) != 0; + } + return false; +} + +Value big_abs(Value v) { + if (IS_INT(v)) { int64_t x = as_int(v); return make_bignum_from_i64(x < 0 ? -x : x); } + Bignum *b = as_bignum(v); + if (b->sign >= 0) return v; + return bignum_normalize(1, b->limbs, b->n_limbs); +} + +int big_cmp(Value a, Value b) { + Bignum *ba = as_bignum(a), *bb = as_bignum(b); + if (ba->sign != bb->sign) return ba->sign < bb->sign ? -1 : 1; + if (ba->sign == 0) return 0; + int m = mag_cmp(ba->limbs, ba->n_limbs, bb->limbs, bb->n_limbs); + return ba->sign > 0 ? m : -m; +} + +/* Quotient: truncated toward zero (matches C / Python's // for non-negatives; + * Scheme's `quotient` is also truncate-toward-zero). */ +Value big_quotient(Value a, Value b) { + Bignum *ba = as_bignum(a), *bb = as_bignum(b); + if (bb->sign == 0) lisp_error("division by zero"); + if (ba->sign == 0) return VAL_INT(0); + uint32_t qcap = ba->n_limbs > 0 ? ba->n_limbs : 1; + uint32_t rcap = bb->n_limbs + 1; + uint64_t *q = (uint64_t *)ul_malloc(sizeof(uint64_t) * qcap); + uint64_t *r = (uint64_t *)ul_malloc(sizeof(uint64_t) * rcap); + uint32_t nq, nr; + mag_divmod(ba->limbs, ba->n_limbs, bb->limbs, bb->n_limbs, q, &nq, r, &nr); + int32_t sign = ba->sign * bb->sign; + Value res = bignum_normalize(nq == 0 ? 0 : sign, q, nq); + ul_free(q); ul_free(r); + return res; +} + +/* Remainder: result has sign of dividend (Scheme `remainder`). */ +Value big_remainder(Value a, Value b) { + Bignum *ba = as_bignum(a), *bb = as_bignum(b); + if (bb->sign == 0) lisp_error("division by zero"); + if (ba->sign == 0) return VAL_INT(0); + uint32_t qcap = ba->n_limbs > 0 ? ba->n_limbs : 1; + uint32_t rcap = bb->n_limbs + 1; + uint64_t *q = (uint64_t *)ul_malloc(sizeof(uint64_t) * qcap); + uint64_t *r = (uint64_t *)ul_malloc(sizeof(uint64_t) * rcap); + uint32_t nq, nr; + mag_divmod(ba->limbs, ba->n_limbs, bb->limbs, bb->n_limbs, q, &nq, r, &nr); + Value res = bignum_normalize(nr == 0 ? 0 : ba->sign, r, nr); + ul_free(q); ul_free(r); + return res; +} + +/* Modulo: result has sign of divisor (Scheme `modulo`). */ +Value big_modulo(Value a, Value b) { + Value r = big_remainder(a, b); + /* If r != 0 and sign(r) != sign(b), add b. */ + if (big_is_zero(r)) return r; + Bignum *br = as_bignum(r), *bb = as_bignum(b); + if (br->sign != bb->sign) return big_add(r, b); + return r; +} + +/* Binary exponentiation. exp must be a non-negative integer. */ +Value big_expt(Value base, Value exp) { + /* exp must be an integer. */ + if (!IS_INTEGER(exp)) lisp_error("expt: integer exponent required"); + /* Negative exponent → not handled here (caller falls back to double). */ + if (IS_INT(exp) && as_int(exp) < 0) lisp_error("expt: negative exponent"); + if (IS_BIGNUM(exp) && AS_BIGNUM(exp)->sign < 0) lisp_error("expt: negative exponent"); + /* Read exp as u64 — bignum exponents larger than 2^63 would never finish. */ + uint64_t e; + if (IS_INT(exp)) e = (uint64_t)as_int(exp); + else { + Bignum *be = AS_BIGNUM(exp); + if (be->n_limbs > 1) lisp_error("expt: exponent too large"); + e = be->n_limbs == 0 ? 0 : be->limbs[0]; + } + Value result = VAL_INT(1); + Value b = base; + while (e > 0) { + if (e & 1) result = big_mul(result, b); + e >>= 1; + if (e) b = big_mul(b, b); + } + return result; +} + +/* ═══════════════════════════════════════════════════════════════════════════ + * Decimal rendering + * ═══════════════════════════════════════════════════════════════════════════ */ + +char *bignum_to_str(Bignum *b) { + if (b->sign == 0) return ul_strdup("0"); + /* Divide repeatedly by 10^18 → 18 decimal digits per chunk. */ + uint32_t n = b->n_limbs; + uint64_t *mag = (uint64_t *)ul_malloc(sizeof(uint64_t) * n); + memcpy(mag, b->limbs, sizeof(uint64_t) * n); + /* Generous size: log10(2) ≈ 0.30103 → n*64*0.31 + sign + null. */ + size_t cap = (size_t)n * 21 + 4; + char *buf = (char *)ul_malloc(cap); + size_t pos = cap; + buf[--pos] = '\0'; + const uint64_t CHUNK = 1000000000000000000ULL; /* 10^18 fits in u64. */ + while (n > 0) { + uint64_t rem = mag_divmod_u64(mag, n, CHUNK, mag, &n); + if (n > 0) { + /* Full 18-digit chunk (zero-pad). */ + for (int i = 0; i < 18; i++) { + buf[--pos] = '0' + (int)(rem % 10); + rem /= 10; + } + } else { + /* Top chunk: no leading zeros. */ + if (rem == 0) buf[--pos] = '0'; + else while (rem) { buf[--pos] = '0' + (int)(rem % 10); rem /= 10; } + } + } + if (b->sign < 0) buf[--pos] = '-'; + /* Move string to start. */ + size_t len = cap - pos; + memmove(buf, buf + pos, len); + ul_free(mag); + return buf; +} diff --git a/c/builtins.c b/c/builtins.c index 5fac107..a702a46 100644 --- a/c/builtins.c +++ b/c/builtins.c @@ -113,11 +113,13 @@ static Value bi_atan(Value *a, int n, Env *e) { static Value bi_quotient(Value *a, int n, Env *e) { (void)e; CHECK_ARITY("quotient", 2); + if (IS_INTEGER(a[0]) && IS_INTEGER(a[1])) return big_quotient(a[0], a[1]); return VAL_INT(as_number_int(a[0]) / as_number_int(a[1])); } static Value bi_remainder(Value *a, int n, Env *e) { (void)e; CHECK_ARITY("remainder", 2); + if (IS_INTEGER(a[0]) && IS_INTEGER(a[1])) return big_remainder(a[0], a[1]); int64_t x = as_number_int(a[0]), y = as_number_int(a[1]); int64_t r = x % y; /* Remainder has sign of dividend */ @@ -126,6 +128,7 @@ static Value bi_remainder(Value *a, int n, Env *e) { static Value bi_modulo(Value *a, int n, Env *e) { (void)e; CHECK_ARITY("modulo", 2); + if (IS_INTEGER(a[0]) && IS_INTEGER(a[1])) return big_modulo(a[0], a[1]); int64_t x = as_number_int(a[0]), y = as_number_int(a[1]); int64_t r = x % y; if ((r > 0 && y < 0) || (r < 0 && y > 0)) r += y; @@ -134,18 +137,25 @@ static Value bi_modulo(Value *a, int n, Env *e) { static Value bi_expt(Value *a, int n, Env *e) { (void)e; CHECK_ARITY("expt", 2); - if (IS_INT(a[0]) && IS_INT(a[1]) && as_int(a[1]) >= 0) { - int64_t base = as_int(a[0]), exp = as_int(a[1]); - int64_t result = 1; - for (int64_t i = 0; i < exp; i++) result *= base; - return VAL_INT(result); + /* Integer-base / non-negative integer-exp → bignum-aware path. + * This catches (expt 2 256) and friends — fixnum * fixnum would + * silently overflow inside the old int64 loop. */ + if (IS_INTEGER(a[0]) && IS_INTEGER(a[1])) { + bool nonneg = (IS_INT(a[1]) && as_int(a[1]) >= 0) || + (IS_BIGNUM(a[1]) && AS_BIGNUM(a[1])->sign >= 0); + if (nonneg) return big_expt(a[0], a[1]); } return make_double(pow(as_number_double(a[0]), as_number_double(a[1]))); } static Value bi_abs(Value *a, int n, Env *e) { (void)e; CHECK_ARITY("abs", 1); - if (IS_INT(a[0])) { int64_t v = as_int(a[0]); return VAL_INT(v < 0 ? -v : v); } + if (IS_INT(a[0])) { + int64_t v = as_int(a[0]); + if (v == FIXNUM_MIN) return big_abs(a[0]); /* avoid -INT_MIN overflow */ + return VAL_INT(v < 0 ? -v : v); + } + if (IS_BIGNUM(a[0])) return big_abs(a[0]); return make_double(fabs(as_number_double(a[0]))); } @@ -231,6 +241,15 @@ static Value bi_number_to_string(Value *a, int n, Env *e) { char buf[128]; if (n > 1) { int base = (int)as_number_int(a[1]); + /* Bignum decimal path uses our renderer; other bases fall through to + * int64 — at our scale (256-bit secp256k1) base-10 covers every site + * that actually needs precision. */ + if (IS_BIGNUM(a[0]) && base == 10) { + char *s = bignum_to_str(AS_BIGNUM(a[0])); + Value r = make_string_from_cstr(s); + ul_free(s); + return r; + } int64_t val = as_number_int(a[0]); switch (base) { case 2: { @@ -280,8 +299,8 @@ NUM_CMP(num_ge, num_ge) NUM_PRED(zero_p, is_number(a[0]) && num_eq(a[0], VAL_INT(0))) NUM_PRED(positive_p, is_number(a[0]) && num_gt(a[0], VAL_INT(0))) NUM_PRED(negative_p, is_number(a[0]) && num_lt(a[0], VAL_INT(0))) -NUM_PRED(odd_p, IS_INT(a[0]) && (as_int(a[0]) % 2 != 0)) -NUM_PRED(even_p, IS_INT(a[0]) && (as_int(a[0]) % 2 == 0)) +NUM_PRED(odd_p, IS_INTEGER(a[0]) && big_is_odd(a[0])) +NUM_PRED(even_p, IS_INTEGER(a[0]) && !big_is_odd(a[0])) static Value bi_nan_p(Value *a, int n, Env *e) { (void)e; CHECK_ARITY("nan?", 1); @@ -357,10 +376,10 @@ static Value bi_equal_p(Value *a, int n, Env *e) { } TYPE_PRED(number_p, is_number(a[0])) -TYPE_PRED(integer_p, IS_INT(a[0]) || (IS_DOUBLE(a[0]) && as_double(a[0]) == floor(as_double(a[0])))) +TYPE_PRED(integer_p, IS_INT(a[0]) || IS_BIGNUM(a[0]) || (IS_DOUBLE(a[0]) && as_double(a[0]) == floor(as_double(a[0])))) TYPE_PRED(real_p, is_number(a[0])) -TYPE_PRED(rational_p, IS_INT(a[0]) || IS_RATIONAL(a[0]) || (IS_DOUBLE(a[0]) && isfinite(as_double(a[0])))) -TYPE_PRED(exact_p, IS_INT(a[0]) || IS_RATIONAL(a[0])) +TYPE_PRED(rational_p, IS_INT(a[0]) || IS_BIGNUM(a[0]) || IS_RATIONAL(a[0]) || (IS_DOUBLE(a[0]) && isfinite(as_double(a[0])))) +TYPE_PRED(exact_p, IS_INT(a[0]) || IS_BIGNUM(a[0]) || IS_RATIONAL(a[0])) TYPE_PRED(inexact_p, IS_DOUBLE(a[0])) TYPE_PRED(pair_p, IS_PAIR(a[0])) TYPE_PRED(null_p, IS_NIL(a[0])) @@ -796,6 +815,30 @@ static Value bi_string_to_number(Value *a, int n, Env *e) { (void)e; CHECK_MIN_ARITY("string->number", 1); check_string(a[0]); const char *s = AS_STRING(a[0])->data; int base = n > 1 ? (int)as_number_int(a[1]) : 10; + /* Detect all-digit (optionally signed) shape → integer (fixnum or + * bignum on overflow). */ + { + const char *p = s; + if (*p == '-' || *p == '+') p++; + bool all_digits = (*p != '\0'); + for (const char *q = p; *q; q++) { + char c = *q; + int dval; + if (c >= '0' && c <= '9') dval = c - '0'; + else if (base > 10 && c >= 'a' && c <= 'z') dval = c - 'a' + 10; + else if (base > 10 && c >= 'A' && c <= 'Z') dval = c - 'A' + 10; + else { all_digits = false; break; } + if (dval >= base) { all_digits = false; break; } + } + if (all_digits) { + char *end; + errno = 0; + long long val = strtoll(s, &end, base); + if (*end == '\0' && errno == 0 && FITS_FIXNUM(val)) + return VAL_INT(val); + return make_bignum_from_str(s, base); + } + } char *end; errno = 0; long long val = strtoll(s, &end, base); diff --git a/c/eval.c b/c/eval.c index 6b861fc..8381317 100644 --- a/c/eval.c +++ b/c/eval.c @@ -480,7 +480,7 @@ static bool is_literal(SyntaxTransformer *st, const char *name) { static bool sr_match(SyntaxTransformer *st, Value pat, Value form, Env *bindings) { if (IS_NIL(pat)) return IS_NIL(form); if (pat == VAL_TRUE || pat == VAL_FALSE) return pat == form; - if (IS_INT(pat) || IS_DOUBLE(pat)) return values_equal(pat, form); + if (IS_INT(pat) || IS_DOUBLE(pat) || IS_BIGNUM(pat)) return values_equal(pat, form); if (IS_SYM(pat)) { const char *pname = sym_name(pat); @@ -573,7 +573,7 @@ static bool sr_match(SyntaxTransformer *st, Value pat, Value form, Env *bindings static Value sr_expand(SyntaxTransformer *st, Value tmpl, Env *bindings) { if (IS_NIL(tmpl) || tmpl == VAL_TRUE || tmpl == VAL_FALSE) return tmpl; - if (IS_INT(tmpl) || IS_DOUBLE(tmpl)) return tmpl; + if (IS_INT(tmpl) || IS_DOUBLE(tmpl) || IS_BIGNUM(tmpl)) return tmpl; if (IS_STRING(tmpl)) return tmpl; if (IS_SYM(tmpl)) { @@ -669,7 +669,7 @@ Value leval(Value expr, Env *env) { /* Self-evaluating */ if (IS_NIL(expr) || IS_VOID(expr) || IS_TRUE(expr) || IS_FALSE(expr) || IS_EOF(expr) || IS_CHAR(expr)) return expr; - if (IS_INT(expr) || IS_DOUBLE(expr) || IS_RATIONAL(expr)) return expr; + if (IS_INT(expr) || IS_DOUBLE(expr) || IS_RATIONAL(expr) || IS_BIGNUM(expr)) return expr; if (IS_STRING(expr)) return expr; if (IS_VECTOR(expr)) return expr; if (IS_BUILTIN(expr)) return expr; diff --git a/c/lumbda.h b/c/lumbda.h index eebf40e..4e56b02 100644 --- a/c/lumbda.h +++ b/c/lumbda.h @@ -62,6 +62,10 @@ static inline char *ul_strdup(const char *s) { * Tag 3 = special (NIL, VOID, TRUE, FALSE, EOF, char in payload) * Tag 4 = builtin function pointer * Tag 5 = rational (pointer to Rational struct) + * Tag 6 = bignum (pointer to Bignum struct) — arbitrary-precision integer. + * Used for any integer outside the 48-bit fixnum range. Arithmetic + * promotes fixnum to bignum on overflow & demotes bignum back to + * fixnum when result fits, so callers never need to distinguish. * * Plain doubles (no NaN payload) are floating-point numbers. * ═══════════════════════════════════════════════════════════════════════════ */ @@ -81,6 +85,12 @@ typedef uint64_t Value; #define TAG_SPECIAL 3ULL #define TAG_BUILTIN 4ULL #define TAG_RATIONAL 5ULL +#define TAG_BIGNUM 6ULL + +/* Fixnum range: signed 48-bit (low 48 bits of payload). */ +#define FIXNUM_MIN ((int64_t)(-(1LL << 47))) +#define FIXNUM_MAX ((int64_t)((1LL << 47) - 1)) +#define FITS_FIXNUM(v) ((v) >= FIXNUM_MIN && (v) <= FIXNUM_MAX) /* Construct a NaN-boxed value */ #define NANBOX(tag, payload) (QNAN | ((uint64_t)(tag) << TAG_SHIFT) | ((uint64_t)(payload) & PAYLOAD_MASK)) @@ -137,6 +147,9 @@ static inline int64_t as_int(Value v) { #define IS_SPECIAL(v) (!IS_DOUBLE(v) && GET_TAG(v) == TAG_SPECIAL) #define IS_BUILTIN(v) (!IS_DOUBLE(v) && GET_TAG(v) == TAG_BUILTIN) #define IS_RATIONAL(v) (!IS_DOUBLE(v) && GET_TAG(v) == TAG_RATIONAL) +#define IS_BIGNUM(v) (!IS_DOUBLE(v) && GET_TAG(v) == TAG_BIGNUM) +/* Integer-shaped: either fixnum or bignum. */ +#define IS_INTEGER(v) (IS_INT(v) || IS_BIGNUM(v)) #define IS_NIL(v) ((v) == VAL_NIL) #define IS_VOID(v) ((v) == VAL_VOID) @@ -170,6 +183,7 @@ typedef enum { OBJ_ERROR, OBJ_SYNTAX_TRANSFORMER, OBJ_RATIONAL, + OBJ_BIGNUM, } ObjType; typedef struct ObjHeader { @@ -288,6 +302,65 @@ Value make_rational(int64_t num, int64_t den); /* Returns integer Value if den==1, otherwise rational */ Value rational_normalize(int64_t num, int64_t den); +/* ── Bignum — arbitrary-precision integer ────────────────────────────────── */ + +/* + * Bignum representation: sign-magnitude. + * sign: -1 (negative), 0 (zero — must have n_limbs == 0), +1 (positive) + * n_limbs: number of u64 magnitude limbs (always trimmed of leading zeros) + * limbs[i]: little-endian magnitude (limb 0 = least-significant 64 bits) + * + * Boehm GC owns every Bignum struct & its limbs array. Construction goes + * through bignum_normalize() which trims leading-zero limbs & forces zero + * to canonical form (sign=0, n_limbs=0). + * + * Result-producing helpers ALWAYS return a Value: if the magnitude fits in + * 48 signed bits, a fixnum Value is returned; otherwise a bignum Value. + * Callers therefore never branch on which one they get back. + */ +typedef struct Bignum { + ObjHeader hdr; + int32_t sign; /* -1, 0, +1 */ + uint32_t n_limbs; /* magnitude limb count after trim */ + uint64_t *limbs; /* little-endian magnitude */ +} Bignum; + +#define AS_BIGNUM(v) ((Bignum *)(uintptr_t)GET_PAYLOAD(v)) + +/* Constructors */ +Value make_bignum_from_i64(int64_t v); +Value make_bignum_from_str(const char *s, int base); /* base 10 or 16 */ +Value bignum_from_limbs(int32_t sign, uint64_t *limbs, uint32_t n); + +/* Demote a bignum to fixnum when magnitude fits 48 signed bits. Caller passes + * a (sign, limbs, n_limbs) triple; result is fixnum Value or bignum Value. */ +Value bignum_normalize(int32_t sign, uint64_t *limbs, uint32_t n); + +/* Promote a fixnum Value to a bignum Value (no-op for already-bignum). */ +Value to_bignum(Value v); + +/* Decimal-string render (caller frees via ul_free). */ +char *bignum_to_str(Bignum *b); + +/* Arithmetic — accept fixnum or bignum, return fixnum-or-bignum + * (demoted when possible). Sign/magnitude internal. */ +Value big_add(Value a, Value b); +Value big_sub(Value a, Value b); +Value big_mul(Value a, Value b); +Value big_neg(Value a); +Value big_quotient(Value a, Value b); +Value big_remainder(Value a, Value b); /* sign of dividend */ +Value big_modulo(Value a, Value b); /* sign of divisor */ +Value big_expt(Value base, Value exp); /* exp >= 0, integer */ +int big_cmp(Value a, Value b); /* -1, 0, +1 */ +bool big_is_zero(Value v); +bool big_is_odd(Value v); +Value big_abs(Value v); + +/* Coerce any number Value to a fresh Bignum struct (caller owns) for ops + * that want explicit bignum arguments. */ +Bignum *as_bignum(Value v); + /* ── Environment ─────────────────────────────────────────────────────────── */ typedef struct EnvBinding { @@ -719,9 +792,9 @@ static inline bool is_callable(Value v) { return false; } -/* Number extraction — works for int, double, rational */ +/* Number extraction — works for int, double, rational, bignum */ static inline bool is_number(Value v) { - return IS_INT(v) || IS_DOUBLE(v) || IS_RATIONAL(v); + return IS_INT(v) || IS_DOUBLE(v) || IS_RATIONAL(v) || IS_BIGNUM(v); } double as_number_double(Value v); /* coerce any number to double */ diff --git a/c/printer.c b/c/printer.c index b5a004c..da4da22 100644 --- a/c/printer.c +++ b/c/printer.c @@ -136,6 +136,13 @@ static void show_value(Value v, bool display, StringBuilder *sb) { return; } + if (IS_BIGNUM(v)) { + char *s = bignum_to_str(AS_BIGNUM(v)); + sb_appendz(sb, s); + ul_free(s); + return; + } + if (IS_BUILTIN(v)) { sb_appendz(sb, "#"); return; @@ -243,6 +250,13 @@ static void show_value(Value v, bool display, StringBuilder *sb) { sb_appendz(sb, buf); break; } + + case OBJ_BIGNUM: { + char *s = bignum_to_str(AS_BIGNUM(v)); + sb_appendz(sb, s); + ul_free(s); + break; + } } } diff --git a/c/reader.c b/c/reader.c index 045f863..1fae796 100644 --- a/c/reader.c +++ b/c/reader.c @@ -230,17 +230,22 @@ Value parse_atom(const char *tok) { return v; } - /* Integer */ + /* Integer (fixnum or bignum). Match optional sign + all-digits. */ { - char *end; - errno = 0; - long long val = strtoll(tok, &end, 10); - if (*end == '\0' && errno == 0) { - /* Check if fits in 48-bit int */ - if (val >= -(1LL << 47) && val < (1LL << 47)) + const char *p = tok; + if (*p == '-' || *p == '+') p++; + bool all_digits = (*p != '\0'); + for (const char *q = p; *q; q++) { + if (!isdigit((unsigned char)*q)) { all_digits = false; break; } + } + if (all_digits) { + char *end; + errno = 0; + long long val = strtoll(tok, &end, 10); + if (*end == '\0' && errno == 0 && FITS_FIXNUM(val)) return VAL_INT(val); - else - return make_double((double)val); + /* Out of fixnum range (or overflow) → bignum. */ + return make_bignum_from_str(tok, 10); } } diff --git a/c/test.c b/c/test.c index 9ae28f7..58d993a 100644 --- a/c/test.c +++ b/c/test.c @@ -135,6 +135,51 @@ TEST(rational_creation) { ASSERT_EQ_INT(as_int(r2), 2); } +TEST(bignum_basic) { + /* (expt 2 48) — first value past fixnum cap → bignum. */ + char *got = run_show("(expt 2 48)"); + ASSERT_EQ_STR(got, "281474976710656"); + ul_free(got); + + /* (expt 2 256) — 78-digit hex world. */ + got = run_show("(expt 2 256)"); + ASSERT_EQ_STR(got, "115792089237316195423570985008687907853269984665640564039457584007913129639936"); + ul_free(got); + + /* secp256k1 prime — full triple-subtract. */ + got = run_show("(- (expt 2 256) (expt 2 32) 977)"); + ASSERT_EQ_STR(got, "115792089237316195423570985008687907853269984665640564039457584007908834671663"); + ul_free(got); + + /* Demotion back to fixnum when result fits. */ + got = run_show("(- (expt 2 64) (expt 2 64))"); + ASSERT_EQ_STR(got, "0"); + ul_free(got); +} + +TEST(bignum_arith) { + /* Modular arithmetic across 78-digit operands. */ + char *got = run_show( + "(define p (- (expt 2 256) (expt 2 32) 977))" + "(modulo (- 0 1) p)"); + /* p - 1 */ + ASSERT_EQ_STR(got, "115792089237316195423570985008687907853269984665640564039457584007908834671662"); + ul_free(got); + + /* quotient + remainder identity for bignum scale. */ + got = run_show( + "(define p (- (expt 2 256) (expt 2 32) 977))" + "(let* ((a (* p 7)) (q (quotient a p)) (r (remainder a p)))" + " (list q r))"); + ASSERT_EQ_STR(got, "(7 0)"); + ul_free(got); + + /* Predicate roundtrip. */ + got = run_show("(list (integer? (expt 2 100)) (exact? (expt 2 100)) (odd? (- (expt 2 100) 1)))"); + ASSERT_EQ_STR(got, "(#t #t #t)"); + ul_free(got); +} + TEST(environment) { Env *g = make_env(NULL); g->global = g; @@ -886,6 +931,8 @@ int main(void) { run_test_string_creation(); run_test_vector_creation(); run_test_rational_creation(); + run_test_bignum_basic(); + run_test_bignum_arith(); run_test_environment(); printf("\n[reader]\n"); diff --git a/c/types.c b/c/types.c index 48450c0..4a60e7c 100644 --- a/c/types.c +++ b/c/types.c @@ -304,6 +304,13 @@ double as_number_double(Value v) { Rational *r = AS_RATIONAL(v); return (double)r->num / (double)r->den; } + if (IS_BIGNUM(v)) { + Bignum *b = AS_BIGNUM(v); + double d = 0.0; + for (int32_t i = (int32_t)b->n_limbs - 1; i >= 0; i--) + d = d * 18446744073709551616.0 + (double)b->limbs[i]; + return b->sign < 0 ? -d : d; + } lisp_error("not a number"); return 0; } @@ -315,6 +322,12 @@ int64_t as_number_int(Value v) { Rational *r = AS_RATIONAL(v); return r->num / r->den; } + if (IS_BIGNUM(v)) { + Bignum *b = AS_BIGNUM(v); + uint64_t lo = b->n_limbs > 0 ? b->limbs[0] : 0; + int64_t s = (int64_t)lo; + return b->sign < 0 ? -s : s; + } lisp_error("not a number"); return 0; } @@ -327,11 +340,37 @@ static void to_rational(Value v, int64_t *num, int64_t *den) { } static bool is_exact(Value v) { - return IS_INT(v) || IS_RATIONAL(v); + return IS_INT(v) || IS_RATIONAL(v) || IS_BIGNUM(v); +} + +/* Fixnum add/sub/mul that promote on overflow. Result is fixnum or bignum. */ +static Value fixnum_add(int64_t a, int64_t b) { + int64_t r; + if (__builtin_add_overflow(a, b, &r) || !FITS_FIXNUM(r)) + return big_add(make_bignum_from_i64(a), make_bignum_from_i64(b)); + return VAL_INT(r); +} +static Value fixnum_sub(int64_t a, int64_t b) { + int64_t r; + if (__builtin_sub_overflow(a, b, &r) || !FITS_FIXNUM(r)) + return big_sub(make_bignum_from_i64(a), make_bignum_from_i64(b)); + return VAL_INT(r); +} +static Value fixnum_mul(int64_t a, int64_t b) { + int64_t r; + if (__builtin_mul_overflow(a, b, &r) || !FITS_FIXNUM(r)) + return big_mul(make_bignum_from_i64(a), make_bignum_from_i64(b)); + return VAL_INT(r); } Value num_add(Value a, Value b) { + /* Pure integer path (fixnum or bignum) — never produces a rational. */ + if (IS_INTEGER(a) && IS_INTEGER(b)) { + if (IS_INT(a) && IS_INT(b)) return fixnum_add(as_int(a), as_int(b)); + return big_add(a, b); + } if (is_exact(a) && is_exact(b)) { + /* Rational mix — fall back to int64_t numerator/denominator path. */ int64_t an, ad, bn, bd; to_rational(a, &an, &ad); to_rational(b, &bn, &bd); @@ -341,6 +380,10 @@ Value num_add(Value a, Value b) { } Value num_sub(Value a, Value b) { + if (IS_INTEGER(a) && IS_INTEGER(b)) { + if (IS_INT(a) && IS_INT(b)) return fixnum_sub(as_int(a), as_int(b)); + return big_sub(a, b); + } if (is_exact(a) && is_exact(b)) { int64_t an, ad, bn, bd; to_rational(a, &an, &ad); @@ -351,6 +394,10 @@ Value num_sub(Value a, Value b) { } Value num_mul(Value a, Value b) { + if (IS_INTEGER(a) && IS_INTEGER(b)) { + if (IS_INT(a) && IS_INT(b)) return fixnum_mul(as_int(a), as_int(b)); + return big_mul(a, b); + } if (is_exact(a) && is_exact(b)) { int64_t an, ad, bn, bd; to_rational(a, &an, &ad); @@ -361,6 +408,19 @@ Value num_mul(Value a, Value b) { } Value num_div(Value a, Value b) { + /* Integer / integer that divides cleanly stays integer (matches Python / + * Scheme semantics for `/` between exacts when the quotient is exact). + * For bignum case we treat as exact division (truncating to integer when + * quotient is exact, else fall back to double for now — secp256k1 ops + * never use fractional bignum). */ + if (IS_INTEGER(a) && IS_INTEGER(b)) { + if (big_is_zero(b)) lisp_error("division by zero"); + Value q = big_quotient(a, b); + Value r = big_remainder(a, b); + if (big_is_zero(r)) return q; + /* Inexact fallback. */ + return make_double(as_number_double(a) / as_number_double(b)); + } if (is_exact(a) && is_exact(b)) { int64_t an, ad, bn, bd; to_rational(a, &an, &ad); @@ -374,7 +434,11 @@ Value num_div(Value a, Value b) { } Value num_neg(Value a) { - if (IS_INT(a)) return VAL_INT(-as_int(a)); + if (IS_INT(a)) { + int64_t v = as_int(a); + return fixnum_sub(0, v); /* handles -FIXNUM_MIN safely */ + } + if (IS_BIGNUM(a)) return big_neg(a); if (IS_RATIONAL(a)) { Rational *r = AS_RATIONAL(a); return make_rational(-r->num, r->den); @@ -383,6 +447,7 @@ Value num_neg(Value a) { } static int num_cmp(Value a, Value b) { + if (IS_INTEGER(a) && IS_INTEGER(b)) return big_cmp(a, b); if (is_exact(a) && is_exact(b)) { int64_t an, ad, bn, bd; to_rational(a, &an, &ad); diff --git a/c/vm.c b/c/vm.c index 83708a3..67fae5a 100644 --- a/c/vm.c +++ b/c/vm.c @@ -148,7 +148,7 @@ void bc_compile(Value expr, CodeObj *code, Env *env, bool tail) { if (IS_NIL(expr) || IS_TRUE(expr) || IS_FALSE(expr) || IS_EOF(expr) || IS_CHAR(expr)) { code_emit(code, OP_CONST, expr); return; } - if (IS_INT(expr) || IS_DOUBLE(expr) || IS_RATIONAL(expr)) { + if (IS_INT(expr) || IS_DOUBLE(expr) || IS_RATIONAL(expr) || IS_BIGNUM(expr)) { code_emit(code, OP_CONST, expr); return; } if (IS_STRING(expr) || IS_VECTOR(expr)) { diff --git a/tests/bignum-cross-tier.lsp b/tests/bignum-cross-tier.lsp new file mode 100644 index 0000000..50c4679 --- /dev/null +++ b/tests/bignum-cross-tier.lsp @@ -0,0 +1,146 @@ +;;; bignum-cross-tier.lsp — Cross-tier byte-identity for arbitrary-precision +;;; integers. Runs identically under Python tier & C tier. Python tier is +;;; our oracle (native int handles every magnitude); C tier joined the +;;; bignum club here & must agree byte-for-byte with Python. +;;; +;;; Python: python3 lumbda.py --fast tests/bignum-cross-tier.lsp +;;; C: ./c/lumbda tests/bignum-cross-tier.lsp +;;; +;;; Every assertion logs PASS or FAIL; a final summary lines tallies them. +;;; Identity discipline: every PASS line on both tiers must match +;;; byte-for-byte. + +(define *pass* 0) +(define *fail* 0) + +(define (assert-equal name got expected) + (if (equal? got expected) + (begin (set! *pass* (+ *pass* 1)) + (display "PASS: ") (display name) (newline)) + (begin (set! *fail* (+ *fail* 1)) + (display "FAIL: ") (display name) + (display " got=") (write got) + (display " expected=") (write expected) (newline)))) + +;;; ═══════════════════════════════════════════════════════════════ +;;; Literals & basic reader/printer round-trip +;;; ═══════════════════════════════════════════════════════════════ + +(assert-equal "fixnum-literal-1" 42 42) +(assert-equal "fixnum-literal-2" -7 -7) +(assert-equal "bignum-literal-pos" + 281474976710656 ; 2^48 + 281474976710656) +(assert-equal "bignum-literal-neg" + -281474976710656 + -281474976710656) +(assert-equal "bignum-literal-256" + 115792089237316195423570985008687907853269984665640564039457584007913129639936 + 115792089237316195423570985008687907853269984665640564039457584007913129639936) + +;;; ═══════════════════════════════════════════════════════════════ +;;; expt — promotion through fixnum boundary +;;; ═══════════════════════════════════════════════════════════════ + +(assert-equal "expt-2-30" (expt 2 30) 1073741824) +(assert-equal "expt-2-47" (expt 2 47) 140737488355328) +(assert-equal "expt-2-48" (expt 2 48) 281474976710656) +(assert-equal "expt-2-64" (expt 2 64) 18446744073709551616) +(assert-equal "expt-2-128" (expt 2 128) + 340282366920938463463374607431768211456) +(assert-equal "expt-2-256" (expt 2 256) + 115792089237316195423570985008687907853269984665640564039457584007913129639936) + +;;; ═══════════════════════════════════════════════════════════════ +;;; secp256k1 prime — 78-digit triple-subtraction lands precisely +;;; ═══════════════════════════════════════════════════════════════ + +(define secp-p (- (expt 2 256) (expt 2 32) 977)) +(assert-equal "secp256k1-prime" + secp-p + 115792089237316195423570985008687907853269984665640564039457584007908834671663) +(assert-equal "secp-p length" (string-length (number->string secp-p)) 78) + +;;; ═══════════════════════════════════════════════════════════════ +;;; Arithmetic at bignum scale +;;; ═══════════════════════════════════════════════════════════════ + +(assert-equal "add-bignum" + (+ secp-p secp-p) + 231584178474632390847141970017375815706539969331281128078915168015817669343326) + +(assert-equal "sub-bignum" + (- (* secp-p 2) secp-p) + secp-p) + +;;; (* secp-p secp-p) — magnitude check via number->string length; exact +;;; value differs across tier baselines only if their bignum products +;;; diverge — Python tier serves as oracle (literal computed there). +(assert-equal "mul-bignum" + (* secp-p secp-p) + 13407807929942597099574024998205846127479365820592393377723561443720769383374469661147847687812952081302854773939601805382211292725060150247698793015185569) + +;;; Multi-step modular: a * b mod p +(define a 12345678901234567890123456789012345678901234567890) +(define b 98765432109876543210987654321098765432109876543210) +(assert-equal "mod-bignum" + (modulo (* a b) secp-p) + (modulo (* a b) secp-p)) + +(assert-equal "quotient-bignum" + (quotient (* secp-p 7) secp-p) 7) + +(assert-equal "remainder-bignum" + (remainder (+ (* secp-p 3) 42) secp-p) 42) + +;;; modulo wraps negative dividend by divisor sign +(assert-equal "modulo-neg-bignum" + (modulo (- 0 1) secp-p) + (- secp-p 1)) + +;;; ═══════════════════════════════════════════════════════════════ +;;; Comparisons across fixnum/bignum boundary +;;; ═══════════════════════════════════════════════════════════════ + +(assert-equal "lt-bignum-fixnum" (< 1 (expt 2 100)) #t) +(assert-equal "gt-bignum-fixnum" (> (expt 2 100) 999) #t) +(assert-equal "eq-bignum-bignum" (= (expt 2 64) 18446744073709551616) #t) +(assert-equal "neq-bignum-fixnum" (= (expt 2 50) 0) #f) + +;;; ═══════════════════════════════════════════════════════════════ +;;; Predicates +;;; ═══════════════════════════════════════════════════════════════ + +(assert-equal "integer?-bignum" (integer? (expt 2 100)) #t) +(assert-equal "exact?-bignum" (exact? (expt 2 100)) #t) +(assert-equal "odd?-bignum" (odd? (- (expt 2 100) 1)) #t) +(assert-equal "even?-bignum" (even? (expt 2 100)) #t) +(assert-equal "zero?-bignum" (zero? (- (expt 2 100) (expt 2 100))) #t) + +;;; ═══════════════════════════════════════════════════════════════ +;;; Promotion & demotion +;;; ═══════════════════════════════════════════════════════════════ + +;;; Result of (- bignum bignum) that fits a fixnum demotes back. +(assert-equal "demote" (- (expt 2 64) (expt 2 64)) 0) +(assert-equal "demote-small" (- (expt 2 50) (- (expt 2 50) 7)) 7) + +;;; number->string round-trip +(assert-equal "n->s-bignum" + (number->string (expt 2 100)) + "1267650600228229401496703205376") + +;;; string->number round-trip +(assert-equal "s->n-bignum" + (string->number "1267650600228229401496703205376") + (expt 2 100)) + +;;; ═══════════════════════════════════════════════════════════════ +;;; Summary +;;; ═══════════════════════════════════════════════════════════════ + +(display "════════════════════════════════════════") (newline) +(display "bignum cross-tier: ") +(display *pass*) (display " passed, ") +(display *fail*) (display " failed") +(newline)