# A Formalisation of a Fast Fourier Transform

Laurent Théry

► To cite this version:

Laurent Théry. A Formalisation of a Fast Fourier Transform. 2022. hal-03807965v4

HAL Id: hal-03807965

https://inria.hal.science/hal-03807965v4

Preprint submitted on 13 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
<figure><img src="image_1.png" /></figure>

HAL
open science

# Revisiting the Fast Fourier Transform in Rocq

Laurent Théry

► To cite this version:

Laurent Théry. Revisiting the Fast Fourier Transform in Rocq. Inria & Université Cote d'Azur, CNRS, I3S, Sophia Antipolis, France. 2025. hal-03807965v5

HAL Id: hal-03807965

https://inria.hal.science/hal-03807965v5

Submitted on 14 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Formalisation of a Fast Fourier Transform

Laurent Théry

Laurent.Thery@sophia.inria.fr

## Abstract

This notes explains how a standard algorithm that constructs the discrete Fourier transform has been formalised and proved correct in the Coq proof assistant using the SSREFLECT extension.

## 1 Introduction

Fast Fourier Transforms are key tools in many areas. In this note, we are going to explain how they have been formalised in a theorem prover like Coq.

There are many ways to motivate this work. As univariate polynomials is at the heart of our formalisation, we take, here, one application : the polynomial multiplication.

If we take two polynomials $p$ $q$ of degree $n$, a naive way of multiplying them requires a quadratic number $n^2$ of multiplications, i.e. we perform the multiplication of each term of $p$ (there are at most $n$ of them) by a term of $q$ (again $n$ of those). It is possible to get a faster algorithm using a fast Fourier transform. The idea is to take an alternative representation for polynomials : a polynomial $p$ of degree $n$ is uniquely defined by its value $p(x_1), ..., p(x_n)$ on $n$ distinct points $x_1, ..., x_n$. This representation is what is called polynomial interpolation and the $x_i$s are called the interpolation points.

A direct way to effectively construct a polynomial $p$ from its evaluation $p(x_1), ..., p(x_n)$ is by using Lagrange polynomial $L_{n,k} = \prod_{i \neq k} (x - x_i) / (x_k - x_i)$. It is easy to check that $L_{n,k}(x_j) = 1$ if $j = k$ or 0 otherwise. Then, we have $p = \sum_{1 \le i \le k} p(x_i) L_{n,i}$. Computing the multiplication in the interpolation world is easy. The polynomial $pq$ is of degree at most $2n$ so if we have the values of $p$ and $q$ at $2n$ points $x_1, ..., x_{2n}$, we simply need to perform $2n$ multiplication to get $p(x_1)q(x_1), ..., p(x_{2n})q(x_{2n})$. Of course, we need algorithms to move from
# Revisiting the Fast Fourier Transform in Rocq

Laurent Théry

Stamp Team - INRIA

Laurent.Thery@inria.fr

## Abstract

This note explains how a standard algorithm that constructs the discrete Fourier transform has been formalised and proved correct in the ROCQ proof assistant using the SSREFLECT extension.

## 1 Introduction

Fast Fourier Transforms are key tools in many areas. In this note, we are going to explain how they have been formalised in a theorem prover like ROCQ. The presentation of this work has three parts. The central part of this work is the formalisation of the algorithm in itself. We have aimed to get a nice and effective formal presentation. Then, two extensions have been developed. The first one builds a dedicated tactic using the Fourier algorithm in order to implement polynomial multiplication. The second one goes toward computer arithmetic and formally proves some standard results about error bounds while computing Fourier transform with floating-point numbers.

Note that this is not the first time this algorithm has been formalised in ROCQ. To our knowledge, the first one was done in 2001 by Venanzio Capretta (see [4]). It is about time to revisit this work and see how concise it can get using existing libraries. Also, this initial effort was mostly interested in the recursive presentation of the algorithm. Here, we also give an iterative version.

The paper is organised as follows. We first describe the formalisation of the main algorithm. Then, we first present the tactic that has been built
the usual polynomial representation to the interpolation one and back. This
is what the fast Fourier transform gives us. The trick is that we are going to
evaluate the polynomials at very specific points of interpolation (the $x_i$s are
roots of unit). What makes it work is that, for those points, the potential $n^2$
values of the different $(x_i)^j$ ($1 \le i, j \le n$) consist of only $n$ distinct values,
so computing the $p(x_1), \dots, p(x_n)$ can be done very efficiently in $n\log(n)$.

Note that this is not the first time fast Fourier transform has been formalised in COQ. To our knowledge, the first one was done in 2001 by Venanzio Capretta (see [2]). It is about time to revisit this work and see how concise it can get using existing libraries. Also, this initial effort was mostly interested in the recursive presentation of the algorithm. Here, we also give an iterative version.

## 2 Formalisation of the recursive algorithm

The algorithm manipulates three kinds of data:

- natural numbers, i.e. elements of type nat;

- elements of an arbitrary integral domain $R$;

- univariate polynomials over $R$, i.e. elements of type {poly $R$}.

The operations we use on natural numbers are the successor and predecessor functions ($i.+1$ and $i.-1$), the doubling and halving functions ($i.*2$, $i./2$ and `uphalf` $i = (i.+1./2)$), the exponentiation ($i.^j$), the division ($i%j$), and the modulo ($i%%j$). For the integral domain, we use the usual ring operations : the addition $x + y$, the multiplication by a scalar $k * : x$, the multiplication $x * y$, and the exponentiation by a natural number $x ^{n}$. Also, a predicate of the library that is useful in our application is the one that indicates that an element $w$ is a $n^{th}$ primitive root of unity. It is written as $n.-primitive_root w$. It means that $w ^{n} = 1$ and that $n$ is actually the smallest non-zero natural number that has this property. In our formalisation, we derive two easy lemmas about primitive roots

