#include "game/sim/colony.h"

#include <cmath>

#include "check.h"
#include "game/sim/economy.h"
#include "game/sim/numeric.h"

using namespace sots::sim;

static TuningTable tuning() {
    TuningTable t;
    t.POPULATION_GROWTH_MOD = 1.2;
    t.POPULATION_GROWTH_EXP = 2.0;
    t.INDSYS_IMPERIAL_POPULATION_MOD = 0.1;
    t.SLAVES_DEATH_RATE = 0.05;
    t.SLAVES_DEATH_RATE_BYHAZARD = 0.5;
    t.SLAVES_DEATH_RATE_BYOUTPUT = 0.1;
    t.SLAVES_MIN_DEATHS = 0;
    t.SLAVES_MAX_DEATHS = -1;
    t.MORALE_INCREASE_OUTPUT = 75;
    t.MORALE_INCREASE_OUTPUT_MOD = 1.1;
    t.MORALE_DECREASE_OUTPUT = 25;
    t.MORALE_DECREASE_OUTPUT_MOD = 0.9;
    t.STATION_BONUS_IMPERIAL_OUTPUT = 0.1;
    t.STATION_BONUS_SHIPCON = 0.25;
    t.ADDICTION_OUTPUT_MOD = 0.5;
    t.ADDICTION_INCOME_MOD = 0.9;
    t.SLAVES_INCOME_MOD = 3.0;
    t.SYSTEMBONUS_MINTURNS = 10;
    t.SYSTEMBONUS_POPBONUS = 0.1;
    t.SYSTEMBONUS_POPBONUS_HOME = 0.2;
    t.SYSTEMBONUS_POPBONUS_INC = 0.01;
    t.SYSTEMBONUS_INFRABONUS = 0.2;
    t.SYSTEMBONUS_INFRABONUS_HOME = 0.5;
    t.SYSTEMBONUS_INFRABONUS_INC = 0.05;
    return t;
}

static void test_capacity() {
    TuningTable t = tuning();
    CapacityInputs c;
    c.planetSize = 5;
    CHECK_EQ(CarryingCapacity(c, t), std::int64_t{500000000});
    c.hazardMod = 0.5;
    CHECK_EQ(CarryingCapacity(c, t), std::int64_t{250000000});
    c.arcologyTech = true;
    CHECK_EQ(CarryingCapacity(c, t), std::int64_t{350000000});
    c.group = PopGroup::Civilian;
    CHECK_EQ(CarryingCapacity(c, t), std::int64_t{450000000});
    c.group = PopGroup::Slaves;
    CHECK_EQ(CarryingCapacity(c, t), std::int64_t{250000000});   // no arcology bonus for slaves
    c.group = PopGroup::Imperial;
    c.groupMaxEnabled = true;
    c.groupMax = 300000000;
    CHECK_EQ(CarryingCapacity(c, t), std::int64_t{300000000});
    c.ownerIsNpc = true;
    CHECK_EQ(CarryingCapacity(c, t), std::int64_t{30000000});
    c.ownerIsNpc = false;
    c.ownerIsDifferentSpecies = true;
    c.crossSpeciesMod = 0.5;
    CHECK_EQ(CarryingCapacity(c, t), std::int64_t{225000000});   // 5e8 x 0.5 x 0.5 + 1e8
    c.species = Species::NPC;
    CHECK_EQ(CarryingCapacity(c, t), std::int64_t{0});
    c.species = Species::Liir;
    c.speciesCanLive = false;
    CHECK_EQ(CarryingCapacity(c, t), std::int64_t{0});
    c.speciesCanLive = true;
    c.planetSize = 0;
    CHECK_EQ(CarryingCapacity(c, t), std::int64_t{100000000});   // arcology alone

    // hazard = clamp01(1 - |suit - ideal| / (tol + 0.1))
    CHECK_NEAR(HazardModifier(0.5, 0.5, 0.2), 1.0, 0.0);
    CHECK_NEAR(HazardModifier(0.6, 0.5, 0.2), 1.0 - 0.1 / 0.3, 1e-12);
    CHECK_NEAR(HazardModifier(0.5, 0.65, 0.2), 0.5, 1e-12);          // symmetric
    CHECK_NEAR(HazardModifier(0.8, 0.5, 0.2), 0.0, 0.0);              // at the band edge
    CHECK_NEAR(HazardModifier(0.9, 0.5, 0.2), 0.0, 0.0);
    CHECK_NEAR(HazardModifier(0.55, 0.5, 0.0), 0.5, 1e-12);           // zero tolerance keeps a 0.1 band
    CHECK_NEAR(HazardModifier(0.7, 0.5, 0.0), 0.0, 0.0);
    // both adaptation techs: 0.2 + 0.75 + 1.5 -> band 2.55
    CHECK_NEAR(HazardModifier(0.8, 0.5, 2.45), 1.0 - 0.3 / 2.55, 1e-12);
}

