lumbda/c/bignum.c

/*
 * 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;
}