Lemma `prim_exp2nS` $n$ ($w : R$) :
(2 ^ {n.+1}).-primitive_root $w \rightarrow w ^{+} (2 ^ {n}) = -1$.

Lemma `prim_sqr` $n$ ($w : R$) :
(2 ^ {n.+1}).-primitive_root $w \rightarrow (2 ^ {n}).-primitive_root (w ^{+} 2)$.
on top of the algorithm. Finally, we mention some results concerning error bounds when computing with floating-point numbers.

# 2 Main algorithms

## The recursive algorithm

The algorithm manipulates three kinds of data:

- natural numbers, i.e. elements of type nat;

- elements of an arbitrary integral domain R;

- univariate polynomials over R, i.e. elements of type {poly R}.

The operations we use on natural numbers are the successor and predecessor functions ($i.+1$ and $i.-1$), the doubling and halving functions ($i.*2$, $i./2$ and $uphalf i = (i.+1./2)$), the exponentiation ($i^j$), the division ($i%j$), and the modulo ($i%%j$). For the integral domain, we use the usual ring operations: the addition $x + y$, the multiplication by a scalar $k * : x$, the multiplication $x * y$, and the exponentiation by a natural number $x ^+ n$. Also, a predicate of the library that is useful in our application is the one that indicates that an element $w$ is a $n^{th}$ primitive root of unity. It is written as $n.-primitive_root w$. It means that $w ^+ n = 1$ and that $n$ is actually the smallest non-zero natural number that has this property. In our formalisation, we derive two easy lemmas about primitive roots

Lemma `prim_exp2nS n (w : R)`: 
(2 ^ n.+1).-primitive_root w -> w ^+ (2 ^ n) = -1.

Lemma `prim_sqr n (w : R)`: 
(2 ^ n.+1).-primitive_root w -> (2 ^ n).-primitive_root (w ^+ 2).

The first lemma is used to simplify expression involving $w^j$. Its proof is the only place where the integrality is used (in order to get from $w^2 = 1$ that either $w = 1$ or $w = -1$). The recursive algorithm of the fast Fourier transform takes as argument a primitive root $w$ and performs some recursive calls with $w ^+ 2$. The lemma `prim_sqr` is then used to prove that the primitive root property of the argument is an invariant of the recursion.

The algorithm is mainly manipulating univariate polynomials. A polynomial is represented by a list whose last element (the leading term), if it exists,
The first lemma is used to simplify expression involving $w^j$. Its proof is the only place where the integrality is used (in order to get from $w^2 = 1$ that either $w = 1$ or $w = -1$). The recursive algorithm of the fast Fourier transform takes as argument a primitive root $w$ and performs some recursive calls with $w \stackrel{\sim}{+} 2$. The lemma *prim_sqr* is then used to prove that the primitive root property of the argument is an invariant of the recursion.

The algorithm is mainly manipulating univariate polynomials. A polynomial is represented by a list whose last element (the leading term), if it exists, is non-zero. The empty list represents the null polynomial. A polynomial $p$ can be automatically converted to a list, so `size p` is understood as the length of the list representing $p$. So, if $p$ is not null, it is the usual degree of the polynomial incremented by one. We can access the $n^{th}$ term of polynomial $p$ by $p'_n$. The leading term is then $p'_-(size p) \cdot (-1)$. The variable of the univariate polynomials is written `$X$` and `$X^n$` its power. Evaluating a polynomial $p$ at point $x$ is written $p. [x]$. Composition two polynomials $p$ and $q$ which consists in lifting the evaluation from points to polynomials $p. [q]$ is written $p \ \text{Po} \ q$. We can turn a function $F$ into a polynomial using `\poly_(i < n) F i` that builds a polynomial of size $n$ whose $i^{th}$ term is $F i$. A special notation is available for constant polynomials where one can write `$c%:P$` in order to build a polynomial that only contains the constant term $c$.

The recursive algorithm is taking two arguments : a polynomial $p$ and a primitive root $w$ of degree $2^n$ and returns the evaluation of $p$ at points $1, w, w^2, ..., w^{2^n-1}$ as the polynomial $p. [1] + p. [w] X + ... + p. [w^{2^n-1}] X^{2^n-1}$. The recursive calls are made on the even and odd parts of the polynomial $p$. If $p = 1 + 2X + 3X^2 + 4X^3 + 5X^4$, its even part is $1 + 3X + 5X^2$ and its even part $2 + 4X$. As these operations on polynomials are not in the library we had to define them.

**Definition even_polymul : {poly R} := \poly_(i < \text{uphalf (size p)}) p'_i.*2.**

**Definition odd_polymul : {poly R} := \poly_(i < (\text{size } p)/2) p'_i.*2.+1.**

It is then easy to derive the key lemma that justifies the decomposition of
the polynomial in the recursive calls