static void test_growth() {
    TuningTable t = tuning();
    // B4: the curve is driven by suitability, not by how full the colony is. `tol` is the
    // owner's SuitTol and doubles as the divisor, so d/tol is the fraction of the habitable
    // band the planet is off by.
    GrowthInputs g;
    g.suitTolerance = 1.0;
    g.idealSuitability = 1.0;
    g.suitability = 0.5;   // half a band off: base 0.5, ^2 = 0.25, x1.2 = 0.3
    g.pop = 500000;
    CHECK_EQ(PopulationGrowthDelta(g, t), std::int64_t{150000});

    g.playerPopMod = 0.5;                         // 75000
    CHECK_EQ(PopulationGrowthDelta(g, t), std::int64_t{75000});
    g.playerPopMod = 1.0;
    g.groupGrowthMult = 2.0;                      // 300000
    CHECK_EQ(PopulationGrowthDelta(g, t), std::int64_t{300000});
    g.groupGrowthMult = 0.0;                      // a zero column is ignored, not applied
    CHECK_EQ(PopulationGrowthDelta(g, t), std::int64_t{150000});
    g.extraFactor = 0.0;                          // so is a zero extra factor
    CHECK_EQ(PopulationGrowthDelta(g, t), std::int64_t{150000});
    g.extraFactor = 1.0;

    g.pop = 0;                                    // an empty group does not grow at all
    CHECK_EQ(PopulationGrowthDelta(g, t), std::int64_t{0});
    g.pop = 1;                                    // trunc(1 x 0.3) == 0 -> forced to 1
    CHECK_EQ(PopulationGrowthDelta(g, t), std::int64_t{1});
    g.pop = 500000;

    g.suitability = 1.0;                          // exactly at the ideal: base 1
    CHECK_EQ(PopulationGrowthDelta(g, t), std::int64_t{600000});   // 1 x 1.2 x 5e5
    g.suitability = 0.0;                          // a whole band off: base 0 -> no growth
    CHECK_EQ(PopulationGrowthDelta(g, t), std::int64_t{0});
    g.suitability = -5.0;                         // clamped up to 0 first, so still a full band
    CHECK_EQ(PopulationGrowthDelta(g, t), std::int64_t{0});
    g.suitability = 0.5;

    g.accommodated = true;                        // suitability ignored entirely: base 1
    CHECK_EQ(PopulationGrowthDelta(g, t), std::int64_t{600000});
    g.accommodated = false;

    g.blockaded = true;
    CHECK_EQ(PopulationGrowthDelta(g, t), std::int64_t{0});
    g.blockaded = false;

    // the suitability distance clamps the planet's own value into [0, 20] before the
    // difference, and is itself capped by the tolerance
    CHECK_NEAR(GrowthSuitabilityDistance(25.0, 0.0, 100.0, false), 20.0, 1e-6);
    CHECK_NEAR(GrowthSuitabilityDistance(-3.0, 5.0, 100.0, false), 5.0, 1e-6);
    CHECK_NEAR(GrowthSuitabilityDistance(0.0, 5.0, 2.0, false), 2.0, 1e-6);
    CHECK_NEAR(GrowthSuitabilityDistance(0.0, 5.0, 2.0, true), 0.0, 0.0);

    // the exponent is clamped into [0.01f, 1000] before pow()
    TuningTable big = t;
    big.POPULATION_GROWTH_EXP = 100000.0;
    TuningTable capped = t;
    capped.POPULATION_GROWTH_EXP = 1000.0;
    CHECK_NEAR(PopulationGrowthFraction(g, big), PopulationGrowthFraction(g, capped), 0.0);
    TuningTable tiny = t;
    tiny.POPULATION_GROWTH_EXP = 1e-9;
    TuningTable floored = t;
    floored.POPULATION_GROWTH_EXP = kGrowthExponentMin;
    CHECK_NEAR(PopulationGrowthFraction(g, tiny), PopulationGrowthFraction(g, floored), 0.0);

    CHECK_EQ(ApplyImperialGrowth(500000, 1000000, 150000), std::int64_t{650000});
    CHECK_EQ(ApplyImperialGrowth(999999, 1000000, 5), std::int64_t{1000000});  // lands on the cap
    CHECK_EQ(ApplyImperialGrowth(2000000, 1000000, 0), std::int64_t{1000000});
    CHECK_EQ(ApplyImperialGrowth(1000000000, 100000000, 0), std::int64_t{950000000});  // shrink capped
    CHECK_EQ(ApplyImperialGrowth(50, 10, 0), std::int64_t{50});     // floor min(pop, 100)
    CHECK_EQ(ApplyImperialGrowth(500, 10, 0), std::int64_t{100});
    // the 50,000,000 cap is on the delta, and it lives in the apply, not in the fraction
    CHECK_EQ(ApplyImperialGrowth(1000, 2000000000, 100000000), std::int64_t{50001000});
    CHECK_EQ(ApplyImperialGrowth(1000, 2000000000, -100), std::int64_t{900});
}

