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