demos / pendulum

A Double Pendulum Swinging in 3D

A cint program swings three double pendulums under gravity in whole-number arithmetic. One leads, and two twins start one and two millionths of a radian from it, yet they still drift away. The consensus tether pulls them back or pushes them apart.

Run the demo

pendulum.ci Emitted C

shell

$ cint run pendulum.ci
after 500 steps (5.0 seconds):
pendulum 0: digest=14995127728425866399
pendulum 1: digest=8897290559184475167
pendulum 2: digest=11575834901304267551
divergence A to B: 1519 um

pendulum.ci

Pendulum A (Leader): Baseline

Pendulum B (Twin 1): +1 μrad Azimuth

Pendulum C (Twin 2): +2 μrad Azimuth

Space Play/Pause · Shift+R Restart Sim · R Reset View

0 / endless · 0.000s

Consensus Tether

Off

◄ Push 0 (off) Pull ►

Twins split Free Twins follow

Tip Separation (Leader A → Twin B): 0 μm (0.0000 m)

Tip Separation (Leader A → Twin C): 0 μm (0.0000 m)

Colors

Kicks

In slingshot mode, drag from any bob to give it a push.

Trail

Keep the whole trail

Trail length 400 ticks (4.0s)

0.5 s 30.0 s

Line weight

0.8px thin

1.5px normal

3.0px thick

Pendulum A (Leader)

Bob 2 Position

(0.000, 0.000, 0.000) m

Polar Angles (θ₁, θ₂)

1.5708 rad

Azimuth Angles (φ₁, φ₂)

0.0000 rad

Angular Speeds (ω₁, ω₂)

0.00 rad/s

64-Bit State Digest

0xCBF29CE484222325

Twins B and C State Digests

Twin B (+1 μrad) Bob 2 / Digest

(0.000, 0.000, 0.000) m

0xCBF29CE484222325

Twin C (+2 μrad) Bob 2 / Digest

(0.000, 0.000, 0.000) m

0xCBF29CE484222325

Tick 500, Page and Program

cintc 0.2.0 compiled pendulum.ci, and cint run printed these digests after 500 steps. The page and the program use the same integer sines, square roots and 64-bit wrapping digests, with no floating point in the physics.

Pendulum	Step 500 State Digest	This page
Leader A (0)	`14995127728425866399`	checking
Twin B (1)	`8897290559184475167`	checking
Twin C (2)	`11575834901304267551`	checking

Loading the engine.

main, run in this page

running

outcome value faults 0 fuel-consumed 41509 stdout-sha256 ebe8373adc18ed8e0b6e6246d295bc15900aec02f393e777011c7977c5ae7683 cintc 0.2.0, GCC: same stdout, 5 of 5 tests pass emitted C: pendulum.c

pendulum.ci

Download .ci

const I64 S = 1_000_000;
const I64 PI = 3_141_593;
const I64 TWO_PI = 6_283_185;
const I64 HALF_PI = 1_570_796;
const I64 G = 9_806_650;
const I64 L = 1_000_000;
const I64 DT = 10_000;
const I64 W_BLEED = 10_000_000;
const I64 W_MAX = 12_000_000;

// Speed governor, the same as the page's: above 10 rad/s an angular speed loses 0.5% a tick,
// and it never passes 12 rad/s
I64 govern(I64 w) {
    I64 v = w;
    if (abs(v) > W_BLEED) {
        v = v - (v * 5) / 1000;
    }
    return clamp(v, -W_MAX, W_MAX);
}

I64 isin(I64 x) {
    I64 v = x % TWO_PI;
    if (v > PI) {
        v = v - TWO_PI;
    }
    if (v < -PI) {
        v = v + TWO_PI;
    }
    if (v > HALF_PI) {
        v = PI - v;
    }
    if (v < -HALF_PI) {
        v = -PI - v;
    }

    I64 x2 = (v * v + S / 2) / S;
    I64 x3 = (x2 * v + S / 2) / S;
    I64 x5 = (x3 * x2 + S / 2) / S;
    I64 x7 = (x5 * x2 + S / 2) / S;
    return v - x3 / 6 + x5 / 120 - x7 / 5040;
}

I64 icos(I64 x) {
    return isin(x + HALF_PI);
}

// The state of the three pendulums: polar angles th1 and th2, azimuths phi1 and phi2 and their
// angular speeds, one element per pendulum. Pendulum 0 leads; 1 and 2 are its twins
I64[3] th1;
I64[3] th2;
I64[3] phi1;
I64[3] phi2;
I64[3] w1;
I64[3] w2;
I64[3] w_phi1;
I64[3] w_phi2;
U64[3] digest;

// Tether strength, -100 to 100: pulling above 0, pushing below
I64 tether = 0;

