// Three bodies under gravity, run three times: a leader and two twins that start one and two
// picounits from it. Units: G = 1, lengths in units of 10^-12 (P), one tick is 10^-4 time units,
// velocities in P per tick. Forces use a softening length of 0.01 so close passes stay finite.
//
// The step is position Verlet in integers: half kick, drift, half kick. A half kick depends only on
// positions, so flipping every velocity and stepping again retraces the path exactly, back to the
// starting bits (Levesque and Verlet, 1993). Each pair force is applied with equal and opposite
// integers, so total momentum is conserved exactly.

const I64 COARSE = 10_000;               // P per force unit u = 10^-8
const I64 EPS2 = 1_000_000_000_000;      // softening 0.01 = 10^6 u, squared
const I64 HALF = 5_000_000_000;          // the half kick's constant over 10^10, see kicks()
const I64 BIG = 2_000_000_000;           // above this many u a square no longer fits
const I64 PUSH_REACH = 100_000_000;      // the push acts within one length unit of the leader

// State, slot copy * 3 + body: copy 0 leads, 1 and 2 are its twins
I64[9] px;
I64[9] py;
I64[9] vx;
I64[9] vy;
I64[3] mass;
U64[3] digest;
I64 tether = 0;

// Half kicks of one copy, filled by kicks(); the kicks of every body at the current positions,
// kept from the end of one tick for the start of the next until load() moves a body
I64[3] hx;
I64[3] hy;
I64[9] kx9;
I64[9] ky9;
Bool kept = false;

// Distance in u between two bodies dx, dy (in u) apart
I64 apart(I64 dx, I64 dy) {
    I64 m = max(abs(dx), abs(dy));
    if (m < BIG) {
        return isqrt(dx * dx + dy * dy);
    }
    I64 s = m / BIG + 1;
    I64 ex = dx / s;
    I64 ey = dy / s;
    return isqrt(ex * ex + ey * ey) * s;
}

// Softened distance in u between two bodies dx, dy (in u) apart
I64 dist(I64 dx, I64 dy) {
    I64 m = max(abs(dx), abs(dy));
    if (m < BIG) {
        return isqrt(dx * dx + dy * dy + EPS2);
    }
    I64 s = m / BIG + 1;
    I64 ex = dx / s;
    I64 ey = dy / s;
    return isqrt(ex * ex + ey * ey) * s;
}

// Half kick per unit mass from a body dx, dy away: G dx / d^3 times dt^2 / 2 in P per tick.
// With lengths in u (10^8 per unit) that is dx / d^3 times 5 * 10^19, applied in three roundings
void kicks(I64 c) {
    for b in 0..3 {
        hx[b] = 0;
        hy[b] = 0;
    }
    for i in 0..3 {
        for j in i + 1..3 {
            I64 dx = (px[c * 3 + j] - px[c * 3 + i]) / COARSE;
            I64 dy = (py[c * 3 + j] - py[c * 3 + i]) / COARSE;
            I64 d = dist(dx, dy);
            I64 kx = muldiv(I64, muldiv(I64, muldiv(I64, HALF, dx, d, half_even), 10_000_000_000, d, half_even), 1, d, half_even);
            I64 ky = muldiv(I64, muldiv(I64, muldiv(I64, HALF, dy, d, half_even), 10_000_000_000, d, half_even), 1, d, half_even);
            hx[i] = hx[i] + mass[j] * kx;
            hy[i] = hy[i] + mass[j] * ky;
            hx[j] = hx[j] - mass[i] * kx;
            hy[j] = hy[j] - mass[i] * ky;
        }
    }
}

void kick(I64 c) {
    for b in 0..3 {
        vx[c * 3 + b] = vx[c * 3 + b] + kx9[c * 3 + b];
        vy[c * 3 + b] = vy[c * 3 + b] + ky9[c * 3 + b];
    }
}

void keep(I64 c) {
    kicks(c);
    for b in 0..3 {
        kx9[c * 3 + b] = hx[b];
        ky9[c * 3 + b] = hy[b];
    }
}

