Rohit's Blog

The last of four, following “Lava: Hardware Design in Haskell” by Bjesse, Claessen, Sheeran and Singh (ICFP 1998). Parts 1, 2 and 3 built the capability hierarchy, the combinator vocabulary, and the symbolic interpretation that turns a circuit into something a solver can chew on.

Everything so far has been demonstrated on a half adder. This post does the paper’s actual demonstration: two genuinely different FFT algorithms, described compactly, run under both interpretations, and proved equivalent.

And the proof needs something the previous post didn’t have — not more gates, but an algebra.

🔗Extending the system with complex numbers

Signal processing needs complex arithmetic, and Lava doesn’t have it. The paper’s answer is its architectural claim in action: adding a datatype and a capability class shouldn’t disturb anything that already works.

Two flavours are needed at once — concrete values for simulation, and abstract variables for verification:

enum CmplxSig:
  case Concrete(re: Double, im: Double)
  case Abstract(v: Var)

And a new capability class, sitting under Arithmetic:

trait CmplxArithmetic[F[_]] extends Arithmetic[F]:
  type Cmplx
  def cplus (a: Cmplx, b: Cmplx): F[Cmplx]
  def csub  (a: Cmplx, b: Cmplx): F[Cmplx]
  def ctimes(a: Cmplx, b: Cmplx): F[Cmplx]
  def w(n: Int, k: Int): F[Cmplx]          // the twiddle factor W_n^k

Extending the existing interpretations to support it is mostly lifting: the complex operations are implemented in terms of the arithmetic operations already there, over concrete doubles in Std and over symbolic terms in Sym.

The w operation is the interesting one, because it means different things to different interpretations. Under simulation it produces an actual complex constant. Under symbolic evaluation it produces an opaque term W(n, k) — a symbol the solver will reason about using laws rather than arithmetic. Same operation, two readings, which is the entire point of the design.

Worth pausing on this in Scala terms. In Composing Contracts part 2 I hit the problem that scala.math.Numeric can’t hold a symbolic quantity, because it demands compare and toDouble on something that has neither. Here the problem doesn’t arise: the arithmetic lives in a capability trait indexed by the interpretation, so ctimes under Sym builds a term and ctimes under Std multiplies doubles, and neither pretends to be a Numeric. Tagless final sidesteps the whole difficulty, and it’s a better answer than the extension methods I settled for there.

🔗What an FFT is, in one paragraph

The Discrete Fourier Transform maps a sequence of N complex numbers to another, where each output depends on every input, weighted by powers of a constant W_N — the Nth root of unity, e^−2πi/N. Those weights are called twiddle factors. Computed directly that’s N² multiplications, which in hardware means area and power you can’t afford.

The Fast Fourier Transforms are algorithms that exploit symmetries in the twiddle factors to do it in N log N. Four laws do the work:

W(n, 0)              = 1
W(n, n)              = 1
W(n, a) · W(n, b)    = W(n, a + b)
W(n, k)              = W(2n, 2k)

Nothing about circuits there. That’s an algebra of roots of unity, and — file this away — it is the entire mathematical content of the correctness argument that comes at the end.

One more fact we’ll need: W₄¹ = −i. Multiplying by it is a rotation, which in hardware is free — you swap the real and imaginary wires and negate one.

🔗The components

A butterfly is the two-point DFT: two inputs to their difference and their sum. It’s the atom every FFT is built from.

def bfly[F[_]](using c: CmplxArithmetic[F]): List[c.Cmplx] => F[List[c.Cmplx]] =
  case List(i1, i2) =>
    for
      o1 <- c.csub (i1, i2)
      o2 <- c.cplus(i1, i2)
    yield List(o1, o2)

A butterfly stage applies a column of these across a bus — pairing wire i with wire i + k/2. Drawn out, that wiring is a thicket of crossing lines. Expressed with the riffle combinator from part 2, it’s three steps: bring the pairs adjacent, run the column, put them back.

def bflys[F[_]: Monad](using c: CmplxArithmetic[F])(n: Int) =
  riffleK >-> raised(n)(two)(bfly) >-> unriffleK

