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: the parsed float64 values are ==, in all four languages, at every point — no tolerance, no "agrees to six decimal places." Each repo replays the same golden files with an exact comparison: go-pflow's stochastic/portable_parity_test.go, pflow-rs's tests/ssa_parity.rs, pflow-xyz's public/petri-ssa_test.ts and pflow-jl's test/test_ssa.jl.
That 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. Getting there exposed a sign error that had been shipping unnoticed for months in three languages' ODE solvers, and two canonicalization bugs in the new Julia port that none of the existing fixtures caught.
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 (Go's copy of the vector is in stochastic/portable_test.go) — 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 Go on a measurable fraction of inputs. The sweep where I found glibc and Go's pure-Go math.log disagreeing on roughly 7% of uniformly-sampled doubles isn't checked in, but the ordinary case that made the point is: log(3.0) is pinned in portable_test.go at 1.0986122886681096, and the test's comment records that glibc returns 1.0986122886681098. Even Go's two implementations only agree "on every normal input tried" — 10,000 random ones in TestPlogMatchesMathLog — and the amd64 assembly returns the wrong value for the smallest subnormal. "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: 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 sqrt already comes free.
The engineering discipline mattered as much as the math, and it went in a specific order. The first goldens came from a from-scratch Python transcription of the spec. go-pflow's stochastic package, already the reference SSA engine, then got the portable path and a generator, cmd/ssa-goldens, and its output had to agree with the Python files on every parsed double before it replaced them — 1,785 doubles over four nets, zero differences, per pflow-jl's test/testdata/ssa/README.md. If the spec's own two implementations can't agree with each other, no downstream port stands a chance. The generator writes every recorded time, every place's mean and standard deviation, and the final marking to JSON for five nets: a linear chain, an SIR model, a dimerization reaction, a net exercising read arcs and inhibitors and self-loops, and a stripped-down coffee shop with combinatorial rate laws. A sixth, timed.json, came later, when the engine gained delayed transitions. Only once the goldens were fixed 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 == (or the language's equivalent, !== in the Deno test, === in Julia), 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 into that README, and the parity test's own doc comment says the same thing shorter — "a difference is a finding about the engine, not about the file" — 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 estimate — Tsit5's Bhat vector, the 5th-order weights minus the embedded 4th-order ones — have to sum to exactly zero, or even a constant derivative produces a nonzero error estimate 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 — 1.7k, 14k, 119k, 994k steps at reltol 1e-3 through 1e-6 on a linear chain — 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 about a million steps on that chain; after the fix the same solve hits the Dtmax cap on every net in the consistency test, about 0.1 ms each. 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 why nobody noticed.
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, using an AccurateOptions solve as the reference. An ODE reference that needed a million steps on a net that small was a smell, not a failure, and chasing it led straight to the sign. (The test's comment still records the before and after.) The fix is a one-character diff, dated 2026-09-02 in solver/tsit5.go. The test that pins it, solver/tsit5_test.go — coefficients sum to zero, the embedded weights satisfy the 4th-order quadrature conditions and not the 5th, the estimate vanishes 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, carried the exact same +1/66. The likely root cause is notation: the reference tableau for Bhat is easy to read as "the embedded weights" when it's actually "the weights minus the embedded weights," and Tsitouras gives the embedded last weight as 1/66 with the 5th-order one at 0, so the difference is -1/66. Whether the ports inherited the sign from Go or each misread the table the same way, the code no longer says. The fix went into all three, each with its own pin: pflow-xyz's public/petri-solver_test.ts asserts the sum and the -1/66, and pflow-rs's ode.rs carries a step-growth test. One canonical implementation with checked ports means a fix travels the same road the bug did.
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 (we found no Julia package that does RDF canonicalization at all, so 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 in the port 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. One was a misreading of hash_related_blank_node's position argument, which means the position of the related node in the quad, not of the node being hashed; the spec prose reads ambiguously enough that the reverse looks plausible, and it produces a valid canonicalization that simply isn't go-pflow's. The other was that every issued blank-node identifier used as hash input needs go-pflow's literal "_:" sigil prepended, not just its bare label — json-gold bakes the sigil into the identifier string, and the Julia port only added it at serialization.
What made these bugs sneaky rather than obvious: URDNA2015's "Hash N-Degree Quads" tie-breaking machinery — the part both bugs lived in — only matters when a document has multiple blank nodes that can't be distinguished by any simpler means, which in a net document means things like several arcs sharing a weight or several places sharing a capacity, each blank node distinguishable only by what owns it. The five pflow-xyz fixtures that had anchored Go/JS CID parity for months, net-a and net-b among them, did not expose either bug. Two new fixtures written for the port did: tie1.jsonld puts three arcs on the same weight: [1], and tie2.jsonld adds three places sharing the same capacity and initial, tripling the tied groups. tie2.jsonld is the one that found both bugs. pflow-jl's test/test_cid.jl now checks all seven fixtures — the five copied from pflow-xyz plus the two tie fixtures — against the content IDs Go's seal.SealJSONLD produces. The CID is a sha256 over the canonical N-Quads, so a match means the N-Quads agreed byte for byte; the test compares CIDs, not the N-Quads text itself. A future implementation of this same algorithm in the next language doesn't get to skip the hard fixture either. The fixture provenance and both bugs are written up in test/testdata/cid/README.md.
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. stochastic.Solve takes one metamodel.Model and dispatches on Options.Method to three engines: the SSA for discrete stochastic sample paths (the default), solver for the continuous ODE relaxation, 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. Along the way it turned up a fifth-order solver quietly running as a first-order one for months — caught by the SSA-versus-ODE consistency gate rather than by the bit comparison itself — and two canonicalization bugs that five existing fixtures had never triggered. Neither would have surfaced from code review; both surfaced from demanding that two things which should agree actually did.
When this post first went up, the chemical Langevin SDE — the third leg of the dispatch — was outside the byte-exact contract. Since 2026-09-10 it is inside it, in the same Go-generates, everyone-replays form: go-pflow's cmd/sde-goldens writes stochastic/testdata/sde/ for the same five nets and options as the SSA set, and pflow-rs (tests/sde_parity.rs), pflow-xyz (parity/sde/, replayed by public/petri-sde_test.ts) and pflow-jl (test/test_sde.jl) replay byte-identical copies with exact comparison. Two of the five, gates and coffeeshop, carry no series at all: the SDE refuses a net gated by read arcs, inhibitors, non-kinetic inputs or a reachable capacity, so those goldens pin the refusal — diverged, reason, caveats — and matching them meant rewriting the Rust and JavaScript refusal messages to Go's exact wording. What the SDE path still lacks is what the SSA path had from the start: an independent oracle. Its goldens come from Go alone, with no Python cross-check, so agreement shows the ports match Go, not that Go matches a written spec.
The code cites the portable SSA spec as ssa-spec.md, by section number, but that document is not in the go-pflow repository alongside it. What a reviewer can check line by line today is stochastic/portable.go (the PRNG, the ported logarithm, and the arithmetic-order rules in its comments), stochastic/portable_test.go (the spec's seeding and logarithm vectors), and stochastic/testdata/README.md (the fixture table, the hashes, and the delayed-transition rule the spec predates).
Revised 2026-11-02 against the repos (tests read, not re-run):
== in all four languages; the missing piece is now only an independent oracle.cmd/ssa-goldens was checked against them and replaced them, not the other way round. There are now six SSA fixtures, not five (timed.json).position bug is a misreading of which node the argument refers to, not an argument-order mistake; the fixture count is five from pflow-xyz plus two tie fixtures; and the Julia test compares CIDs, not N-Quads text.stochastic.Solve, not metamodel.Model; the "never fixed by regenerating" sentence is in testdata/README.md, not the test file; the 1e-5 tolerance loosening is not recorded in the code, so it was replaced with the step counts the test does record; the 7% figure for glibc is flagged as an unrecorded measurement; ssa-spec.md is not in the repo.Related: JSON-LD as Declarative Infrastructure · Small Models > LLMs · The Category Settle