Take a Petri net — a linear decay chain, an SIR epidemic model, a dimerization reaction, a coffee shop with combinatorial rate laws on its ingredient arcs. Simulate it stochastically with Gillespie's direct method, seed 42, in Go. Now simulate the identical net with the identical seed in Rust. In JavaScript. In Julia. Compare every value in the output — the ensemble mean at every recorded time, the standard deviation, the final marking.
They match. Not "close within tolerance." Not "agrees to six decimal places." The parsed float64 bits are ==, in all four languages, at every point.
That sentence is the whole claim, and it is a stronger claim than this ecosystem has made before. The Go↔JS state-root parity that has anchored the pflow ecosystem since the beginning compares hashes — a canonicalized document either produces the same content ID or it doesn't, and a single differing byte anywhere in the pipeline is invisible until it changes the hash. This is a different, harder thing: a numerical simulation, run independently by four different language runtimes, agreeing on the actual floating-point values a random process produced. Randomness that isn't random across implementations. That's the trick, and getting there exposed two real bugs that had been shipping, unnoticed, for months.
A stochastic simulation touches three things that vary silently across languages and runtimes, and every one of them has to be nailed down or the outputs diverge on step one.
The random stream. math/rand in Go, rand in Rust's rand crate, Math.random() in JS, and Julia's default RNG are four different algorithms with four different internal states. "Same seed" across them is a category error — there's no shared meaning of "seed 42" to agree on. The fix is to stop using any of them: implement one PRNG, by hand, identically in all four languages. I picked SplitMix64 to seed xoshiro256**, both simple enough to write out as an unambiguous spec and fast enough not to matter. SplitMix64 turns one 64-bit seed into the four 64-bit words xoshiro256** needs as initial state:
x = S
for i in 0..3:
x = x + 0x9E3779B97F4A7C15
z = x
z = (z XOR (z >> 30)) * 0xBF58476D1CE4E5B9
z = (z XOR (z >> 27)) * 0x94D049BB133111EB
z = z XOR (z >> 31)
s[i] = z
Seed 42 pins down to a specific, checkable state — s[0] = 0xBDD732262FEB6E95, s[1] = 0x28EFE333B266F103, and so on — and every port's unit test asserts that exact vector before it's allowed to touch a Petri net. Get the seeding wrong and every downstream draw is wrong, but wrong in a way that still looks like a plausible simulation, which is the worst kind of bug to catch by eye.
The logarithm. Gillespie's algorithm draws a waiting time as -log(uniform). log is a transcendental function, and IEEE 754 does not require transcendentals to be correctly rounded the way it requires +, -, *, /, and sqrt to be. Go's own runtime math.Log on amd64 is hand-written assembly; on s390x it's a third, different implementation; glibc's log disagrees with both on a measurable fraction of inputs — I found glibc and Go's pure-Go math.log disagreeing on roughly 7% of uniformly-sampled doubles, including something as ordinary as log(3.0). "The runtime logs are close enough" is a comfortable belief and it's false. The fix, again, is to stop delegating: stochastic/portable.go carries an explicit, byte-for-byte port of Go's own pure-Go math.log — not the assembly, the portable fallback — and every language ports that function, not their own libm. Grep the file and the runtime logarithm is nowhere in it; the sampler that calls it lives elsewhere, deliberately, so a reviewer doesn't have to trust a comment.
Arithmetic order. Modern compilers are allowed to fuse a multiply and an add into a single fused-multiply-add instruction when it doesn't change the mathematical result — but FMA rounds once instead of twice, which does change the bit pattern. This happens silently on arm64 and not on amd64 for the exact same source line. The fix is unglamorous: every expression in the spec is written in the order it must evaluate, with explicit intermediate float64 conversions wherever the compiler would otherwise be free to contract it. No Kahan summation, no pairwise reduction, no SIMD — anything that changes the order of floating-point operations changes the bits, so nothing is allowed to.
That's the whole list, and it's shorter than it sounds because of one fact worth sitting with: sqrt needs none of this. IEEE 754 requires sqrt to be correctly rounded — the standard doesn't budge on it the way it does for log, sin, exp. So when the chemical Langevin SDE variant needed Gaussian noise, the sampler could use sqrt freely and it would already agree across every conformant runtime with zero porting work. I chose the Marsaglia polar method for exactly this reason — it needs sqrt and the already-pinned log, nothing else — over Box-Muller, which needs sin and cos and would have meant porting a second transcendental for a feature that didn't strictly need one. Three pinned primitives turned out to be enough for two solver types, because one of the four basic operations already comes free.
The engineering discipline mattered as much as the math, and it went in a specific order. Go generated the fixtures first, since it was already the reference SSA engine promoted from petri-pilot — cmd/ssa-goldens runs five of them (a linear chain, an SIR model, a dimerization reaction, a net exercising read arcs and inhibitors and self-loops, a stripped-down coffee shop with combinatorial rate laws) and writes every recorded time, every place's mean and standard deviation, and the final marking to JSON. Before any port touched those fixtures, a from-scratch Python transcription of the spec had to agree bit-for-bit with a separate Go cross-check program — if the spec's own two implementations can't agree with each other, no downstream port stands a chance. Only then did Rust, JavaScript and Julia each implement the SplitMix64/xoshiro256**/portable-log contract from the written specification alone, without looking at each other's code; agreement on every double, once all three were done, is evidence the spec said enough, not evidence that everyone quietly copied the same mistake.
Every comparison in the resulting tests is ==, never a tolerance — portable_parity_test.go in Go, and its counterpart in each other language, replay every golden file and assert bit-equality on every recorded value, because a tolerance would let a real divergence hide under "close enough," which is the bug class this exercise exists to catch. And a parity test that would still pass after a real regression isn't safety, it's the appearance of it: flipping the last bit of one golden value and watching the suite go red, on purpose, once, is how you find out the gate actually bites.
The goldens themselves carry their own provenance chain: each JSON file's _comment names the exact go-pflow commit that generated it, stochastic/testdata/README.md carries a table of sha256 hashes for every file, and the rule for all three downstream repos is cp plus sha256sum verification — never regenerate independently, because a divergent regeneration is indistinguishable from a real bug until someone diffs the bytes. "A difference means the engine's sample path changed; it is never fixed by regenerating" is written directly into the test file's own doc comment, so the rule survives whoever reads the code next.
This project also carries an ordinary differential equation solver — Tsit5, a fifth-order Runge-Kutta method used whenever a net is analyzed as a continuous relaxation instead of a discrete stochastic process. Fifth-order methods estimate their own error using an embedded lower-order solution: compute both a 5th-order and a 4th-order answer from the same function evaluations, and the difference between them estimates how wrong the 5th-order answer might be. The coefficients that produce that 4th-order embedded answer — Tsit5's Bhat vector — have to sum to exactly zero when the underlying function is constant, or the estimator has a built-in bias that scales with step size regardless of how smooth the solution actually is.
Go's Bhat vector had +1/66 in its last slot. It should have been -1/66. The seven coefficients summed to 2/66 instead of 0.
The consequence wasn't a crash or an obviously wrong answer — it was quietly, uniformly bad step control. An error estimate that's supposed to shrink with the fifth power of step size was instead shrinking with the first power, so the adaptive solver treated a perfectly smooth trajectory as if it needed constant refinement. Accepted step counts grew roughly 10x for every decade you tightened the relative tolerance, instead of the roughly 1.6x a genuine fifth-order estimate gives. AccurateOptions — the "give me a trustworthy answer" preset at reltol = 1e-6 — was costing on the order of a million steps on nets simple enough that a hundred should have sufficed. It worked. It was just burning three or four orders of magnitude more compute than it needed to, silently, and the answers it eventually produced were fine — which is exactly why nobody noticed for as long as they didn't.
It surfaced because of the discipline the SSA work had just put in place: stochastic/consistency_test.go checks that the SSA's long-run ensemble average converges to the ODE relaxation's trajectory on the same net, and that test needed its tolerance relaxed to 1e-5 to pass — a smell, not a failure. Chasing why an ODE solver needed loosening to agree with a stochastic simulation that should have been the noisier of the two led straight to the sign. The fix is a one-character diff. The test that pins it — coefficients sum to zero, the estimate vanishes identically for a constant derivative, accepted-step growth stays under 3x per decade of tolerance from 1e-3 to 1e-6 — is 129 lines, because a bug that shipped for months needs more than a spot check to make sure it can't come back quietly.
The two ports that mattered most, JavaScript and Rust, had independently carried the exact same +1/66. Not a coincidence — both were transcribed from the same Butcher tableau reference, and both transcriptions made the identical sign error, because the reference notation for Bhat is easy to misread as "the embedded weights" when it's actually "the weights minus the embedded weights." Three languages, one bug, one shared root cause. Fixing it once and re-deriving the port in the other two, rather than three people independently rediscovering it, is the entire reason this ecosystem insists on one canonical implementation with byte-checked ports instead of three teams each writing their own Tsit5 from the paper.
The second bug lived somewhere stranger: not in a numerical method but in RDF Dataset Canonicalization — URDNA2015, the algorithm that turns a JSON-LD graph into a deterministic, order-independent string so that two semantically identical documents hash to the same content ID regardless of what order their fields happened to be written in. Porting this to Julia (no existing package does RDF canonicalization at all — it had to be written from scratch, alongside a hand-rolled JSON-LD expander scoped to exactly this schema's context) meant reproducing go-pflow's actual output, not just the URDNA2015 spec's prose, because a specification this fiddly has room for a compliant implementation to still disagree with a specific one on edge cases the spec text doesn't make unambiguous.
Two bugs turned up, and neither was findable by reading the spec harder. Both were found by instrumenting go-pflow's own piprate/json-gold dependency with debug prints and diffing what it actually did against what the Julia port did on the same input: an argument ordering mistake in hash_related_blank_node's position parameter, and a missing detail in how every issued blank-node identifier needs to be hashed together with go-pflow's literal "_:" sigil prefix, not just its numeric suffix.
Here's the part that made these bugs sneaky rather than obvious: URDNA2015's "Hash N-Degree Quads" tie-breaking machinery — the part both bugs lived in — only activates when a document has multiple blank nodes that can't be distinguished by any simpler means, which in practice means multiple arcs of the same weight connecting the same structural pattern. Every simple test fixture, every net with distinguishable arc weights, sails through canonicalization without ever exercising that code path. The bug was invisible on net-a and net-b, the two fixtures that had anchored Go/JS CID parity for months — both are hand-built nets where the arcs happen not to tie. It took a deliberately adversarial fixture — tie2.jsonld, built specifically to force several same-weight arcs into contention — to make the divergence show up at all. Verification against seven existing parity fixtures plus two new tie-breaking ones now checks not just the resulting content ID but the full canonical N-Quads text, byte for byte, so a future implementation of this same algorithm in the next language doesn't get to skip the hard fixture either.
The lesson generalizes past RDF: a fixture set built entirely from "the models I already had" will systematically miss whatever code path only your hardest inputs exercise, and a spec's prose is not a substitute for diffing against a real, instrumented, independently-generated reference — go-pflow's own dependency, questioned with print statements, not trusted on the strength of matching the written algorithm.
The pflow ecosystem's whole premise is that a Petri net is a declarative schema — places, transitions, arcs, rates — and everything downstream should be derived from that declaration, not independently reimplemented and hoped into agreement. metamodel.Model dispatches to three different engines from one declared net: solver for the continuous ODE relaxation, the SSA promoted from petri-pilot for discrete stochastic sample paths, and the chemical Langevin SDE for the continuous-noise middle ground between them. That dispatch existing in one place, off one input type, is half the thesis.
The other half is what this post is actually about: the same declared net, run through four independently written language implementations, has to produce the same numbers — not similar numbers, the same numbers — or "declare once, derive everywhere" is marketing rather than a property you can check. Four implementations agreeing to the last bit of a float64, on a genuinely random process, is that property made checkable. It found a fifth-order solver quietly running as a first-order one for months, and a canonicalization bug two independent test fixtures had never once triggered. Neither would have surfaced from code review, and both surfaced from insisting on bytes instead of vibes.
The chemical Langevin SDE — the third leg of the dispatch trio — isn't part of this byte-exact contract yet. There's no external reference to generate cross-language goldens from independent of Go's own implementation, the way SSA has ssa-spec.md and an independent Python cross-check. What exists today is each of the four languages checking its own Gaussian sampler bit-exact against Go's seed-42 output, and each running the same three-way consistency gate (SDE's mean tracks the ODE, SDE's variance at scale tracks the SSA's) against its own SSA rather than a shared reference. Promoting that to the same standard as the SSA path — an independent spec, an independent oracle, byte-identical goldens copied and hash-verified into all four repos — is the obvious next capability to close, and the playbook above is exactly how to close it.
The full spec — every PRNG state transition, the portable logarithm's exact port, and the arithmetic-order rules a reviewer can check line by line — lives in go-pflow's ssa-spec.md.
Related: JSON-LD as Declarative Infrastructure · Small Models > LLMs · The Category Settle