#include <stdckdint.h>
#include "bresenham-distribution.h"

/* (a + b) mod m. Require a < m, b < m, and m > 0. */
static uintmax_t bd_addmod(uintmax_t a, uintmax_t b, uintmax_t m)
{
    uintmax_t sum;
    if (ckd_add(&sum, a, b) || sum >= m)
        sum -= m;  /* Unsigned wraparound is intentional when addition overflowed. */
    return sum;
}

/* (a * b) mod m, with the same preconditions as bd_addmod. */
static uintmax_t bd_mulmod(uintmax_t a, uintmax_t b, uintmax_t m)
{
    uintmax_t product;
    if (!ckd_mul(&product, a, b))
        return product % m;

    /* The wrapped product alone cannot recover the remainder: use doubling. */
    product = 0;
    while (b) {
        if (b & 1)
            product = bd_addmod(product, a, m);
        b >>= 1;
        if (b)
            a = bd_addmod(a, a, m);
    }
    return product;
}

uintmax_t bd_count(uintmax_t N, uintmax_t B, uintmax_t b, uintmax_t shift)
{
    if (N == UINTMAX_MAX || B == 0 || b >= B || shift >= B)
        return UINTMAX_MAX;

    uintmax_t base = N / B;
    uintmax_t R = N % B;
    if (R == 0)
        return base;

    /* (b - shift) mod B, without an overflowing intermediate sum. */
    uintmax_t k = b >= shift ? b - shift : B - (shift - b);

    /* The floor difference is one exactly when (k * R) mod B >= B - R. */
    return base + (bd_mulmod(k, R, B) >= B - R);
}

uintmax_t bd_count_0(uintmax_t N, uintmax_t B, uintmax_t b)
{
    return bd_count(N, B, b, 0);
}

uintmax_t bd_count_c(uintmax_t N, uintmax_t B, uintmax_t b)
{
    if (N == UINTMAX_MAX || B == 0 || b >= B)
        return UINTMAX_MAX;

    uintmax_t base = N / B;
    uintmax_t R = N % B;
    if (R == 0)
        return base;

    /* Center heavy buckets when sparse (or tied), light buckets when dense. */
    bool dense = R > B / 2;
    uintmax_t minority = dense ? B - R : R;

    /* Split a longest run between minority buckets across the endpoints. */
    uintmax_t left_rotation = ((B - 1) / minority) / 2;
    uintmax_t shift = left_rotation == 0 ? 0 : B - left_rotation;
    uintmax_t marked = bd_count(minority, B, b, shift);

    /* marked is 0 or 1; complement it when it marks a light bucket. */
    return base + (dense ? 1 - marked : marked);
}