// Gravitational wave chirp: 3.5PN binary inspiral in checked 1024-bit integers (I1024).
// Evaluates the TaylorT4 post-Newtonian phasing of compact binary coalescences
// down to the Innermost Stable Circular Orbit (ISCO) at v = 1/sqrt(6), followed
// by spliced quasinormal mode (QNM) ringdown of the remnant Kerr black hole.
// Three copies run concurrently: Leader A and Twins B and C with small mass
// offsets, linked by a relativistic consensus tether on orbital phase.

const I1024 S = 1_0000000000_0000000000_0000000000_0000000000_0000000000_0000000000;
const I1024 PI_S = 3_1415926535_8979323846_2643383279_5028841971_6939937510_5820974945;
const I1024 TWO_PI_S = 6_2831853071_7958647692_5286766559_0057683943_3879875021_1641949889;
const I1024 HALF_PI_S = 1_5707963267_9489661923_1321691639_7514420985_8469968755_2910487472;
const I1024 GAMMA_E_S = 5772156649_0153286060_6512090082_4024310421_5933593992_3598805767;
const I1024 LN2_S = 6931471805_5994530941_7232121458_1765680755_0013436025_5254120680;

// 3.5PN TaylorT4 coefficients for equal masses (eta = 1/4)
// dv/dt = (8/5 M) * v^9 * [ 1 + c2*v^2 + c3*v^3 + c4*v^4 + c5*v^5 + (c6_const + c6_log*ln(4v))*v^6 + c7*v^7 ]
const I1024 C2_S = -2_8988095238_0952380952_3809523809_5238095238_0952380952_3809523810;
const I1024 C3_S = 12_5663706143_5917295385_0573533118_0115367886_7759750042_3283899778;
const I1024 C4_S = 3_7785080467_3721340388_0070546737_2134038800_7054673721_3403880071;
const I1024 C5_S = -37_9983111434_1940393188_1496159666_3682183848_1083053699_4215601711;
const I1024 C6_CONST_S = 96_4276677848_7411383988_2700892954_6386460020_2111406056_1272711614;
const I1024 C6_LOG_S = -16_3047619047_6190476190_4761904761_9047619047_6190476190_4761904762;
const I1024 C7_S = 55_0194268962_0219904397_5394701879_1245856752_9477212486_9139825749;

// ISCO velocity v_ISCO = 1 / sqrt(6)
const I1024 V_ISCO_S = 4082482904_6386301636_6214012450_9818986609_9124677611_1688072115;

// Bound limits for range reduction in ln(u)
const I1024 LIMIT_LO_S = 7071067811_8654752440_0844362104_8490392848_3593768847_4036588340;
const I1024 LIMIT_HI_S = 1_4142135623_7309504880_1688724209_6980785696_7187537694_8073176680;

// Integration step dt = 25 M
const I1024 DT_S = 25 * S;

// QNM ringdown frequency M * omega_QNM = 0.55 and decay factor per tick (0.96)
const I1024 OMEGA_QNM_S = 55 * S / 100;
const I1024 QNM_DECAY_S = 96 * S / 100;

// Output display scaling (10^6 for micrometers / micro-radians)
const I1024 S_MICRONS = 1_000_000;
const I1024 S_RATIO = S / S_MICRONS;

// State across three copies: 0: Leader, 1: Twin B, 2: Twin C
I1024[3] v_state;
I1024[3] phi_state;
I1024[3] r_state;
I1024[3] hp_state;
I1024[3] hx_state;
I1024[3] f_state;
I1024[3] m_state;
I1024[3] amp_state;
I64[3] mode_state;
U64[3] digest_state;

I64[6] body_x;
I64[6] body_y;

I64 tether = 0;

I1024 qmul(I1024 a, I1024 b) {
    return muldiv(I1024, a, b, S, half_even);
}

I1024 qdiv(I1024 a, I1024 b) {
    return muldiv(I1024, a, S, b, half_even);
}

I1024 qsqrt(I1024 a) {
    return isqrt(a * S);
}

