I did a bit of research earlier in my life[0] to study the handling of loops and without using Monte Carlo simulation. The result was actually workable if incredibly resource intensive to the point of being impractical. If I had chosen to do it again, I might’ve accepted using Monte Carlo simulations while still supporting loops.
[0]: https://github.com/kccqzy/probabilistic-program-inference/bl... Shameless self promotion I know! I put quite a bit of effort into that README and the code.
I remember encountering this idea written in a book written by Ed Catmull of Pixar fame (can't find the title sorry, but it was written in the 80s), but generally comes from signal processing as a way of avoiding aliasing artifacts..
The core idea is to make programming, which is a discrete and discontinuous domain, into a well-behaved band limited signal. Otherwise you get aliasing (or jaggies), which can happen even INSIDE a surface, if the shader's like that.
The code idea for this is the step function which is the integral of the dirac delta. step(x) returns 1 for all x >0 and 0 otherwise. Step is not a well-behaved function in the sense, that it changes infinitely quickly at x=0. But once we know what we want, we can replace it with something like that, that's well behaved.
Consider the example pseudocode
color = x> 5? green:blue;
can be rewritten as
color = blue + step(x-5)(green-blue)With the two being equivalent.
Now if we put the code into a shader, we get jaggies. So to combat the value changing infinitely fast, we go for a function that's like step, but changes smoothly* from 0 to 1 around x=0. Enter smoothstep: color = blue + smoothstep(x-(5+EPSILION),(x-EPSILON), x)*(green-blue)
And so we defined a 'transition zone' of +-EPSILON(an arbitrary number). While any smooth function can work, smoothstep is chosen because it has a smooth first and second derivative (meaning even if you want to get the rate of change, something that often pops up in computer graphics, the result will be still well behaved).
Pixar's Renderman shading language (which is remarkably similar to GLSL/HLSL/C), used to do this automatically for you. Essentially it could take arbitrary code peppered with if statements, and turn it into a continuous function.
Which is kinda cool imo.
It's also a cool trick in the age of AI. Since you have a function that's well-behaved, you can do things like gradient descent to train an AI to synthetize a function for you. You can even say, that you don't need exact results, you can accept some error.
In this case your program optimization problem can be reframed from doing idempotent transformations on the list of instructions, to getting a program that generates a target function whose error is no greater than some (mathematical) reference function.
The engine is Rust, the JIT is built on Cranelift, there is also a WASM backend so everything runs in the browser too.
Full disclosure, I could only finish it now because of AI agents. In my experience they are amazing at the runtime and the numerical code, but pretty bad at language design, so I kept that part for myself.
It's a toy language. Ask me anything!
Seems worth an investigation and maybe mention on the article.
applied to the step function, you would get a smooth cutoff function
https://en.wikipedia.org/wiki/Mollifier#Smooth_cutoff_functi...
this is also related somewhat to the notion of differentiable programming. RELU is (roughly) the same as x * step(x). In differentiable programming one can replace it with smooth approximations, cf "softplus"
https://arxiv.org/pdf/2403.14606
That book also has a chapter on control flow, which is very similar to what you're talking about.
Unrolling an if statement into x = b (result of one branch) + (1-b) (result of the other branch) is also incredibly common in cryptography. If `b` is a "secret" variable, an if statement may leak the value of it via the branch predictor/speculative execution. The way around this is to compute both branches, and then select them with the above arithmetic expression. This mostly works, though compilers are tediously smart, and so one often has to be careful how with how you precisely do it.
Stan and PyMC beat Noise at the thing they’re built for, fitting a posterior to lots of continuous data with their HMC/NUTS samplers, and NumPy beats it at raw array crunching. Conditioning in Noise is rejection-based, so it works great for a handful of discrete observations but becomes useless for ten thousand continuous measurements, and there is no stateful simulation yet (no Markov chains yet). Where Noise wins when you have a probability question and you wanna know the answer without much hassle.
So use Noise for the whiteboard stage of a problem, when you want to run the math you just wrote, and move to Stan or PyMC when you need a real posterior, or to NumPy and JAX when you need to go to production.
Did you ever look at symbolic or exact operators instead of purely monte-Carlo?
I remember reading the paper for distr, they used a mix of Fourier transforms for convolution and symbolic reduction to build a probabilistic computing library in R. I attempted building a small python library for this, but for my problems the CLT ended up sufficient to approximate the results faster, so I went with that.
You may enjoy reading the paper, it’s not groundbreaking but is a nice presentation of relationships/operations. Maybe it’ll inspire some features for you.
1. It's obvious that you are a big fan of Cranelift. I'd be interested to hear more about your experience in practical terms. For example could you share any insight about use cases where it is best suited, and where it might be better to look elsewhere? Did you hit any pain points? What was its killer feature for NoiseLang?
2. You wrote: "My favorite trick is in the RNG. Generating random numbers is a serial dependency chain, so instead of fighting that, the kernel runs four independent streams at once and lets the out-of-order core overlap them. This trick ended up beating a hand-written SIMD kernel!" What does this mean exactly? you just ran a scalar kernel 4-wide using SIMD instructions? or you interleaved 4 scalar copies of the same algorithm? did you generate the code or just duplicate the streams by hand?
I know MCMC isn’t your goal, but seems like this could be used for ABC-MCMC (as is?)
Would also be nice to have an option to plot using a KDE vs histograms.
(Also your FM example seems to be technically PM)
There are multiple versions of the list. The authoritative site appears to be https://github.com/hagezi/dns-blocklists, and making a fairly random choice, I used the "medium" version of the "Threat intelligence feed", and specifically the one marked "Link" for AdBlock. That took me to https://cdn.jsdelivr.net/gh/hagezi/dns-blocklists@latest/adb..., and manualmeida.dev does appear in that list.
The software I'm using is Little Snitch on a Mac, but since the entry is in the list, that's not the problem.
What noiselang does it parse, convert to a execution graph then carefully use the Cranelift api to describe the program, and under the hood, it will find different opetimization and generate byte code that runs directly in the host CPU.
The reason for it to be killer is that it allows NoiseLang to run as native speeds, with very little compiler/optimization work. it's a very simple repo.
For the RNG, this was a discovery myself, when profiling, i found the many benchmarks were limited by the speed of the RNG itself, ie, if i could genenrate random numbers faster, the simulation would be faster. xoshiro's next number is computed from its current state. So to get number N+1 you must have finished number N. It's a chain: A → B → C → D. Your CPU can run maybe 6 integer operations per cycle, but a chain only ever offers it one to run. Five of the six lanes sit empty.
I tried to use SIMD to speed this up, but still hit the limit, even if using SIMD, it still had to wait for the next number, a massive speed up came from realizing that i can keep four independent xoshiro256++ states and emits four samples per loop iteration, i += 4. Since the four state-update chains share no registers, the out-of-order core issues them in the same cycles instead of stalling on one serial chain.
SIMD gives you more work per instruction. But I wasn't short on work, I was waiting. A 2-wide xoshiro still needs state N before it can compute state N+1, so the chain is the same length and I wait at every link, I just get two numbers per link instead of one.
And each link costs more. xoshiro rotates a 64-bit word every round, and NEON has no 64-bit rotate, so that becomes three instructions instead of one. Twice the numbers, three times the wait.
Four streams wins because it leaves the chain alone. It just runs four of them at once, and the CPU was already idle enough to overlap them for free.
Fair! My thinking was that PM of a single tone signal (the one i use in the demo is equivalent to FM, but shifted a bit). And implementing real FM for decoding is a lot more noisy, but I will add some callout in the article.
Truth be told, you motivated me to write the exact FM with the differenciation, maybe. Could be interesting to simulate PM vs FM for non single tone signals, to see how FM does even better!
For the scope of the language it never even comes up, because Noise is a simulator, it does not evaluate densities, it draws samples.
The point is that every value goes through the same operators. Add them, compare them, pass them to a function, put one in the condition of an if. You can even use a random variable to define another random variable:
bias ~ unif(0, 1) flips ~[10] bernoulli(bias) // bernoulli just took a distribution where a number normally goes.
and in if-stataments:
DistributionC = if DistributionA < DistributionB { 0 } else { 1 }
But you right, dirac only applies to continuous functions, in Noise is only refers to the dirac measure. I found this article a fun/nerd to make my point that everything "acts" as a distribution from the DX perspective, but under the hood 5 is just 5.
And a constant collapses back to a plain integer in the graph anyway, so 5 costs nothing.
Since the advent of AI I have been delighted by the number of old project ideas I've been able to execute. There are just too many ideas to implement them all by yourself - but AIs don't seem to mind in the least.
It's a brave new world.
[1] And feels philosophically like the unification in the underlying maths between discrete and continuous probability that you get when you apply measure theory
like this:
D ~ unif_int(1, 6); Print("P(rolled a 6 | rolled > 3) =", P(D == 6 | D > 3));
or:
loss ~ unif(0, 1000); claim = if loss > 200 { loss - 200 } else { 0 }; p = P(claim > 0); Print("P(insurer pays a claim) =", p)
notice that "claim" is also a random variable! result of a if expression
N = 10**6
D = np.random.randint(1, 6, N)
print("P(rolled a 6 | rolled > 3) =", ((D == 6) | (D > 3)).mean())
loss = np.random.uniform(0, 1000, N)
claim = np.where(loss > 200, loss - 200, 0)
p = (claim > 0).mean()
print("P(insurer pays a claim) =", p)
It’s concise enough that people generally wouldn’t bother writing a library. Unless they really want their custom syntax, then perhaps they write a parser.During my telecommunications degree I took a course on signals and noise, I spent a lot of evenings writing probability by hand: expectations, variances, the odds of two random variables landing in some region. It always sucked, when I tried to run it on a computer, so much boilerplate.
That wish became NoiseLang. I started it about nine years ago, however, I never finished it. Only recently, I brought it back thanks to AI tools and something far more ambitious than what I could have built alone the first time.
The whole language hangs on one idea, that every value is a probability distribution. A plain number is a Dirac spike, a distribution with all its weight on a single value. Since constants and random variables are the same kind of object, every operator in the language maps distributions to distributions.
A name always refers to one fixed node, the same way X is the same X across a whole page of math. So X + X is 2X and X - X is exactly 0,. If you want variable independence you write separate draws, ie, using ~ multiple times, or ~[N] to draw N independent variables into a vector.
X ~ unif_int(1, 6)
Y ~ unif_int(1, 6)
X + Y # two independent dice, a real 2d6 distribution
Nothing runs until you ask for some results, for example, P(X + Y < 10), at that moment it forces the runtime to run millions of simulations (across all cores, if available) and return an estimate with a standard error attached.
Bday = unif_int(1, 365)
days ~[23] Bday # 23 people in a room
P(has_duplicates(days)) # the birthday paradox, about 0.507
This is much easier to watch than to describe, so I built some cool demos!
A = 5D1 ~ unif_int(1, 6)D2 ~ unif_int(1, 6)S = A + D1 + D2E(S)
draws 0 A ≈ —
Figure 1. The same histogram through four steps: a constant is a single spike, the tilde spreads it into a die, adding two dice curves it toward a bell, and a query reads the mean back.
A = 5 looks like a plain constant, and it is, but Noise sees it as a probability distribution where every draw lands on 5. Statisticians call this a Dirac delta, a single infinitely-thin spike.
D1 ~ unif_int(1, 6) binds a fair die, a discrete uniform where all six faces are equally likely. The spike fans out into six flat bars, so D1 is now uncertain, but you still write it like any other variable.
Add the constant and two independent dice with S = A + D1 + D2. The flat shapes convolve into a triangle peaked at 12, already curving toward a bell, and if you join a few more they sharpen into a true Gaussian. That's the central limit theorem showing up for free!
Nothing runs until you ask. E(S), the expected value, fires a Monte Carlo pass that averages millions of rolls into one number, and by the law of large numbers the running mean settles on 12. You wrote a few lines of math, and an expert-level kernel ran underneath.
The design was never the hard part, because a parser and a tree-walking interpreter for this language is a weekend of work. The problem was everything else, writing a efficient Monte Carlo runtime, instead of a naive interpreter, conditional bayesian inference, and more.
Current version is a compiler, a JIT (using the amazing Cranelift), a WASM backend, and a pile of careful numerical code, so for a cute-weekend project it stayed permanently out of reach.
At my day job and side projects, I am experimenting with the boundaries of what today’s AI agents can do. For example, I am also porting a game I built 15 years ago for iOs, in archaic Objective-C to a modern game engine (with relative success).
With NoiseLand, I realized, AI is great at building the JIT parts, the runtime parts, the numerical parts, but it sucks at coming up with good language design ideas, many times overriding existing language features for different purposes, or coming with with different syntax for non-orthogonal features.
Under the hood, ~ and the distribution constructors build an append-only DAG called the RvGraph. This graph is the single source of truth, which later are converted into three different code paths:
fallback
fallback
X ~ unif(1, 6)
RvGraph
batch interpreter
Cranelift JIT
WASM emitter
One shared module defines what the graph means, so the two code generators stay thin and cannot drift apart. Anything a backend can’t successfully compile falls back to the interpreter, and the results stay identical across backends and core counts. All tests run in all three code paths and compared to be bit-identical.
All the performance work is about one loop: draw a few million samples, evaluate the expression on each, and reduce the results. A handful of techniques carry most of it, while keeping the results deterministic (that was the hard part).
Kernel fusion keeps every intermediate value in registers, so an arithmetic-heavy expression stay on registers. The PRNG (xoshiro256++) compiles into the kernel, and the ln, sin, and cos become inline polynomial approximations, speeding up the kernel by a factor of 2.
My favorite trick is in the RNG. Generating random numbers is a serial dependency chain, so instead of fighting that, the kernel runs four independent streams at once and lets the out-of-order core overlap them. This trick ended up beating a hand-written SIMD kernel!
On my 14-core M4 Pro, a one-line P(...) sustains around 5.8 billion samples per second and scales about 9.6× from one core to all of them. Per core, the generated kernel runs within about 1.15× of hand-written Rust compiled by LLVM. The same fused loop, emitted as WASM, runs at roughly half to three-quarters of native speed inside V8.
X ~ rand::unif(-1, 1)Y ~ rand::unif(-1, 1) C = X^2 + Y^2 < 1 4 * P(C)
darts 0 inside 0 π ≈ —
Figure 4. Monte Carlo π. The circle covers π/4 of the square, so four times the teal share converges on π as the darts accumulate.
This is the same ~ as before, used twice. X ~ unif(-1, 1) and Y ~ unif(-1, 1) each draw a number anywhere in [-1, 1], and together they pick a random point in the square.
C = X² + Y² < 1 is true exactly when the point falls inside the unit circle. In Noise that comparison is itself a random variable, a Bernoulli that lands true or false on each draw.
Each dart is one draw of (X, Y), teal when it lands inside the circle and grey when it doesn't. They scatter evenly, because the uniform draws don't favor any point over another.
The circle covers π/4 of the square's area, so the teal share hovers around π/4, and that makes 4 * P(C) an estimate of π. The more darts you throw, the sharper the estimate gets.
NoiseLang is a toy language, you probably should not use it for anything serious, however I wish this language existed during my university days.
For a language nerd, it’s a small, static random-variable algebra with forward Monte Carlo, expression-based, rejection-based conditioning, language.
You might ask, how does it compare to NumPy or Stan? NumPy makes you write the simulation yourself, and Stan makes you declare a model and wait for a sampler. Noise lets you write the probability as math while running Monte Carlo under the hood to get the answer.
| Noise | NumPy | Stan | PyMC | |
|---|---|---|---|---|
| What it is | RV algebra + Monte Carlo | array numerics | declarative Bayesian modeling | Bayesian modeling in Python |
| Core abstraction | every value is a distribution | n-d arrays | data / parameters / model blocks | with Model() context |
| Getting an answer | P / E / Var / Q force sampling | you write the loop | compile → NUTS → summarize | pm.sample() → arviz |
| Inference | forward MC + rejection conditioning | none (you build it) | HMC / NUTS at scale | HMC / NUTS, VI, SMC |
| Setup to first answer | zero — type in a browser | import numpy | install toolchain, compile | install stack, build a context |
| Visualization | built into the language | matplotlib (separate) | bayesplot (external) | arviz (separate) |
| Runs in the browser | yes — WASM, no install | no | no | no |
Table 1. Where Noise sits. Stan and PyMC beat it at fitting a posterior to lots of data, and NumPy beats it at raw array crunching, but Noise gets you from a probability question to a visible answer faster than any of them.
Stan and PyMC beat Noise at the thing they’re built for, fitting a posterior to lots of continuous data with their HMC/NUTS samplers, and NumPy beats it at raw array crunching. Conditioning in Noise is rejection-based, so it works great for a handful of discrete observations but becomes useless for ten thousand continuous measurements, and there is no stateful simulation yet (no Markov chains yet). Where Noise wins when you have a probability question and you wanna know the answer without much hassle.
So use Noise for the whiteboard stage of a problem, when you want to run the math you just wrote, and move to Stan or PyMC when you need a real posterior, or to NumPy and JAX when you need to go to production.
Going back to my university days, there was this subject called “Señales Aleatorias Y Ruido”, which is a spanish translation of “Random Signals and Noise”.

The textbook, in all its glory.
This subject was the nightmare of many students, including myself, in fact, I failed it. Truth be told, I didn’t put enough effort into it during the first year, but it changed when I had the take the same subject again in the second year. The professor was great, and made the subject interesting. The things that blew my mind was how he could model why FM survives a noisy channel when AM doesn’t.
So, here is my tribute to the subject, a one-screen Noise program that models why FM survives a noisy channel when AM doesn’t.
msg = 0.3 * signal::sine(3); am_modulate(m) = 1 + m;fm_modulate(m, swing) = exp(i * swing * m); static ~ signal::noise_white_complex(0.4);rx_am = am_modulate(msg) + static;rx_fm = fm_modulate(msg, 3) + static; am_demodulate(c) = abs(c) - 1;fm_demodulate(c, swing) = arg(c) / swing;rec_am = am_demodulate(rx_am);rec_fm = fm_demodulate(rx_fm, 3); Print("AM error", E(mse(rec_am, msg)));Print("FM error", E(mse(rec_fm, msg)));
scroll through the steps →
Figure 3. AM vs FM end to end. The carrier is a phasor, AM writes the message into its length and FM into its angle, so the same static corrupts what AM reads while barely touching what FM reads.
In Noise a value can carry a whole signal, and here it's a gentle tone, a slow three-cycle sine. This is what we want to send through a noisy channel and get back intact.
The carrier is a complex number, a spinning arrow. AM's modulator puts the message in the arrow's length, so the tip slides in and out, crossing the unit circle as the message rises and falls.
FM puts the same message in the angle instead, so the tip rides along the unit circle while its length stays fixed. The information is identical, just encoded in rotation.
Now identical white noise hits both tips. It smears AM along exactly the radius it reads from, but it only nudges FM around the circle, barely changing the angle. You can see it in the shape of the two clouds.
Demodulate: read the length back for AM and the angle back for FM, laid over the original message. You can already see it, the AM trace comes back ragged while the FM one hugs the original.
mse averages the error over the waveform, and E averages that over many draws of the static. Same signal energy, same noise energy, yet FM comes back about 5× cleaner, and in the small-noise limit the advantage approaches swing² = 9. That is the bandwidth it spent on the angle paying off.
Everything above runs on @noiselang/core, thanks to the Rust engine compiled to WebAssembly.
npm install @noiselang/core
import { run } from "@noiselang/core";
const result = await run(`
X ~ rand::unif(-1, 1);
Y ~ rand::unif(-1, 1);
4 * P(X^2 + Y^2 < 1)
`);
console.log(result.value); // "3.1415…" — the last statement's value
console.log(result.output); // everything Print(...) emitted
run never throws, failures come back on result.error with a source span. There is also runWithIntrospection, the API behind the variable inspector at noiselang.com.
NoiseLang is playable in the browser at noiselang.com. Open it, type X ~ unif(-1, 1); Y ~ unif(-1, 1); 4 * P(X^2 + Y^2 < 1), and watch a few million draws estimate π from your browser tab.