// Colonies: carrying capacity, population growth, infrastructure and terraforming
// point conversion, slave deaths, output split, build-queue consumption.
#pragma once

#include <cstddef>
#include <cstdint>
#include <vector>

#include "game/sim/species.h"
#include "game/sim/tuning.h"

namespace sots::sim {

// Population group kinds (rows of the per-type population table).
enum class PopGroup : int { Imperial = 0, Civilian = 1, Slaves = 2 };

// Xenotech level a player has reached against one species: bit k of the per-species
// flag word is set when the k-th xenotech for that species is researched. The bit order
// is the order of the xenotech families (translation 1/2/3, incorporate, addict,
// temperance, subjugate, accommodate, proliferate). CONFIDENCE: high on the order.
struct SpeciesTechFlags {
    bool translation1 = false;   // bit 0: slave death rate x0.8
    bool translation2 = false;   // bit 1: slave death rate -0.2
    bool translation3 = false;   // bit 2: slave death rate -0.2
    bool incorporate = false;    // bit 3
    bool addict = false;         // bit 4
    bool temperance = false;     // bit 5: addiction to this species is cured
    bool subjugate = false;      // bit 6
    bool accommodate = false;    // bit 7: no suitability penalty on capacity (hazard = 1)
    bool proliferate = false;    // bit 8