// Horner series for ln(u) reduced by ln 2
I1024 qln(I1024 u_in) {
    I1024 u = u_in;
    I1024 p = 0;

    while (u < LIMIT_LO_S) {
        u = u * 2;
        p = p - 1;
    }
    while (u > LIMIT_HI_S) {
        u = u / 2;
        p = p + 1;
    }

    I1024 z = qdiv(u - S, u + S);
    I1024 z2 = qmul(z, z);
    I1024 term = z;
    I1024 total = 0;

    for k in 0..45 {
        I1024 denom = ((2 * k + 1) as I1024) * S;
        total = total + qdiv(term, denom);
        term = qmul(term, z2);
    }

    return 2 * total + p * LN2_S;
}

// Fixed-point sine via Taylor series on [-pi/2, pi/2]
I1024 qsin(I1024 x) {
    I1024 v = x % TWO_PI_S;
    if (v > PI_S) {
        v = v - TWO_PI_S;
    } else if (v < -PI_S) {
        v = v + TWO_PI_S;
    }
    if (v > HALF_PI_S) {
        v = PI_S - v;
    } else if (v < -HALF_PI_S) {
        v = -PI_S - v;
    }

    I1024 v2 = qmul(v, v);
    I1024 term = v;
    I1024 total = 0;
    I1024 sgn = 1;

    for k in 0..16 {
        I1024 term_val = term;
        if (sgn > 0) {
            total = total + term_val;
        } else {
            total = total - term_val;
        }
        sgn = -sgn;
        I1024 denom = (((2 * k + 2) * (2 * k + 3)) as I1024) * S;
        term = qdiv(qmul(term, v2), denom);
    }
    return total;
}

I1024 qcos(I1024 x) {
    return qsin(x + HALF_PI_S);
}

// 3.5PN acceleration dv/dt = (8/5 M) * v^9 * A(v)
I1024 taylor_t4_dvdt(I1024 v, I1024 m) {
    I1024 v2 = qmul(v, v);
    I1024 v3 = qmul(v2, v);
    I1024 v4 = qmul(v2, v2);
    I1024 v5 = qmul(v4, v);
    I1024 v6 = qmul(v3, v3);
    I1024 v7 = qmul(v4, v3);
    I1024 v9 = qmul(v5, v4);

    I1024 ln_4v = qln(4 * v);
    I1024 c6_term = C6_CONST_S + qmul(C6_LOG_S, ln_4v);

    I1024 a = S + qmul(C2_S, v2) + qmul(C3_S, v3) + qmul(C4_S, v4) + qmul(C5_S, v5) + qmul(c6_term, v6) + qmul(C7_S, v7);
    I1024 v9_a = qmul(v9, a);
    I1024 prefactor = qdiv(8 * S / 5, m);
    return qmul(prefactor, v9_a);
}