// Consensus tether on the twins' velocities. Pulling (k > 0) is a critically damped spring toward
// the leader's body with stiffness k^4 / 10^11 per tick squared: below k = 10 it barely holds
// against close passes, from 10 to 40 it draws the twins in over a few seconds of play, and at 100
// it settles in about 150 ticks and tears apart any pair the softening has bound. Pushing (k < 0)
// is a weak spring away from the leader, undamped, that lets go one length unit out, so it starts
// the split and gravity does the rest
void pull_twins() {
    I64 k = tether;
    if (k == 0) {
        return;
    }
    for c in 1..3 {
        for b in 0..3 {
            I64 l = b;
            I64 t = c * 3 + b;
            if (k > 0) {
                vx[t] = vx[t] + muldiv(I64, px[l] - px[t], k * k * k * k, 100_000_000_000, half_even) + muldiv(I64, vx[l] - vx[t], 2 * k * k, 316_228, half_even);
                vy[t] = vy[t] + muldiv(I64, py[l] - py[t], k * k * k * k, 100_000_000_000, half_even) + muldiv(I64, vy[l] - vy[t], 2 * k * k, 316_228, half_even);
            } else if (apart((px[l] - px[t]) / COARSE, (py[l] - py[t]) / COARSE) < PUSH_REACH) {
                vx[t] = vx[t] - muldiv(I64, px[l] - px[t], k * k, 250_000_000_000, half_even);
                vy[t] = vy[t] - muldiv(I64, py[l] - py[t], k * k, 250_000_000_000, half_even);
            }
        }
    }
}

// One tick of all three copies
void step() {
    if (!kept) {
        for c in 0..3 {
            keep(c);
        }
        kept = true;
    }
    for c in 0..3 {
        kick(c);
        for b in 0..3 {
            px[c * 3 + b] = px[c * 3 + b] + vx[c * 3 + b];
            py[c * 3 + b] = py[c * 3 + b] + vy[c * 3 + b];
        }
        keep(c);
        kick(c);
        I64 o = c * 3;
        digest[c] = (digest[c] *% 1_099_511_628_211) ^ (px[o] as% U64) ^ ((py[o] as% U64) <<% 13) ^ ((px[o + 1] as% U64) <<% 26) ^ ((py[o + 1] as% U64) <<% 39) ^ ((px[o + 2] as% U64) <<% 52) ^ ((py[o + 2] as% U64) <<% 7);
    }
    pull_twins();
}

void put(I64 b, I64 m, I64 x, I64 y, I64 u, I64 v) {
    mass[b] = m;
    for c in 0..3 {
        px[c * 3 + b] = x;
        py[c * 3 + b] = y;
        vx[c * 3 + b] = u;
        vy[c * 3 + b] = v;
    }
}

// Loads a preset into all three copies, then moves the twins' first body 1 and 2 P along x.
// 0: the Pythagorean problem (masses 3, 4, 5 at rest on a 3-4-5 triangle), chaotic.
// 1: the figure eight of Chenciner and Montgomery (three equal masses), stable
export I64 start(I64 preset) {
    if (preset == 1) {
        put(0, 1, 970_004_360_000, -243_087_530_000, 46_620_368, 43_236_573);
        put(1, 1, -970_004_360_000, 243_087_530_000, 46_620_369, 43_236_573);
        put(2, 1, 0, 0, -93_240_737, -86_473_146);
    } else {
        put(0, 3, 1_000_000_000_000, 3_000_000_000_000, 0, 0);
        put(1, 4, -2_000_000_000_000, -1_000_000_000_000, 0, 0);
        put(2, 5, 1_000_000_000_000, -1_000_000_000_000, 0, 0);
    }
    px[3] = px[3] + 1_000_000;
    px[6] = px[6] + 2_000_000;
    for c in 0..3 {
        digest[c] = 14_695_981_039_346_656_037;
    }
    kept = false;
    return preset;
}

