#include "game/sim/colony.h"

#include <algorithm>
#include <cmath>

#include "game/sim/numeric.h"

namespace sots::sim {

namespace {
constexpr std::int64_t kMaxPopStep = 50000000;
}

double HazardModifier(double suitability, double idealSuitability, double suitTolerance) {
    // Every step is double: the 0.1 in the image is a true double (0x3fb999999999999a),
    // not a widened float, and nothing here is narrowed. The clamp is applied low end
    // first, then high end, which is what lets a NaN fall through unchanged.
    const double band = suitTolerance + 0.1;
    const double v = 1.0 - std::fabs(suitability - idealSuitability) / band;
    if (!(v == v)) return v;  // 0/0 when the band is zero and the planet is at the ideal
    const double lo = v > 0.0 ? v : 0.0;
    return lo < 1.0 ? lo : 1.0;
}

std::int64_t QuantiseCapacity(std::int64_t v) {
    // The capacity helper rounds down to a whole ten -- but only above ten, so small values
    // and negatives pass through untouched.
    return v > 10 ? (v / 10) * 10 : v;
}

std::int64_t CarryingCapacity(const CapacityInputs& in, const TuningTable& t) {
    if (IsNpcSpecies(in.species)) return 0;
    if (!in.speciesCanLive) return 0;

    // `Size x 1e8` is an exact 64-bit *integer* product (the 1e8 is the immediate
    // 0x05f5e100), and the whole floating-point factor is formed first, in this order.
    const std::int64_t base = static_cast<std::int64_t>(in.planetSize) * 100000000LL;
    const double factor = in.hazardMod * (in.groupCapacityMult * in.speciesGrowthFactor *
                                          (in.ownerIsDifferentSpecies ? in.crossSpeciesMod : 1.0));
    std::int64_t result = QuantiseCapacity(Ftoi64(static_cast<double>(base) * factor));
    if (in.arcologyTech) {
        if (in.group == PopGroup::Imperial) result += 100000000;
        else if (in.group == PopGroup::Civilian) result += 200000000;
        // Slaves get nothing: the helper falls through both decrements and returns 0.
    }
    if (in.groupMaxEnabled) result = std::min(result, in.groupMax);
    if (in.ownerIsNpc) {
        result = QuantiseCapacity(
            Ftoi64(static_cast<double>(result) * t.INDSYS_IMPERIAL_POPULATION_MOD));
    }
    return result;
}

double GrowthSuitabilityDistance(double suit, double ideal, double tolerance,
                                 bool accommodated) {
    if (accommodated) return 0.0;
    const double clamped = F32(ClampT(suit, 0.0, 20.0));
    const double d = F32(std::fabs(F32(ideal - clamped)));
    return tolerance < d ? tolerance : d;
}

double PopulationGrowthFraction(const GrowthInputs& in, const TuningTable& t) {
    const double d = GrowthSuitabilityDistance(in.suitability, in.idealSuitability,
                                               in.suitTolerance, in.accommodated);
    const double q = F32(d / in.suitTolerance);  // 0/0 when the tolerance is zero: NaN
    const double base = F32(1.0 - Clamp01(q));
    const double exponent = ClampT(t.POPULATION_GROWTH_EXP, kGrowthExponentMin, kGrowthExponentMax);
    double g = F32(std::pow(base, exponent));
    g = F32(Clamp01(g));
    if (!(g > 0)) return g;  // zero, or a NaN that must survive
    if (t.POPULATION_GROWTH_MOD > 0) g = F32(g * t.POPULATION_GROWTH_MOD);
    if (in.playerPopMod > 0) g = F32(g * in.playerPopMod);
    if (in.extraFactor > 0) g = F32(g * in.extraFactor);
    if (in.groupGrowthMult > 0) g = F32(g * in.groupGrowthMult);
    return g;
}

std::int64_t PopulationGrowthDelta(const GrowthInputs& in, const TuningTable& t) {
    if (in.blockaded) return 0;
    if (in.pop <= 0) return 0;
    const double g = PopulationGrowthFraction(in, t);
    if (g == 0.0) return 0;
    const std::int64_t delta = Ftoi64(static_cast<double>(in.pop) * g);
    if (delta != 0) return delta;
    // A truncated-to-zero delta becomes exactly 1, but only for a strictly positive
    // fraction -- a NaN fraction leaves the colony alone.
    return g > 0 ? 1 : 0;
}

std::int64_t ApplyImperialGrowth(std::int64_t pop, std::int64_t capacity, std::int64_t delta) {
    if (delta >= 0 && delta > kMaxPopulationStep) delta = kMaxPopulationStep;
    const std::int64_t grown = pop + delta;
    std::int64_t result;
    if (grown <= capacity) {
        result = grown;
    } else if (pop <= capacity) {
        result = capacity;  // growing past the cap simply lands on it: no shrink
    } else {
        std::int64_t over = capacity > pop ? capacity - pop : pop - capacity;  // |cap - pop|
        if (over >= kMaxPopulationStep) over = kMaxPopulationStep;
        const std::int64_t floorValue = pop >= 100 ? 100 : pop;
        const std::int64_t shrunk = pop - over;
        result = shrunk > floorValue ? shrunk : floorValue;
    }
    return result < 0 ? 0 : result;
}

double InfrastructurePointsNeeded(double infra) {
    // A real ceil() on the double quotient; the result stays a double and is compared
    // against the (also double) point pool, so nothing is truncated on the way.
    return std::ceil((1.0 - infra) / kInfraPerPoint);
}

double InfrastructureGain(double points) {
    // Three separate 80-bit steps in this exact order, then a clamp at zero and one
    // narrowing to float32 on the store. Folding them into one x3.3e-5 multiply is not
    // bit-identical.
    const double v = points / 500.0 * 0.01 * 1.65;
    return F32(v > 0.0 ? v : 0.0);
}

double ApplyInfrastructurePoints(double infra, double pool, double* pointsUnused) {
    const double spend = std::min(pool, InfrastructurePointsNeeded(infra));
    if (pointsUnused) *pointsUnused = std::max(pool - spend, 0.0);
    return InfrastructureGain(spend);
}

double ApplyInfrastructureDelta(double infra, double delta) {
    if (infra >= 1.0) return infra;  // the apply is a no-op once the colony is built out
    const double v = delta + infra;
    return F32(v > 1.0 ? 1.0 : v);
}

double DecayUnownedInfrastructure(double infra) {
    // The image holds `(double)0.02f`, not the decimal 0.02, and `Infra` is a 4-byte float, so
    // the difference is rounded to single precision before the floor test. The floor itself is
    // `result <= 0 -> 0`, i.e. an exact zero also takes the zero branch.
    const double v = F32(infra - kUnownedInfraDecay);
    return v > 0.0 ? v : 0.0;
}

double TerraformPointsNeeded(double suit, double ideal, double terraMod) {
    // The per-point yield -- terraforming modifier included -- is folded into the *need*,
    // so a player with a better modifier needs proportionally fewer points. The gap is
    // narrowed to float32 before the fabs.
    const double gap = std::fabs(F32(ideal - suit));
    const double per = std::fabs(terraMod * kTerraformNeedFactor / 20000.0);
    return std::ceil(gap / per);
}

double TerraformDelta(double points, double terraMod, double suit, double ideal) {
    const double sign = suit > ideal ? -1.0 : 1.0;  // strictly greater; at the ideal, +1
    return F32(points * 1.5 * kTerraform12 * terraMod * sign / 20000.0);
}

double ApplyTerraformDelta(double suit, double delta, double ideal) {
    const double v = suit + delta;
    // The apply clamps on whichever side it approached the ideal from, so suitability
    // never overshoots.
    if (suit < ideal) return F32(v > ideal ? ideal : v);
    if (suit > ideal) return F32(v < ideal ? ideal : v);
    return suit;
}

double SlaveDeathRate(double slaveOutputRate, double suit, double ideal,
                      const SpeciesTechFlags& flags, const TuningTable& t, bool owned) {
    if (!owned) return 1.0;  // an unowned system reports the full rate, not zero
    // The mod chain: 0.8f replaces 1 outright for the first translation tech, then each of
    // the next two subtracts (double)0.2f. There is no clamp at zero.
    double mod = flags.translation1 ? kSlaveModBase : 1.0;
    if (flags.translation2) mod = F32(mod - kSlaveModStep);
    if (flags.translation3) mod = F32(mod - kSlaveModStep);
    // Term order matters because every step is stored back to a float32: the hazard term is
    // added to the base first, and only then the output term.
    const double base = F32(t.SLAVES_DEATH_RATE);
    const double hazard = F32(std::fabs(F32(ideal - suit)));
    double v = F32(F32(hazard * t.SLAVES_DEATH_RATE_BYHAZARD) + base);
    v = F32(v + F32(slaveOutputRate * t.SLAVES_DEATH_RATE_BYOUTPUT));
    return F32(v * mod);
}

std::int64_t SlaveDeaths(std::int64_t slaves, double rate, double plagueRate,
                         const TuningTable& t) {
    // The worst plague at the system contributes an *additive* term to the rate, not a
    // multiplier -- easy to miss, and it is the only path by which a plague kills slaves.
    std::int64_t d = Ftoi64(static_cast<double>(slaves) * (rate + plagueRate));
    // Either bound is disabled by ANY negative value, not specifically by -1.
    if (t.SLAVES_MIN_DEATHS >= 0 && d < t.SLAVES_MIN_DEATHS) d = t.SLAVES_MIN_DEATHS;
    if (t.SLAVES_MAX_DEATHS >= 0 && d > t.SLAVES_MAX_DEATHS) d = t.SLAVES_MAX_DEATHS;
    if (d < 0) d = 0;
    if (d > slaves) d = slaves;
    return d;
}

OutputRates NormaliseOutputRates(const OutputRates& raw, bool suitAtIdeal, bool infraFull) {
    OutputRates r = raw;
    // Suppression first: a planet already at its ideal cannot terraform, and a colony
    // whose infrastructure plus pending bonus has reached 1 cannot build more.
    if (suitAtIdeal) r.terraform = 0.0;
    if (infraFull) r.infra = 0.0;

    // A slider at or below the threshold counts as off. The threshold is the widened
    // float literal (double)1e-4f.
    if (!(r.trade > kOutputRateThreshold)) r.trade = 0.0;
    if (!(r.construction > kOutputRateThreshold)) r.construction = 0.0;
    if (!(r.infra > kOutputRateThreshold)) r.infra = 0.0;
    if (!(r.terraform > kOutputRateThreshold)) r.terraform = 0.0;

    // Only the pinned channel is clamped into [0, 1]; the other three are left alone
    // because they are about to be rescaled anyway.
    r.trade = Clamp01(r.trade);

    // The sum excludes the pinned channel, and it is rounded to float32 at every step.
    double sum = 0;
    sum = F32(sum + r.construction);
    sum = F32(sum + r.infra);
    sum = F32(sum + r.terraform);
    if (sum == 0.0) {
        // The all-zero fallback seeds the three unpinned channels with the threshold value
        // itself -- respecting the two suppressions -- and re-sums.
        r.construction = kOutputRateThreshold;
        r.infra = infraFull ? 0.0 : kOutputRateThreshold;
        r.terraform = suitAtIdeal ? 0.0 : kOutputRateThreshold;
        sum = 0;
        sum = F32(sum + r.construction);
        sum = F32(sum + r.infra);
        sum = F32(sum + r.terraform);
        if (sum == 0.0) return r;  // everything suppressed: nothing left to rescale
    }
    const double share = 1.0 - r.trade;
    r.construction = F32(F32(r.construction / sum) * share);
    r.infra = F32(F32(r.infra / sum) * share);
    r.terraform = F32(F32(r.terraform / sum) * share);
    return r;
}

PopTypeConstants PopTypeOf(PopGroup g, const TuningTable& t) {
    // The three rows are built in code from x87 literals; only the slave row reads the
    // data files. `SLAVES_OUTPUT_MOD` / `SLAVES_INCOME_MOD` are not yet fields of
    // TuningTable, so the slave row's modifiers arrive as zero until a loader supplies
    // them -- which is the honest state, not a silent 1.0.
    switch (g) {
        case PopGroup::Imperial:
            return PopTypeConstants{1.0, 1.0, 50000000};
        case PopGroup::Civilian:
            // 0.33 is a float literal in the image, so it widens to 0.33000001311302185.
            return PopTypeConstants{F32(0.33), F32(0.33), 20000000};
        case PopGroup::Slaves:
            return PopTypeConstants{t.SLAVES_OUTPUT_MOD, t.SLAVES_INCOME_MOD, 0};
    }
    return PopTypeConstants{};
}

double MoraleOutputMultiplier(int morale, const TuningTable& t) {
    // An entry of exactly zero means "no morale record for this species": the original
    // returns 1 without consulting the thresholds.
    if (morale == 0) return 1.0;
    if (morale >= t.MORALE_INCREASE_OUTPUT) {
        // The modifier is used only when strictly positive; otherwise the multiplier is 1.
        return t.MORALE_INCREASE_OUTPUT_MOD > 0.0 ? t.MORALE_INCREASE_OUTPUT_MOD : 1.0;
    }
    if (morale <= t.MORALE_DECREASE_OUTPUT) {
        return t.MORALE_DECREASE_OUTPUT_MOD > 0.0 ? t.MORALE_DECREASE_OUTPUT_MOD : 1.0;
    }
    return 1.0;
}

double GroupOutput(const GroupOutputInputs& in, const TuningTable& t) {
    const double count = static_cast<double>(in.count);
    if (!(count > 0.0)) return 0.0;

    const double q = count / kOutputPopulationDivisor;

    double stationFactor = 1.0;
    if (in.owned && in.group == PopGroup::Imperial) {
        const double b =
            t.STATION_BONUS_IMPERIAL_OUTPUT > 0.0 ? t.STATION_BONUS_IMPERIAL_OUTPUT : 0.0;
        stationFactor = 1.0 + static_cast<double>(in.stations) * b;
    }

    double morale = 1.0;
    if (in.group == PopGroup::Civilian && in.owned && !in.independent) {
        morale = MoraleOutputMultiplier(in.morale, t);
    }

    // The original's association, and each step rounded to double the way its x87 does.
    const double sf18 = Narrow(stationFactor * kOutputPopulationFactor);
    const double a = Narrow(PopTypeOf(in.group, t).outputMod * sf18);
    const double b = Narrow(a * morale);
    const double v = Narrow(b * q);
    return v > 0.0 ? v : 0.0;
}

float StripMineFraction(const StripMineInputs& in) {
    const double fi = F32(static_cast<double>(in.infraBonus) + static_cast<double>(in.infra));
    const double pop = static_cast<double>(in.population) / 100.0;
    // The original's helper is an odd-symmetric cube root: pow(|x|, 1/3) with the sign
    // carried through, using the double 1/3 rather than std::cbrt.
    const double root = pop >= 0.0 ? std::pow(pop, 1.0 / 3.0) : -std::pow(-pop, 1.0 / 3.0);
    double r = Clamp01(root * 0.01);
    if (0.0001 + r >= 1.0) r = fi;
    return static_cast<float>(r < fi ? r : fi);
}

double OverHarvestDemand(const OverHarvestInputs& in) {
    const double avail = static_cast<double>(in.resourcesAvailable);
    double b = 0.0;
    if (in.overHarvestRate > 0.0) {
        // The population sum is an int32 add in the original and the product is formed as
        // rate x available x scale, in that order.
        const double scale = Clamp01(static_cast<double>(in.population) * 1e-05);
        const double v = in.overHarvestRate * avail * scale;
        b = v > 1.0 ? v : 1.0;
    }
    const double base = in.owned ? static_cast<double>(in.speciesBaseDemand) : 0.0;
    const double t = base + b;
    const double lo = t > 0.0 ? t : 0.0;
    return avail < lo ? avail : lo;
}

double SystemBaseOutput(const BaseOutputInputs& in, const TuningTable& t) {
    OverHarvestInputs oh;
    oh.overHarvestRate = in.overHarvestRate;
    oh.resourcesAvailable = in.resourcesAvailable;
    oh.population = in.imperialPopulation;
    oh.speciesBaseDemand = in.speciesBaseDemand;
    const double harvestTerm = Narrow(OverHarvestDemand(oh) * F32(in.speciesResourceOutput));

    StripMineInputs sm;
    sm.population = in.imperialPopulation;
    sm.infra = in.infra;
    sm.infraBonus = in.infraBonus;
    const double resourceTerm =
        Narrow(Narrow(static_cast<double>(in.transitResources + in.resourcesAvailable) *
                      static_cast<double>(StripMineFraction(sm))) *
               0.9);

    GroupOutputInputs g;
    g.stations = in.stations;
    g.independent = in.independent;

    g.group = PopGroup::Imperial;
    g.count = in.imperialPopulation;
    g.morale = in.imperialMorale;
    const double imperial = GroupOutput(g, t);

    g.group = PopGroup::Civilian;
    g.count = in.civilianPopulation;
    g.morale = in.civilianMorale;
    const double civilian = GroupOutput(g, t);

    g.group = PopGroup::Slaves;
    g.count = in.slavePopulation;
    g.morale = 0;
    const double slaves = GroupOutput(g, t);

    // The original's association: ((slaves + (civilian + (imperial + 0.0))) + (harvest + resource)).
    const double popTerm = Narrow(slaves + Narrow(civilian + Narrow(imperial + 0.0)));
    return Narrow(popTerm + Narrow(harvestTerm + resourceTerm));
}

double TotalSystemOutputRaw(const OutputModifiers& m, const TuningTable& t) {
    if (!m.owned || m.rebelling) return 0.0;
    const double addiction = m.addictionPhase3 ? t.ADDICTION_OUTPUT_MOD : 1.0;
    double v = m.baseOutput;
    v = Narrow(v * m.playerOutMod);
    v = Narrow(v * m.systemOutMod);
    v = Narrow(v * m.techOutMod);
    v = Narrow(v * m.rebOutMod);
    v = Narrow(v * m.scOutMod);
    return Narrow(addiction * v);
}

double TotalSystemOutput(const OutputModifiers& m, const TuningTable& t) {
    // Rounded half-to-even and kept as a double: the channel splits multiply this value,
    // and only the reported slot 0 truncates it to an int.
    return RoundHalfEven(TotalSystemOutputRaw(m, t));
}

OutputSplit SplitOutput(double total, const OutputRates& rates) {
    // Each channel is rounded half-to-even independently and stays a double.
    OutputSplit s;
    s.trade = RoundHalfEven(total * rates.trade);
    s.construction = RoundHalfEven(rates.construction * total);
    s.terraform = RoundHalfEven(rates.terraform * total);
    s.infra = RoundHalfEven(rates.infra * total);
    return s;
}

int ConstructionPoints(double constructionShare, int stations, const TuningTable& t) {
    // Truncating, not rounding -- this slot goes through the float-to-int helper.
    // C3 correction, from the instruction stream of 0x00746830: the bonus is ignored unless
    // it is STRICTLY positive (the same unloaded-table guard the output term carries), and
    // the association is `k x (b x cons) + cons`, not `cons x (1 + b x k)`. Both differences
    // are invisible while no system has a shipyard station, which is the whole corpus.
    const double b = t.STATION_BONUS_SHIPCON > 0.0 ? t.STATION_BONUS_SHIPCON : 0.0;
    return Ftol(static_cast<double>(stations) * (b * constructionShare) + constructionShare);
}

OutputSplit SplitLeftover(double leftover, const OutputRates& rates, bool suitAtIdeal,
                          bool infraFull) {
    double wt, wf, wi;
    // The "construction was the only slider" case is an exact equality against 1, not a
    // >= test: after normalisation a pure-construction colony has exactly 1 there.
    if (rates.construction == 1.0) {
        wt = 1.0;
        wf = suitAtIdeal ? 0.0 : 1.0;
        wi = infraFull ? 0.0 : 1.0;
    } else {
        wt = rates.trade;
        wf = rates.terraform;
        wi = rates.infra;
    }
    // C3 correction: the original accumulates `wi + (wf + wt)` on the x87 stack, in that
    // association. Reordering it is not free in floating point.
    const double sum = wi + (wf + wt);
    OutputSplit s;
    if (sum <= 0 || leftover <= 0) {
        s.trade = std::max(0.0, leftover);
        return s;
    }
    // Independently rounded, so the three need not add back up to the leftover.
    s.trade = RoundHalfEven(wt * leftover / sum);
    s.terraform = RoundHalfEven(wf * leftover / sum);
    s.infra = RoundHalfEven(wi * leftover / sum);
    return s;
}

int GroupIncome(PopGroup group, std::int64_t count, const TuningTable& t) {
    // `fld DWORD [row+0x14]` -- the income column is a float32 in the table.
    const double mod = F32(PopTypeOf(group, t).incomeMod);
    return Ftol(Narrow(mod * (static_cast<double>(count) / kIncomePopulationDivisor)));
}

double PopulationIncome(PopGroup group, const PopIncomeRow (&rows)[kSpeciesCount],
                        bool owned, bool independent, const TuningTable& t) {
    double sum = 0.0;
    for (int sp = 0; sp < kSpeciesCount; ++sp) {
        const PopIncomeRow& r = rows[sp];
        if (!(r.count > 0)) continue;
        // Only the civilian row takes morale, and it takes it through the same helper the
        // output term uses -- including the "no owner / independent / no record" bypass.
        double morale = 1.0;
        if (group == PopGroup::Civilian && owned && !independent) {
            morale = MoraleOutputMultiplier(r.morale, t);
        }
        const double addiction = r.addicted ? t.ADDICTION_INCOME_MOD : 1.0;
        const double base = static_cast<double>(GroupIncome(group, r.count, t));
        // The SECOND truncation: per species, after both factors.
        sum += static_cast<double>(Ftol(Narrow(Narrow(base * morale) * addiction)));
    }
    return sum;
}

double SuitabilityCostMod(double suitability, double idealSuitability, double suitTolerance,
                          bool rebelAI, bool owned, bool vonNeumann) {
    if (vonNeumann) return 0.0;   // the original's very first test
    if (!owned) return 20.0;      // ... and it logs a warning
    if (rebelAI) return 0.0;
    // All three operands are 4-byte floats in the original (`server->IdealSuit[sp]`,
    // `sys->Suit`, `owner->SuitTol`); the subtraction and the compare are then done on the
    // x87 with no store back, so nothing is narrowed here. The compare is `<=`: a distance
    // exactly at the tolerance is charged as itself, not as the cap.
    const double d = std::fabs(idealSuitability - suitability);
    return d <= suitTolerance ? d : suitTolerance;
}

double SystemMoneyIncomeRaw(const SystemMoneyInputs& in) {
    // Whole blocks of five trade points; the same literal is the modulus and the multiplier.
    // The `+ 0.0` is a real instruction (an .rdata zero) and the association below is the
    // original's -- x87 addition is not associative, so neither is reorderable.
    const double blocks = (in.tradePoints - std::fmod(in.tradePoints, 5.0)) * 5.0;
    double t = Narrow(Narrow(blocks + 0.0) + in.popIncomeImperial);
    t = Narrow(in.popIncomeCivilian + t);
    t = Narrow(in.slaveIncome + t);

    // Every per-player multiplier is stored back through a 4-byte float before it is used.
    t = Narrow(F32(in.speciesIncomeFactor) * t);
    const double diff = F32(F32(in.difficultyIncomeMult) * F32(in.serverIncomeMod));
    t = Narrow(diff * Narrow(F32(in.playerIncMod) * t));

    const double cost =
        Narrow(F32(in.speciesCostFactor) * Narrow(Narrow(in.suitCostMod * 10000.0) * 1.5));
    return t - cost;
}

int SystemMoneyIncome(const SystemMoneyInputs& in) { return Ftol(SystemMoneyIncomeRaw(in)); }

int SystemMaxIncome(double totalOutput, const SystemMoneyInputs& in) {
    SystemMoneyInputs m = in;
    // The max-income rate vector is trade = 1 and every other channel 0, so the trade
    // points are the rounded total and the two cascade channels contribute nothing.
    m.tradePoints = RoundHalfEven(totalOutput);
    const int money = SystemMoneyIncome(m);
    return money > 0 ? money : 0;
}

double IdealSuitability(const IdealSuitabilityInputs& in) {
    if (!in.owned) return in.systemSuitability;
    double v = in.ownerIdealSuitability;
    if (in.independent) v = in.serverIdealSuitability;
    // An `!=` against the sentinel, so a NaN override would also win. Nothing in the corpus
    // exercises either side of that.
    if (in.systemOverride != kIdealSuitabilityNoOverride) v = in.systemOverride;
    return v;
}

RepairPassResult RepairShipsInOrbit(int points, int repairDemand) {
    RepairPassResult r;
    if (points <= 0 || repairDemand <= 0) {
        r.left = points;
        return r;
    }
    r.spent = points < repairDemand ? points : repairDemand;
    r.left = points - r.spent;
    return r;
}

SystemOutput ComputeSystemOutput(const SystemOutputInputs& in, const TuningTable& t) {
    SystemOutput o;
    const OutputRates r = NormaliseOutputRates(in.rates, in.suitAtIdeal, in.infraFull);
    o.normalisedRates = r;

    // One rounding of the total, then one rounding per channel off that same value.
    const double total = RoundHalfEven(in.totalOutputRaw);
    o.totalOutput = Ftol(total);
    const OutputSplit split = SplitOutput(total, r);
    o.tradePoints = split.trade;

    // --- construction: the queue first, then the repair pass -----------------------------
    o.construction = ConstructionPoints(split.construction, in.shipyardStations, t);
    o.constructionToQueue =
        in.buildQueueDemand < o.construction ? in.buildQueueDemand : o.construction;
    int rem = o.construction - o.constructionToQueue;
    if (rem < 0) rem = 0;
    if (rem > 0) {
        const RepairPassResult rep = RepairShipsInOrbit(rem, in.repairDemand);
        o.constructionToRepair = rep.spent;
        rem = rep.left;
    }

    // --- the leftover redistribution -----------------------------------------------------
    // `SplitLeftover`'s `infraFull` argument is the RAW `Infra == 1` test, not the
    // `Infra + ibon >= 1` one the normaliser used.
    const OutputSplit left =
        rem > 0 ? SplitLeftover(static_cast<double>(rem), r, in.suitAtIdeal, in.infraExactlyOne)
                : OutputSplit{};
    o.leftoverToTrade = left.trade;

    // --- infrastructure ------------------------------------------------------------------
    const double infraNeed = std::ceil((1.0 - in.infra) / 3.3e-5);
    const double poolInfra = left.infra + split.infra;
    const double spendInfra = poolInfra < infraNeed ? poolInfra : infraNeed;
    double leftInfra = poolInfra - spendInfra;
    if (!(leftInfra > 0.0)) leftInfra = 0.0;
    // Three separate 80-bit steps, not one x3.3e-5.
    const double infraGain = spendInfra / 500.0 * 0.01 * 1.65;
    o.infraDelta = F32(infraGain > 0.0 ? infraGain : 0.0);

    // --- terraforming: the infrastructure leftover lands in THIS pool ---------------------
    const double terraNeed = std::ceil(in.terraformPointsNeeded);
    const double poolTerra = (left.terraform + split.terraform) + leftInfra;
    const double spendTerra = poolTerra < terraNeed ? poolTerra : terraNeed;
    double leftTerra = poolTerra - spendTerra;
    if (!(leftTerra > 0.0)) leftTerra = 0.0;
    o.leftoverToMoney = leftTerra;
    // The same helper the colony pass uses; its sign test is `suit > ideal`, which is the
    // original's `IdealSuitability() < Suit` with the operands swapped.
    o.suitabilityDelta = TerraformDelta(spendTerra, in.terraformMod,
                                        in.terraformDown ? 1.0 : 0.0, 0.0);

    // --- money ----------------------------------------------------------------------------
    SystemMoneyInputs m = in.money;
    m.tradePoints = (o.leftoverToTrade + o.tradePoints) + leftTerra;
    o.money = SystemMoneyIncome(m);
    return o;
}

BonusApplyResult ApplyPopulationBonus(std::int64_t& pop, std::int64_t capacity,
                                      std::int64_t& pendingBonus, bool owned, bool homeSystem) {
    BonusApplyResult r;
    if (pendingBonus <= 0) return r;
    if (!owned) {  // an unowned system keeps no pool at all
        pendingBonus = 0;
        return r;
    }
    if (pop >= capacity) return r;
    r.resetTurnsDeveloping = !homeSystem;
    const std::int64_t applied = std::min(capacity - pop, pendingBonus);
    pop += applied;
    pendingBonus -= applied;
    r.applied = applied != 0;
    return r;
}

BonusApplyResult ApplyInfrastructureBonus(double& infra, double& pendingBonus, bool homeSystem) {
    BonusApplyResult r;
    if (!(pendingBonus > 0)) return r;
    if (infra >= 1.0) return r;
    r.resetTurnsDeveloping = !homeSystem;
    const double room = F32(1.0 - infra);
    const double applied = std::min(pendingBonus, room);
    // The original stores the literal 1.0 rather than the sum when the pool covers the whole
    // remainder, so a colony topped up this way lands exactly on 1 with no rounding residue.
    infra = applied == room ? 1.0 : F32(infra + applied);
    pendingBonus = F32(pendingBonus - applied);
    r.applied = true;
    return r;
}

void AccrueSystemBonus(const SystemBonusInputs& in, std::int64_t& popBonus, double& infraBonus,
                       const TuningTable& t) {
    if (!in.stable) return;
    if (in.turnsOwned <= t.SYSTEMBONUS_MINTURNS) return;
    if (in.turnsDeveloping <= t.SYSTEMBONUS_MINTURNS) return;

    const double cap = static_cast<double>(in.capacity);
    const std::int64_t popTarget =
        in.ownerSpeciesEligible ? Ftol(std::max(t.SYSTEMBONUS_POPBONUS, 0.0) * cap) : 0;
    const std::int64_t popInc = Ftol(t.SYSTEMBONUS_POPBONUS_INC * cap);
    popBonus += std::min(std::max<std::int64_t>(popInc, 0), std::max<std::int64_t>(popTarget - popBonus, 0));

    // The original reads the infrastructure target through a helper that already returns 0
    // for an ineligible owner, and clamps it at 0 before differencing.
    const double infraTarget =
        std::max(0.0, in.ownerSpeciesEligible ? t.SYSTEMBONUS_INFRABONUS : 0.0);
    infraBonus += std::min(std::max(t.SYSTEMBONUS_INFRABONUS_INC, 0.0), std::max(infraTarget - infraBonus, 0.0));
}

BuildQueueResult ProcessBuildQueue(std::vector<BuildOrder>& queue, int points) {
    BuildQueueResult r;
    if (points > 0) {
        for (BuildOrder& o : queue) {
            if (o.constructionLeft > points) {
                o.constructionLeft -= points;
                points = 0;
                break;
            }
            if (o.moneyCost > 0) {
                if (!o.moneyAvailable) continue;  // refused: skip this order, keep going
                r.moneyCharged = SaturatingAdd(r.moneyCharged, o.moneyCost);
            }
            points -= o.constructionLeft;
            o.constructionLeft = 0;
            r.completedOrderIds.push_back(o.orderId);
        }
    }
    // The removal sweep is separate and unconditional: every order at or below zero goes,
    // including one that was already finished before this turn.
    queue.erase(std::remove_if(queue.begin(), queue.end(),
                               [](const BuildOrder& o) { return o.constructionLeft <= 0; }),
                queue.end());
    r.pointsLeft = points;
    return r;
}

// ---------------------------------------------------------------------------------------
// Countdown nibbles
// ---------------------------------------------------------------------------------------

int CountdownFor(const ColonyCountdowns& c, int player) {
    if (player < 0 || player >= kMaxCountdownPlayers) return 0;
    return static_cast<int>((c.counters >> (4 * player)) & 0xfULL);
}

void SetCountdown(ColonyCountdowns& c, int player, int value) {
    if (player < 0 || player >= kMaxCountdownPlayers) return;  // the original refuses too
    const std::uint64_t v = static_cast<std::uint64_t>(ClampT(value, 0, 15));
    const std::uint64_t mask = 0xfULL << (4 * player);
    c.counters = (c.counters & ~mask) | (v << (4 * player));
}

void TickCountdowns(ColonyCountdowns& c, int playerCount) {
    if (c.counters == 0) return;  // the whole sweep is skipped when nothing is counting
    for (int i = 0; i < playerCount; ++i) {
        if (i >= kMaxCountdownPlayers) continue;  // no nibble exists; neither word is touched
        const int v = CountdownFor(c, i);
        if (v != 0) SetCountdown(c, i, v - 1);
        else c.active &= ~(1u << i);
    }
}

// ---------------------------------------------------------------------------------------
// Addiction
// ---------------------------------------------------------------------------------------

AddictionPhase AddictionPhaseOf(int startTurn, int currentTurn, int phase2Start, int phase3Start) {
    if (startTurn == 0) return AddictionPhase::None;
    const int elapsed = currentTurn - startTurn;
    if (elapsed > phase3Start) return AddictionPhase::Terminal;
    if (elapsed > phase2Start) return AddictionPhase::Established;
    return AddictionPhase::Onset;
}

// ---------------------------------------------------------------------------------------
// The per-system turn pass
// ---------------------------------------------------------------------------------------

ColonyTurnResult ProcessColonyTurn(ColonyTurnState& s, const ColonyTurnInputs& in,
                                   const TuningTable& t) {
    ColonyTurnResult r;

    // 1. An unowned system's infrastructure rots. An owned one is left alone here; its
    //    infrastructure moves in the output pass, which is a different phase of the turn.
    if (!in.owned) s.infra = DecayUnownedInfrastructure(s.infra);

    // 2. The two pending bonus pools are drained into the colony, owned or not. Either one
    //    resets the turns-developing counter when it fires on a colony that is not the
    //    owner's home system -- so a colony that is still absorbing a bonus never reaches
    //    the system-bonus gate, which is a real feedback loop and easy to miss.
    const BonusApplyResult ib = ApplyInfrastructureBonus(s.infra, s.infraBonus, in.homeSystem);
    const BonusApplyResult pb =
        ApplyPopulationBonus(s.pop, in.imperialCapacity, s.popBonus, in.owned, in.homeSystem);
    if (in.owned && (ib.resetTurnsDeveloping || pb.resetTurnsDeveloping)) s.turnsDeveloping = 0;

    // (an independent colony's imperial<->civilian drift runs here; not modelled)

    // 3. Turns-developing: incremented while the colony is stable, reset the moment it is not.
    if (in.stable) ++s.turnsDeveloping;
    else s.turnsDeveloping = 0;

    // 4. Long-stability bonus accrual, which reads the ntdev just written above.
    if (in.owned) {
        SystemBonusInputs b;
        b.stable = in.stable;
        b.turnsOwned = in.turnsOwned;
        b.turnsDeveloping = s.turnsDeveloping;
        b.ownerSpeciesEligible = in.ownerSpeciesEligible;
        b.capacity = in.imperialCapacity;
        AccrueSystemBonus(b, s.popBonus, s.infraBonus, t);
    }

    // (plague, the build queue, imperial and civilian growth, the resource debit and the
    //  in-orbit refuel run here; all of them are the declared input boundary)

    // 5. The turn's resource total is consumed and reset.
    s.totalResources = 0;

    // 6. Growth halts expire every turn: a blockade has to be re-asserted.
    s.growthHalted[0] = s.growthHalted[1] = s.growthHalted[2] = false;

    // 7. The two per-player countdowns tick.
    TickCountdowns(s.battles, in.playerCount);
    TickCountdowns(s.recon, in.playerCount);

    // (slaves and rebellion run here; not modelled)

    // 8. Addiction. Only an owned, non-independent colony sweeps it, and only species that
    //    actually have civilians here. Temperance suppresses the addiction outright; otherwise
    //    the onset and terminal phases each raise their own morale event and the established
    //    phase (2) raises none at all.
    if (in.owned && !in.independent) {
        for (int i = 0; i < kSpeciesCount; ++i) {
            if (!in.civilianPresent[i]) continue;
            if (in.temperance[i]) {
                r.moraleEvents.push_back({i, kMoraleEventAddictionSuppressed, -1});
                continue;
            }
            switch (AddictionPhaseOf(in.addictionStart[i], in.currentTurn, in.addictionPhase2Start,
                                     in.addictionPhase3Start)) {
                case AddictionPhase::Onset:
                    r.moraleEvents.push_back({i, kMoraleEventAddictionOnset, +1});
                    break;
                case AddictionPhase::Terminal:
                    r.moraleEvents.push_back({i, kMoraleEventAddictionTerminal, -2});
                    break;
                case AddictionPhase::None:
                case AddictionPhase::Established:
                    break;
            }
        }
    }

    return r;
}

}  // namespace sots::sim