void step() {
    I1024 leader_phi = phi_state[0];
    for c in 0..3 {
        I1024 v = v_state[c];
        I1024 phi = phi_state[c];
        I1024 m = m_state[c];
        I1024 amp = amp_state[c];
        I64 mode = mode_state[c];

        if (mode == 0) {
            // 3.5PN Inspiral
            I1024 dv = qmul(taylor_t4_dvdt(v, m), DT_S);
            I1024 v_next = v + dv;

            if (v_next >= V_ISCO_S) {
                v_next = V_ISCO_S;
                mode_state[c] = 1; // Transition to QNM ringdown
                amp = qmul(v_next, v_next);
            }

            // Orbital angular velocity omega = v^3 / M
            I1024 v3 = qmul(v_next, qmul(v_next, v_next));
            I1024 omega = qdiv(v3, m);

            // Relativistic consensus tether acting directly on phase
            if (c > 0 && tether != 0) {
                I1024 dphi = phi - leader_phi;
                I1024 torque = muldiv(I1024, dphi, (tether as I1024) * S, 25000 * S, half_even);
                omega = omega - torque;
            }

            I1024 phi_next = phi + qmul(omega, DT_S);
            I1024 r = qdiv(m, qmul(v_next, v_next));
            amp = qmul(v_next, v_next);

            v_state[c] = v_next;
            phi_state[c] = phi_next;
            r_state[c] = r;
            amp_state[c] = amp;
            hp_state[c] = qmul(amp, qcos(2 * phi_next));
            hx_state[c] = qmul(amp, qsin(2 * phi_next));
            f_state[c] = qdiv(omega, PI_S);
        } else {
            // Post-ISCO QNM Ringdown
            I1024 amp_next = qmul(amp, QNM_DECAY_S);
            I1024 phi_next = phi + qmul(OMEGA_QNM_S, DT_S);
            I1024 r = qdiv(r_state[c] + 17 * S / 10, 2 * S); // Horizon contraction

            amp_state[c] = amp_next;
            phi_state[c] = phi_next;
            r_state[c] = r;
            hp_state[c] = qmul(amp_next, qcos(2 * phi_next));
            hx_state[c] = qmul(amp_next, qsin(2 * phi_next));
            f_state[c] = qdiv(OMEGA_QNM_S, PI_S);
        }

        // Coordinates for the 2 bodies in orbital plane
        I1024 sep = r_state[c];
        I1024 current_phi = phi_state[c];
        I1024 half_r = sep / 2;
        I1024 x1_wide = qmul(half_r, qcos(current_phi));
        I1024 y1_wide = qmul(half_r, qsin(current_phi));

        I64 x1 = (x1_wide / S_RATIO) as I64;
        I64 y1 = (y1_wide / S_RATIO) as I64;

        body_x[c * 2] = x1;
        body_y[c * 2] = y1;
        body_x[c * 2 + 1] = -x1;
        body_y[c * 2 + 1] = -y1;

        // FNV-1a 64-bit rolling state digest
        I64 v_scaled = (v_state[c] / S_RATIO) as I64;
        I64 phi_scaled = (phi_state[c] / S_RATIO) as I64;
        I64 hp_scaled = (hp_state[c] / S_RATIO) as I64;
        I64 hx_scaled = (hx_state[c] / S_RATIO) as I64;
        U64 phi_limb1 = (phi_state[c] >> 64) as% U64;
        U64 phi_limb2 = (phi_state[c] >> 128) as% U64;

        U64 d = digest_state[c];
        d = (d ^ (v_scaled as% U64)) *% 1_099_511_628_211;
        d = (d ^ (phi_scaled as% U64)) *% 1_099_511_628_211;
        d = (d ^ phi_limb1) *% 1_099_511_628_211;
        d = (d ^ phi_limb2) *% 1_099_511_628_211;
        d = (d ^ (hp_scaled as% U64)) *% 1_099_511_628_211;
        d = (d ^ (hx_scaled as% U64)) *% 1_099_511_628_211;
        digest_state[c] = d;
    }
}

export I64 start(I64 preset) {
    I1024 v0 = 2 * S / 10; // v0 = 0.20
    if (preset == 1) {
        v0 = 25 * S / 100; // NS-NS preset: v0 = 0.25
    } else if (preset == 2) {
        v0 = 32 * S / 100; // Late inspiral plunge: v0 = 0.32
    }

    for c in 0..3 {
        v_state[c] = v0;
        phi_state[c] = 0;
        m_state[c] = S;
        amp_state[c] = qmul(v0, v0);
        mode_state[c] = 0;
        digest_state[c] = 14_695_981_039_346_656_037;
    }

    // Mass perturbation for twins: +10^-18 M and +2*10^-18 M
    m_state[1] = S + 1_0000000000_0000000000_0000000000_0000000000_00;
    m_state[2] = S + 2_0000000000_0000000000_0000000000_0000000000_00;

    for c in 0..3 {
        I1024 r = qdiv(m_state[c], qmul(v0, v0));
        r_state[c] = r;
        hp_state[c] = amp_state[c];
        hx_state[c] = 0;
        f_state[c] = qdiv(qmul(v0, qmul(v0, v0)), PI_S);

        I1024 half_r = r / 2;
        I64 x1 = (half_r / S_RATIO) as I64;
        body_x[c * 2] = x1;
        body_y[c * 2] = 0;
        body_x[c * 2 + 1] = -x1;
        body_y[c * 2 + 1] = 0;
    }

    tether = 0;
    return preset;
}