    static SpeciesTechFlags FromBits(unsigned bits) {
        SpeciesTechFlags f;
        f.translation1 = bits & 0x001;
        f.translation2 = bits & 0x002;
        f.translation3 = bits & 0x004;
        f.incorporate  = bits & 0x008;
        f.addict       = bits & 0x010;
        f.temperance   = bits & 0x020;
        f.subjugate    = bits & 0x040;
        f.accommodate  = bits & 0x080;
        f.proliferate  = bits & 0x100;
        return f;
    }
};

// ---------------------------------------------------------------------------------------
// Capacity
// ---------------------------------------------------------------------------------------

// Hazard modifier on carrying capacity from planet suitability vs the species' ideal:
//   clamp01(1 - |suit - ideal| / (SuitTol + 0.1))
// Linear, no exponent; the +0.1 gives every species a habitable band even at zero
// tolerance. SuitTol starts at the species value and is raised by the atmospheric and
// gravitational adaptation techs. Skipped (treated as 1) when the player has the
// accommodate xenotech for the species or is the rebel AI -- the caller decides that.
// CONFIDENCE: high -- read with its constant.
double HazardModifier(double suitability, double idealSuitability, double suitTolerance);

struct CapacityInputs {
    int planetSize = 0;                  // Size
    Species species = Species::Human;    // species of the population group
    bool speciesCanLive = true;          // species can survive on this planet class
    PopGroup group = PopGroup::Imperial;
    double groupCapacityMult = 1.0;      // per-group-type capacity column
    double speciesGrowthFactor = 1.0;    // per-species factor
    bool ownerIsDifferentSpecies = false;
    double crossSpeciesMod = 1.0;        // applied when the owner is another species
    double hazardMod = 1.0;              // from HazardModifier (or 1 when accommodated)
    bool arcologyTech = false;           // adds a flat 1e8 (imperial) / 2e8 (civilian)
    bool groupMaxEnabled = false;        // per-group hard cap present
    std::int64_t groupMax = 0;
    bool ownerIsNpc = false;             // owner species is the independent race
};

// Round a capacity down to a whole ten -- but only above ten. B4 finding: every carrying
// capacity goes through this, so a colony's cap is always a multiple of 10 in practice.
// CONFIDENCE: high -- read from the helper's truncate/divide/multiply sequence.
std::int64_t QuantiseCapacity(std::int64_t v);

// cap = quantise(ftoi64((Size x 1e8) x (hazard x (groupMult x speciesFactor x crossSpecies))))
//       + arcology flat bonus (1e8 imperial / 2e8 civilian, **0 for slaves**); clamped to
//       the group max; x INDSYS_IMPERIAL_POPULATION_MOD and re-quantised when the owner is
//       the NPC race. The NPC species itself, and a species that cannot live there, get 0.
// `Size x 1e8` is an exact 64-bit integer product; only the modifier chain is floating
// point. CONFIDENCE: high on the shape, the order and the special cases.
std::int64_t CarryingCapacity(const CapacityInputs& in, const TuningTable& t);

// ---------------------------------------------------------------------------------------
// Growth
// ---------------------------------------------------------------------------------------

struct GrowthInputs {
    std::int64_t pop = 0;
    bool blockaded = false;          // growth halted at this system for this group
    double suitability = 0;          // the planet's Suit
    double idealSuitability = 0;     // the growing species' ideal
    double suitTolerance = 0;        // the owner's SuitTol
    bool accommodated = false;       // the owner ignores suitability entirely (RebAI etc.)
    double playerPopMod = 1.0;       // PopMod, applied only when > 0
    double extraFactor = 1.0;        // the caller's per-call factor, applied only when > 0
    double groupGrowthMult = 0.0;    // per-group-type growth column; applied only if > 0
};

// The clamp the engine puts on the growth exponent before calling pow().
constexpr double kGrowthExponentMin = 0.009999999776482582;  // (double)0.01f
constexpr double kGrowthExponentMax = 1000.0;

// Suitability distance the growth curve is built from:
//   min( |ideal - clamp(suit, 0, 20)| , tolerance )
// clamped through float32 at every step, and 0 when the owner is accommodated. The 0..20
// clamp on the planet's suitability is easy to miss and matters at the extremes.
// CONFIDENCE: high.
double GrowthSuitabilityDistance(double suit, double ideal, double tolerance, bool accommodated);

// Growth fraction for one population group.
//
// B4 CORRECTION -- the curve is **not** logistic in population. There is no `pop / capacity`
// term anywhere in the chain; the capacity is never even passed in. The base of the power
// is a *suitability* term:
//   base = 1 - clamp01( GrowthSuitabilityDistance(...) / tolerance )
//   g    = clamp01( pow(base, clamp(POPULATION_GROWTH_EXP, 0.01f, 1000)) )
//   if g > 0: g x= POPULATION_GROWTH_MOD (if > 0) x PopMod (if > 0)
//                x extraFactor (if > 0) x groupGrowthMult (if > 0)
// Every guard is a strict `> 0` and every intermediate is stored back to a float32 slot.
// A zero tolerance makes the quotient 0/0, which propagates NaN through the whole chain --
// a real hazard, reproduced rather than papered over. CONFIDENCE: high.
double PopulationGrowthFraction(const GrowthInputs& in, const TuningTable& t);

// The per-turn delta: `trunc(pop x g)`, or exactly 1 when that truncates to 0 while g is
// strictly positive. A group with no population, or a blockaded one, yields 0. There is no
// 50,000,000 cap here -- that lives in the apply. CONFIDENCE: high.
std::int64_t PopulationGrowthDelta(const GrowthInputs& in, const TuningTable& t);

// The largest population step a colony can take in one turn, up or down.
constexpr std::int64_t kMaxPopulationStep = 50000000;

// Apply growth to an imperial population:
//   delta = min(delta, 50,000,000) when delta >= 0
//   new   = pop + delta
//   new <= cap                 -> new
//   new >  cap and pop <= cap  -> exactly cap        (no shrink at all)
//   new >  cap and pop >  cap  -> max(pop - min(|cap - pop|, 50,000,000), min(pop, 100))
//   finally, floored at 0
// B4 correction: the shrink is computed from the *old* population and only runs when the
// colony was already over the cap; a colony that merely grows past the cap simply lands on
// it. CONFIDENCE: high -- read branch by branch.
std::int64_t ApplyImperialGrowth(std::int64_t pop, std::int64_t capacity, std::int64_t delta);

// ---------------------------------------------------------------------------------------
// Infrastructure and terraforming
// ---------------------------------------------------------------------------------------

// The infrastructure divisor is a TRUE double 3.3e-05 in the image (not a widened float).
constexpr double kInfraPerPoint = 3.3e-5;

// The two terraforming factors are widened *float* literals, and they differ: the "points
// needed" side uses 1.8000000715255737 (= 1.5 x (double)1.2f, stored as one constant),
// the yield side multiplies by the exact double 1.5 and then by (double)1.2f. The two
// therefore agree exactly, which is what keeps need and yield consistent.
constexpr double kTerraform12 = 1.2000000476837158;      // (double)1.2f
constexpr double kTerraformNeedFactor = 1.8000000715255737;  // 1.5 x (double)1.2f

// Points needed to bring infrastructure from `infra` to 1.0: ceil((1 - infra) / 3.3e-5).
// A real ceil on the double quotient. CONFIDENCE: high.
double InfrastructurePointsNeeded(double infra);

// Infrastructure gained from spending `points`: `points / 500 x 0.01 x 1.65` evaluated in
// that order on the x87, clamped at zero and narrowed to float32 once on the store. The
// arithmetic is deliberately not folded into a single x3.3e-5. CONFIDENCE: high.
double InfrastructureGain(double points);

// Spend from a pool toward full infrastructure; returns the *delta* and reports the
// points the pool has left (which cascade into the terraforming pool). CONFIDENCE: high.
double ApplyInfrastructurePoints(double infra, double pool, double* pointsUnused);

// Add an infrastructure delta to a colony: a no-op once infrastructure is at 1, clamped
// to 1 otherwise, stored as float32. CONFIDENCE: high.
double ApplyInfrastructureDelta(double infra, double delta);

// Per-turn infrastructure loss of an unowned system. The image's constant is the widened
// float literal `(double)0.02f` = 0.019999999552965164, not the decimal 0.02 (B4).
constexpr double kUnownedInfraDecay = 0.019999999552965164;

// Unowned systems lose `kUnownedInfraDecay` infrastructure per turn: the difference is stored
// back into a 4-byte float and then floored, `result <= 0 -> 0`. CONFIDENCE: high -- read with
// its constant's exact bits and the fcomp against the .rdata zero.
double DecayUnownedInfrastructure(double infra);

// Points needed to terraform from `suit` to `ideal`:
//   ceil( |float32(ideal - suit)| / |TerraMod x 1.8000000715255737 / 20000| )
// B4 correction: the player's terraforming modifier is inside the *need*, not only the
// yield -- a better modifier needs proportionally fewer points, and that is what keeps
// need and yield consistent. CONFIDENCE: high.
double TerraformPointsNeeded(double suit, double ideal, double terraMod);

// Suitability change from spending `points`:
//   float32( points x 1.5 x (double)1.2f x TerraMod x sign / 20000 )
// with `sign` = -1 only when `suit > ideal` strictly (a planet already at the ideal takes
// +1). CONFIDENCE: high.
double TerraformDelta(double points, double terraMod, double suit, double ideal);

// Add a suitability delta, clamped so the planet stops at the ideal from whichever side
// it approached; stored as float32. CONFIDENCE: high.
double ApplyTerraformDelta(double suit, double delta, double ideal);

// ---------------------------------------------------------------------------------------
// Slaves
// ---------------------------------------------------------------------------------------

// The two xenotech factors of the slave-death modifier, as the image holds them.
constexpr double kSlaveModBase = 0.800000011920929;   // float32 0.8f, replaces the 1
constexpr double kSlaveModStep = 0.20000000298023224;  // (double)0.2f, subtracted

// Per-turn slave death rate:
//   mod  = (translation1 ? 0.8f : 1) - 0.2 x translation2 - 0.2 x translation3   [no clamp]
//   rate = ((|ideal - suit| x BYHAZARD + DEATH_RATE) + SRs x BYOUTPUT) x mod
// B4 correction: the hazard term is folded into the base *before* the output term, and every
// step is stored back to a float32; an unowned system reports 1.0, not 0. CONFIDENCE: high.
double SlaveDeathRate(double slaveOutputRate, double suit, double ideal,
                      const SpeciesTechFlags& flags, const TuningTable& t, bool owned = true);

// Deaths this turn:
//   d = ftoi64(slaves x (rate + plagueRate))
//   then MIN_DEATHS/MAX_DEATHS, each **disabled by any negative value**, then clamped into
//   [0, slaves].
// B4 correction: the worst plague at the system contributes an *additive* rate term -- that
// is the only path by which a plague kills slaves -- and both bounds are clamps, not just
// the maximum. CONFIDENCE: high.
std::int64_t SlaveDeaths(std::int64_t slaves, double rate, double plagueRate,
                         const TuningTable& t);

// ---------------------------------------------------------------------------------------
// Output split
// ---------------------------------------------------------------------------------------

struct OutputRates {
    double trade = 0;         // SRt
    double construction = 0;  // SRsc
    double terraform = 0;     // SRtf
    double infra = 0;         // SRi
};

// A slider at or below this counts as off, and it is also the value the all-zero fallback
// seeds the unpinned channels with. The image holds the widened float literal
// `(double)1e-4f`.
constexpr double kOutputRateThreshold = 9.999999747378752e-05;

// Normalise the player's output sliders.
//
// B4 CORRECTION -- this is not a symmetric "rescale all four to sum 1". The routine takes a
// *pinned* channel, and every caller in the game passes the trade slider (explicitly or by
// leaving the argument null, which selects it). So:
//   * terraform -> 0 when suitability is exactly at the ideal (an exact ==, no epsilon);
//     infra -> 0 when float32(infra + pending bonus) >= 1;
//   * any channel at or below (double)1e-4f -> 0;
//   * **only the trade slider** is clamped into [0, 1];
//   * the sum covers construction + infra + terraform only, accumulated in float32;
//   * an exactly-zero sum seeds those three with 1e-4f (still honouring the two
//     suppressions) and re-sums -- an equal split over three channels, not four;
//   * each of the three becomes `(rate / sum) x (1 - trade)`, float32 at both stores.
// Trade comes out unchanged. CONFIDENCE: high -- read step by step, including the pinned
// argument's default and the float32 accumulation.
OutputRates NormaliseOutputRates(const OutputRates& raw, bool suitAtIdeal, bool infraFull);

// ---------------------------------------------------------------------------------------
// Base output: the population, resource and over-harvest terms
//
// LANE N CORRECTION (2026-09-08). The previous model here treated the whole of a system's
// output as one multiplicative chain, `base x morale x stationFactor x ...`, with `baseOutput`
// an unresolved input. That is the wrong shape. The original sums **three independent terms**
// and applies morale and the station bonus to only some of them; the five player/system
// multipliers at the end are the only genuinely global factors.
//
// The population term itself is linear and is carried entirely by the executable: output
// points per head are `typeOutputMod x 1.8 / 500000`. See
// sots-re findings/subsystems/output-term.md.
// ---------------------------------------------------------------------------------------

// The per-population-type constants. The original builds this three-row table **in code**
// from x87 literals; only the slave row reads the data files. CONFIDENCE: high (read out of
// the initialiser, including its x87 register rotation).
struct PopTypeConstants {
    double outputMod = 0;                   // multiplies the per-capita output rate
    double incomeMod = 0;                   // multiplies the per-capita income rate
    std::int64_t maxPopulation = 0;         // the row's population cap
};

// Slave-row values come from the data files, so they are taken from the tuning table.
PopTypeConstants PopTypeOf(PopGroup g, const TuningTable& t);

// Output points contributed per head, before the type modifier: `1.8 / 500000`. Both are
// .rdata literals, i.e. facts about the algorithm rather than about the shipped data.
constexpr double kOutputPopulationFactor = 1.7999999999999998;  // the image's (double)1.8
constexpr double kOutputPopulationDivisor = 500000.0;

struct GroupOutputInputs {
    PopGroup group = PopGroup::Imperial;
    std::int64_t count = 0;      // heads in this group (imperial: Pop + pbon)
    int morale = 0;              // the system's Morale int[7] entry for this species
    int stations = 0;            // stations at the system
    bool owned = true;           // the system has an owner
    bool independent = false;    // sys.indi != null: morale is bypassed
};

// One (group, species) row's contribution:
//   count <= 0                        -> 0
//   q  = count / 500000
//   sf = 1 + stations x STATION_BONUS_IMPERIAL_OUTPUT, imperial groups of an owned system only
//        (and only while the constant is > 0)
//   mo = morale multiplier, civilian groups only
//   max(0, typeOutputMod x (sf x 1.8) x mo x q)
// CONFIDENCE: high -- read instruction by instruction, including the association of the
// multiplies, which is not free in 80-bit x87.
double GroupOutput(const GroupOutputInputs& in, const TuningTable& t);

// Morale effect on output. A morale entry of exactly 0 means "no record", and the original
// returns 1 rather than consulting the thresholds; each modifier is also ignored unless it
// is strictly positive, which matters because an unloaded tuning table has them at 0.
// CONFIDENCE: high (both guards read from the branch).
double MoraleOutputMultiplier(int morale, const TuningTable& t);

// Resource extraction efficiency: `min( clamp01(cbrt((Pop + pbon)/100) x 0.01), Infra + ibon )`,
// narrowed to float32 on the way in and on the way out. The `>= 1 - 1e-4` branch substitutes
// the infrastructure term outright. CONFIDENCE: high.
struct StripMineInputs {
    std::int64_t population = 0;   // Pop + pbon
    float infra = 0;               // Infra
    float infraBonus = 0;          // ibon
};
float StripMineFraction(const StripMineInputs& in);

// The over-harvest resource demand (0x007483b0). Also the term the resource ledger charges
// against the stock, and -- multiplied by the species' resource-output factor -- one of the
// three summands of a system's output.
//   B = overHarvestRate > 0 ? max(rate x resourcesAvailable x clamp01((Pop+pbon) x 1e-5), 1) : 0
//   return min(resourcesAvailable, max(speciesBaseDemand + B, 0))
// CONFIDENCE: high on the structure; the `max(.., 1)` floor is UNVERIFIED behaviourally
// because every save in the corpus carries SRoh = 0.
struct OverHarvestInputs {
    double overHarvestRate = 0;       // the normalised SRoh slider
    std::int64_t resourcesAvailable = 0;  // Res, plus MRes + ARes2 when the owner strip-mines
    std::int64_t population = 0;      // Pop + pbon
    int speciesBaseDemand = 0;        // SpeciesDef +0x4c, from the data files
    bool owned = true;
};
double OverHarvestDemand(const OverHarvestInputs& in);

struct BaseOutputInputs {
    // population
    std::int64_t imperialPopulation = 0;      // Pop + pbon, owner species only
    std::int64_t civilianPopulation = 0;      // Pop2 + pbon2 of the owner species, incl. surplus
    std::int64_t slavePopulation = 0;         // summed over species
    int imperialMorale = 0;                   // unused: imperial output ignores morale
    int civilianMorale = 0;
    int stations = 0;
    bool independent = false;
    // resources
    std::int64_t transitResources = 0;        // TRes
    std::int64_t resourcesAvailable = 0;      // Res (+ MRes + ARes2 when strip-mining)
    float infra = 0;
    float infraBonus = 0;
    // over-harvest
    double overHarvestRate = 0;               // SRoh
    int speciesBaseDemand = 0;                // SpeciesDef +0x4c
    double speciesResourceOutput = 0;         // SpeciesDef +0x50 (a float, widened)
};

// The sum of the three terms, before the global multipliers:
//   overHarvestDemand x speciesResourceOutput
// + (TRes + resourcesAvailable) x stripMineFraction x 0.9
// + imperial + civilian + slave population output
// CONFIDENCE: high on the terms; the summation order below reproduces the original's,
// which matters because x87 addition is not associative.
double SystemBaseOutput(const BaseOutputInputs& in, const TuningTable& t);

struct OutputModifiers {
    double baseOutput = 0;             // SystemBaseOutput
    bool addictionPhase3 = false;
    double scOutMod = 1.0;             // ScOutMod    (player +0x12c)
    double rebOutMod = 1.0;            // RebOutMod   (player +0x128)
    double techOutMod = 1.0;           // the game-setup output multiplier (player +0x224)
    double systemOutMod = 1.0;         // sys.OutMod  (+0x7c)
    double playerOutMod = 1.0;         // OutMod      (player +0x124)
    bool owned = true;                 // no owner -> 0
    bool rebelling = false;            // rbfl != 0 -> 0
};

// total = roundHalfEven( addiction x base x OutMod x sys.OutMod x techOut x RebOutMod x ScOutMod )
// The five multipliers are applied in exactly that order in one uninterrupted 80-bit chain.
// B4 correction retained: the engine's "round" is `fistp`/`fild` -- round to nearest, ties to
// EVEN -- and the result stays a double that the channel splits multiply; only the reported
// total slot truncates it.
// CONFIDENCE: high. Note the return of the original is the *unrounded* double; the rounding
// is done by its caller before the channel split, and is folded in here because every caller
// in the game does it.
double TotalSystemOutput(const OutputModifiers& m, const TuningTable& t);

// The unrounded value the original returns, for a hook that compares its return bit for bit.
double TotalSystemOutputRaw(const OutputModifiers& m, const TuningTable& t);

struct OutputSplit {
    double trade = 0;
    double construction = 0;
    double terraform = 0;
    double infra = 0;
};

// roundHalfEven(total x rate) for each channel, each kept as a double. CONFIDENCE: high.
OutputSplit SplitOutput(double total, const OutputRates& rates);

// Construction points after the shipyard station bonus. **Truncating**, not rounding.
// CONFIDENCE: high.
int ConstructionPoints(double constructionShare, int stations, const TuningTable& t);

// Redistribute unspent construction points L over trade/terraform/infra. Weights are the
// normalised rates of those three channels, or 1 / (suit != ideal) / (infra != 1) when the
// construction rate is exactly 1. Each share is rounded half-to-even independently, so the
// three need not add back up to L. CONFIDENCE: medium on the weights, high on the rounding.
OutputSplit SplitLeftover(double leftover, const OutputRates& rates, bool suitAtIdeal,
                          bool infraFull);

// ---------------------------------------------------------------------------------------
// Money: the population income law and a system's money output
// ---------------------------------------------------------------------------------------
//
// The income chain is a SECOND chain off the same three-row population table, and it is not
// the output chain with a different constant: income per head is `typeIncomeMod / 14000`,
// with neither the 1.8 factor nor the 500000 divisor the output law uses, and it truncates
// TWICE per (group, species) row -- once inside the per-row term and once after the morale
// and addiction factors. Summing the species first and truncating once is wrong on any
// colony carrying more than one species or a non-unit morale/addiction factor.
// See sots-re findings/subsystems/income-term.md.

// Money contributed per head, before the type modifier. An .rdata literal.
constexpr double kIncomePopulationDivisor = 14000.0;

// One (group, count) row's money: `ftol( typeIncomeMod x (count / 14000) )`. The type
// modifier is a float32 in the table, so it is narrowed before the multiply.
// CONFIDENCE: high -- read instruction by instruction (the whole function is five
// instructions and a tail-jump into the float-to-int helper).
int GroupIncome(PopGroup group, std::int64_t count, const TuningTable& t);

struct PopIncomeRow {
    std::int64_t count = 0;   // heads of this species in this group
    int morale = 0;           // the system's Morale entry for this species
    bool addicted = false;    // the system's addiction entry for this species is non-zero
};

// The population income of one group type, summed over the seven species rows:
//   for each species with count > 0:
//       mo  = civilian groups only: the morale multiplier (the SAME helper output uses)
//       ad  = addicted ? ADDICTION_INCOME_MOD : 1
//       sum += (double) ftol( (double)GroupIncome(group, count) x mo x ad )
// Slave groups take neither morale nor the civilian capacity surplus.
// CONFIDENCE: high on the loop and both truncations; the civilian capacity-surplus term
// (which the original adds to the owner species' civilian count when the colony is at its
// cap) is the CALLER's business -- fold it into that row's `count`. No colony in the
// corpus is at its cap, so that term is UNEXERCISED either way.
double PopulationIncome(PopGroup group, const PopIncomeRow (&rows)[kSpeciesCount],
                        bool owned, bool independent, const TuningTable& t);

// Suitability term of the money cost: how far the planet is from the owner species'
// ideal, capped by the owner's SuitTol (so the tolerance techs also cap the cost).
// The rebel AI pays nothing, and so does a system carrying a von Neumann machine
// (`vnh`, the very first test in the original); an unowned system is charged as if 20 away.
//
// `idealSuitability` is the SERVER's per-species baseline (`server->IdealSuit[species]`),
// not the owner's own `IdealSuit` field, and the species is the system's population species
// (`indi->indsp` on an independent colony) rather than the owner's. In the 11-save corpus
// the two IdealSuit sources agree for every player, so **which one the original reads is
// not falsifiable here** -- it is read from the instruction stream only.
// CONFIDENCE: high.
double SuitabilityCostMod(double suitability, double idealSuitability, double suitTolerance,
                          bool rebelAI, bool owned, bool vonNeumann = false);

struct SystemMoneyInputs {
    double tradePoints = 0;              // the system's trade-channel output this turn
    double popIncomeImperial = 0;        // PopulationIncome(Imperial, ...)
    double popIncomeCivilian = 0;        // PopulationIncome(Civilian, ...), incl. surplus
    double slaveIncome = 0;              // PopulationIncome(Slaves, ...)
    double speciesIncomeFactor = 1.0;    // ConstantsOf(owner).incomeFactor (1 if unowned)
    double playerIncMod = 1.0;           // IncMod (1 if unowned)
    double serverIncomeMod = 1.0;        // game-option income modifier (Sim tag `IncMod`)
    double difficultyIncomeMult = 1.0;   // DifficultyIncomeMod's table entry
    double speciesCostFactor = 1.0;      // ConstantsOf(owner).hazardCostFactor
    double suitCostMod = 0;              // from SuitabilityCostMod
};

// Money a system contributes this turn (`TradePointsToMoney`, returning the DOUBLE the
// original returns -- its caller is what truncates):
//   t = slaves + (civilian + (((trade - fmod(trade, 5)) x 5 + 0) + imperial))
//   t = f32(speciesIncomeFactor) x t
//   t = f32(f32(difficultyIncomeMult) x f32(serverIncomeMod)) x (f32(IncMod) x t)
//   cost = f32(speciesCostFactor) x (suitCostMod x 10000 x 1.5)
//   return t - cost
// Whole five-point blocks of trade are worth FIVE money each after being multiplied by
// five -- i.e. 25 money per block of five points; the same literal is the modulus and the
// multiplier. Every per-player multiplier is narrowed to float32 before it is applied, and
// the additive association above is the original's, which x87 does not let us reorder.
// CONFIDENCE: high -- read instruction by instruction.
double SystemMoneyIncomeRaw(const SystemMoneyInputs& in);

// The truncated form its callers store. CONFIDENCE: high.
int SystemMoneyIncome(const SystemMoneyInputs& in);

// `ServerSystem::ComputeMaxIncome` -- what `ComputeBudget` and `UpdateBankruptcyLimits`
// actually sum. It runs the output total through `ComputeOutputFromRates` with the rate
// vector "all output to trade" and takes the money channel:
//   trade = roundHalfEven(totalOutput)          (rate 1.0, so the second rounding is a no-op)
//   return max( ftol(SystemMoneyIncomeRaw(...)), 0 )
// The `max(.., 0)` is the original's `jg`: a negative-income colony contributes zero to the
// empire total rather than reducing it.
//
// The two cascade terms `ComputeOutputFromRates` can add to the money channel -- unspent
// industry and unspent terraforming points -- are provably ZERO under this rate vector:
// both are `0 - min(0, need)` with `need >= 0` (infra is clamped to 1, and the terraform
// point count takes fabs after its sign multiply). They are therefore not inputs here.
// CONFIDENCE: high; the science cascade being zero rests on the repair helper returning 0
// for zero science points, which is INFERRED rather than read.
int SystemMaxIncome(double totalOutput, const SystemMoneyInputs& in);

// ---------------------------------------------------------------------------------------
// The turn path: `ServerSystem::ComputeOutput`
// ---------------------------------------------------------------------------------------
//
// `ComputeBudget` has two modes and they take a system's money from DIFFERENT functions.
// Projected mode calls `ComputeMaxIncome` (SystemMaxIncome above). The turn calls
// `ComputeOutput`, which runs `ComputeOutputFromRates` with the system's OWN stored sliders,
// so the construction, terraform and infrastructure channels are funded and three edges that
// are provably dead under the max-income vector are live:
//
//   * the build queue and the ship-repair pass consume construction points;
//   * whatever they leave over is redistributed across trade / terraform / infrastructure
//     and the TRADE share is added to the money channel;
//   * unspent infrastructure points cascade into the terraform pool, and unspent terraform
//     points cascade into the money channel -- two hops, not one.
//
// See sots-re findings/subsystems/output-turn-path.md. Note that the field the earlier
// income-term note called "science" is the SHIP CONSTRUCTION slider; there is no science
// channel in this function (research is bought with money at the empire level).

// `ServerSystem::IdealSuitability` (0x00745d60) -- the suitability the terraform channel
// aims at, and the `==` the rate normaliser tests against.
//   unowned                       -> the system's own suitability (so it is "at its ideal")
//   independent colony            -> the SERVER's per-species baseline for `indi->indsp`
//   otherwise                     -> the OWNER's own `IdealSuit` field
//   and a system-level `dsu` override wins over all three when it is not the sentinel.
// CONFIDENCE: high on the branch order. The sentinel is FLT_MAX: that is INFERRED from the
// corpus (every system carries exactly FLT_MAX there) rather than read out of the data files.
constexpr double kIdealSuitabilityNoOverride = 3.4028234663852886e+38;  // FLT_MAX

struct IdealSuitabilityInputs {
    bool owned = true;
    double systemSuitability = 0;      // sys.Suit
    double ownerIdealSuitability = 0;  // owner's IdealSuit field
    bool independent = false;          // sys.hindi
    double serverIdealSuitability = 0; // server->IdealSuit[indi.indsp], independent only
    double systemOverride = kIdealSuitabilityNoOverride;  // sys.dsu
};
double IdealSuitability(const IdealSuitabilityInputs& in);

// C3 note on `TerraformPointsNeeded` above (0x00746890): the original does NOT round -- the
// `ceil` belongs to `ComputeOutputFromRates`, which applies it to the returned double. Our
// version folds the `ceil` in, which is harmless because `ceil` is idempotent and every
// caller applies it, but the boundary is worth stating. The sign multiply inside the divisor
// is cancelled by a `fabs`, so the result is always >= 0, and a zero `TerraMod` yields +inf,
// which makes the terraform channel absorb its whole pool with nothing cascading to money.
// That branch is UNEXERCISED -- no corpus player carries TerraMod 0.

// `ServerSystem::RepairShipsInOrbit` (0x00751590) -- **the side effect** that makes
// `ComputeOutputFromRates` unsafe to call for its value. It hands each damaged ship of the
// owner's fleets at the system a share of the construction points left over after the build
// queue, round-robin, and returns what is left.
//
// The round robin is EQUIVALENT to `points - min(points, demand)` and this is a proof rather
// than an observation: the per-pass share is `max(points / shipCount, 1)`, so every ship with
// a positive remaining cost takes at least one point per pass, and the loop's only early exit
// requires every remaining cost to be zero. So it ends either with the points exhausted or
// with the demand met. CONFIDENCE: high; the equivalence is pinned by a test.
struct RepairPassResult {
    int spent = 0;
    int left = 0;
};
RepairPassResult RepairShipsInOrbit(int points, int repairDemand);

struct SystemOutputInputs {
    // --- the rate vector, exactly as the system stores it (NOT normalised) ---
    OutputRates rates;
    // The two suppressions the normaliser applies. `suitAtIdeal` is an exact `==` against
    // IdealSuitability(); `infraFull` is `float32(Infra + ibon) >= 1`.
    bool suitAtIdeal = false;
    bool infraFull = false;
    // The leftover-weight test reads the RAW `Infra` against 1.0 and does NOT add the pending
    // bonus, so it is a different predicate from `infraFull` and is carried separately.
    bool infraExactlyOne = false;