static void test_infra_terraform() {
    CHECK_NEAR(InfrastructurePointsNeeded(0.0), 30304.0, 0.0);   // ceil(30303.03)
    CHECK_NEAR(InfrastructurePointsNeeded(1.0), 0.0, 0.0);
    CHECK_NEAR(InfrastructureGain(500), 0.0165, 1e-7);
    CHECK_NEAR(InfrastructureGain(1000), 0.033, 1e-7);
    CHECK_NEAR(InfrastructureGain(-100), 0.0, 0.0);              // clamped at zero
    double unused = -1;
    double delta = ApplyInfrastructurePoints(0.5, 100000, &unused);
    CHECK_NEAR(unused, 100000.0 - 15152.0, 0.0);                 // ceil(0.5 / 3.3e-5)
    CHECK_NEAR(ApplyInfrastructureDelta(0.5, delta), 1.0, 0.0);  // clamped to exactly 1
    delta = ApplyInfrastructurePoints(0.5, 1000, &unused);
    CHECK_NEAR(unused, 0.0, 0.0);
    CHECK_NEAR(ApplyInfrastructureDelta(0.5, delta), 0.533, 1e-7);
    CHECK_NEAR(ApplyInfrastructureDelta(1.0, 0.5), 1.0, 0.0);    // already built out: no-op
    CHECK_NEAR(DecayUnownedInfrastructure(0.5), 0.48, 1e-7);
    CHECK_NEAR(DecayUnownedInfrastructure(0.01), 0.0, 0.0);
    // the decay constant is the widened float literal, not the decimal 0.02
    CHECK(DecayUnownedInfrastructure(0.5) != 0.5 - 0.02);

    // the terraforming modifier is inside the point count, and the result is a ceil
    CHECK_NEAR(TerraformPointsNeeded(0.5, 0.8, 1.0), 3334.0, 0.0);
    CHECK_NEAR(TerraformPointsNeeded(0.8, 0.5, 1.0), 3334.0, 0.0);
    CHECK_NEAR(TerraformPointsNeeded(0.5, 0.8, 2.0), 1667.0, 0.0);  // twice the modifier, half the points
    CHECK_NEAR(TerraformPointsNeeded(0.5, 0.5, 1.0), 0.0, 0.0);
    CHECK_NEAR(TerraformDelta(1000, 1.0, 0.5, 0.8), 0.09, 1e-7);
    CHECK_NEAR(TerraformDelta(1000, 1.0, 0.8, 0.5), -0.09, 1e-7);
    CHECK_NEAR(TerraformDelta(1000, 1.0, 0.5, 0.5), 0.09, 1e-7);   // at the ideal the sign is +1
    CHECK_NEAR(TerraformDelta(1000, 2.0, 0.5, 0.8), 0.18, 1e-7);
    CHECK_NEAR(TerraformDelta(0, 2.0, 0.5, 0.8), 0.0, 0.0);
    // suitability stops at the ideal from whichever side it came
    CHECK_NEAR(ApplyTerraformDelta(0.5, 0.09, 0.8), 0.59, 1e-7);
    CHECK_NEAR(ApplyTerraformDelta(0.75, 0.09, 0.8), 0.8, 1e-7);
    CHECK_NEAR(ApplyTerraformDelta(0.85, -0.09, 0.8), 0.8, 1e-7);
    CHECK_NEAR(ApplyTerraformDelta(0.8, 0.09, 0.8), 0.8, 0.0);
}

static void test_slaves() {
    TuningTable t = tuning();
    SpeciesTechFlags f;
    // 0.5 x 0.1 + 0.2 x 0.5 + 0.05 = 0.2
    CHECK_NEAR(SlaveDeathRate(0.5, 0.3, 0.5, f, t), 0.2, 1e-6);
    f.translation1 = true;
    CHECK_NEAR(SlaveDeathRate(0.5, 0.3, 0.5, f, t), 0.16, 1e-6);
    f.translation2 = f.translation3 = true;
    CHECK_NEAR(SlaveDeathRate(0.5, 0.3, 0.5, f, t), 0.08, 1e-6);
    SpeciesTechFlags none;
    CHECK_NEAR(SlaveDeathRate(0.0, 0.5, 0.5, none, t), 0.05, 1e-6);   // base rate only
    // an unowned system short-circuits to a rate of 1, not 0
    CHECK_NEAR(SlaveDeathRate(0.5, 0.3, 0.5, none, t, false), 1.0, 0.0);

    SpeciesTechFlags bits = SpeciesTechFlags::FromBits(0x087);          // bits 0,1,2,7
    CHECK(bits.translation1 && bits.translation2 && bits.translation3 && bits.accommodate);
    CHECK(!bits.incorporate && !bits.addict && !bits.temperance && !bits.subjugate && !bits.proliferate);
    CHECK_NEAR(SlaveDeathRate(0.5, 0.3, 0.5, bits, t), 0.08, 1e-6);
    CHECK(SpeciesTechFlags::FromBits(0x100).proliferate);
    CHECK(SpeciesTechFlags::FromBits(0x020).temperance);

    CHECK_EQ(SlaveDeaths(1000, 0.2, 0.0, t), std::int64_t{200});
    CHECK_EQ(SlaveDeaths(0, 0.2, 0.0, t), std::int64_t{0});
    CHECK_EQ(SlaveDeaths(1000, 0.2, 0.1, t), std::int64_t{300});  // the plague rate ADDS
    t.SLAVES_MIN_DEATHS = 300;
    CHECK_EQ(SlaveDeaths(1000, 0.2, 0.0, t), std::int64_t{300});
    CHECK_EQ(SlaveDeaths(100, 0.2, 0.0, t), std::int64_t{100});   // never more than present
    t.SLAVES_MAX_DEATHS = 150;
    CHECK_EQ(SlaveDeaths(1000, 0.2, 0.0, t), std::int64_t{150});
    t.SLAVES_MIN_DEATHS = -1;                                      // any negative disables it
    t.SLAVES_MAX_DEATHS = -7;
    CHECK_EQ(SlaveDeaths(1000, 0.2, 0.0, t), std::int64_t{200});
    t.SLAVES_MIN_DEATHS = 0;
    t.SLAVES_MAX_DEATHS = -1;
    CHECK_EQ(SlaveDeaths(7, 0.2, 0.0, t), std::int64_t{1});       // 1.4 truncates
}