// The tips after the last step, and the unit vector that unit() leaves
I64[3] step_x2;
I64[3] step_y2;
I64[3] step_z2;
I64 ux = 0;
I64 uy = 0;
I64 uz = 0;

// The joints and tips of the current state, which locate() fills: x1 y1 z1 x2 y2 z2 per pendulum
I64[18] places;

// Unit vector of a link from the sine and cosine of its polar angle and its azimuth
void unit(I64 s, I64 c, I64 phi) {
    I64 cphi = icos(phi);
    I64 sphi = isin(phi);
    I64 x = (s * cphi + S / 2) / S;
    I64 y = (s * sphi + S / 2) / S;
    I64 z = c;
    I64 mag = isqrt(x * x + y * y + z * z);
    if (mag == 0) {
        mag = 1;
    }
    ux = (x * S + mag / 2) / mag;
    uy = (y * S + mag / 2) / mag;
    uz = (z * S + mag / 2) / mag;
}

// The nominal start: both links level, the twins one and two millionths of a radian round from
// the leader, everything at rest
void start() {
    for p in 0..3 {
        th1[p] = HALF_PI;
        th2[p] = HALF_PI;
        phi1[p] = p;
        phi2[p] = 0;
        w1[p] = 0;
        w2[p] = 0;
        w_phi1[p] = 0;
        w_phi2[p] = 0;
        digest[p] = 14_695_981_039_346_656_037;
    }
}

// One tick of 10 ms for all three pendulums
void step() {
    I64[3] alpha1;
    I64[3] alpha2;
    I64[3] alpha_phi1;
    I64[3] alpha_phi2;
    I64[3] s1_arr;
    I64[3] c1_arr;
    I64[3] s2_arr;
    I64[3] c2_arr;

    for p in 0..3 {
        I64 dth = th1[p] - th2[p];
        I64 s1 = isin(th1[p]);
        I64 c1 = icos(th1[p]);
        I64 s2 = isin(th2[p]);
        I64 c2 = icos(th2[p]);
        s1_arr[p] = s1;
        c1_arr[p] = c1;
        s2_arr[p] = s2;
        c2_arr[p] = c2;

        I64 sd = isin(dth);
        I64 cd = icos(dth);
        I64 c2d = icos(2 * dth);

        I64 den = 3 * S - c2d;
        if (den == 0) {
            den = 1;
        }

        I64 w1_sq = (w1[p] * w1[p] + S / 2) / S;
        I64 w2_sq = (w2[p] * w2[p] + S / 2) / S;

        I64 num1 = -3 * (G * s1 / S) - (G * isin(th1[p] - 2 * th2[p]) / S) - 2 * (sd * (w2_sq + (w1_sq * cd / S)) / S);
        alpha1[p] = (num1 * S) / den;

        I64 g_c1 = (G * c1) / S;
        I64 num2 = 2 * (sd * (2 * w1_sq + 2 * g_c1 + (w2_sq * cd / S)) / S);
        alpha2[p] = (num2 * S) / den;

        I64 cent1 = ((w_phi1[p] * w_phi1[p] / S) * (s1 * c1 / S)) / S;
        I64 cent2 = ((w_phi2[p] * w_phi2[p] / S) * (s2 * c2 / S)) / S;
        alpha1[p] = alpha1[p] + cent1;
        alpha2[p] = alpha2[p] + cent2;

        I64 dphi = phi1[p] - phi2[p];
        I64 sd_phi = isin(dphi);
        I64 tau = (G * sd_phi) / S;
        alpha_phi1[p] = -tau;
        alpha_phi2[p] = tau;
    }

    // Consensus tether. Pulling (k > 0): a spring toward the leader on each angle, stiffness up to
    // gravity's own (G / L), critically damped. Pushing (k < 0): a weak spring away from the leader
    // with negative damping, so the separation grows slowly until the speed governor holds it
    I64 coupling_k = tether;
    if (coupling_k > 0) {
        I64 kp = (G * coupling_k * coupling_k) / 10_000;
        I64 kd = 2 * isqrt(kp * S);
        for p in 1..3 {
            alpha1[p] = alpha1[p] + (isin(th1[0] - th1[p]) * kp) / S + ((w1[0] - w1[p]) * kd) / S;
            alpha2[p] = alpha2[p] + (isin(th2[0] - th2[p]) * kp) / S + ((w2[0] - w2[p]) * kd) / S;
            alpha_phi1[p] = alpha_phi1[p] + (isin(phi1[0] - phi1[p]) * kp) / S + ((w_phi1[0] - w_phi1[p]) * kd) / S;
            alpha_phi2[p] = alpha_phi2[p] + (isin(phi2[0] - phi2[p]) * kp) / S + ((w_phi2[0] - w_phi2[p]) * kd) / S;
        }
    }
    if (coupling_k < 0) {
        I64 c_ang = (coupling_k * 100_000) / 100;
        I64 c_vel = (coupling_k * 2_000) / 100;
        for p in 1..3 {
            alpha1[p] = alpha1[p] + (isin(th1[0] - th1[p]) * c_ang) / S + ((w1[0] - w1[p]) * c_vel) / 1000;
            alpha2[p] = alpha2[p] + (isin(th2[0] - th2[p]) * c_ang) / S + ((w2[0] - w2[p]) * c_vel) / 1000;
            alpha_phi1[p] = alpha_phi1[p] + (isin(phi1[0] - phi1[p]) * c_ang) / S + ((w_phi1[0] - w_phi1[p]) * c_vel) / 1000;
            alpha_phi2[p] = alpha_phi2[p] + (isin(phi2[0] - phi2[p]) * c_ang) / S + ((w_phi2[0] - w_phi2[p]) * c_vel) / 1000;
        }
    }

    for p in 0..3 {
        w1[p] = w1[p] + (alpha1[p] * DT) / S;
        w2[p] = w2[p] + (alpha2[p] * DT) / S;
        w_phi1[p] = w_phi1[p] + (alpha_phi1[p] * DT) / S;
        w_phi2[p] = w_phi2[p] + (alpha_phi2[p] * DT) / S;
        w1[p] = govern(w1[p]);
        w2[p] = govern(w2[p]);
        w_phi1[p] = govern(w_phi1[p]);
        w_phi2[p] = govern(w_phi2[p]);

        th1[p] = th1[p] + (w1[p] * DT) / S;
        th2[p] = th2[p] + (w2[p] * DT) / S;
        phi1[p] = phi1[p] + (w_phi1[p] * DT) / S;
        phi2[p] = phi2[p] + (w_phi2[p] * DT) / S;

        // The polar sines and cosines are from before this tick's update; the azimuths are updated
        unit(s1_arr[p], c1_arr[p], phi1[p]);
        I64 x1 = (L * ux) / S;
        I64 y1 = (L * uy) / S;
        I64 z1 = (L * uz) / S;
        unit(s2_arr[p], c2_arr[p], phi2[p]);
        I64 x2 = x1 + (L * ux) / S;
        I64 y2 = y1 + (L * uy) / S;
        I64 z2 = z1 + (L * uz) / S;

        step_x2[p] = x2;
        step_y2[p] = y2;
        step_z2[p] = z2;

        digest[p] = (digest[p] *% 1_099_511_628_211) ^ (th1[p] as% U64) ^ ((phi1[p] as% U64) <<% 11) ^ ((th2[p] as% U64) <<% 22) ^ ((phi2[p] as% U64) <<% 33) ^ ((x2 as% U64) <<% 44) ^ ((y2 as% U64) <<% 27) ^ ((z2 as% U64) <<% 7);
    }
}