That is the same line as in the paper, modulo syntax, and it’s the one I’d point at if someone asked what combinators buy you. The permutation and its inverse bracket the computation, so the crossings never appear.

Twiddle multiplication scales one wire by a constant:

def wMult[F[_]](using c: CmplxArithmetic[F])(n: Int, k: Int)(a: c.Cmplx): F[c.Cmplx] =
  c.w(n, k).flatMap(twd => c.ctimes(twd, a))

and multiplying a whole bus by −i is that, at W₄¹, mapped across:

def minusI[F[_]](using c: CmplxArithmetic[F]) = (xs: List[c.Cmplx]) =>
  xs.traverse(wMult(4, 1))

Plus bit reversal, which part 2 derived from repeated riffling.

🔗Two algorithms

The radix-2 FFT is the standard one. Input length a power of two, halved repeatedly until you reach length-2 DFTs — butterflies — then recombined stage by stage. In Lava it’s a bit reversal composed with n stages, where each stage is a twiddle layer feeding a butterfly layer, replicated across the bus by raised (n-i) two:

type Fft[F[_]] = Int => List[Cmplx] => F[List[Cmplx]]

def radix2[F[_]](using c: CmplxArithmetic[F]): Fft[F] = n =>
  bitRev(n) >-> compose((1 to n).toList.map(stage(n, _)))

def stage[F[_]](using c: CmplxArithmetic[F])(n: Int, i: Int) =
  raised(n - i)(two)(twid(i) >-> bflys(i - 1))

def twid[F[_]](using c: CmplxArithmetic[F])(i: Int) =
  one(decmap(1 << (i - 1))((k, x) => wMult(1 << i, k)(x)))

The radix-2² FFT is the less well-known one, and at a glance you might mistake it for a reversed radix-2. It isn’t. The paper is careful about how it differs: it assumes the input length is a power of four; each stage uses two different butterfly networks rather than one; the twiddle multiplications are modified; and −i multiplication stages are inserted between the butterfly layers.

Structurally, a stage becomes: a butterfly layer, then −i applied to one quarter of the bus, then a second butterfly layer on both halves, then the twiddle layer — the whole thing replicated by raised (m-i) (two ∘ two), since we’re now quartering rather than halving. The stages run in decreasing order and the bit reversal comes at the end rather than the start.

I’ve given radix-2² structurally rather than as a line-by-line port. The scan of the paper I worked from garbles several operators — exponents come out as subtraction, list literals lose their commas — and its twiddle layer uses a column combinator that isn’t defined in the text I have. I’d rather describe the shape accurately than transcribe something I can’t read with confidence. The radix-2 above I’m confident in.

What matters for the punchline is that the two circuits look nothing alike. Different butterfly networks, different twiddle placement, extra multiplication stages in one and not the other, bit reversal at opposite ends. No local inspection tells you they compute the same function.

The paper notes the corresponding VHDL descriptions would be several times longer, which I believe: every one of those stage descriptions is a generate loop with index arithmetic attached to instances.

🔗Simulating

The standard interpretation runs them on concrete inputs — supply an exponent and a list of complex values, get a list of complex values:

val input = List(Concrete(1, 4), Concrete(2, -2), Concrete(3, 2), Concrete(1, 2))

scala> simulate(radix2(using Std)(2)(input))

Ordinary testing, on the same description. Which is worth stating plainly because it’s the thing VHDL forced you to give up: you are not simulating a model of the circuit that might have drifted from the thing you’ll fabricate. There is one description.

🔗Proving them equivalent

Now the demonstration. Following the pattern from part 3 — the property is itself a circuit — we build a description that feeds symbolic complex variables to both algorithms and asks whether the outputs agree:

def fftSame[F[_]](using c: CmplxArithmetic[F] & Symbolic[F])(n: Int): F[c.Bit] =
  for
    inp  <- newCmplxVector(1 << (2 * n))
    out1 <- radix2 (2 * n)(inp)
    out2 <- radix22(n)(inp)
    eq   <- equalLists(out1, out2)
  yield eq

