feat(stdio): printf family

src/stdio/vfprintf_float.c

@@ -0,0 +1,842 @@
+/* internal unsigned 2^32-limb bignum core for printf float formatting (todo 16)
+ * clean-room implementation; nothing here is exported */
+/*
+ * Invariants: limbs are little-endian (value = sum d[i] * 2^(32*i)); n == 0
+ * means value 0; no trailing-zero-limb invariant is kept -- callers normalize
+ * via vfpn_bn_norm. Shifts and rounding operate in place; allocation failure
+ * is reported as -1 and leaves errno to __libc_malloc.
+ */
+#include "../internal/malloc.h"
+#include <stddef.h>
+#include <stdint.h>
+
+struct vfpn_bn
+{
+    uint32_t fixed[640];
+    uint32_t *d;
+    size_t n;
+    size_t cap;
+    int heap;
+};
+
+static __attribute__((unused)) void
+vfpn_bn_init(struct vfpn_bn *b)
+{
+    b->d = b->fixed;
+    b->n = 0;
+    b->cap = 640;
+    b->heap = 0;
+}
+
+static __attribute__((unused)) void
+vfpn_bn_free(struct vfpn_bn *b)
+{
+    if (b->heap)
+    {
+        __libc_free(b->d);
+    }
+}
+
+static __attribute__((unused)) int
+vfpn_bn_reserve(struct vfpn_bn *b, size_t need)
+{
+    if (need <= b->cap)
+    {
+        return 0;
+    }
+
+    uint32_t *new = __libc_malloc(need * sizeof(uint32_t));
+    if (!new)
+    {
+        return -1;
+    }
+
+    size_t copy = (b->cap < need) ? b->cap : need;
+    __builtin_memcpy(new, b->d, copy * sizeof(uint32_t));
+    if (b->heap)
+    {
+        __libc_free(b->d);
+    }
+
+    b->d = new;
+    b->cap = need;
+    b->heap = 1;
+    return 0;
+}
+
+static __attribute__((unused)) void
+vfpn_bn_norm(struct vfpn_bn *b)
+{
+    while (b->n > 0 && b->d[b->n - 1] == 0)
+    {
+        b->n--;
+    }
+}
+
+static __attribute__((unused)) int
+vfpn_bn_set_u64(struct vfpn_bn *b, uint64_t v)
+{
+    size_t need = 0;
+    for (uint64_t t = v; t != 0; t >>= 32)
+    {
+        need++;
+    }
+
+    if (vfpn_bn_reserve(b, need) != 0)
+    {
+        return -1;
+    }
+
+    b->n = 0;
+    while (v != 0)
+    {
+        b->d[b->n++] = (uint32_t)v;
+        v >>= 32;
+    }
+
+    return 0;
+}
+
+static __attribute__((unused)) int
+vfpn_bn_mul_small(struct vfpn_bn *b, uint32_t m)
+{
+    if (m == 0)
+    {
+        b->n = 0;
+        return 0;
+    }
+
+    if (vfpn_bn_reserve(b, b->n + 1) != 0)
+    {
+        return -1;
+    }
+
+    uint64_t carry = 0;
+    for (size_t i = 0; i < b->n; i++)
+    {
+        uint64_t cur = (uint64_t)b->d[i] * m + carry;
+        b->d[i] = (uint32_t)cur;
+        carry = cur >> 32;
+    }
+
+    if (carry != 0)
+    {
+        b->d[b->n++] = (uint32_t)carry;
+    }
+
+    return 0;
+}
+
+static __attribute__((unused)) int
+vfpn_bn_shl(struct vfpn_bn *b, size_t bits)
+{
+    if (bits == 0)
+    {
+        return 0;
+    }
+
+    size_t wl = bits / 32;
+    size_t s = bits % 32;
+    size_t newn = b->n + wl + ((s != 0) ? 1 : 0);
+    if (vfpn_bn_reserve(b, newn) != 0)
+    {
+        return -1;
+    }
+
+    if (s == 0)
+    {
+        for (size_t i = b->n; i-- > 0;)
+        {
+            b->d[i + wl] = b->d[i];
+        }
+
+        for (size_t i = 0; i < wl; i++)
+        {
+            b->d[i] = 0;
+        }
+    }
+    else
+    {
+        for (size_t j = newn; j-- > 0;)
+        {
+            uint64_t cur = 0;
+            if (j >= wl && j - wl < b->n)
+            {
+                cur |= (uint64_t)b->d[j - wl] << s;
+            }
+
+            if (j > wl && j - wl - 1 < b->n)
+            {
+                cur |= (uint64_t)b->d[j - wl - 1] >> (32 - s);
+            }
+
+            b->d[j] = (uint32_t)cur;
+        }
+    }
+
+    b->n = newn;
+    return 0;
+}
+
+static __attribute__((unused)) int
+vfpn_bn_add_one(struct vfpn_bn *b)
+{
+    size_t i = 0;
+    while (i < b->n)
+    {
+        b->d[i]++;
+        if (b->d[i] != 0)
+        {
+            return 0;
+        }
+
+        i++;
+    }
+
+    if (vfpn_bn_reserve(b, b->n + 1) != 0)
+    {
+        return -1;
+    }
+
+    b->d[b->n] = 1;
+    b->n++;
+    return 0;
+}
+
+static __attribute__((unused)) int
+vfpn_bn_is_zero(const struct vfpn_bn *b)
+{
+    return b->n == 0;
+}
+
+static __attribute__((unused)) int
+vfpn_bn_cmp(const struct vfpn_bn *a, const struct vfpn_bn *b)
+{
+    if (a->n < b->n)
+    {
+        return -1;
+    }
+
+    if (a->n > b->n)
+    {
+        return 1;
+    }
+
+    for (size_t i = a->n; i-- > 0;)
+    {
+        if (a->d[i] < b->d[i])
+        {
+            return -1;
+        }
+
+        if (a->d[i] > b->d[i])
+        {
+            return 1;
+        }
+    }
+
+    return 0;
+}
+
+static __attribute__((unused)) int
+vfpn_bn_shr_round(struct vfpn_bn *b, size_t bits)
+{
+    if (bits == 0)
+    {
+        return 0;
+    }
+
+    size_t w = bits / 32;
+    size_t s = bits % 32;
+
+    if (w >= b->n)
+    {
+        /* Whole value discarded: result is 0, rounded up iff the original
+         * value exceeds 2^(bits-1). That can only happen when bits == 32*n
+         * (s == 0 and w == n), so compare the top limb against half. */
+        int round = 0;
+        if (s == 0 && w == b->n && b->n > 0)
+        {
+            uint32_t top = b->d[b->n - 1];
+            if (top > 0x80000000u)
+            {
+                round = 1;
+            }
+            else if (top == 0x80000000u)
+            {
+                for (size_t i = b->n - 1; i-- > 0;)
+                {
+                    if (b->d[i] != 0)
+                    {
+                        round = 1;
+                        break;
+                    }
+                }
+            }
+        }
+
+        b->n = 0;
+        if (round)
+        {
+            return vfpn_bn_add_one(b);
+        }
+
+        return 0;
+    }
+
+    /* Round bit and sticky are captured from the discarded low `bits` bits
+     * of the original value before the in-place shift. */
+    int round_bit;
+    int sticky;
+    if (s > 0)
+    {
+        round_bit = (int)((b->d[w] >> (s - 1)) & 1u);
+        sticky = (b->d[w] & ((1u << (s - 1)) - 1u)) != 0;
+        for (size_t i = 0; i < w; i++)
+        {
+            if (b->d[i] != 0)
+            {
+                sticky = 1;
+                break;
+            }
+        }
+    }
+    else
+    {
+        round_bit = (int)((b->d[w - 1] >> 31) & 1u);
+        sticky = (b->d[w - 1] & 0x7FFFFFFFu) != 0;
+        for (size_t i = 0; i + 1 < w; i++)
+        {
+            if (b->d[i] != 0)
+            {
+                sticky = 1;
+                break;
+            }
+        }
+    }
+
+    /* In-place right shift by `bits`; ascending so every source limb is
+     * still readable when its destination is written. */
+    size_t old_n = b->n;
+    if (s == 0)
+    {
+        for (size_t j = 0; j + w < old_n; j++)
+        {
+            b->d[j] = b->d[j + w];
+        }
+    }
+    else
+    {
+        for (size_t j = 0; j + w < old_n; j++)
+        {
+            uint32_t lo = b->d[j + w] >> s;
+            uint32_t hi = (j + w + 1 < old_n) ? (b->d[j + w + 1] << (32 - s)) : 0;
+            b->d[j] = lo | hi;
+        }
+    }
+
+    b->n = old_n - w;
+    vfpn_bn_norm(b);
+
+    int odd = (b->n > 0 && (b->d[0] & 1u)) ? 1 : 0;
+    if (round_bit && (sticky || odd))
+    {
+        return vfpn_bn_add_one(b);
+    }
+
+    return 0;
+}
+
+static __attribute__((unused)) uint32_t
+vfpn_bn_divmod_small_1e9(struct vfpn_bn *b)
+{
+    uint64_t r = 0;
+    for (size_t i = b->n; i-- > 0;)
+    {
+        uint64_t cur = (r << 32) | b->d[i];
+        b->d[i] = (uint32_t)(cur / 1000000000ULL);
+        r = cur % 1000000000ULL;
+    }
+
+    vfpn_bn_norm(b);
+    return (uint32_t)r;
+}
+
+static __attribute__((unused)) int
+vfpn_bn_pow5(struct vfpn_bn *b, size_t k)
+{
+    int r = vfpn_bn_set_u64(b, 1);
+    if (r != 0)
+    {
+        return r;
+    }
+
+    for (size_t i = 0; i < k; i++)
+    {
+        r = vfpn_bn_mul_small(b, 5);
+        if (r != 0)
+        {
+            return r;
+        }
+    }
+
+    return 0;
+}
+
+/* --- exact decimal rounding layer ------------------------------------ */
+
+/*
+ * Number of bits needed to hold |b| (index of the highest set bit plus one);
+ * 0 for the value zero. The caller may leave a trailing zero limb in place,
+ * so the scan starts from the recorded length and skips down.
+ */
+static __attribute__((unused)) size_t
+vfpn_bn_bitlen(const struct vfpn_bn *b)
+{
+    if (b->n == 0)
+    {
+        return 0;
+    }
+
+    size_t i = b->n - 1;
+    while (i > 0 && b->d[i] == 0)
+    {
+        i--;
+    }
+
+    if (b->d[i] == 0)
+    {
+        return 0;
+    }
+
+    return i * 32 + (size_t)(32 - __builtin_clz(b->d[i]));
+}
+
+static __attribute__((unused)) int
+vfpn_bn_get_bit(const struct vfpn_bn *b, size_t i)
+{
+    return (int)((b->d[i / 32] >> (i % 32)) & 1u);
+}
+
+/*
+ * Set bit i of b. Zero-extends the limb array when the bit lies beyond the
+ * current length; ORing never clears, so the value only grows.
+ */
+static __attribute__((unused)) int
+vfpn_bn_set_bit(struct vfpn_bn *b, size_t i, int bit)
+{
+    if (!bit)
+    {
+        return 0;
+    }
+
+    if (b->n <= i / 32)
+    {
+        size_t need = i / 32 + 1;
+        if (vfpn_bn_reserve(b, need) != 0)
+        {
+            return -1;
+        }
+
+        for (size_t j = b->n; j < need; j++)
+        {
+            b->d[j] = 0;
+        }
+
+        b->n = need;
+    }
+
+    b->d[i / 32] |= (1u << (i % 32));
+    return 0;
+}
+
+/*
+ * Shift |b| left by one bit in place (limb ripple). The caller's cap must
+ * already make room for the possible extra top limb; the reserve here is
+ * only a safety net and normally returns immediately.
+ */
+static __attribute__((unused)) int
+vfpn_bn_shl1(struct vfpn_bn *b)
+{
+    if (b->n == 0)
+    {
+        return 0;
+    }
+
+    if (vfpn_bn_reserve(b, b->n + 1) != 0)
+    {
+        return -1;
+    }
+
+    uint32_t carry = 0;
+    for (size_t i = 0; i < b->n; i++)
+    {
+        uint32_t cur = b->d[i];
+        b->d[i] = (cur << 1) | carry;
+        carry = cur >> 31;
+    }
+
+    if (carry != 0)
+    {
+        b->d[b->n++] = carry;
+    }
+
+    return 0;
+}
+
+/*
+ * In-place limb-wise subtraction a -= b. The caller guarantees a >= b, so
+ * the borrow always clears before the top of a; a is renormalized.
+ */
+static __attribute__((unused)) void
+vfpn_bn_sub_inplace(struct vfpn_bn *a, const struct vfpn_bn *b)
+{
+    uint32_t borrow = 0;
+    size_t i = 0;
+    for (; i < b->n; i++)
+    {
+        uint64_t diff = (uint64_t)a->d[i] - (uint64_t)b->d[i] - borrow;
+        a->d[i] = (uint32_t)diff;
+        borrow = (uint32_t)((diff >> 32) & 1u);
+    }
+
+    while (i < a->n && borrow != 0)
+    {
+        a->d[i]--;
+        borrow = (a->d[i] == 0xFFFFFFFFu) ? 1u : 0u;
+        i++;
+    }
+
+    vfpn_bn_norm(a);
+}
+
+/*
+ * Q = N / D, R = N mod D by restoring binary long division. N may be zero;
+ * D must be nonzero. Q and R are written over whatever they held. The
+ * quotients bits are produced most-significant first; the loop runs over the
+ * dividend's bit positions from bitlen(N)-1 down to 0, pulling one new
+ * dividend bit into R and subtracting D whenever R >= D.
+ */
+static __attribute__((unused)) int
+vfpn_divmod_binary(struct vfpn_bn *Q, struct vfpn_bn *R, const struct vfpn_bn *N,
+                   const struct vfpn_bn *D)
+{
+    if (vfpn_bn_is_zero(D))
+    {
+        return -1;
+    }
+
+    if (vfpn_bn_reserve(Q, N->n + 2) != 0)
+    {
+        return -1;
+    }
+
+    if (vfpn_bn_reserve(R, N->n + 2) != 0)
+    {
+        return -1;
+    }
+
+    Q->n = 0;
+    R->n = 0;
+
+    size_t nbits = vfpn_bn_bitlen(N);
+    for (size_t i = nbits; i-- > 0;)
+    {
+        if (vfpn_bn_shl1(R) != 0)
+        {
+            return -1;
+        }
+
+        if (vfpn_bn_set_bit(R, 0, vfpn_bn_get_bit(N, i)) != 0)
+        {
+            return -1;
+        }
+
+        if (vfpn_bn_cmp(R, D) >= 0)
+        {
+            vfpn_bn_sub_inplace(R, D);
+            if (vfpn_bn_set_bit(Q, i, 1) != 0)
+            {
+                return -1;
+            }
+        }
+    }
+
+    vfpn_bn_norm(Q);
+    vfpn_bn_norm(R);
+    return 0;
+}
+
+/*
+ * Compare 2*a against b (both canonical base-2^32 digit strings) without
+ * building the doubled value: doubling a limb's carry is the next limb's
+ * (d[i] >> 31), which the loop below folds in, and an extra top digit
+ * appears when the original top limb had its high bit set. All digits of 2*a
+ * stay below 2^32, so a plain digitwise comparison is exact.
+ */
+static __attribute__((unused)) int
+vfpn_bn_cmp_doubled(const struct vfpn_bn *a, const struct vfpn_bn *b)
+{
+    if (a->n == 0)
+    {
+        return vfpn_bn_is_zero(b) ? 0 : -1;
+    }
+
+    size_t extra = (a->d[a->n - 1] >> 31) ? 1u : 0u;
+    size_t la = a->n + extra;
+    if (la != b->n)
+    {
+        return (la > b->n) ? 1 : -1;
+    }
+
+    for (size_t i = la; i-- > 0;)
+    {
+        uint32_t da;
+        if (extra != 0 && i == a->n)
+        {
+            da = 1; /* carry out of the doubled top limb */
+        }
+        else
+        {
+            da = (uint32_t)(a->d[i] << 1);
+            if (i > 0)
+            {
+                da += (a->d[i - 1] >> 31);
+            }
+        }
+
+        if (da != b->d[i])
+        {
+            return (da < b->d[i]) ? -1 : 1;
+        }
+    }
+
+    return 0;
+}
+
+/*
+ * R = round_half_even(mant * 2^exp2 * 10^s), computed exactly as an integer
+ * bignum (mant > 0; s may be negative). s >= 0 scales R by 5^s and by a
+ * single power of two, rounding once at the end (nearest-even). s < 0 builds
+ * D = 5^(-s) and shifts either N or D so that value = N / D, divides exactly,
+ * and rounds the quotient half-even off the remainder: 2*rem > D rounds up,
+ * 2*rem == D rounds up only when the quotient is odd. Returns 0 or -1 on
+ * allocation failure (R is then unspecified).
+ */
+static __attribute__((unused)) int
+vfpn_round_scale10(struct vfpn_bn *R, uint64_t mant, int exp2, long s)
+{
+    struct vfpn_bn N;
+    struct vfpn_bn D;
+    struct vfpn_bn Q;
+    struct vfpn_bn rem;
+    vfpn_bn_init(&N);
+    vfpn_bn_init(&D);
+    vfpn_bn_init(&Q);
+    vfpn_bn_init(&rem);
+
+    if (s >= 0)
+    {
+        if (vfpn_bn_set_u64(R, mant) != 0)
+        {
+            goto fail;
+        }
+
+        for (long i = 0; i < s; i++)
+        {
+            if (vfpn_bn_mul_small(R, 5) != 0)
+            {
+                goto fail;
+            }
+        }
+
+        long long sh = (long long)exp2 + (long long)s;
+        if (sh >= 0)
+        {
+            if (vfpn_bn_shl(R, (size_t)sh) != 0)
+            {
+                goto fail;
+            }
+        }
+        else
+        {
+            if (vfpn_bn_shr_round(R, (size_t)(-sh)) != 0)
+            {
+                goto fail;
+            }
+        }
+
+        vfpn_bn_free(&N);
+        vfpn_bn_free(&D);
+        vfpn_bn_free(&Q);
+        vfpn_bn_free(&rem);
+        return 0;
+    }
+
+    /* value = mant * 2^exp2 / (2^t * 5^t) with t = -s: make it N / D. */
+    {
+        long long t = -(long long)s;
+        long long sh = (long long)exp2 - t;
+
+        if (vfpn_bn_pow5(&D, (size_t)t) != 0)
+        {
+            goto fail;
+        }
+
+        if (vfpn_bn_set_u64(&N, mant) != 0)
+        {
+            goto fail;
+        }
+
+        if (sh >= 0)
+        {
+            if (vfpn_bn_shl(&N, (size_t)sh) != 0)
+            {
+                goto fail;
+            }
+
+            vfpn_bn_norm(&N);
+        }
+        else
+        {
+            if (vfpn_bn_shl(&D, (size_t)(-sh)) != 0)
+            {
+                goto fail;
+            }
+
+            vfpn_bn_norm(&D);
+        }
+
+        if (vfpn_divmod_binary(&Q, &rem, &N, &D) != 0)
+        {
+            goto fail;
+        }
+
+        int c = vfpn_bn_cmp_doubled(&rem, &D);
+        if (c > 0 || (c == 0 && Q.n > 0 && (Q.d[0] & 1u) != 0))
+        {
+            if (vfpn_bn_add_one(&Q) != 0)
+            {
+                goto fail;
+            }
+        }
+
+        if (vfpn_bn_reserve(R, Q.n) != 0)
+        {
+            goto fail;
+        }
+
+        if (Q.n > 0)
+        {
+            __builtin_memcpy(R->d, Q.d, Q.n * sizeof(uint32_t));
+        }
+
+        R->n = Q.n;
+    }
+
+    vfpn_bn_free(&N);
+    vfpn_bn_free(&D);
+    vfpn_bn_free(&Q);
+    vfpn_bn_free(&rem);
+    return 0;
+
+fail:
+    vfpn_bn_free(&N);
+    vfpn_bn_free(&D);
+    vfpn_bn_free(&Q);
+    vfpn_bn_free(&rem);
+    return -1;
+}
+
+/*
+ * Decimal digits of |v| into out (most significant first, no terminating
+ * NUL); returns the digit count or -1 when outcap is too small. The value
+ * zero yields a single '0'. Groups of nine decimal digits are stripped from
+ * the low end into a fixed stack (a 640-limb value caps at ~6165 digits, far
+ * inside 700 groups), then re-emitted top group unpadded, lower groups
+ * zero-padded to nine.
+ */
+static __attribute__((unused)) int
+vfpn_bn_to_dec(const struct vfpn_bn *v, char *out, size_t outcap)
+{
+    if (vfpn_bn_is_zero(v))
+    {
+        if (outcap < 1)
+        {
+            return -1;
+        }
+
+        out[0] = '0';
+        return 1;
+    }
+
+    struct vfpn_bn w;
+    vfpn_bn_init(&w);
+    if (vfpn_bn_reserve(&w, v->n) != 0)
+    {
+        return -1;
+    }
+
+    w.n = v->n;
+    if (v->n > 0)
+    {
+        __builtin_memcpy(w.d, v->d, v->n * sizeof(uint32_t));
+    }
+
+    uint32_t groups[700];
+    size_t ngroups = 0;
+    while (!vfpn_bn_is_zero(&w))
+    {
+        if (ngroups >= 700)
+        {
+            vfpn_bn_free(&w);
+            return -1;
+        }
+
+        groups[ngroups++] = vfpn_bn_divmod_small_1e9(&w);
+    }
+
+    vfpn_bn_free(&w);
+
+    size_t top = ngroups - 1;
+    uint32_t t = groups[top];
+    size_t tdig = 1;
+    for (uint32_t x = t / 10; x != 0; x /= 10)
+    {
+        tdig++;
+    }
+
+    size_t ndigits = tdig + (ngroups - 1) * 9;
+    if (outcap < ndigits)
+    {
+        return -1;
+    }
+
+    size_t pos = 0;
+    char tmp[10];
+    size_t ntmp = 0;
+    do
+    {
+        tmp[ntmp++] = (char)('0' + t % 10);
+        t /= 10;
+    } while (t != 0);
+
+    while (ntmp > 0)
+    {
+        out[pos++] = tmp[--ntmp];
+    }
+
+    for (size_t gi = ngroups - 1; gi-- > 0;)
+    {
+        uint32_t val = groups[gi];
+        for (int k = 8; k >= 0; k--)
+        {
+            out[pos + (size_t)k] = (char)('0' + val % 10);
+            val /= 10;
+        }
+
+        pos += 9;
+    }
+
+    return (int)ndigits;
+}