// Triplet simulation run for n steps with coupling parameter k (-100 to 100)
// Returns final tip divergence between Leader A (0) and Twin B (1) in um
I64 run_triplet(I64 steps, I64 coupling_k, inout U64[3] out_digests) {
    start();
    tether = coupling_k;
    for i in 0..steps {
        step();
    }
    for p in 0..3 {
        out_digests[p] = digest[p];
    }

    I64 dx = step_x2[1] - step_x2[0];
    I64 dy = step_y2[1] - step_y2[0];
    I64 dz = step_z2[1] - step_z2[0];
    return isqrt(dx * dx + dy * dy + dz * dz);
}

// The page drives the program through these exports. A state slot is field * 3 + pendulum, with
// the fields in the order th1, th2, phi1, phi2, w1, w2, w_phi1, w_phi2
export I64 load(I64 slot, I64 value) {
    I64 p = slot % 3;
    I64 f = slot / 3;
    if (f == 0) { th1[p] = value; }
    if (f == 1) { th2[p] = value; }
    if (f == 2) { phi1[p] = value; }
    if (f == 3) { phi2[p] = value; }
    if (f == 4) { w1[p] = value; }
    if (f == 5) { w2[p] = value; }
    if (f == 6) { w_phi1[p] = value; }
    if (f == 7) { w_phi2[p] = value; }
    return value;
}

export I64 read(I64 slot) {
    I64 p = slot % 3;
    I64 f = slot / 3;
    if (f == 0) { return th1[p]; }
    if (f == 1) { return th2[p]; }
    if (f == 2) { return phi1[p]; }
    if (f == 3) { return phi2[p]; }
    if (f == 4) { return w1[p]; }
    if (f == 5) { return w2[p]; }
    if (f == 6) { return w_phi1[p]; }
    return w_phi2[p];
}

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

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

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

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