static void test_output() {
    TuningTable t = tuning();
    // B4: the trade slider is PINNED. Everything else is rescaled to what is left of 1.
    OutputRates r = NormaliseOutputRates({1, 1, 1, 1}, false, false);
    CHECK_NEAR(r.trade, 1.0, 0.0);
    CHECK_NEAR(r.construction, 0.0, 0.0);   // nothing left over for the other three
    CHECK_NEAR(r.terraform, 0.0, 0.0);
    CHECK_NEAR(r.infra, 0.0, 0.0);

    r = NormaliseOutputRates({0.25, 0.25, 0.25, 0.25}, false, false);
    CHECK_NEAR(r.trade, 0.25, 0.0);         // untouched
    CHECK_NEAR(r.construction, 0.25, 1e-7);
    CHECK_NEAR(r.terraform, 0.25, 1e-7);
    CHECK_NEAR(r.infra, 0.25, 1e-7);

    r = NormaliseOutputRates({0, 1, 1, 1}, true, false);
    CHECK_NEAR(r.terraform, 0.0, 0.0);
    CHECK_NEAR(r.construction, 0.5, 1e-7);
    CHECK_NEAR(r.infra, 0.5, 1e-7);

    // all-zero: the three unpinned channels split evenly, trade stays at zero
    r = NormaliseOutputRates({0, 0, 0, 0}, false, false);
    CHECK_NEAR(r.trade, 0.0, 0.0);
    CHECK_NEAR(r.construction, 1.0 / 3.0, 1e-6);
    CHECK_NEAR(r.terraform, 1.0 / 3.0, 1e-6);
    CHECK_NEAR(r.infra, 1.0 / 3.0, 1e-6);
    // ... and the suppressions still apply inside the fallback
    r = NormaliseOutputRates({0, 0, 0, 0}, true, true);
    CHECK_NEAR(r.construction, 1.0, 1e-7);
    CHECK_NEAR(r.terraform, 0.0, 0.0);
    CHECK_NEAR(r.infra, 0.0, 0.0);

    // a slider at or below the threshold counts as off
    r = NormaliseOutputRates({kOutputRateThreshold, 1, 0, 0}, false, false);
    CHECK_NEAR(r.trade, 0.0, 0.0);
    CHECK_NEAR(r.construction, 1.0, 1e-7);

    r = NormaliseOutputRates({-1, 3, 0, 1}, false, true);
    CHECK_NEAR(r.trade, 0.0, 0.0);
    CHECK_NEAR(r.construction, 1.0, 1e-7);
    CHECK_NEAR(r.infra, 0.0, 0.0);

    CHECK_NEAR(MoraleOutputMultiplier(80, t), 1.1, 0.0);
    CHECK_NEAR(MoraleOutputMultiplier(75, t), 1.1, 0.0);
    CHECK_NEAR(MoraleOutputMultiplier(50, t), 1.0, 0.0);
    CHECK_NEAR(MoraleOutputMultiplier(25, t), 0.9, 0.0);
    // A morale entry of exactly 0 means "no record": the thresholds are not consulted,
    // which matters because 0 <= MORALE_DECREASE_OUTPUT would otherwise apply the penalty.
    CHECK_NEAR(MoraleOutputMultiplier(0, t), 1.0, 0.0);
    // A modifier that is not strictly positive is ignored (an unloaded tuning table has
    // every field at zero, and a zero multiplier would silently wipe the term).
    TuningTable zero;
    zero.MORALE_INCREASE_OUTPUT = 60;
    CHECK_NEAR(MoraleOutputMultiplier(80, zero), 1.0, 0.0);

    OutputModifiers m;
    m.baseOutput = 1000;
    CHECK_NEAR(TotalSystemOutput(m, t), 1000.0, 0.0);
    m.addictionPhase3 = true;
    CHECK_NEAR(TotalSystemOutput(m, t), 500.0, 0.0);
    m.addictionPhase3 = false;
    m.playerOutMod = 0.5;
    m.systemOutMod = 0.5;
    CHECK_NEAR(TotalSystemOutput(m, t), 250.0, 0.0);
    // B4: the engine's round is ties-to-EVEN, so 302.5 goes DOWN to 302 (it used to be 303)
    m.baseOutput = 1210;
    CHECK_NEAR(TotalSystemOutput(m, t), 302.0, 0.0);
    m.baseOutput = 0;
    CHECK_NEAR(TotalSystemOutput(m, t), 0.0, 0.0);
    // No owner and a rebelling system both return zero before any multiplier runs.
    m.baseOutput = 1000;
    m.playerOutMod = 1.0;
    m.systemOutMod = 1.0;
    m.owned = false;
    CHECK_NEAR(TotalSystemOutput(m, t), 0.0, 0.0);
    m.owned = true;
    m.rebelling = true;
    CHECK_NEAR(TotalSystemOutput(m, t), 0.0, 0.0);
    m.rebelling = false;

    // ---- the population -> output term (lane N) ----------------------------------------
    // Output per head is typeOutputMod x 1.8 / 500000; the imperial row's modifier is 1.
    GroupOutputInputs g;
    g.group = PopGroup::Imperial;
    g.count = 2000000000LL;
    CHECK_NEAR(GroupOutput(g, t), 7200.0, 1e-9);             // 2e9 / 5e5 x 1.8
    g.stations = 2;                                          // 1 + 2 x 0.1
    CHECK_NEAR(GroupOutput(g, t), 8640.0, 1e-9);
    g.stations = 0;
    g.count = 0;
    CHECK_NEAR(GroupOutput(g, t), 0.0, 0.0);
    g.count = -5;
    CHECK_NEAR(GroupOutput(g, t), 0.0, 0.0);
    // The station bonus is imperial-only, and civilians carry the morale multiplier.
    g.group = PopGroup::Civilian;
    g.count = 500000000LL;
    g.stations = 4;
    g.morale = 0;
    CHECK_NEAR(GroupOutput(g, t), 500000000.0 / 500000.0 * F32(0.33) * 1.8, 1e-9);
    g.morale = 80;                                           // above MORALE_INCREASE_OUTPUT
    CHECK_NEAR(GroupOutput(g, t), 500000000.0 / 500000.0 * F32(0.33) * 1.8 * 1.1, 1e-9);
    g.independent = true;                                    // an independent colony has no morale
    CHECK_NEAR(GroupOutput(g, t), 500000000.0 / 500000.0 * F32(0.33) * 1.8, 1e-9);

    // The whole base-output sum on the reference save's human homeworld, with the two
    // data-file species fields left at zero so only the terms the executable carries move.
    BaseOutputInputs b;
    b.imperialPopulation = 2000000000LL;   // Pop 1e9 + pbon 1e9
    b.civilianPopulation = 500000000LL;
    b.civilianMorale = 75;
    b.transitResources = 0;
    b.resourcesAvailable = 5000;
    b.infra = 1.0f;
    b.infraBonus = 1.0f;
    b.overHarvestRate = 0.0;
    // cbrt(2e9/100) x 0.01 = 2.71 -> clamps to 1, and a clamped value at or above 1 - 1e-4
    // is *substituted* by the infrastructure term rather than capping it. Infra + ibon = 2
    // here, so the fraction is 2, not 1 -- the branch is a substitution, not a min, and a
    // pending infrastructure bonus can push a colony's extraction above unity.
    CHECK_NEAR(StripMineFraction({b.imperialPopulation, b.infra, b.infraBonus}), 2.0f, 0.0);
    // ... and a colony whose population term has not saturated is capped by it as usual.
    CHECK_NEAR(StripMineFraction({1000000, 0.4f, 0.0f}), 0.21544346f, 1e-6f);
    CHECK_NEAR(StripMineFraction({100000000LL, 0.4f, 0.0f}), 0.4f, 0.0);
    const double expected = 7200.0 + 500000000.0 / 500000.0 * F32(0.33) * 1.8 * 1.1 + 9000.0;
    CHECK_NEAR(SystemBaseOutput(b, t), expected, 1e-6);
    // Linear in population: a tenth of the imperial pop is a tenth of that term.
    b.imperialPopulation = 200000000LL;
    b.civilianPopulation = 0;
    b.resourcesAvailable = 0;
    CHECK_NEAR(SystemBaseOutput(b, t), 720.0, 1e-9);
    // With SRoh = 0 the over-harvest demand degenerates to min(available, speciesBaseDemand).
    OverHarvestInputs oh;
    oh.resourcesAvailable = 5000;
    oh.population = 2000000000LL;
    oh.speciesBaseDemand = 120;
    CHECK_NEAR(OverHarvestDemand(oh), 120.0, 0.0);
    oh.speciesBaseDemand = 9000;
    CHECK_NEAR(OverHarvestDemand(oh), 5000.0, 0.0);
    // ... and with SRoh > 0 it adds rate x available x clamp01(pop x 1e-5), floored at 1.
    oh.speciesBaseDemand = 0;
    oh.overHarvestRate = 0.5;
    CHECK_NEAR(OverHarvestDemand(oh), 2500.0, 1e-9);         // clamp01(2e9 x 1e-5) = 1
    oh.population = 10000;                                   // clamp01(0.1)
    CHECK_NEAR(OverHarvestDemand(oh), 250.0, 1e-9);
    oh.population = 0;                                       // the floor of 1, not 0
    CHECK_NEAR(OverHarvestDemand(oh), 1.0, 0.0);

    // The population-type table the executable builds in code.
    CHECK_NEAR(PopTypeOf(PopGroup::Imperial, t).outputMod, 1.0, 0.0);
    CHECK_NEAR(PopTypeOf(PopGroup::Civilian, t).outputMod, F32(0.33), 0.0);
    CHECK_EQ(static_cast<int>(PopTypeOf(PopGroup::Imperial, t).maxPopulation), 50000000);
    CHECK_EQ(static_cast<int>(PopTypeOf(PopGroup::Civilian, t).maxPopulation), 20000000);
    CHECK_NEAR(RoundHalfEven(0.5), 0.0, 0.0);
    CHECK_NEAR(RoundHalfEven(1.5), 2.0, 0.0);
    CHECK_NEAR(RoundHalfEven(-2.5), -2.0, 0.0);

    OutputSplit s = SplitOutput(1000, {0.5, 0.25, 0.125, 0.125});
    CHECK_NEAR(s.trade, 500.0, 0.0);
    CHECK_NEAR(s.construction, 250.0, 0.0);
    CHECK_NEAR(s.terraform, 125.0, 0.0);
    CHECK_NEAR(s.infra, 125.0, 0.0);

    CHECK_EQ(ConstructionPoints(250, 2, t), 375);
    CHECK_EQ(ConstructionPoints(250, 0, t), 250);
    CHECK_EQ(ConstructionPoints(3, 1, t), 3);                // 3.75 TRUNCATES, it does not round

    OutputSplit l = SplitLeftover(100, {0.5, 0.25, 0.125, 0.125}, false, false);
    CHECK_NEAR(l.trade, 67.0, 0.0);
    CHECK_NEAR(l.terraform, 17.0, 0.0);
    CHECK_NEAR(l.infra, 17.0, 0.0);                          // and 67+17+17 != 100
    l = SplitLeftover(100, {0, 1, 0, 0}, true, false);        // construction rate exactly 1
    CHECK_NEAR(l.trade, 50.0, 0.0);
    CHECK_NEAR(l.terraform, 0.0, 0.0);
    CHECK_NEAR(l.infra, 50.0, 0.0);
    l = SplitLeftover(100, {0, 1, 0, 0}, true, true);
    CHECK_NEAR(l.trade, 100.0, 0.0);
    l = SplitLeftover(0, {0.5, 0.25, 0.125, 0.125}, false, false);
    CHECK_NEAR(l.trade, 0.0, 0.0);
}

