Skip to content

Add product DFTs to gr_dft; make acb_dft a thin wrapper #2774

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
fredrik-johansson wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Add product DFTs to gr_dft; make acb_dft a thin wrapper #2774

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
41881b3
Add product DFTs to gr_dft; make acb_dft a thin wrapper around gr_dft…
fredrik-johansson Jul 9, 2026
ea153f8
Fixes
fredrik-johansson Jul 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -262,7 +262,7 @@ set(_BUILD_DIRS

gr gr_generic gr_vec gr_mat
gr_poly gr_mpoly gr_ore_poly gr_series
gr_special
gr_special gr_dft

calcium
fmpz_mpoly_q
Expand Down
234 changes: 16 additions & 218 deletions doc/source/acb_dft.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,224 +81,22 @@ precomputed schemes then return results accurate to about *prec* bits.

Compute the inverse DFT of *v* into *w*.

DFT on products
These functions are thin wrappers around the :doc:`gr_dft <gr_dft>`
module, which computes DFTs over generic rings. For complex balls
the transform is carried out in fixed-point arithmetic with
rigorous error bounds whenever the input permits, falling back to
ball arithmetic otherwise, and the precomputation object is the
:type:`gr_dft_acb_pre_t` plan itself.

Obsolete functions
-------------------------------------------------------------------------------

A finite abelian group is isomorphic to a product of cyclic components

.. math::

G = \bigoplus_{i=1}^r \mathbb Z/n_i\mathbb Z

Characters are product of component characters and the DFT reads

.. math::

\hat f(x_1,\dots x_r) = \sum_{y_1\dots y_r} f(y_1,\dots y_r)
e^{-2i \pi \sum\frac{x_i y_i}{n_i}}

We assume that `f` is given by a vector of length `\prod n_i` corresponding
to a lexicographic ordering of the values `y_1,\dots y_r`, and the computation
returns the same indexing for values of `\hat f`.

.. function:: void acb_dft_prod(acb_ptr w, acb_srcptr v, slong * cyc, slong num, slong prec)

Computes the DFT on the group product of *num* cyclic components of sizes *cyc*. Assume the entries
of *v* are indexed according to lexicographic ordering of the cyclic
components.

.. type:: acb_dft_prod_struct

.. type:: acb_dft_prod_t

Stores a fast DFT scheme on a product of cyclic groups.

An *acb_dft_prod_t* is defined as an array of *acb_dft_prod_struct*
of length 1, permitting it to be passed by reference.

.. function:: void acb_dft_prod_init(acb_dft_prod_t t, slong * cyc, slong num, slong prec)

Stores in *t* a DFT scheme for the product of *num* cyclic components whose sizes are given in the array *cyc*.

.. function:: void acb_dft_prod_clear(acb_dft_prod_t t)

Clears *t*.

.. function:: void acb_dft_prod_precomp(acb_ptr w, acb_srcptr v, const acb_dft_prod_t prod, slong prec)

Sets *w* to the DFT of *v*. Assume the entries are lexicographically
ordered according to the product of cyclic groups initialized in *prod*.

Convolution
-------------------------------------------------------------------------------

For functions `f` and `g` on `G` we consider the convolution

.. math::

(f \star g)(x) = \sum_{y\in G} f(x-y)g(y)

.. function:: void acb_dft_convol_naive(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec)

.. function:: void acb_dft_convol_rad2(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec)

.. function:: void acb_dft_convol(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec)

Sets *w* to the convolution of *f* and *g* of length *len*.

The *naive* version simply uses the definition.

The *rad2* version embeds the sequence into a power of 2 length and
uses the formula

.. math::

\widehat{f \star g}(\chi) = \hat f(\chi)\hat g(\chi)

to compute it using three radix 2 FFT.

The default version uses radix 2 FFT unless *len* is a product of small
primes where a non padded FFT is faster.

FFT algorithms
-------------------------------------------------------------------------------

Fast Fourier transform techniques allow to compute efficiently
all values `\hat f(\chi)` by reusing common computations.

Specifically, if `H\triangleleft G` is a subgroup of size `M` and index
`[G:H]=m`, then writing `f_x(h)=f(xh)` the translate of `f` by representatives
`x` of `G/H`, one has a decomposition

.. math::

\hat f(\chi) = \sum_{x\in G/H} \overline{\chi(x)} \hat{f_x}(\chi_{H})

so that the DFT on `G` can be computed using `m` DFT on `H` (of
appropriate translates of `f`), then `M` DFT on `G/H`, one for
each restriction `\chi_{H}`.

This decomposition can be done recursively.

Naive algorithm
...............................................................................

.. function:: void acb_dft_naive(acb_ptr w, acb_srcptr v, slong n, slong prec)

Computes the DFT of *v* into *w*, where *v* and *w* have size *n*,
using the naive `O(n^2)` algorithm.

.. type:: acb_dft_naive_struct

.. type:: acb_dft_naive_t

.. function:: void acb_dft_naive_init(acb_dft_naive_t t, slong len, slong prec)

.. function:: void acb_dft_naive_clear(acb_dft_naive_t t)

Stores a table of roots of unity in *t*.
The length *len* is stored as *t->n*.

.. function:: void acb_dft_naive_precomp(acb_ptr w, acb_srcptr v, const acb_dft_naive_t t, slong prec)

Sets *w* to the DFT of *v* of size *t->n*, using the naive algorithm data *t*.

CRT decomposition
...............................................................................

.. function:: void acb_dft_crt(acb_ptr w, acb_srcptr v, slong n, slong prec)

Computes the DFT of *v* into *w*, where *v* and *w* have size *len*,
using CRT to express `\mathbb Z/n\mathbb Z` as a product of cyclic groups.

.. type:: acb_dft_crt_struct

.. type:: acb_dft_crt_t

.. function:: void acb_dft_crt_init(acb_dft_crt_t t, slong len, slong prec)

.. function:: void acb_dft_crt_clear(acb_dft_crt_t t)

Initialize a CRT decomposition of `\mathbb Z/n\mathbb Z` as a direct product
of cyclic groups.
The length *len* is stored as *t->n*.

.. function:: void acb_dft_crt_precomp(acb_ptr w, acb_srcptr v, const acb_dft_crt_t t, slong prec)

Sets *w* to the DFT of *v* of size *t->n*, using the CRT decomposition scheme *t*.

Cooley-Tukey decomposition
...............................................................................

.. function:: void acb_dft_cyc(acb_ptr w, acb_srcptr v, slong n, slong prec)

Computes the DFT of *v* into *w*, where *v* and *w* have size *n*,
using each prime factor of `m` of `n` to decompose with
the subgroup `H=m\mathbb Z/n\mathbb Z`.

.. type:: acb_dft_cyc_struct

.. type:: acb_dft_cyc_t

.. function:: void acb_dft_cyc_init(acb_dft_cyc_t t, slong len, slong prec)

.. function:: void acb_dft_cyc_clear(acb_dft_cyc_t t)

Initialize a decomposition of `\mathbb Z/n\mathbb Z` into cyclic subgroups.
The length *len* is stored as *t->n*.

.. function:: void acb_dft_cyc_precomp(acb_ptr w, acb_srcptr v, const acb_dft_cyc_t t, slong prec)

Sets *w* to the DFT of *v* of size *t->n*, using the cyclic decomposition scheme *t*.

Radix 2 decomposition
...............................................................................

.. function:: void acb_dft_rad2(acb_ptr w, acb_srcptr v, int e, slong prec)

Computes the DFT of *v* into *w*, where *v* and *w* have size `2^e`,
using a radix 2 FFT.

.. function:: void acb_dft_inverse_rad2_precomp_inplace(acb_ptr v, const acb_dft_rad2_t t, slong prec)

Computes the inverse DFT of *v* in place, of size *t->n*, using the
precomputed radix 2 scheme *t*.

.. type:: acb_dft_rad2_struct

.. type:: acb_dft_rad2_t

.. function:: void acb_dft_rad2_init(acb_dft_rad2_t t, int e, slong prec)

.. function:: void acb_dft_rad2_clear(acb_dft_rad2_t t)

Initialize and clear a radix 2 FFT of size `2^e`, stored as *t->n*.

.. function:: void acb_dft_rad2_precomp(acb_ptr w, acb_srcptr v, const acb_dft_rad2_t t, slong prec)

Sets *w* to the DFT of *v* of size *t->n*, using the precomputed radix 2 scheme *t*.

Bluestein transform
...............................................................................

.. function:: void acb_dft_bluestein(acb_ptr w, acb_srcptr v, slong n, slong prec)

Computes the DFT of *v* into *w*, where *v* and *w* have size *n*,
by conversion to a radix 2 one using Bluestein's convolution trick.

.. type:: acb_dft_bluestein_struct

.. type:: acb_dft_bluestein_t

Stores a Bluestein scheme for some length *n* : that is a :type:`acb_dft_rad2_t` of size
`2^e \geq 2n-1` and a size *n* array of convolution factors.

.. function:: void acb_dft_bluestein_init(acb_dft_bluestein_t t, slong len, slong prec)
void acb_dft_bluestein_clear(acb_dft_bluestein_t t)

Initialize and clear a Bluestein scheme to compute DFT of size *len*.

.. function:: void acb_dft_bluestein_precomp(acb_ptr w, acb_srcptr v, const acb_dft_bluestein_t t, slong prec)

Sets *w* to the DFT of *v* of size *t->n*, using the precomputed Bluestein scheme *t*.
The remaining functions of this module (product DFTs, convolutions,
and direct access to the naive, CRT, cyclic, radix-2 and Bluestein
algorithms) have been removed. Product DFTs with complex ball input
and output are provided by :func:`gr_dft_acb_prod`; the individual
algorithms, transforms over other rings and further functionality are
available in the :doc:`gr_dft <gr_dft>` module, and cyclic
convolutions are easily expressed through forward and inverse
transforms.

10 changes: 8 additions & 2 deletions doc/source/acb_dirichlet.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -471,7 +471,7 @@ to consider the Fourier transform on Conrey labels as
g(a) = \sum_{b\bmod q}\overline{\chi_q(a,b)}f(b)


.. function:: void acb_dirichlet_dft_conrey(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec)
.. function:: void acb_dirichlet_dft_index(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec)

Compute the DFT of *v* using Conrey indices.
This function assumes *v* and *w* are vectors
Expand Down Expand Up @@ -503,6 +503,12 @@ to consider the Fourier transform on Conrey labels as
This function assumes *v* and *w* are vectors of size *G->q*.
All values at index not coprime to *G->q* are ignored.

Both functions are thin wrappers around the :doc:`gr_dft <gr_dft>`
module (see :func:`gr_dft_dirichlet` and
:func:`gr_dft_dirichlet_index`), which computes the product DFT
over the Conrey components in fixed-point arithmetic with rigorous
error bounds whenever possible.

Dirichlet L-functions
-------------------------------------------------------------------------------

Expand Down Expand Up @@ -581,7 +587,7 @@ Dirichlet L-functions
Hurwitz zeta function and a discrete Fourier transform.
The output *res* is assumed to have length *G->phi_q* and values
are stored by lexicographically ordered
Conrey logs. See :func:`acb_dirichlet_dft_conrey`.
Conrey logs. See :func:`acb_dirichlet_dft_index`.

If *precomp* is *NULL*, each Hurwitz zeta function value is computed
directly. If a pre-initialized *precomp* object is provided, this will be
Expand Down
Loading
Loading