// Fills places from the current angles and returns the distance in um from the leader's tip to
// twin B's
export I64 locate() {
    for p in 0..3 {
        unit(isin(th1[p]), icos(th1[p]), phi1[p]);
        I64 x1 = (L * ux) / S;
        I64 y1 = (L * uy) / S;
        I64 z1 = (L * uz) / S;
        unit(isin(th2[p]), icos(th2[p]), phi2[p]);
        places[6 * p] = x1;
        places[6 * p + 1] = y1;
        places[6 * p + 2] = z1;
        places[6 * p + 3] = x1 + (L * ux) / S;
        places[6 * p + 4] = y1 + (L * uy) / S;
        places[6 * p + 5] = z1 + (L * uz) / S;
    }
    return spread(1);
}

export I64 place(I64 slot) {
    return places[slot];
}

// Distance in um from the leader's tip to the tip of pendulum p, from places
export I64 spread(I64 p) {
    I64 dx = places[6 * p + 3] - places[3];
    I64 dy = places[6 * p + 4] - places[4];
    I64 dz = places[6 * p + 5] - places[5];
    return isqrt(dx * dx + dy * dy + dz * dz);
}

void main() {
    U64[3] digests;
    I64 dist = run_triplet(500, 0, digests);
    "after 500 steps (5.0 seconds):\n";
    "pendulum 0: digest={digests[0]}\n";
    "pendulum 1: digest={digests[1]}\n";
    "pendulum 2: digest={digests[2]}\n";
    "divergence A to B: {dist} um\n";
}

test "zero force preserves bit-for-bit cint digests at step 500" {
    U64[3] digests;
    _ = run_triplet(500, 0, digests);
    assert(digests[0] == 14_995_127_728_425_866_399);
    assert(digests[1] == 8_897_290_559_184_475_167);
    assert(digests[2] == 11_575_834_901_304_267_551);
}

test "positive consensus tether reduces chaotic divergence" {
    U64[3] d_free;
    I64 dist_free = run_triplet(500, 0, d_free);

    U64[3] d_bound;
    I64 dist_bound = run_triplet(500, 50, d_bound);

    // positive coupling strictly collapses separation
    assert(dist_bound < dist_free);
}

test "negative coupling amplifies chaotic divergence" {
    U64[3] d_free;
    I64 dist_free = run_triplet(500, 0, d_free);

    U64[3] d_repulse;
    I64 dist_repulse = run_triplet(500, -50, d_repulse);

    // negative coupling pushes trajectories apart faster
    assert(dist_repulse > dist_free);
}

I64 safe_isqrt(I64 val) {
    return isqrt(val);
}

I64 checked_kinetic_energy(I64 w) {
    return w * w;
}

test "negative radicand triggers domain fault" expect_fault E_DOMAIN {
    I64 invalid_len_sq = -100;
    _ = safe_isqrt(invalid_len_sq);
}

test "supercritical velocity triggers overflow fault" expect_fault E_OVERFLOW {
    I64 runaway_w = 4_000_000_000;
    _ = checked_kinetic_energy(runaway_w);
}

What the Page Runs

The page runs pendulum.c, the C that cintc 0.2.0 emits from pendulum.ci, compiled to 64-bit WebAssembly in a background worker. A browser without 64-bit WebAssembly runs the same arithmetic rewritten in JavaScript integers, and the two agree on every value of 3,500 ticks at five tether settings.

Sines come from a seventh-order series and lengths from an integer square root, in millionths.
Every tick folds each pendulum's angles and tip position into a 64-bit digest. At tick 500 the page shows the digests that cint run prints.
Above 10 rad/s an angular speed loses 0.5% a tick, and it never passes 12 rad/s. The program and the page both apply this governor.

The Consensus Tether

The tether acts on each of the four angles of the two twins, toward the leader.

Pulling, the tether is a spring toward the leader. At full strength it is as stiff as gravity's own restoring pull, and it is damped so it settles without overshoot.
Pushing, the tether is a weak spring away from the leader with negative damping. The separation grows slowly until the speed governor holds it.
A change takes hold at the tick on screen, as the slider moves, and every earlier tick stays as it was. Copy link records each change, so the link it copies replays the same run.

Tip distance between the leader and the first twin after 500 ticks, in micrometers, from `cint run`
Tether	Push 100	Push 50	Off	Pull 25	Pull 50	Pull 100
Distance	416,045	11,402	1,519	253	63	16

The Receipt

outcome value
run.program pendulum.ci
run.steps 500
fuel-consumed 41509
stdout-sha256 ebe8373adc18ed8e0b6e6246d295bc15900aec02f393e777011c7977c5ae7683
tests 5/5

The reference interpreter and cintc 0.2.0 print the same bytes. The file's tests pin the three digests at tick 500 and check that the tether pulls and pushes. Two more stop on a negative square root and on an overflowing speed.

The emitted C is the output of cint emit-c on the same file.