static void test_system_money() {
    CHECK_NEAR(SuitabilityCostMod(0.3, 0.5, 0.15, false, true), 0.15, 0.0);   // capped by SuitTol
    CHECK_NEAR(SuitabilityCostMod(0.45, 0.5, 0.15, false, true), 0.05, 1e-12);
    CHECK_NEAR(SuitabilityCostMod(0.3, 0.5, 0.15, true, true), 0.0, 0.0);     // rebel AI pays nothing
    CHECK_NEAR(SuitabilityCostMod(0.5, 0.5, 0.15, false, false), 20.0, 0.0);  // unowned

    SystemMoneyInputs m;
    m.tradePoints = 20;                                   // 4 blocks x 5 = 100
    CHECK_EQ(SystemMoneyIncome(m), 100);
    m.tradePoints = 4.9;                                  // no whole block
    CHECK_EQ(SystemMoneyIncome(m), 0);
    m.tradePoints = 123;                                  // 24 blocks -> 600
    CHECK_EQ(SystemMoneyIncome(m), 600);

    m.popIncomeImperial = 100;
    m.popIncomeCivilian = 50;
    m.slaveIncome = 10;
    m.speciesIncomeFactor = 1.1;                          // Zuul
    m.speciesCostFactor = 0.7;
    m.suitCostMod = 0.2;
    // (600 + 160) x 1.1 = 836.0000000000001; cost 0.7 x 0.2 x 15000 = 2100;
    // -1263.9999999999998 truncates toward zero
    CHECK_EQ(SystemMoneyIncome(m), -1263);
    m.speciesIncomeFactor = 1.0;                          // 760 - 0.25 x 15000 = -2990 exactly
    m.speciesCostFactor = 1.0;
    m.suitCostMod = 0.25;
    CHECK_EQ(SystemMoneyIncome(m), -2990);

    SystemMoneyInputs n;
    n.tradePoints = 10;
    n.speciesIncomeFactor = 0.8;                          // Morrigi
    CHECK_EQ(SystemMoneyIncome(n), 40);

    n.speciesIncomeFactor = 1.0;
    n.tradePoints = 100;                                  // 500
    n.playerIncMod = 1.2;                                 // 600
    n.serverIncomeMod = 0.5;
    n.difficultyIncomeMult = 2.0;                         // x1 net
    CHECK_EQ(SystemMoneyIncome(n), 600);

    SystemMoneyInputs unowned;
    unowned.suitCostMod = 20.0;                           // 20 x 15000
    CHECK_EQ(SystemMoneyIncome(unowned), -300000);

    // the cost is not scaled by the income multipliers
    SystemMoneyInputs c;
    c.tradePoints = 10;                                   // 50
    c.playerIncMod = 3.0;                                 // 150
    c.suitCostMod = 0.01;                                 // cost 150
    CHECK_EQ(SystemMoneyIncome(c), 0);

    CHECK_NEAR(ConstantsOf(Species::Zuul).incomeFactor, 1.1, 0.0);
    CHECK_NEAR(ConstantsOf(Species::Zuul).hazardCostFactor, 0.7, 0.0);
    CHECK_NEAR(ConstantsOf(Species::Morrigi).incomeFactor, 0.8, 0.0);
    CHECK_NEAR(ConstantsOf(Species::Human).incomeFactor, 1.0, 0.0);
    CHECK(!ConstantsOf(Species::Zuul).systemBonusEligible);
    CHECK(ConstantsOf(Species::Hiver).systemBonusEligible);

    // A von Neumann machine at the system zeroes the cost before anything else is looked at.
    CHECK_NEAR(SuitabilityCostMod(0.3, 0.5, 0.15, false, true, true), 0.0, 0.0);
    CHECK_NEAR(SuitabilityCostMod(0.5, 0.5, 0.15, false, false, true), 0.0, 0.0);
}

