This post has interactive code editors and animations on the original version where you can run and modify every Rust snippet in the browser.
Roots of Unity
Welcome back, fellow zk-enjoyer.
In the last post, we used polynomials as a data structure and verified them with a single number (Schwartz-Zippel, remember?). But there's a question that should be bugging you: "how do we do this efficiently?"
The answer is called "roots of unity." I promise the idea is simpler than the name suggests, and it's what makes ZK proofs better (about 50,000 times better, for a typical ZK circuit, terms and conditions may apply [TM]).
The Road Map
- Polynomials over finite fields
- Roots of unity <- you are here
- Execution trace
- Polynomial interpolation via FFT/NTT
- Applying constraints to polynomials
- Packing into a single polynomial
- Commitment scheme (starting with KZG)
- Proof AND Verification
Why Are There Two Ways to Store a Polynomial?
Before we get to roots of unity, you need to see the problem they solve.
You can represent a polynomial in two different ways.
Coefficient form is what you already know: store the coefficients as a vector. P(x) = 1 + 2x + 3x^2 becomes [1, 2, 3]. Adding two polynomials is cheap: add the vectors element by element, O(N). But multiplying them? Every coefficient of P has to multiply with every coefficient of Q. That's O(N^2).
Evaluation form stores what the polynomial outputs at N chosen points. Instead of [1, 2, 3], you store [P(x_0), P(x_1), ..., P(x_{n-1})]. Now multiplying two polynomials is trivial: just multiply the values at each point, O(N).
Think of it like a coin. Coefficient form and evaluation form are two views of the same polynomial. Each view makes a different operation cheap: coefficients for addition, evaluations for multiplication. The catch is flipping the coin.
Flipping from coefficients to points is evaluation: plug N points into the polynomial, O(N^2). Flipping back is interpolation: find the polynomial that passes through N given points, also O(N^2).
For N = 2^20 coefficients (a common ZK domain size), each flip costs (2^20)^2 = 2^40: about a trillion multiplications. Evaluation form is where the heavy operations happen, but the flips are the bottleneck.
Unless you pick the right evaluation points.
So What Are Roots of Unity?
For the 8 roots of unity, you start with a number omega. Then omega^0 = 1, omega^1, omega^2, ..., omega^7 gives you 8 distinct values, and omega^8 = 1 again.
On the complex plane, roots of unity sit on a circle and multiplying by omega is a visible rotation. That's the intuition.
In a finite field, there's no circle: omega is a single integer (a big one), and "multiply by omega" is just modular multiplication. Nothing rotates. But the algebra is the same: N distinct values that cycle back to 1, with identical symmetries. The circle helps you see the structure. The code below is how roots of unity actually look in ZK.
But in many zkSNARK frameworks and libraries you will see a thing named omega and operations called rotation. This makes no sense for integers. The name comes from the complex plane, where it does make sense, and that's why those operations and primitives are named that way.
use ark_bn254::Fr;
use ark_ff::FftField;
fn main() {
// The primitive 8th root of unity for BN254.
// "Primitive" means: omega^8 = 1,
// and omega^k != 1 for 0 < k < 8.
let omega = Fr::get_root_of_unity(8 as u64).expect("Field has this root");
println!("All 8th roots of unity for BN254:\n");
let mut current = Fr::from(1); // start at omega^0 = 1
for i in 0..=8 {
if i == 8 {
println!(" omega^{} = {:?} <-- back to 1!", i, current);
} else {
println!(" omega^{} = {:?}", i, current);
}
current *= omega;
}
}
Those numbers are huge (welcome to finite fields), but the punchline is at the bottom: omega^0 and omega^8 are the same value. The cycle is the same as in the complex plane, and it also has exactly 8 distinct values.
There's one more thing to see here. The "root" in "root of unity" means the same as in "root of an equation": a number that makes the equation equal zero. Every root of unity satisfies omega^8 = 1, or equivalently omega^8 - 1 = 0. Not just one omega, but every root:
(omega^0)^8 = 1
(omega^1)^8 = 1
(omega^2)^8 = 1
...
(omega^7)^8 = 1
That's the definition. Why? Because you can reorder the exponents:
(omega^2)^8 = (omega^8)^2 = 1^2 = 1
use ark_bn254::Fr;
use ark_ff::{FftField, Field};
fn main() {
let n: u64 = 8;
let omega = Fr::get_root_of_unity(n).expect("Field has this root");
let one = Fr::from(1);
println!("Every 8th root of unity, raised to the 8th power:\n");
for k in 0..n {
let omega_k = omega.pow([k]); // first: omega^k
let result = omega_k.pow([n]); // then: (omega^k)^8
println!(" (omega^{})^8 = {:?}", k, result);
assert_eq!(result, one);
}
println!("\nAll equal 1. They're all roots of z^8 - 1 = 0.");
}
This will matter later: the polynomial z^8 - 1 has all eight roots of unity as its roots, which means that no matter what omega you put into that equation, you will always get 0.
That fact gives us a way to check constraints across all roots at once, and of course we can use roots of unity as the x coordinate for encoding our data, instead of one by one (post #5).
How Does This Speed Things Up?
Take a polynomial P(x) and split it by even and odd coefficients:
P(x) = P_even(x^2) + x * P_odd(x^2)
For example, take a degree-7 polynomial with 8 coefficients:
P(x) = 3 + 1x + 4x^2 + 1x^3 + 5x^4 + 9x^5 + 2x^6 + 6x^7
Pull out the even-index coefficients (positions 0, 2, 4, 6) and the odd-index coefficients (positions 1, 3, 5, 7):
P_even(y) = 3 + 4y + 5y^2 + 2y^3
P_odd(y) = 1 + 1y + 9y^2 + 6y^3
Plug in y = x^2 and recombine:
P(x) = (3 + 4x^2 + 5x^4 + 2x^6) + x * (1 + 1x^2 + 9x^4 + 6x^6)
One degree-7 polynomial just became two degree-3 polynomials. Each of those can be split the same way. Split P_even by its even/odd coefficients:
P_even_even(y) = 3 + 5y
P_even_odd(y) = 4 + 2y
And split P_odd the same way:
P_odd_even(y) = 1 + 9y
P_odd_odd(y) = 1 + 6y
Four degree-1 polynomials. Split once more, each has two coefficients, so the split produces eight degree-0 polynomials (just the constants: 3, 5, 4, 2, 1, 9, 1, 6).
Eight constants. "Evaluation" of a degree-0 polynomial is just returning the number. That's the base case.
The full cascade: 8 -> 4 -> 2 -> 1 coefficients per sub-polynomial, log2(8) = 3 levels of splitting.
Now the computation runs in reverse, from the bottom up. At each level, the butterfly combines two values a and b into a + omega^k * b and a - omega^k * b.
The omega^k multiplier is called the twiddle factor: it's just a power of omega that changes at each level. At the first level it's 1, so the butterfly simplifies to plain a + b and a - b:
Level 1: evaluate degree-1 polynomials at 2nd roots (omega^k = 1)
P_even_even(omega^0) = 3 + 1*5 = 8, P_even_even(omega^4) = 3 - 1*5 = -2
P_even_odd(omega^0) = 4 + 1*2 = 6, P_even_odd(omega^4) = 4 - 1*2 = 2
P_odd_even(omega^0) = 1 + 1*9 = 10, P_odd_even(omega^4) = 1 - 1*9 = -8
P_odd_odd(omega^0) = 1 + 1*6 = 7, P_odd_odd(omega^4) = 1 - 1*6 = -5
All additions and subtractions, the twiddle factor at this level is just 1.
Level 2: combine Level 1 results into degree-3 polynomials (twiddle factors: 1 and omega^2)
P_even(omega^0) = 8 + 1*6 = 14, P_even(omega^4) = 8 - 1*6 = 2
P_even(omega^2) = -2 + omega^2 * 2, P_even(omega^6) = -2 - omega^2 * 2
P_odd(omega^0) = 10 + 1*7 = 17, P_odd(omega^4) = 10 - 1*7 = 3
P_odd(omega^2) = -8 + omega^2 * (-5), P_odd(omega^6) = -8 - omega^2 * (-5)
Same butterfly, but now the twiddle factor alternates between 1 and omega^2. Where it's 1, we get clean numbers. Where it's omega^2, the expression stays symbolic, but in a finite field omega is a concrete integer, so these all resolve to specific values.
Level 3: combine Level 2 results into the full degree-7 polynomial (twiddle factors: 1, omega, omega^2, omega^3)
P(omega^0) = 14 + 1*17 = 31, P(omega^4) = 14 - 1*17 = -3
P(omega^1) = P_even(omega^2) + omega * P_odd(omega^2)
P(omega^2) = 2 + omega^2 * 3, P(omega^6) = 2 - omega^2 * 3
All eight evaluations, four twiddle factors. Each line: left and right share the same computation, the right just flips the sign. That's the butterfly: same work, two outputs.
Why "butterfly"? Look at each crossing in the diagram: two lines cross and diverge, like a pair of wings. If you squint hard enough, you can see it. It would be a shame to waste such a cool name.
(The original post has an interactive SVG butterfly diagram showing the full 8-point FFT with all three stages.)
The pattern widens at each stage: stage 1 combines adjacent pairs, stage 2 reaches across two, stage 3 across four. At every crossing, two inputs produce two outputs, same computation, just + and -.
For our 8-coefficient polynomial, that's 3 stages of 8 operations each, N * log2(N) total instead of N^2. That's O(N log N) vs O(N^2).
Why does the halving work? omega^1 and omega^5 sit opposite each other on the circle, so squaring both gives the same result: (omega^1)^2 = (omega^5)^2 = omega^2. Eight roots collapse to four distinct squared values, four to two, two to one.
Here's the full cascade: all 8 roots of unity, squaring down level by level. Each level, paired roots produce the same squared value, halving the distinct values.
use ark_bn254::Fr;
use ark_ff::{FftField, Field};
fn main() {
let omega = Fr::get_root_of_unity(8 as u64).expect("8th root exists");
let roots: Vec<Fr> = (0..8).map(|k| omega.pow([k])).collect();
// Square every root: opposite pairs collapse.
// 8 inputs -> only 4 distinct squared values.
println!("8 roots -> square -> 4 distinct values:");
for k in 0..4 {
let sq_k = roots[k] * roots[k];
let sq_opp = roots[k + 4] * roots[k + 4];
assert_eq!(sq_k, sq_opp);
println!(" (omega^{})^2 == (omega^{})^2 check", k, k + 4);
}
// Square again: 4 -> 2 distinct values.
let level1: Vec<Fr> = (0..4).map(|k| roots[k] * roots[k]).collect();
println!("\n4 values -> square -> 2 distinct values:");
for k in 0..2 {
assert_eq!(level1[k] * level1[k], level1[k + 2] * level1[k + 2]);
println!(" level1[{}]^2 == level1[{}]^2 check", k, k + 2);
}
// Once more: 2 -> 1.
let level2: Vec<Fr> = (0..2).map(|k| level1[k] * level1[k]).collect();
assert_eq!(level2[0] * level2[0], level2[1] * level2[1]);
println!("\n2 values -> square -> 1 value check");
}
Three levels of halving, so half the computation at each level is shared with the level above.
Here's the math for a real ZK domain:
N = 2^20 (1,048,576 points, a common circuit size):
- Naive: N^2 = 2^40 ~ 1.1 trillion field multiplications
- NTT: N * log2(N) = 2^20 * 20 ~ 21 million
- Speedup: ~50,000x
N = 2^25 (~33 million points, large production circuits):
- Naive: N^2 = 2^50 ~ 10^15
- NTT: N * log2(N) = 2^25 * 25 ~ 839 million
- Speedup: > 1,000,000x
Here are both approaches side by side, step by step, counting every field multiplication. First, the naive way: for each root of unity, walk all coefficients and compute P(x) = a_0 + a_1*x + a_2*x^2 + ...
use ark_bn254::Fr;
use ark_ff::FftField;
fn main() {
let n = 8usize;
let omega = Fr::get_root_of_unity(n as u64).expect("root exists");
let coeffs: Vec<Fr> = vec![10, 20, 30, 40, 0, 0, 0, 0]
.into_iter().map(|v| Fr::from(v)).collect();
let mut total_muls = 0u64;
let mut results = Vec::new();
println!("Naive: evaluate P(x) at each root, one by one");
println!("(each * = one field multiplication)\n");
let mut point = Fr::from(1);
for k in 0..n {
let mut muls_here = 0u64;
let mut result = Fr::from(0);
let mut x_pow = Fr::from(1);
for j in 0..n {
result += coeffs[j] * x_pow;
muls_here += 1;
if j < n - 1 {
x_pow *= point;
muls_here += 1;
}
}
if k < n - 1 {
point *= omega;
muls_here += 1;
}
total_muls += muls_here;
results.push(result);
println!(" P(w^{}): {} ({} muls)", k, "*".repeat(muls_here as usize), muls_here);
}
println!();
for k in 0..n {
println!(" P(omega^{}) = {:?}", k, results[k]);
}
println!("\nTotal field multiplications: {}", total_muls);
}
Now the same polynomial, same roots, but using the butterfly. Split into even/odd, recurse, combine with + and -. Count the multiplications.
use ark_bn254::Fr;
use ark_ff::FftField;
fn ntt(coeffs: &[Fr], omega: Fr, depth: usize, muls_per_depth: &mut Vec<u64>) -> Vec<Fr> {
let n = coeffs.len();
if n == 1 {
return coeffs.to_vec();
}
if depth >= muls_per_depth.len() {
muls_per_depth.resize(depth + 1, 0);
}
let even: Vec<Fr> = coeffs.iter().step_by(2).cloned().collect();
let odd: Vec<Fr> = coeffs.iter().skip(1).step_by(2).cloned().collect();
let omega_sq = omega * omega;
let y_even = ntt(&even, omega_sq, depth + 1, muls_per_depth);
let y_odd = ntt(&odd, omega_sq, depth + 1, muls_per_depth);
let mut result = vec![Fr::from(0); n];
let mut twiddle = Fr::from(1);
for k in 0..n / 2 {
let t = twiddle * y_odd[k];
muls_per_depth[depth] += 1;
result[k] = y_even[k] + t;
result[k + n / 2] = y_even[k] - t;
if k < n / 2 - 1 {
twiddle *= omega;
muls_per_depth[depth] += 1;
}
}
result
}
fn main() {
let n = 8usize;
let omega = Fr::get_root_of_unity(n as u64).expect("root exists");
let coeffs: Vec<Fr> = vec![10, 20, 30, 40, 0, 0, 0, 0]
.into_iter().map(|v| Fr::from(v)).collect();
let mut muls_per_depth: Vec<u64> = Vec::new();
let result = ntt(&coeffs, omega, 0, &mut muls_per_depth);
println!("NTT butterfly: all evaluations at once");
println!("(each * = one field multiplication)\n");
let mut total = 0u64;
for (d, &m) in muls_per_depth.iter().enumerate() {
total += m;
let size = n >> d;
println!(" depth {} (n={:>2}): {} ({} muls)", d, size, "*".repeat(m as usize), m);
}
println!();
for k in 0..n {
println!(" P(omega^{}) = {:?}", k, result[k]);
}
println!("\nTotal field multiplications: {}", total);
}
Same results, different cost. For N = 8 the naive approach uses 127 field multiplications, the butterfly uses 17. That's already 7x fewer, and the ratio keeps growing with N.
This algorithm is the Fast Fourier Transform (FFT). When it runs over finite fields instead of complex numbers, it's called the Number Theoretic Transform (NTT). Same structure, same recursion, different arithmetic.
One more thing: the recursive halving needs to split cleanly at every level (N -> N/2 -> N/4 -> ... -> 1). That only works if N is a power of two.
This is why every ZK framework uses domain sizes like 2^16 (65,536) or 2^20 (1,048,576). Not an arbitrary design choice: a direct consequence of the algorithm. When a ZK tutorial says "pad your trace to a power of two," now you know why. We'll talk about traces in the next post.
The field itself also needs to support this: p - 1 must be divisible by large powers of 2 so the recursive halving has room to work. This is called "2-adicity."
BN254 has a 2-adicity of 28, meaning it supports NTT domains up to 2^28 (about 268 million points). Plenty for most ZK circuits.
Let's See It Work
Enough theory. Let's use roots of unity for what they're actually for in ZK.
Remember post #1? We stored data inside a polynomial by choosing x values and evaluating.
Roots of unity give us those x values for free: omega^0, omega^1, ..., omega^{n-1}. You assign each data value to a root, run the inverse NTT, and out come the polynomial's coefficients. That's interpolation in O(N log N).
This is what a ZK prover does with the execution trace to get the data in polynomial form (next post).
use ark_bn254::Fr;
use ark_poly::{
EvaluationDomain, Radix2EvaluationDomain,
DenseUVPolynomial, Polynomial, univariate::DensePolynomial,
};
fn main() {
// Build a domain of 4 roots of unity
let domain = Radix2EvaluationDomain::<Fr>::new(4).unwrap();
let omega = domain.group_gen();
println!("Domain size: {}", domain.size());
println!("Generator (omega): {:?}\n", omega);
// Encode the word "cool" as ASCII values:
// 'c'=99, 'o'=111, 'o'=111, 'l'=108
// Each value is assigned to a root of unity:
// data[0] at omega^0, data[1] at omega^1, ..., data[3] at omega^3
let mut data = vec![
Fr::from(99), Fr::from(111), Fr::from(111), Fr::from(108),
];
// Inverse NTT: find the polynomial whose evaluations at roots of unity
// give back this data. This is interpolation. O(N log N).
domain.ifft_in_place(&mut data);
// data now holds coefficients. Build the polynomial.
let poly = DensePolynomial::from_coefficients_vec(data);
// Verify by evaluating at roots of unity, just like post #1.
// If interpolation worked, P(omega^i) gives back our original data.
println!("Evaluating the interpolated polynomial:\n");
let mut point = Fr::from(1); // omega^0 = 1
for i in 0..4 {
let value = poly.evaluate(&point);
println!(" P(omega^{}) = {:?}", i, value);
point *= omega;
}
println!("\nThat's [99, 111, 111, 108] = \"cool\" in ASCII.");
println!("The inverse NTT found the right polynomial.");
}
The inverse NTT took four ASCII values and found a polynomial that encodes them. We verified it the old-fashioned way: by calling evaluate() at omega^0, omega^1, omega^2, omega^3 and getting back 99, 111, 111, 108. Same evaluate() from post #1, but now the x values are roots of unity instead of arbitrary numbers.
The EvaluationDomain trait covers more operations beyond FFT/iFFT, but these two are the core.
What You Just Did
You generated roots of unity over a finite field, saw why their symmetry enables O(N log N) evaluation, and used the inverse NTT to interpolate a polynomial from raw data.
"Roots of unity"? It is just another name. "Unity" means 1 (the multiplicative identity), "roots" means solutions to an equation. So the whole thing literally means "solutions to z^n = 1."
Leonhard Euler studied them in 1751 in Recherches sur les racines imaginaires des equations. Carl Friedrich Gauss wrote 74 pages formalizing the theory in 1801 in Disquisitiones Arithmeticae, Section VII. The name comes from their work. Now you've used the tool behind it.
Next Up: Execution Trace
We have our fast polynomial engine. Now: how do we represent an actual program, with all its steps and state, as polynomials?
That's the execution trace. See you there.
If you want to discuss this with other Rust devs learning ZK, join the Discord.
What stuck with me the most: when I first saw omega and "rotation" in ZK code, I had no idea what was going on. It felt like geometry, not finite fields.
In geometry, Greek letters usually name angles, and rotating by angle omega makes perfect sense. In integers modulo a prime? Not so much.
But then I started looking for roots of unity in other places in mathematics and found them in FFT. And I figured out that (yet again) these are just names. They convey some meaning, but they're also misleading if you take them literally. When you know the origin story, it kinda makes sense why they're called roots of unity even for finite fields.
This is post #2 in my "first principles" ZK series for Rust devs. Interactive version with runnable code and animations on rustarians.com.
I'm running a free live session on March 17: From Polynomials to Proof in 60 Minutes. Range check circuit built step by step, with a real use case.
Questions? Found a bug? linkedin.com/in/adam-smolarek
Top comments (0)