    // --- the output total, unrounded (lane N's TotalSystemOutputRaw) ---
    double totalOutputRaw = 0;

    // --- construction ---
    int shipyardStations = 0;   // StationCount(sys, owner, 1)
    int buildQueueDemand = 0;   // sum of `conleft` over the system's build queue
    // Sum of `Ship::RepairCost` over the owner's damaged ships in orbit. NOT modelled from
    // the wire anywhere yet; a caller that cannot compute it must leave it 0 and say so.
    int repairDemand = 0;

    // --- infrastructure ---
    double infra = 0;   // sys.Infra, for `ceil((1 - Infra) / 3.3e-5)`

    // --- terraforming ---
    double terraformPointsNeeded = 0;   // TerraformPointsNeeded(...)
    bool terraformDown = false;         // IdealSuitability() < sys.Suit
    double terraformMod = 1.0;          // owner's TerraMod

    // --- money: every field except `tradePoints`, which this function computes ---
    SystemMoneyInputs money;
};

struct SystemOutput {
    int totalOutput = 0;            // out[0], truncated
    int money = 0;                  // out[3]  <- the ONLY slot ComputeBudget reads
    int construction = 0;           // out[7]
    int constructionToQueue = 0;    // out[8]
    int constructionToRepair = 0;   // out[9]
    double infraDelta = 0;          // out[10], a float32
    double suitabilityDelta = 0;    // out[11], a float32
    // Reported so a caller can see which edges actually carried anything.
    double tradePoints = 0;         // round(total x SRt)
    double leftoverToTrade = 0;     // the construction leftover's trade share
    double leftoverToMoney = 0;     // the terraform leftover that reached the money channel
    OutputRates normalisedRates;
};

// `ServerSystem::ComputeOutput` restricted to the channels the campaign has models for.
// out[1], out[2], out[4], out[5] and out[6] -- the resource ledger, the trade-route income
// pair and the repair demand -- are NOT produced here: none of them feeds `out[3]`, they have
// their own inputs, and inventing them would be coverage theatre.
// CONFIDENCE: high on the channel algebra and the rounding sites (every one read off the
// instruction stream). The repair spend is only as good as `repairDemand`.
SystemOutput ComputeSystemOutput(const SystemOutputInputs& in, const TuningTable& t);

// ---------------------------------------------------------------------------------------
// System bonus and build queue
// ---------------------------------------------------------------------------------------

// Both bonus-pool helpers reset the system's turns-developing counter when they fire on a
// colony that is not the owner's home system. That is reported rather than written, so the
// helpers stay pure; the turn pass applies it.
struct BonusApplyResult {
    bool applied = false;             // the pool actually moved
    bool resetTurnsDeveloping = false;  // the non-home `ntdev = 0` the original performs
};

// Apply a pending population bonus. Nothing happens unless `bonus > 0`; an unowned system
// simply drops the whole pool. Otherwise `applied = min(bonus, cap - pop)` with the colony
// already at or over the cap short-circuiting, `pop += applied`, `bonus -= applied`.
// The population and the capacity are 32-bit in the original (it takes only the low word of
// the capacity), which matters only for colonies past 2^31. CONFIDENCE: high -- read
// instruction by instruction, including the unowned drop and the home-system test.
BonusApplyResult ApplyPopulationBonus(std::int64_t& pop, std::int64_t capacity,
                                      std::int64_t& pendingBonus, bool owned = true,
                                      bool homeSystem = true);

// Apply a pending infrastructure bonus. Nothing happens unless `bonus > 0` and
// `infra < 1`; then `applied = min(bonus, 1 - infra)` and `infra` becomes **exactly 1.0**
// when the pool covers the whole remainder (the original stores the literal, not the sum),
// otherwise `infra + applied`. Every value is rounded to float32. `bonus -= applied`.
// CONFIDENCE: high -- read instruction by instruction, including the exact-1.0 snap.
BonusApplyResult ApplyInfrastructureBonus(double& infra, double& pendingBonus,
                                          bool homeSystem = true);

struct SystemBonusInputs {
    bool stable = false;
    int turnsOwned = 0;                 // current turn - acquisition turn
    // Turns the colony has been *developing* (the system's `ntdev` counter, bumped by this
    // same pass just before the accrual and reset to 0 whenever the system is not stable).
    // B4 correction: the second gate is `ntdev`, not the rebellion turn `rbtn` -- read off
    // the two `cmp ... , SYSTEMBONUS_MINTURNS` at the head of the original.
    int turnsDeveloping = 0;
    bool ownerSpeciesEligible = true;   // ConstantsOf(owner).systemBonusEligible (Zuul: false)
    std::int64_t capacity = 0;          // imperial carrying capacity
};

// Accrue the long-stability bonus when stable, owned for more than SYSTEMBONUS_MINTURNS
// and developing for as long:
//   target      = eligible ? ftol(max(POPBONUS, 0) x cap) : 0
//   popBonus   += min(max(ftol(POPBONUS_INC x cap), 0), max(target - popBonus, 0))
//   infraBonus += min(max(INFRABONUS_INC, 0), max((eligible ? INFRABONUS : 0) - infraBonus, 0))
// The *_HOME keys are not used by the per-turn accrual (only when a home system's
// bonus is initialised, which is not modelled here). Applied next turn by the Apply*
// functions above. CONFIDENCE: high -- increment, target and gating read with constants.
void AccrueSystemBonus(const SystemBonusInputs& in, std::int64_t& popBonus, double& infraBonus,
                       const TuningTable& t);

struct BuildOrder {
    int designId = 0;
    int orderId = 0;
    int constructionCost = 0;   // con
    int constructionLeft = 0;   // conleft
    int moneyCost = 0;          // charged to the owner on completion (0 = free)
    bool moneyAvailable = true; // the owner's answer to the charge; false skips the order
};

struct BuildQueueResult {
    std::vector<int> completedOrderIds;
    int moneyCharged = 0;
    int pointsLeft = 0;   // the function's return value: points not consumed by any order
};

// FIFO consumption of construction points.
//   points <= 0                -> nothing happens, the points are returned unchanged
//   conleft > points           -> `conleft -= points`, `points = 0`, **stop**
//   otherwise                  -> if the design costs money, ask the owner to pay; a refusal
//                                 SKIPS this order and moves to the next one rather than
//                                 stopping the pass. On success build the ship,
//                                 `points -= conleft`, `conleft = 0`, continue.
// Removal is a separate sweep afterwards that unlinks **every** order with `conleft <= 0`,
// including ones that were already at zero before this turn. More than one order can
// complete in a turn. CONFIDENCE: high -- B4 corrected the money-refusal path (continue, not
// stop), the removal sweep's predicate, and that the leftover is the return value.
//
// The `points <= 0` case skips **everything**, the removal sweep included: the entry test
// branches to the epilogue, which returns the argument. Corrected by lane B6 from the
// instruction stream, and carried as a LABELLED HYPOTHESIS rather than a result because no
// save can exercise it: an order can only reach `conleft <= 0` inside this pass, and this
// pass erases it before returning, so a queue never *starts* a turn with one -- unless a
// design with zero construction cost is ever queued, which no corpus save has done.
BuildQueueResult ProcessBuildQueue(std::vector<BuildOrder>& queue, int points);

// ---------------------------------------------------------------------------------------
// Per-player countdown nibbles (`Bats2`, `rcex`)
// ---------------------------------------------------------------------------------------
//
// A system carries two 64-bit words in which each player owns a 4-bit counter: `Bats2` (a
// battle cooldown) and `rcex` (a recon cooldown). Alongside each sits a 32-bit "someone is
// counting" mask with one bit per player. The per-turn pass ticks every counter down by one
// and clears the mask bit of a player whose counter has already reached zero.
//
// Only the first 15 players have a nibble: the original skips index >= 15 entirely (it neither
// ticks nor clears), which is why the mask is a separate word. A counter is clamped to 15 on
// write. CONFIDENCE: high -- read instruction by instruction, including the >= 15 guard.

constexpr int kMaxCountdownPlayers = 15;

struct ColonyCountdowns {
    std::uint64_t counters = 0;  // player i's counter occupies bits [4i, 4i+4)
    std::uint32_t active = 0;    // the companion mask, bit i per player
};

int CountdownFor(const ColonyCountdowns& c, int player);
void SetCountdown(ColonyCountdowns& c, int player, int value);  // clamped to [0, 15]

// One turn of decay over `playerCount` players. A zero counter clears the player's mask bit;
// a non-zero one is written back one lower. Players at or past index 15 are untouched, and
// the whole sweep is skipped when the counter word is already zero. CONFIDENCE: high.
void TickCountdowns(ColonyCountdowns& c, int playerCount);

// ---------------------------------------------------------------------------------------
// Addiction phases and their morale events
// ---------------------------------------------------------------------------------------

// Morale-event ids raised by the addiction sweep, with the per-species morale delta each
// carries (the deltas come from the event-id table the original inlines).
constexpr int kMoraleEventAddictionSuppressed = 0x1b;  // delta -1: the owner has temperance
constexpr int kMoraleEventAddictionOnset = 0x1c;       // delta +1: phase 1
constexpr int kMoraleEventAddictionTerminal = 0x1d;    // delta -2: phase 3

enum class AddictionPhase : int { None = 0, Onset = 1, Established = 2, Terminal = 3 };

// Phase of one species' addiction at this system. `startTurn` is the system's `adt[species]`
// stamp; a zero stamp means the species was never addicted here.
//   0            -> None
//   elapsed > PHASE3_START -> Terminal, elapsed > PHASE2_START -> Established, else Onset
// Both comparisons are strict. CONFIDENCE: high.
AddictionPhase AddictionPhaseOf(int startTurn, int currentTurn, int phase2Start, int phase3Start);

struct MoraleEventRaised {
    int species = 0;
    int eventId = 0;
    int moraleDelta = 0;
};

// ---------------------------------------------------------------------------------------
// The per-system turn pass (`ServerSystem::ProcessTurn`)
// ---------------------------------------------------------------------------------------
//
// The original is a dispatcher: most of a colony turn happens in callees (plague, imperial
// and civilian growth, resources, refuelling, slaves, rebellion, the build queue). What this
// function reproduces is the set of words the dispatcher writes *itself*, plus the two
// bonus-pool helpers whose arithmetic is already modelled above. Everything else is an input
// boundary -- see docs/B4.md for the list and why each entry is there.

struct ColonyTurnState {
    double infra = 0;             // Infra
    double infraBonus = 0;        // ibon, the pending infrastructure bonus
    std::int64_t pop = 0;         // Pop (imperial); overwritten later in the turn by growth
    std::int64_t popBonus = 0;    // pbon, the pending population bonus
    int turnsDeveloping = 0;      // ntdev
    int totalResources = 0;       // TRes, zeroed every turn
    bool growthHalted[3] = {false, false, false};  // haltv, cleared every turn
    ColonyCountdowns battles;     // Bats2 + its mask
    ColonyCountdowns recon;       // rcex + its mask
};

struct ColonyTurnInputs {
    bool owned = false;           // PID != null
    bool independent = false;     // indi != null (an independent colony record exists)
    bool homeSystem = true;       // this is the owner's home system
    bool stable = false;          // IsStable()
    int playerCount = 0;          // players in the game (bounds the countdown sweep)
    int currentTurn = 0;          // StrategyServer ModCount
    int turnsOwned = 0;           // ModCount - TAcq
    bool ownerSpeciesEligible = true;  // the owner's species may accrue system bonuses
    std::int64_t imperialCapacity = 0;  // MaxPop for the imperial group