static void test_population_income() {
    TuningTable t = tuning();

    // Income per head is typeIncomeMod / 14000 -- no 1.8, no 500000. An imperial billion
    // is 1e9/14000 = 71428.57..., truncated.
    CHECK_EQ(GroupIncome(PopGroup::Imperial, 1000000000, t), 71428);
    CHECK_EQ(GroupIncome(PopGroup::Imperial, 13999, t), 0);
    CHECK_EQ(GroupIncome(PopGroup::Imperial, 14000, t), 1);
    // The civilian row's modifier is the float32 0.33, so half a billion civilians give
    // ftol(0.33000001311302185 x 35714.2857...) = 11785.
    CHECK_EQ(GroupIncome(PopGroup::Civilian, 500000000, t), 11785);
    CHECK_EQ(GroupIncome(PopGroup::Slaves, 14000, t), 3);   // SLAVES_INCOME_MOD = 3

    PopIncomeRow rows[kSpeciesCount] = {};
    rows[0].count = 1000000000;
    CHECK_NEAR(PopulationIncome(PopGroup::Imperial, rows, true, false, t), 71428.0, 0.0);

    // Two species truncate SEPARATELY, so the sum is not the truncation of the sum.
    PopIncomeRow two[kSpeciesCount] = {};
    two[0].count = 20999;   // -> 1
    two[2].count = 20999;   // -> 1
    CHECK_NEAR(PopulationIncome(PopGroup::Imperial, two, true, false, t), 2.0, 0.0);
    CHECK_EQ(GroupIncome(PopGroup::Imperial, 41998, t), 2);   // ... which happens to agree here
    PopIncomeRow three[kSpeciesCount] = {};
    three[0].count = 13999;  // -> 0
    three[2].count = 13999;  // -> 0
    CHECK_NEAR(PopulationIncome(PopGroup::Imperial, three, true, false, t), 0.0, 0.0);
    CHECK_EQ(GroupIncome(PopGroup::Imperial, 27998, t), 1);   // ... and here it does NOT

    // Morale applies to the civilian row only, and the product truncates again.
    PopIncomeRow mor[kSpeciesCount] = {};
    mor[0].count = 500000000;
    mor[0].morale = 80;                                        // >= MORALE_INCREASE_OUTPUT (75)
    CHECK_NEAR(PopulationIncome(PopGroup::Civilian, mor, true, false, t),
               std::floor(11785.0 * 1.1), 0.0);
    // ... but not to the imperial row.
    CHECK_NEAR(PopulationIncome(PopGroup::Imperial, mor, true, false, t),
               static_cast<double>(GroupIncome(PopGroup::Imperial, 500000000, t)), 0.0);
    // ... and an independent colony bypasses morale entirely.
    CHECK_NEAR(PopulationIncome(PopGroup::Civilian, mor, true, true, t), 11785.0, 0.0);

    // Addiction multiplies every row, imperial included.
    PopIncomeRow add[kSpeciesCount] = {};
    add[0].count = 1000000000;
    add[0].addicted = true;
    CHECK_NEAR(PopulationIncome(PopGroup::Imperial, add, true, false, t),
               std::floor(71428.0 * 0.9), 0.0);
}