**Lemma poly_even_odd p : (even_polymul p \ \text{Po} 'X^2) + (odd_polymul p \ \text{Po} 'X^2) * 'X = p.**
is non-zero. The empty list represents the null polynomial. A polynomial $p$ can be automatically converted to a list, so `size p` is understood as the length of the list representing $p$. So, if $p$ is not null, it is the usual degree of the polynomial incremented by one. We can access the $n^{th}$ term of polynomial $p$ by $p'_n$. The leading term is then $p'_-(size p).-1$. The unique variable of the univariate polynomials is written `$X$` and `$X^n$` its power. Evaluating a polynomial $p$ at point $x$ is written $p. [x]$. Composition two polynomials $p$ and $q$ which consists in lifting the evaluation from points to polynomials $p. [q]$ is written $p \ \text{Po} q$. We can turn a function $F$ into a polynomial using `\poly_(i < n) F i$ that builds a polynomial of size at most $n$ whose $i^{th}$ term is $F i$. A special notation is available for constant polynomials where one can write `$c%:P$` in order to build a polynomial that only contains the constant term $c$.

Our recursive algorithm is taking two arguments : a polynomial $p$ and a primitive root $w$ of degree $2^n$ and returns the evaluation of $p$ at points $1, w, w^2, ..., w^{2^n-1}$ as the polynomial $p. [1] + p. [w] X + ... + p. [w^{2^n-1}] X^{2^n-1}$. So, the return type is also a polynomial. The recursive calls are made on the even and odd parts of the polynomial $p$. If $p = 1 + 2X + 3X^2 + 4X^3 + 5X^4$, its even part is $1 + 3X + 5X^2$ and its even part $2 + 4X$. These operations on polynomials are in the library.

**Definition even polynomial $p$ : {poly $R$} := $\text{poly}_-(i < \text{uphalf}(\text{size } p)) p'_i.*2.$**

**Definition odd polynomial $p$ : {poly $R$} := $\text{poly}_-(i < (\text{size } p)/2) p'_i.*2.+1.$**

There is also the key lemma that justifies the decomposition of the polynomial in the recursive calls

**Lemma polynomial $p$ : $(even\_pol y p \ \text{Po} 'X^2) + (odd\_pol y p \ \text{Po} 'X^2)) * 'X = p.$**

We are now ready to present the algorithm we have proved correct:

**Fixpoint fft (n : nat) (w : R) (p : {poly R}) : {poly R} :=**
**if n is $n_1$.+1 then**
**let $ev :=$ `fft $n_1$ (`w` + 2) (`even polynomial $p$`) in**
**let $ov :=$ `fft $n_1$ (`w` + 2) (`odd polynomial $p$`) in**
**$\text{poly}_-(i < 2^-(n_1.+1))$ let $j := i \%% (2^j n_1)$ in $ev' - j + ov' - j * w^j + i$**
**else ($p'_0$)%:P.**
We are now ready to present the algorithm we have proved correct:

Fixpoint fft (n : nat) (w : R) (p : {poly R}) : {poly R} :=
  if n is n₁.+1 then
    let ev := fft n₁ (w ^+ 2) (even_polyp) in
    let ov := fft n₁ (w ^+ 2) (odd_polyp) in
    \poly_(i < 2 ^ (n₁.+1)) let j := i %% (2 ^ n₁) in ev'_j + ov'_j * w ^+ i
  else (p'_0)%;P.

It takes a natural number *n*, a root of unity *w* and a polynomial *p* and returns a polynomial whose coefficients are the value at the interpolation point *w*¹.
More formally, the correctness lemma is the following :

Lemma fftE n (w : R) p :
    size p ≤ 2 ^ n → (2 ^ n) . -primitive_root w →
    fft n w p = \poly_(i < 2 ^ n) p . [w ^+ i].

The proof is straightforward. It is done by induction. We have to prove the
equality of two polynomials, so we show that their *i*th coefficients are equal

(fft n (w ^+ 2) (even_polyp))'_(i %% 2 ^ n) +
(fft n (w ^+ 2) (odd_polyp))'_(i %% 2 ^ n) * w ^+ i = p . [w ^+ i]

with the assumption that *w* is a primitive root of order 2<sup>*n*+1</sup>. Using the
induction hypothesis on the left of the equality and the decomposition lemma
*odd\_even\_polyE* on the right, we get

(even_polyp) . [(w ^+ 2) ^+ (i %% 2 ^ n)] +
(odd_polyp) . [(w ^+ 2) ^+ (i %% 2 ^ n)] * w ^+ i =
(even_polyp) . [(w ^+ i) ^+ 2] + (odd_polyp) . [(w ^+ i) ^+ 2] * w ^+ i

This means we are left with proving the following equality

$$
(w ^{+ 2}) ^{+} (i ^{+} 2) \hat{n} = (w ^{+ i}) ^{+} 2
$$
It takes a natural number $n$, a root of unity $w$ and a polynomial $p$ and returns a polynomial whose coefficients are the values at the interpolation point $w^i$. More formally, the correctness lemma is the following:

**Lemma** `fftE n (w : R) p :`
size $p \le 2$ $\hat{n} \to (2 \hat{n})$. `-primitive_root w` $\to$
`fft n w p = \poly_{i < 2 \hat{n})} p. [w^{\hat{i}}].`

The proof is straightforward. It is done by induction. We have to prove the equality of two polynomials, so we show that their $i^{th}$ coefficients are equal

($fft n (w^{\hat{+}2}) (even\_poly p))'_{-(i \% \% 2 \hat{n})} +$
($fft n (w^{\hat{+}2}) (odd\_poly p))'_{-(i \% \% 2 \hat{n})} * w^{\hat{+}i} = p. [w^{\hat{+}i}]$

with the assumption that $w$ is a primitive root of order $2^{n+1}$. Using the induction hypothesis on the left of the equality and the decomposition lemma `odd_even_polyp` on the right, we get

($even\_poly p$). $[(w^{\hat{+}2})^{\hat{+}} (i \% \% 2 \hat{n})] +$
($odd\_poly p$). $[(w^{\hat{+}2})^{\hat{+}} (i \% \% 2 \hat{n})] * w^{\hat{+}i} =$
($even\_poly p$). $[(w^{\hat{+}i})^{\hat{+}2}] + (odd\_poly p). [(w^{\hat{+}i})^{\hat{+}2}] * w^{\hat{+}i}$

This means we are left with proving the following equality

($w^{\hat{+}2})^{\hat{+}} (i \% \% 2 \hat{n}) = (w^{\hat{+}i})^{\hat{+}2}$

that directly follows from the fact that $w^{2^{n+1}} = 1$.

Finally, we also prove an alternative version of the algorithm that more explicitly exhibits the data path, the so-called butterfly.

Fixpoint `fft1 n w p : {poly R} :=`
if $n$ is $n_1$. +1 then
let $ev := \text{fft1 } n_1 (w^{\hat{+}2}) (\text{even\_poly } p)$ in
let $ov := \text{fft1 } n_1 (w^{\hat{+}2}) (\text{odd\_poly } p)$ in
$\sum_{j < 2 \hat{n}_1}$
$((ev'_j + ov'_j * w^{\hat{+}j}) * 'X^{\hat{j}} + (ev'_j - ov'_j * w^{\hat{+}j}) * 'X^{\hat{j}}(j + 2 \hat{n}_1))$
else $(p'_0)\%:P.$
that directly follows from the fact that $w^{2n+1} = 1$.

Finally, we also prove an alternative version of the algorithm that more explicitly exhibits the data path, the so-called butterfly.

Fixpoint fft1 n w p : {poly R} :=
  if n is n₁.+1 then
    let ev := fft1 n₁ (w ^+ 2) (evenosi p) in
    let ov := fft1 n₁ (w ^+ 2) (odds p) in
    \sum_(j < 2 ^ n₁)
      ((ev'_j + ov'_j * w ^+ j) *: 'X^j +
       (ev'_j - ov'_j * w ^+ j) *: 'X^ (j + 2 ^ n₁))
  else (p'_0)%;P.

It is straightforward to prove that both recursive versions compute the same thing.

Lemma fft1E n (w : R) p : (2 ^ n) . -primitive_root w Ñ fft1 n w p = fft n w p.

3 Formalisation of the iterative algorithm

In our formalisation, we are going to derive an iterative version from the
recursive in a very straightforward way. Let us explain it on an example with
a polynomial of degree 7 (2³ - 1). The depth of the recursion is 3 and the
binary tree of the recursive calls looks like:

<figure><img src="image_1.png" /></figure>
It is straightforward to prove that both recursive versions compute the same thing.

Lemma `fft1E n (w : R) p : (2 ^ n). -primitive_root w -> fft1 n w p = fft n w p.`

## The iterative algorithm

In our formalisation, we are going to derive an iterative version from the recursive in a very straightforward way. Let us explain it on an example with a polynomial of degree 7 (2³ - 1). The depth of the recursion is 3 and the binary tree of the recursive calls looks like:

<figure><img src="image_3.png" /></figure>

The idea is to put all the results at depth *i* in a single polynomial. Here, at depth 3 it is a polynomial containing 8 sub polynomials of degree 0, at depth 2, 4 polynomials of degree 1, at depth 1, 2 polynomials of degree 3 and finally one polynomial of degree 7. The final result is built bottom up. Initially we start with the polynomials that contains all the leaves. Then, one step of the iteration simply take the results at depth *i* and returns the results at depth *i* - 1.

Let us first concentrate on the initial value, the values of all the leaves. If we remember that we use an even/odd partition, putting the even part on the left and the odd part on the right, this means that if we look at the
The idea is to put all the results at depth $i$ in a single polynomial. Here, at depth 3 it is a polynomial containing 8 sub polynomials of degree 0, at depth 2, 4 polynomials of degree 1, at depth 1, 2 polynomials of degree 3 and finally one polynomial of degree 7. The final result is built bottom up. Initially we start with the polynomials that contains all the leaves. Then, one step of the iteration simply take the results at depth $i$ and returns the results at depth $i-1$.

Let us first concentrate on the initial value, the values of all the leaves. If we remember that we use an even/odd partition, putting the odd part on the left and the right part on the right, this means that if we look at the binary representation, the bits from the right to the left let us sort. If we take our example, we have to sort

$$[0; 1; 2; 3; 4; 5; 6; 7]$$

With their binary representation, it gives

$$[0 \to 000; 1 \to 001; 2 \to 010 \to; 3 \to 011; 4 \to 100; 5 \to 101; 6 \to 110; 7 \to 111]$$

Reversing them, we get

$$[0 \to 000; 1 \to 100; 2 \to 010; 3 \to 110; 4 \to 001; 5 \to 101; 6 \to 011; 7 \to 111]$$

Translating them back to decimal numbers, we have

$$[0 \to 0 \quad ; 1 \to 4 \quad ; 2 \to 2 \quad ; 3 \to 6 \quad ; 4 \to 1 \quad ; 5 \to 5 \quad ; 6 \to 3 \quad ; 7 \to 7 \quad ]$$

To build this initial polynomial, we first define digitn b n m that computes the m^th digit of n in base b. We then use it to reverse a number : rdigitn b n m reverses the n bits in base b of m. Finally, the initial polynomial with 2^n terms is created by reverse_polym n p using an appropriate permutation of the coefficient of p.

Definition digitn b n m := (n %/ b ^ m) %% b.

Definition rdigitn b n m := \sum_(i < n) digitn b m (n.-1 - i) * b ^ i.

Definition reverse_polym n (p : {poly R}) := \poly_(i < 2 ^ n) p'_(rdigitn 2 n i).

On our example, reverse_polym 3 p returns

$$p_0 + p_4X + p_2X^2 + p_6X^3 + p_1X^4 + p_5X^5 + p_3X^6 + p_7X^7$$
binary representation, the bits in reverse order from the right to the left gives us the sorting order. If we take our example, we have to sort

$$[0; 1; 2; 3; 4; 5; 6; 7]$$

With their binary representation, it gives

$$[0 \to 000; 1 \to 001; 2 \to 010 \to 3 \to 011; 4 \to 100; 5 \to 101; 6 \to 110; 7 \to 111]$$

Reversing them, we get

$$[0 \to 000; 1 \to 100; 2 \to 010; 3 \to 110; 4 \to 001; 5 \to 101; 6 \to 011; 7 \to 111]$$

Translating them back to natural numbers, we have

$$[0 \to 0 \quad ; \quad 1 \to 4 \quad ; \quad 2 \to 2 \quad ; \quad 3 \to 6 \quad ; \quad 4 \to 1 \quad ; \quad 5 \to 5 \quad ; \quad 6 \to 3 \quad ; \quad 7 \to 7 \quad ]$$

To build this initial polynomial, we first define digitn b n m that computes the m^th digit of n in base b. We then use it to reverse a number : rdigitn b n m reverses the n first bits in base b of m. Finally, the initial polynomial with 2^n terms is created by reverse_polym n p using an appropriate permutation of the coefficient of p.

Definition digitn b n m := (n %/ b ^ m) %% b.

Definition rdigitn b n m := \sum_(i < n) digitn b m (n.-1 - i) * b ^ i.

Definition reverse_polym n (p : {poly R}) := \poly_(i < 2 ^ n) p'_(rdigitn 2 n i).

On our example, reverse_polym 3 p returns a polynomial with 2³ terms.

$$p_0 + p_4X + p_2X^2 + p_6X^3 + p_1X^4 + p_5X^5 + p_3X^6 + p_7X^7$$

Now, we want to express that after each step of the iteration we get all the results at depth i. In the recursion, the polynomial is split in two using the even and odd part, but the results are glued together using the left to right concatenation. So, we need some operations to get the low terms or the high term of a polynomial.
Definition takepoly m (p : {poly R}) := \poly_(i < m) p'_i.
Definition drop_ply m (p : {poly R}) := \poly_(i < \text{size } p - m) p'_(i + m).
Lemma poly_take_drop m p : take_ply m p + drop_ply m p * 'X^m = p.

`take_ply m p` returns the polynomial that has the `m` low terms of `p` while
`drop_ply m p` returns the polynomial with the high terms of `p` skipping the
`m` low terms.

We want to express that we have in a polynomial $q$ all the results of
calling the recursive algorithm $p$ cutting at depth $n$, so every leaf is a call
of $fft_1$ with the appropriate part of the polynomial $p$. This is done using a
recursive predicate. If $n$ is not zero, we split $p$ in two with even and odd part
and $q$ with left part and right part. If $n$ is zero, we are at a leaf so the result
must be a call to $fft_1$.

Fixpoint all_resultsff't1 n m w p q :=
    if n is n1.+1 then
        all_resultsff't1 n1 m w (even_ply p) (take_ply (2 ^ (m + n1))) q) /
        all_resultsFFT1 n1 m w (odd_ply p) (drop_ply (2 ^ (m + n1))) q
    else q = fFt1 m w p.

For the initial polynomial, we prove that, given a polynomial $p$ of size less
than $2^n$ the reverse polynomial has all the results of $fft_1$ 0 at depth $n$

Lemma all_results ``fft1 ``reverse ``poly p n w :
size p ≤ 2 ^ n → all_results ``fft1 n 0 w p (reverse ``poly n p ).

This lemma is proved by induction on $n$. The key lemma for proving it is
the fact `reverse_ply` decomposes nicely using odd and even parts.

Lemma reverse_plyS n p :
reverse_ply n.+1 p =
reverse_ply n (even_ply p) + reverse_ply n (odd_ply p) * 'X^ (2 ^ n).

Now, we can define a step of the iteration. If we are at depth $m$, we have
$2^{m+1}$ results of size $2^n$
**Lemma reverse_polys n p:**
reverse_polys n.+1 p =
reverse_polys n (even_polys p) + reverse_polys n (odd_polys p) * 'X^(2 ^ n).

Now, we can define a step of the iteration. If we are at depth $m$, we have $2^{m+1}$ results of size $2^n$

**Definition step m n w (p : {poly R}):**
\sum_(l < 2^m)
let ev := \poly_(i < 2^m) p'_(i + l * 2^m + 1) in
let ov := \poly_(i < 2^m) p'_(i + l * 2^m + 1 + 2^m) in
\sum_(j < 2^m)
((ev'_j + ov'_j * w^+ + j) *: 'X^(j + l * 2^m + 1) +
(ev'_j - ov'_j * w^+ + j) *: 'X^(j + l * 2^m + 1 + 2^m)).

and the correctness is that applying a step decreases the depth $m$ while increasing the size $n$.

**Lemma all_results_ftt1_step m n w (p q : {poly R}):**
size p ≤ 2 ^ (m + n).+1 →
size q ≤ 2 ^ (m + n).+1 →
all_results_ftt1 m.+1 n (w^+ 2) p q →
all_results_ftt1 m n.+1 w p (step m n w q).

Again, the key lemmas are the ones that show that step behaves well with `take_polys` and `drop_polys`.

**Lemma take_step m n w (p : {poly R}):**
size p ≤ 2 ^ (m + n).+2 →
take_polys (2 ^ (m + n).+1) (step m.+1 n w p) =
step m n w (take_polys (2 ^ (m + n).+1) p).

**Lemma drop_step m n w (p : {poly R}):**
size p ≤ 2 ^ (m + n).+2 →
drop_polys (2 ^ (m + n).+1) (step m.+1 n w p) =
step m n w (drop_polys (2 ^ (m + n).+1) p).

Now, we code the iteration of step, it is straightforward to prove that we get the same result as the recursive algorithm.
Definition step m n w ( p : {poly R } ) :=
\sum_(l < 2 ^ m)
let ev := \poly_(i < 2 ^ n) p'_(i + l * 2 ^ n.+1) in
let ov := \poly_(i < 2 ^ n) p'_(i + l * 2 ^ n.+1 + 2 ^ n) in
\sum_(j < 2 ^ n)
((ev'_j + ov'_j * w ^+ j) * : 'X^(j + l * 2 ^ n.+1) +
(ev'_j - ov'_j * w ^+ j) * : 'X^(j + l * 2 ^ n.+1 + 2 ^ n)).

and the correctness is that applying a step decreases the depth *m* while
increasing the size *n*.

Lemma all_results_ftt1_step m n w (p q : {poly R}):
size p ≤ 2 ^ (m + n).+1 →
size q ≤ 2 ^ (m + n).+1 →
all_results_ftt1 m.+1 n (w ^+ 2) p q →
all_results_ftt1 m n.+1 w p (step m n w q).

Again, the key lemmas are the ones that show that *step* behaves well with *takeosi* and *droposi*.

Lemma take_step m n w (p : {poly R}):
size p ≤ 2 ^ (m + n).+2 →
takeosi (2 ^ (m + n).+1) (step m.+1 n w p) =
step m n w (takeosi (2 ^ (m + n).+1) p).

Lemma drop_step m n w (p : {poly R}):
size p ≤ 2 ^ (m + n).+2 →
droposi (2 ^ (m + n).+1) (step m.+1 n w p) =
step m n w (droposi (2 ^ (m + n).+1) p).

Now, we code the iteration of *step*, it is straightforward to prove that we get
the same result as the recursive algorithm.

Fixpoint istep_aux m n w p :=
  if m is m_1.+1 then istep_aux m_1 n.+1 w (step m_1 n (w ^+ (2 ^ m_1)) p) else p.
Definition istep n w p := istep_aux n 0 w (reverseosi poly n p).
Lemma istep_ftt1 n p w : size p ≤ 2 ^ n → istep n w p = fft1 n w p.

Similarly we prove the correctness of an alternative and more direct iterative
version.
Fixpoint istep_aux m n w p :=
  if m is m₁.+1 then istep_aux m₁ n.+1 w (step m₁ n (w ^+ (2 ^ m₁)) p) else p.

Definition istep n w p := istep_aux n 0 w (reverse_pol y n p).

Lemma istep_ft t₁ n p w : size p ≤ 2 ^ n → istep n w p = fft t₁ n w p.

Similarly we prove the correctness of an alternative and more direct iterative
version.

Definition step₁ m n w (p : {poly R}) :=
\poly_(i < 2 ^ (m + n).+1)

let j := i % 2 ^ n.+1 in

if j < 2 ^ n then

$$
p'_i + p'_{-(i + 2 ^ n) * w ^+ j}
$$

else

$$
p'_{-(i - 2 ^ n) - p'_i * w ^+ (j - 2 ^ n)}
$$

Lemma step₁E m n w p : step₁ m n w p = step m n w p.

4 Inverse algorithm

We have seen how to go from polynomial to interpolation. What about the
other way around. In fact, we can use the same algorithm. From the direct
direction $R$ was an integral ring. Here, we need a field $F$. The notation for
inverse is $x^{-1}$. The idea of the inverse algorithm is to use $1/w$ instead of $w$.

Definition ifft n w p : {poly F} := (2 ^ n) % : R ^ -1 % : P * (ft t n w ^ -1 p).

where $n%:R$ is the coercion for natural number $n$ into $F$. Its correctness
follows if the characteristic of $F$ is not 2-.

Lemma fftK n (w : F) p :
2%:R != 0 → size p ≤ 2 ^ n → (2 ^ n).-primitive_root w →

ifft n w (fft n w p) = p.
Definition step₁ m n w (p : {poly R}) :=
\poly_(i < 2 ^ (m + n).+1)
let j := i %% 2 ^ n.+1 in
if j < 2 ^ n then
  p'_i + p'_(i + 2 ^ n) * w ^+ j
else
  p'_(i - 2 ^ n) - p'_i * w ^+ (j - 2 ^ n).
Lemma step₁ E m n w p : step₁ m n w p = step m n w p.

## The inverse algorithm

We have seen how to go from polynomial to interpolation. What about the other way around. In fact, we can use the same algorithm. From the direct direction $R$ was an integral ring. Here, we need a field $F$. The notation for inverse is $x^{-1}$. The idea of the inverse algorithm is to use $1/w$ instead of $w$.

Definition ifft n w p : {poly F} := (2 ^ n) % : R ^ -1 % : P * (ftt n w ^-1 p).

where $n%:R$ is the coercion for natural number $n$ into $F$. Its correctness follows if the characteristic of $F$ is not 2-.

Lemma fftK n (w : F) p :
2%:R != 0 → size p ≤ 2 ^ n → (2 ^ n).-primitive_root w →
ifft n w (ftt n w p) = p.

# 3 Deriving a tactic

One application of Fourier algorithm is to implement polynomial multiplication. If we take two polynomials $p$ $q$ of degree $n$, a naive way of multiplying them requires a quadratic number $n^2$ of multiplications, i.e. we perform the multiplication of each term of $p$ (there are at most $n$ of them) by a term of $q$ (again $n$ of those). It is possible to get a faster algorithm using a fast Fourier transform. The idea is to take an alternative representation for polynomials : a polynomial $p$ of degree $n$ is uniquely defined by its value $p(x_1), \dots, p(x_n)$ on $n$
# 5 Conclusion

We have presented our formalisation of the fast Fourier transform. It makes used intensively of the polynomial library and the big operators [1]. The complete source code is about 600 lines. It is available at

https://github.com/thery/mathcomp-extra/blob/master/fft.v

# References

[1] Yves Bertot, Georges Gonthier, Sidi Ould Biha, and Ioana Paşca. Canonical Big Operators. In *TPHOLs*, volume 5170 of *LNCS*, Montreal, Canada, 2008.

[2] Venanzio Capretta. Certifying the Fast Fourier Transform with Coq. In *TPHOLs'01*, number 2152 in LNCS, pages 169–184, Edinburgh, Scotland, 2001.
distinct points $x_1, \dots, x_n$. This representation is what is called polynomial interpolation and the $x_i$s are called the interpolation points. A direct way to effectively construct a polynomial $p$ from its evaluation $p(x_1), \dots, p(x_n)$ is by using Lagrange polynomial $L_{n,k} = \prod_{i \neq k} (x-x_i)/(x_k-x_i)$. It is easy to check that $L_{n,k}(x_j) = 1$ if $j=k$ or $j=0$ otherwise. Then, we have $p = \sum_{1 \le i \le k} p(x_i)L_{n,i}$.

Computing the multiplication in the interpolation world is easy. The polynomial $pq$ is of degree at most $2n$ so if we have the values of $p$ and $q$ at $2n$ points $x_1, \dots, x_{2n}$, we simply need to perform $2n$ multiplication to get $p(x_1)q(x_1), \dots, p(x_{2n})q(x_{2n})$. Of course, we need algorithms to move from the usual polynomial representation to the interpolation one and back. This is what the fast Fourier transform gives us. The trick is that we are going to evaluate the polynomials at very specific points of interpolation (the $x_i$s are roots of unit). What makes it work is that, for those points, the potential $n^2$ values of the different $(x_i)^j$ ($1 \le i, j \le n$) consist of only $n$ distinct values, so computing the $p(x_1), \dots, p(x_n)$ can be done very efficiently in $n\log(n)$.

Formally this gives

**Definition poly\_mul (p₁ p₂ : {poly R}) :=**
($\text{\texttt{poly\_}}(i < \text{minn (size p₁) (size p₂)}) (p₁' - i * p₂' - i)$).

**Definition fft\_mul n w (p₁ p₂ : {poly F}) :=**
let $p'_1 := \text{fft n w p₁ in}$
let $p'_2 := \text{fft n w p₂ in}$
$\text{ifft n w (poly\_mul } p'_1 p'_2)$.

with its associated correctness theorem.

**Lemma fft\_mul\_correct n (p₁ p₂ : {poly F}) w :**
$2\%:R != 0 \Rightarrow (\text{size } p₁ + \text{size } p₂) \cdot -1 \le 2 \sim n \Rightarrow (2 \sim n) \cdot \text{-primitive\_root } w \Rightarrow \text{fft\_mul n w } p₁ p₂ = p₁ * p₂$.

In RocQ there is a reflexive tactic ring tactic that solves equations in a commutative ring structure [6]. For example, let us consider the equation $(a-b)*(a+b) = a^2-b^2$. The tactic first subtracts the two components of the equation $(a-b)*(a+b)-(a^2-b^2)$ and computes a normal form by developing products and simplifying terms. If the normal form is 0, the initial equation holds. We have developed an alternative tactic that works only for univariate polynomials over an arbitrary ring. It can then be used to prove equations
over univariate polynomials. It works similarly by computing a norm form
but using the Fourier transform when developing products. More precisely,
we use the following datatype to represent object on which to compute the
normal form.

Inductive pexpr :=
  | Pcons (a : nat) | Pxn (n : nat) | Px
  | Popp (p : pexpr) | Padd (p₁ p₂ : pexpr) | Pmult (p₁ p₂ : pexpr).

with the interpretation defined as follows.

Fixpoint pexpr2poly (e : pexpr) : {poly R} :=
match e with
  | Pcons a ⇒ a%:R
  | Pxn n ⇒ 'X ^ n
  | Px ⇒ 'X
  | Popp p ⇒ — (pexpr2poly p)
  | Padd p₁ p₂ ⇒ pexpr2poly p₁ + pexpr2poly p₂
  | Pmult p₁ p₂ ⇒ pexpr2poly p₁ * pexpr2poly p₂
end.

Note that we limit coefficients to be natural numbers. This means that
we can inject an element of *pexpr* over any ring *R* into a polynomial with
integer coefficients. In order to use Fourier multiplication, we need a field
with primitive roots. So, we are actually going to perform the computation
on *C*. More precisely, we use the capability of the library *COQINTERVAL* [5]
to compute with real numbers. Real numbers are axiomatised in *ROCQ*. So
we cannot directly compute with them as we do for example with integer
numbers. For example, *PI* is defined as 2*p* where *p* is the value between 7/8
and 7/4 such that cos *p* = 0. The *COQINTERVAL* library provides for any real
expression *e* an interval of rational numbers (*r*₁, *r*₂) such that *r*₁ ≤ *e* ≤ *r*₂.
If we apply it with the constant *PI* and the default precision, it gives the
following interval (7074237752028440 * 2*⁻⁵¹, 7074237752028441 * 2*⁻⁵¹). For
example, it can then be used to prove automatically standard approximations
of *PI* like the following.
## 4 The algorithm and floating-point numbers

The previous section was a toy application of the Fourier transform. Here, we are interested in more realistic applications of the Fourier transform. A typical way of computing Fourier transform is by using floating-point numbers. As the computation is not exact, it is important to quantify the error. In this section, we explain how this has been done formally. For this, we follow the pen-and-paper proof given in [3]. Our formalisation is based on the library FLOCQ[1]. In this library, the floating-point numbers are a subset of the real numbers. There is a rounding function RN that is parametrized by a rounding mode r, a basis β, and a format f. In the following, when not stating, the rounding mode is “round to the nearest”, the basis 2, and a format corresponding to the IEEE 754 standard but with no upper bound on the exponent (no overflow). This corresponds to the setting of our reference paper [3].

Prior to tackling the algorithm, we first need to quantify the error made by complex multiplications. Two types of multiplication are used one without fused multiply add (FMA) and one without. Here is the code of the basic one without FMA.

Definition brent_imult ($z_1$ $z_2$ : C) : C :=
let ($x_1, y_1$) := $z_1$ in
let ($x_2, y_2$) := $z_2$ in
(RN (RN($x_1 * z_2$) - RN($y_1 * y_2$)), RN (RN ($x_1 * y_2$) + RN ($y_1 * z_2$))).

In ROCQ, complex numbers are implemented as a pair of real numbers. In order to multiply to complex numbers $z_1$ and $z_2$, we first extract their real parts and imaginary parts and then compose the result adding a rounding operation after each arithmetic operation. To establish a bound on the error we follow the pen-and-paper proof given in [2].

Lemma error_brent_imult $z_1$ $z_2$ :
$2 \hat{=} 5 * u \le 1 \Rightarrow C \mod (brent\_imult $z_1$ $z_2$ - $z_1$ * $z_2$) \le \sqrt{5} * u * C \mod ($z_1$ * $z_2$).

where $p$ is the precision (the number of bits of the mantissa), $u$ the unit in the last place ($u = \beta^{1-p}/2$), and $C \mod$ the usual modulus for complex numbers.
For the multiplication using FMA, it is derived from an algorithm given by Kahan to compute the 2x2 determinant [7]. Here is the algorithm that uses 2 FMAs :

Definition kahan a b c d :=
    let w := RN (b * c) in
    let e := RN (w - b * c) in
    let f := RN (a * d - w) in
    RN (f + e).

and the corresponding error as given in [7].

Lemma kahan_bound :
even beta ⇒ Rabs (kahan a b c d - (a * d - b * c)) ≤ 2 * u * Rabs (a * d - b * c).

Note that this theorem is valid only if the basis $\beta$ is even.

Once we get the elementary block proved, it is straightforward to follow
the proof given in [3].

Lemma cnorm_floatDexact fma n (w : C) (l : list C) :
let omega k := u + gk n.-1 w fma k * (1 + u) in
size l = (2 ^ n) ⇒
(forall i : nat, (i < 2 ^ n) ⇒ iformat (get l i)) limp
(w ^ 2 ^ n.-1) = 1 ⇒
(fma = false ⇒ 2 ^ 5 * u ≤ 1) ⇒
cnorm n (csub n (fft_float fma n w l) (fft_exact n w l)) ≤
cnorm n (fft_exact n w l) * ((\prod_(i < n) (1 + omega i) - 1)).

Our theorem is very close to the one of the pen-and-paper proof [3].
**Theorem 7.** Assume radix-2, precision-p arithmetic, with rounding unit $u = 2^{-p}$. Let $\hat{Z}$ be the computed $2^n$-point FFT of $z \in \mathbb{C}^{2^n}$, and let Z be the exact value. Then

$$ \|\hat{Z} - Z\|_2 \le \|Z\|_2 \cdot \left( \prod_{i=1}^{n} (1 + \Omega_i) - 1 \right), $$

with

$$ \Omega_k = u + g_k(1+u), $$

$$ g_k = \begin{cases} 0 & \text{if } k = 1, 2, \\ \Delta_k^{\max} + \rho_{\times} \cdot (1 + \Delta_k^{\max}) & \text{otherwise,} \end{cases} $$

$$ \Delta_k^{\max} = \max_{\{\omega \ 2^k\text{-th root of } 1\}} \Delta_{\omega}, $$

$$ \Delta_{\omega} = |\hat{\omega} - \omega|, $$

$$ \rho_{\times} = \begin{cases} u\sqrt{5} & \text{if (10) is used (naive multiplication),} \\ 2u & \text{if (11) is used (multiplication with FMA).} \end{cases} $$

$$ \rho_{\times} = \begin{cases} u\sqrt{5} & \text{if (10) is used (naive multiplication),} \\ 2u & \text{if (11) is used (multiplication with FMA).} \end{cases} $$

$$ \rho_{\times} = \begin{cases} u\sqrt{5} & \text{if (10) is used (naive multiplication),} \\ 2u & \text{if (11) is used (multiplication with FMA).} \end{cases} $$

Only the condition $2^5 * u \le 1$ is mentioned explicitly since it comes from [2].

# 5 Conclusion

We have presented our formalisation of the Fast Fourier transform. The complete source code is about 11000 lines. It is available at

https://github.com/thery/mathcomp-extra/blob/master/fft.v

We hope that what we have presented shows the versatility of what can be done in a proof assistant. In the first part, the formalisation of the algorithm itself, we believe it illustrates the benefit of having the Mathematical component library that contains algebraic structures, big operators and polynomials. This makes the formalisation much easier compared to [4]. The second part with the toy tactic was mostly motivated by our eagerness to execute our algorithm on actual values. The tactic per se is not really of practical use. The multiplication using Fourier transform outperforms the naive multiplication only when the number of coefficients is quite high. This is very unlikely to occur while verifying equation over univariate polynomials. Finally, the application for floating-point numbers opens the opportunity to verify actual code that makes use of the FFT.
References

[1] Sylvie Boldo and Guillaume Melquiond. Flocq: A Unified Library for Proving Floating-Point Algorithms in Coq. In *20th IEEE Symposium on Computer Arithmetic*, pages 243–252. IEEE Computer Society, 2011.

[2] Richard P. Brent, Colin Percival, and Paul Zimmermann. Error bounds on complex floating-point multiplication. *Mathematics of Computation*, 76(259):1469–1481, 2007.

[3] Nicolas Brisebarre, Mioara Joldes, Jean-Michel Muller, Ana-Maria Nanes, and Joris Picot. Error analysis of some operations involved in the cooley-tukey fast fourier transform. *ACM TOMS*, 46(2):11:1–11:27, 2020.

[4] Venanzio Capretta. Certifying the Fast Fourier Transform with Coq. In *TPHOLs’01*, number 2152 in LNCS, pages 169–184, Edinburgh, Scotland, 2001.

[5] Paul Geneau de Lamarlière and Guillaume Melquiond. Coqinterval, November 2024.

[6] Benjamin Grégoire and Assia Mahboubi. Proving Equalities in a Commutative Ring Done Right in Coq. In *TPHOLs*, volume 3603 of *Lecture Notes in Computer Science*, pages 98–113. Springer, 2005.

[7] Claude-Pierre Jeannerod, Nicolas Louvet, and Jean-Michel Muller. Further analysis of kahan's algorithm for the accurate computation of $2 \times 2$ determinants. *Mathematics of Computation*, 82(284):2245–2264, 2013.