Note the intersection type CmplxArithmetic[F] & Symbolic[F] — this circuit needs both capabilities, and Scala 3 says so directly where Haskell would list two constraints. A small thing, but it’s the kind of place where Scala’s type system is doing the same work with less ceremony.

Run the symbolic interpretation and you get a large first-order term for each output, built from plus, times, sub and opaque W(n, k) symbols.

And here is where Boolean reasoning runs out. In part 3 the half-adder equivalence was propositional: gates over bits, a SAT solver’s home ground. This is not that. These are terms over complex numbers, and the two circuits are equal only because of algebraic facts about roots of unity. No amount of case-splitting on bits will get you there.

So the solver has to be given the algebra. The paper supplies it as theories — user-defined, handed to the prover as an option — of two kinds:

arithmetic laws: multiplication distributes over addition and subtraction, multiplication is associative, 1 is the identity;
twiddle-factor laws: the four listed earlier.

With those, the paper’s prover shows the size-4 circuits equivalent, and they report having proved size 16 and size 64 as well.

The modern spelling is SMT-LIB with uninterpreted functions and quantified axioms — declare W as an uninterpreted function, assert the four laws as universally quantified equalities, declare plus/times/sub likewise with the ring axioms, emit both circuits as terms, and ask whether their difference is satisfiable. Same structure, a solver you can brew install.

As in part 3: I’ve written this and not run it. Equational reasoning over quantified axioms is exactly the regime where solvers get slow or time out, and I’d expect real fiddling with triggers and encodings before size 16 goes through. I’ll report what actually happens rather than assume it works.

🔗What the whole series was about

Four posts, and the arc is worth stating.

One description, four interpretations. Simulation under Id, symbolic evaluation under State, and from the symbolic netlist both a solver input and structural VHDL. The description never changed. That’s the property VHDL can’t give you, and it’s why the authors resisted designing another hardware description language and embedded one instead.

Combinators are what make it readable. riffle, two, raised, decmap — none of which mention Circuit, all of which therefore work under every interpretation. The FFTs fit on a page because the regularity had somewhere to live.

The proof needed an algebra, not a solver. This is the part I keep turning over. The solver is a proof engine; the actual content is four equations about roots of unity. The same shape as the tile algebra, where the laws did the design work and the tests merely checked the transcription — except here the equational reasoning is performed by machine over a term with thousands of nodes, which is not something you’d do by hand.

And Scala 3 held up well. Type members gave each interpretation its own datatypes — the generalisation the paper’s future work asked for. Intersection types expressed multiple capabilities directly. Tagless final is idiomatic enough that the central trick needed no explanation.

The costs were real but small: path-dependent types mean using clauses come first and signatures read oddly; every gate is a monadic bind, which the paper itself concedes against Hydra’s monad-free directness; and cats asks for tailRecM in places the paper never had to think about.

Nothing here is a reason to write Haskell. As ever, it’s a reason to know which of my Scala habits are choices and which are consequences — and this time, one of the consequences went Scala’s way.

🔗Sources

Bjesse, Per, Koen Claessen, Mary Sheeran & Satnam Singh. Lava: Hardware Design in Haskell. ICFP 1998. The complex-number extension and CmplxArithmetic class, the butterfly, the riffle-based butterfly stage, twiddle multiplication, the −i stage, the radix-2 and radix-2² descriptions, the equivalence question, and the use of user-supplied arithmetic and twiddle-factor theories are all theirs.
Jones, Geraint. Deriving the fast Fourier algorithm by calculation. 1990 — the Ruby work the FFT derivation builds on.
Bjesse, Per. Specification of signal processing programs in a pure functional language. Master’s thesis, Chalmers, 1997.
He, Shousheng. Radix-2² work (1995), cited by the paper for the second algorithm.

The Scala 3 translation, the intersection-type encoding of multiple capabilities, the SMT-LIB modernisation, and the observation that tagless final sidesteps the Numeric problem from the contracts series are mine. The radix-2² description is structural rather than a line-by-line port, for the reason noted above.