    // Addiction sweep, one entry per species slot.
    bool civilianPresent[kSpeciesCount] = {};  // civilians (incl. pending bonus) > 0
    bool temperance[kSpeciesCount] = {};       // the owner holds the temperance xenotech
    int addictionStart[kSpeciesCount] = {};    // adt[species]
    int addictionPhase2Start = 0;              // ADDICTION_PHASE2_START
    int addictionPhase3Start = 0;              // ADDICTION_PHASE3_START
};

struct ColonyTurnResult {
    std::vector<MoraleEventRaised> moraleEvents;  // from the addiction sweep, in species order
    int rngDraws = 0;  // words this pass consumes from the server generator
};

// One `ServerSystem::ProcessTurn`, restricted to what the function body writes itself:
//   1. unowned system: `Infra -= UNOWNED_INFRA_DECAY`, floored at 0 (`Infra <= 0` -> 0);
//   2. `ApplyInfrastructureBonus` then `ApplyPopulationBonus` (both run even when unowned);
//   2b. either bonus firing on a colony that is not the owner's home system resets `ntdev`;
//   3. `ntdev` ++ when stable, else reset to 0;
//   4. `AccrueSystemBonus` (gated on owned/stable/turnsOwned/ntdev);
//   5. `TRes = 0`;
//   6. `haltv[0..2] = false`;
//   7. the two countdown sweeps;
//   8. the addiction sweep, which raises one morale event per addicted species.
// The order above is the original's. Steps 3-4 read `stable`, which is a callee's verdict and
// therefore an input here. The pass consumes no RNG at all -- and neither do ProcessPlague,
// the civilian growth pass or ProcessSlaves, all three swept to call depth one and provably
// draw-free. **The only RNG consumer in a colony turn is ProcessRebellion**, and its draw
// count is data dependent (one per million rebels in the counting pass, plus a
// short-circuiting per-species roll), so a system with no rebellion must leave the generator
// untouched -- which is what makes the generator a usable declared region here.
// CONFIDENCE: high on every step listed; the sub-passes themselves are not modelled.
ColonyTurnResult ProcessColonyTurn(ColonyTurnState& s, const ColonyTurnInputs& in,
                                   const TuningTable& t);

}  // namespace sots::sim