static void test_max_income() {
    // The whole chain, with the numbers a level-1 AI Zuul colony produces: the trade points
    // are the rounded output total, blocks of five are worth five each after the x5, the
    // species factor is 1.1 and the difficulty income modifier another 1.1.
    SystemMoneyInputs m;
    m.popIncomeImperial = 71428;
    m.speciesIncomeFactor = ConstantsOf(Species::Zuul).incomeFactor;
    m.speciesCostFactor = ConstantsOf(Species::Zuul).hazardCostFactor;
    m.suitCostMod = 0.0;                       // a homeworld sits exactly at its ideal
    m.difficultyIncomeMult = DifficultyModsFor(1, /*isAI=*/true, /*isNpc=*/false).incomeMult;
    const int ai = SystemMaxIncome(12345.0, m);
    m.difficultyIncomeMult = DifficultyModsFor(1, /*isAI=*/false, /*isNpc=*/false).incomeMult;
    const int human = SystemMaxIncome(12345.0, m);
    // The AI's advantage on this row is exactly the 1.1 in the difficulty table.
    CHECK(ai > human);
    CHECK_NEAR(static_cast<double>(ai) / static_cast<double>(human), 1.1, 1e-5);

    // A colony whose money comes out negative contributes ZERO to the empire total rather
    // than reducing it -- the `jg` at the end of ComputeMaxIncome.
    SystemMoneyInputs bad;
    bad.suitCostMod = 20.0;                    // 20 x 15000 of cost against no income
    CHECK_EQ(SystemMoneyIncome(bad), -300000);
    CHECK_EQ(SystemMaxIncome(0.0, bad), 0);

    // The rate vector is trade = 1, so the trade points are the half-to-even rounded total.
    SystemMoneyInputs r;
    CHECK_EQ(SystemMaxIncome(20.5, r), SystemMoneyIncome([] {
                 SystemMoneyInputs x;
                 x.tradePoints = 20.0;         // 20.5 ties to the even neighbour
                 return x;
             }()));
    CHECK_EQ(SystemMaxIncome(21.5, r), SystemMoneyIncome([] {
                 SystemMoneyInputs x;
                 x.tradePoints = 22.0;
                 return x;
             }()));
}

static void test_difficulty_table() {
    // Level 0 gives the break to the human; levels 1 and 2 give it to the AI.
    const DifficultyMods e_ai = DifficultyModsFor(0, true, false);
    const DifficultyMods e_pl = DifficultyModsFor(0, false, false);
    CHECK_NEAR(e_ai.maintenanceDivisor, 1.0, 0.0);
    CHECK_NEAR(e_ai.incomeMult, 1.0, 0.0);
    CHECK_NEAR(e_pl.maintenanceDivisor, 1.5, 0.0);
    CHECK_NEAR(e_pl.incomeMult, 1.5, 0.0);

    const DifficultyMods n_ai = DifficultyModsFor(1, true, false);
    CHECK_NEAR(n_ai.maintenanceDivisor, 3.0, 0.0);
    CHECK_NEAR(n_ai.incomeMult, 1.1, 1e-7);
    CHECK_NEAR(n_ai.researchMult, 1.5, 0.0);
    CHECK_NEAR(DifficultyModsFor(1, false, false).incomeMult, 1.0, 0.0);

    const DifficultyMods h_ai = DifficultyModsFor(2, true, false);
    CHECK_NEAR(h_ai.maintenanceDivisor, 1000000.0, 0.0);
    CHECK_NEAR(h_ai.incomeMult, 1.7, 1e-7);
    CHECK_NEAR(h_ai.researchMult, 2.0, 0.0);
    // Hard maintenance really is "divided by a million", i.e. free.
    CHECK_EQ(MaintenanceCost(999999, h_ai.maintenanceDivisor), 0);

    // An NPC player takes the non-AI triple whatever its AI flag says ...
    CHECK_NEAR(DifficultyModsFor(1, true, true).incomeMult, 1.0, 0.0);
    // ... and an out-of-range level falls back to all ones rather than failing.
    CHECK_NEAR(DifficultyModsFor(-1, true, false).incomeMult, 1.0, 0.0);
    CHECK_NEAR(DifficultyModsFor(3, true, false).incomeMult, 1.0, 0.0);
    CHECK_NEAR(DifficultyModsFor(3, true, false).maintenanceDivisor, 1.0, 0.0);
}