// The page drives the program through these exports. A state slot is field * 9 + copy * 3 + body,
// fields px, py, vx, vy
export I64 load(I64 slot, I64 value) {
    I64 s = slot % 9;
    I64 f = slot / 9;
    kept = false;
    if (f == 0) { px[s] = value; }
    if (f == 1) { py[s] = value; }
    if (f == 2) { vx[s] = value; }
    if (f == 3) { vy[s] = value; }
    return value;
}

export I64 read(I64 slot) {
    I64 s = slot % 9;
    I64 f = slot / 9;
    if (f == 0) { return px[s]; }
    if (f == 1) { return py[s]; }
    if (f == 2) { return vx[s]; }
    return vy[s];
}

export U64 seed(I64 c, U64 value) {
    digest[c] = value;
    return value;
}

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

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

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

// Flips every velocity. With the tether at 0, as many ticks again bring every body back to its start
export I64 flip() {
    for s in 0..9 {
        vx[s] = -vx[s];
        vy[s] = -vy[s];
    }
    return 0;
}

// The largest distance in u from a body of copy c to the same body of the leader
export I64 spread(I64 c) {
    I64 most = 0;
    for b in 0..3 {
        most = max(most, apart((px[c * 3 + b] - px[b]) / COARSE, (py[c * 3 + b] - py[b]) / COARSE));
    }
    return most;
}

// Total momentum along x and y of copy c, in mass times P per tick
export I64 momentum(I64 c, I64 axis) {
    I64 total = 0;
    for b in 0..3 {
        if (axis == 0) {
            total = total + mass[b] * vx[c * 3 + b];
        } else {
            total = total + mass[b] * vy[c * 3 + b];
        }
    }
    return total;
}

// The largest coordinate difference in P between the leader now and the leader at slot values s0
I64 offset(I64[12] s0) {
    I64 most = 0;
    for b in 0..3 {
        most = max(most, abs(px[b] - s0[b]));
        most = max(most, abs(py[b] - s0[3 + b]));
        most = max(most, abs(vx[b] + s0[6 + b]));
        most = max(most, abs(vy[b] + s0[9 + b]));
    }
    return most;
}

I64 there_and_back(I64 preset, I64 ticks) {
    _ = start(preset);
    I64[12] s0;
    for b in 0..3 {
        s0[b] = px[b];
        s0[3 + b] = py[b];
        s0[6 + b] = vx[b];
        s0[9 + b] = vy[b];
    }
    _ = advance(ticks);
    _ = flip();
    _ = advance(ticks);
    return offset(s0);
}

void main() {
    _ = start(0);
    _ = advance(20_000);
    "after 20000 ticks (2.0 time units), Pythagorean problem:\n";
    "leader digest={digest[0]}\n";
    "twin 1 digest={digest[1]}\n";
    "twin 2 digest={digest[2]}\n";
    "spread twin 1={spread(1)} u, twin 2={spread(2)} u\n";
    "momentum x={momentum(0, 0)}, y={momentum(0, 1)}\n";
    "20000 ticks out and 20000 back, largest miss: {there_and_back(0, 20_000)} P\n";
}

test "time reversal returns to the starting bits" {
    assert(there_and_back(0, 3_000) == 0);
    assert(there_and_back(1, 3_000) == 0);
}

test "pair forces conserve momentum exactly" {
    _ = start(1);
    I64 mx = momentum(0, 0);
    I64 my = momentum(0, 1);
    _ = advance(2_000);
    assert(momentum(0, 0) == mx && momentum(0, 1) == my);
}

test "the tether pulls the twins back onto the leader" {
    _ = start(0);
    _ = advance(300_000);
    I64 wide = spread(1);
    _ = pull(100);
    _ = advance(20_000);
    assert(wide > 10_000_000 && spread(1) * 100 < wide);
}

test "the digests at tick 20000 match the receipt" {
    _ = start(0);
    _ = advance(20_000);
    assert(digest[0] == 15_110_298_595_637_020_077);
    assert(digest[1] == 9_072_904_302_982_474_941);
    assert(digest[2] == 13_663_343_393_609_666_067);
}