export I64 advance(I64 ticks) {
    for i in 0..ticks {
        step();
    }
    return ticks;
}

export I64 pull(I64 k) {
    tether = k;
    return k;
}

export I64 read_slot(I64 c, I64 slot) {
    if (slot == 0) { return body_x[c * 2]; }
    if (slot == 1) { return body_y[c * 2]; }
    if (slot == 2) { return body_x[c * 2 + 1]; }
    if (slot == 3) { return body_y[c * 2 + 1]; }
    if (slot == 4) { return (r_state[c] / S_RATIO) as I64; }
    if (slot == 5) { return (v_state[c] / S_RATIO) as I64; }
    if (slot == 6) { return (phi_state[c] / S_RATIO) as I64; }
    if (slot == 7) { return (f_state[c] / S_RATIO) as I64; }
    if (slot == 8) { return (hp_state[c] / S_RATIO) as I64; }
    if (slot == 9) { return (hx_state[c] / S_RATIO) as I64; }
    if (slot == 10) { return mode_state[c]; }
    if (slot == 11) { return (amp_state[c] / S_RATIO) as I64; }
    if (slot == 12) { return ((phi_state[c] - phi_state[0]) / S_RATIO) as I64; }
    return 0;
}

export U64 digest_of(I64 c) {
    return digest_state[c];
}

// -----------------------------------------------------------------------------
// Tests
// -----------------------------------------------------------------------------

test "3.5PN acceleration agrees with analytical expansion at v=0.2" {
    I1024 v_02 = 2 * S / 10;
    I1024 dv = taylor_t4_dvdt(v_02, S);
    // dv/dt(0.2) = 8.07381572020032153471378043462343847605407087544984810e-7 in scale S
    I1024 expected = 8073_8157202003_2153471378_0434623438_4760540708_7544984810;
    I1024 diff = dv - expected;
    if (diff < 0) { diff = -diff; }
    // Agreement to within 1 part in 10^30
    assert(diff < 1_0000000000_0000000000_0000000000);
}

test "quadrupole strain matches analytical scaling" {
    I1024 v1 = S / 10;
    I1024 v2 = 2 * S / 10;
    I1024 amp1 = qmul(v1, v1);
    I1024 amp2 = qmul(v2, v2);
    // Quadrupole amplitude scales as v^2, so doubling v quadruples amplitude
    assert(amp2 == 4 * amp1);
}

test "ISCO transition halts circular inspiral at 1/sqrt(6)" {
    _ = start(0);
    // At start, mode is 0 (inspiral)
    assert(mode_state[0] == 0);
    // Run past ISCO
    _ = advance(1300);
    // Must be in QNM ringdown mode (1)
    assert(mode_state[0] == 1);
    // Separation has contracted below ISCO separation (r_ISCO = 6)
    I1024 r_isco = 6 * S;
    assert(r_state[0] < r_isco);
}

test "consensus tether pulls twins into relativistic phase lock" {
    _ = start(0);
    _ = pull(100);
    _ = advance(800);
    // Under strong tether pull, twin dephasing remains tightly locked (|dphi| < 10^-15)
    I1024 dphi_b = phi_state[1] - phi_state[0];
    if (dphi_b < 0) { dphi_b = -dphi_b; }
    assert(dphi_b < 1_0000000000_0000000000_0000000000_0000000000_00000);
}

test "the digests at tick 1000 match the receipt" {
    _ = start(0);
    _ = advance(1000);
    assert(digest_of(0) == 14_576_398_344_542_445_859);
    assert(digest_of(1) == 16_148_053_973_575_091_900);
    assert(digest_of(2) == 6_304_118_998_254_214_534);
}

void main() {
    _ = start(0);
    _ = advance(1000);
    U64 d0 = digest_of(0);
    U64 d1 = digest_of(1);
    U64 d2 = digest_of(2);
    "chirp receipt: tick=1000 d0={d0} d1={d1} d2={d2}\n";
}