static void test_bonuses() {
    TuningTable t = tuning();
    std::int64_t pop = 900, bonus = 500;
    ApplyPopulationBonus(pop, 1000, bonus);
    CHECK_EQ(pop, std::int64_t{1000});
    CHECK_EQ(bonus, std::int64_t{400});
    pop = 1200;                                    // over cap: nothing applied
    ApplyPopulationBonus(pop, 1000, bonus);
    CHECK_EQ(pop, std::int64_t{1200});
    CHECK_EQ(bonus, std::int64_t{400});

    double infra = 0.95, ibon = 0.1;
    BonusApplyResult br = ApplyInfrastructureBonus(infra, ibon);
    CHECK_NEAR(infra, 1.0, 0.0);      // the pool covered the remainder: exactly 1, not 0.999...
    CHECK_NEAR(ibon, 0.05, 1e-7);
    CHECK(br.applied);
    CHECK(!br.resetTurnsDeveloping);  // a home system is not reset
    br = ApplyInfrastructureBonus(infra, ibon, /*homeSystem=*/false);
    CHECK(!br.applied);               // already at 1: nothing happens, no reset either
    infra = 0.5;
    ibon = 0.1;
    br = ApplyInfrastructureBonus(infra, ibon, /*homeSystem=*/false);
    CHECK_NEAR(infra, 0.6, 1e-7);
    CHECK_NEAR(ibon, 0.0, 1e-7);
    CHECK(br.resetTurnsDeveloping);   // a non-home colony absorbing a bonus resets ntdev

    std::int64_t up = 10, ub = 500;
    ApplyPopulationBonus(up, 1000, ub, /*owned=*/false);
    CHECK_EQ(up, std::int64_t{10});   // an unowned system drops the whole pool
    CHECK_EQ(ub, std::int64_t{0});

    SystemBonusInputs in;
    in.stable = true;
    in.turnsOwned = 11;
    in.turnsDeveloping = 11;
    in.capacity = 1000000;
    std::int64_t pbon = 0;
    double ibonus = 0;
    AccrueSystemBonus(in, pbon, ibonus, t);
    CHECK_EQ(pbon, std::int64_t{10000});           // 1e6 x 0.01
    CHECK_NEAR(ibonus, 0.05, 1e-12);
    for (int i = 0; i < 20; ++i) AccrueSystemBonus(in, pbon, ibonus, t);
    CHECK_EQ(pbon, std::int64_t{100000});          // capped at 1e6 x 0.1
    CHECK_NEAR(ibonus, 0.2, 1e-12);                // capped at INFRABONUS

    // the increment truncates: 12345 x 0.005 = 61.725 -> 61; target 1234.5 -> 1234
    TuningTable small = t;
    small.SYSTEMBONUS_POPBONUS_INC = 0.005;
    in.capacity = 12345;
    std::int64_t p3 = 0;
    double i3 = 0;
    AccrueSystemBonus(in, p3, i3, small);
    CHECK_EQ(p3, std::int64_t{61});
    for (int i = 0; i < 30; ++i) AccrueSystemBonus(in, p3, i3, small);
    CHECK_EQ(p3, std::int64_t{1234});
    p3 = 5000;                                     // already above the target: untouched
    AccrueSystemBonus(in, p3, i3, small);
    CHECK_EQ(p3, std::int64_t{5000});
    in.capacity = 1000000;

    // Zuul never accrue either bonus
    std::int64_t pz = 0;
    double iz = 0;
    in.ownerSpeciesEligible = false;
    AccrueSystemBonus(in, pz, iz, t);
    CHECK_EQ(pz, std::int64_t{0});
    CHECK_NEAR(iz, 0.0, 0.0);
    in.ownerSpeciesEligible = true;

    std::int64_t p2 = 0;
    double i2 = 0;
    in.turnsOwned = 10;                            // not strictly more than MINTURNS
    AccrueSystemBonus(in, p2, i2, t);
    CHECK_EQ(p2, std::int64_t{0});
    in.turnsOwned = 11;
    in.stable = false;
    AccrueSystemBonus(in, p2, i2, t);
    CHECK_EQ(p2, std::int64_t{0});
    in.stable = true;
    in.turnsDeveloping = 10;
    AccrueSystemBonus(in, p2, i2, t);
    CHECK_EQ(p2, std::int64_t{0});
}

static void test_build_queue() {
    std::vector<BuildOrder> q = {{1, 11, 100, 100, 50}, {2, 12, 200, 200, 0}, {3, 13, 300, 300, 70}};
    BuildQueueResult r = ProcessBuildQueue(q, 250);
    CHECK_EQ(r.completedOrderIds.size(), std::size_t{1});
    CHECK_EQ(r.completedOrderIds[0], 11);
    CHECK_EQ(r.moneyCharged, 50);
    CHECK_EQ(r.pointsLeft, 0);
    CHECK_EQ(q.size(), std::size_t{2});
    CHECK_EQ(q[0].orderId, 12);
    CHECK_EQ(q[0].constructionLeft, 50);
    CHECK_EQ(q[1].constructionLeft, 300);

    r = ProcessBuildQueue(q, 700);
    CHECK_EQ(r.completedOrderIds.size(), std::size_t{2});
    CHECK_EQ(r.moneyCharged, 70);
    CHECK_EQ(r.pointsLeft, 350);
    CHECK(q.empty());

    r = ProcessBuildQueue(q, 100);                 // empty queue: points pass through
    CHECK_EQ(r.pointsLeft, 100);

    std::vector<BuildOrder> exact = {{1, 21, 100, 100, 0}};
    r = ProcessBuildQueue(exact, 100);             // exactly enough completes
    CHECK_EQ(r.completedOrderIds.size(), std::size_t{1});
    CHECK(exact.empty());

    std::vector<BuildOrder> zero = {{1, 31, 100, 100, 0}};
    r = ProcessBuildQueue(zero, 0);
    CHECK(r.completedOrderIds.empty());
    CHECK_EQ(zero[0].constructionLeft, 100);
}

int main() {
    test_capacity();
    test_growth();
    test_infra_terraform();
    test_slaves();
    test_output();
    test_system_money();
    test_population_income();
    test_max_income();
    test_difficulty_table();
    test_bonuses();
    test_build_queue();
    return simtest::finish("test_colony");
}