From 41881b3fa2f64af7c343f23152d71673298ed141 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Thu, 9 Jul 2026 17:58:06 +0200 Subject: [PATCH 1/2] Add product DFTs to gr_dft; make acb_dft a thin wrapper around gr_dft and remove obsolete algorithm functions --- doc/source/acb_dft.rst | 234 +------------- doc/source/acb_dirichlet.rst | 10 +- doc/source/gr_dft.rst | 73 +++++ src/acb_dft.h | 297 +----------------- src/acb_dft/bluestein.c | 154 --------- src/acb_dft/convol.c | 35 --- src/acb_dft/convol_dft.c | 40 --- src/acb_dft/convol_mullow.c | 31 -- src/acb_dft/convol_naive.c | 30 -- src/acb_dft/convol_rad2.c | 87 ------ src/acb_dft/crt.c | 210 ------------- src/acb_dft/cyc.c | 121 -------- src/acb_dft/impl.h | 19 -- src/acb_dft/naive.c | 95 ------ src/acb_dft/precomp.c | 157 +++------- src/acb_dft/prod.c | 65 ---- src/acb_dft/profile/p-convol.c | 136 -------- src/acb_dft/profile/p-dft.c | 145 --------- src/acb_dft/rad2.c | 121 -------- src/acb_dft/rad2_threaded.c | 129 -------- src/acb_dft/step.c | 84 ----- src/acb_dft/test/main.c | 2 - src/acb_dft/test/t-accuracy.c | 256 +++++---------- src/acb_dft/test/t-convol.c | 135 -------- src/acb_dft/test/t-dft.c | 220 +++++-------- src/acb_dirichlet/dft.c | 59 +--- src/acb_dirichlet/platt_multieval.c | 9 + src/acb_dirichlet/test/t-dft.c | 126 ++++---- src/acb_theta/ql_jet_fd.c | 8 +- src/gr_dft.h | 74 +++++ src/gr_dft/acb.c | 413 +++++++++++++++++++++--- src/gr_dft/dirichlet.c | 149 +++++++++ src/gr_dft/prod.c | 466 ++++++++++++++++++++++++++++ src/gr_dft/profile/p-gr_dft_acb.c | 80 +++-- src/gr_dft/test/main.c | 4 + src/gr_dft/test/t-acb_dft.c | 11 +- src/gr_dft/test/t-dirichlet_dft.c | 224 +++++++++++++ src/gr_dft/test/t-prod.c | 244 +++++++++++++++ 38 files changed, 1985 insertions(+), 2768 deletions(-) delete mode 100644 src/acb_dft/bluestein.c delete mode 100644 src/acb_dft/convol.c delete mode 100644 src/acb_dft/convol_dft.c delete mode 100644 src/acb_dft/convol_mullow.c delete mode 100644 src/acb_dft/convol_naive.c delete mode 100644 src/acb_dft/convol_rad2.c delete mode 100644 src/acb_dft/crt.c delete mode 100644 src/acb_dft/cyc.c delete mode 100644 src/acb_dft/impl.h delete mode 100644 src/acb_dft/naive.c delete mode 100644 src/acb_dft/prod.c delete mode 100644 src/acb_dft/profile/p-convol.c delete mode 100644 src/acb_dft/profile/p-dft.c delete mode 100644 src/acb_dft/rad2.c delete mode 100644 src/acb_dft/rad2_threaded.c delete mode 100644 src/acb_dft/step.c delete mode 100644 src/acb_dft/test/t-convol.c create mode 100644 src/gr_dft/dirichlet.c create mode 100644 src/gr_dft/prod.c create mode 100644 src/gr_dft/test/t-dirichlet_dft.c create mode 100644 src/gr_dft/test/t-prod.c diff --git a/doc/source/acb_dft.rst b/doc/source/acb_dft.rst index 3f42b12b5a..2ca224065b 100644 --- a/doc/source/acb_dft.rst +++ b/doc/source/acb_dft.rst @@ -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 ` + 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 ` module, and cyclic +convolutions are easily expressed through forward and inverse +transforms. diff --git a/doc/source/acb_dirichlet.rst b/doc/source/acb_dirichlet.rst index b7ce603a53..eeac807d83 100644 --- a/doc/source/acb_dirichlet.rst +++ b/doc/source/acb_dirichlet.rst @@ -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 @@ -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 ` + 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 ------------------------------------------------------------------------------- @@ -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 diff --git a/doc/source/gr_dft.rst b/doc/source/gr_dft.rst index 0d145fbdf7..f4a81c72d5 100644 --- a/doc/source/gr_dft.rst +++ b/doc/source/gr_dft.rst @@ -704,3 +704,76 @@ Internal functions Bluestein chirp-z transform (``GR_DFT_ALG_BLUESTEIN``). *res* must not alias *vec*. The inverse transform is computed as the forward transform of the cyclically reversed input. + +DFT on products of cyclic groups +------------------------------------------------------------------------------- + +The following functions compute the DFT on a product group +`\mathbb{Z}/c_0 \times \cdots \times \mathbb{Z}/c_{k-1}` with row-major +indexing (the first component varies slowest), following the same +algorithm as :func:`acb_dft_prod`: one cyclic DFT along each axis, +with no twiddle factors in between. Component plans are shared +between axes of equal length. The lines of one axis are independent +and are transformed in parallel when multiple threads are available +and the ring is marked thread-safe. + +.. function:: int gr_dft_prod_precomp_init(gr_dft_prod_pre_t P, const ulong * cyc, slong num, int flags, gr_ctx_t ctx) + int gr_dft_prod_precomp_init_root(gr_dft_prod_pre_t P, gr_srcptr w, ulong order, const ulong * cyc, slong num, int flags, gr_ctx_t ctx) + + Initializes a plan for the product DFT with component lengths + *cyc*. The first version uses canonical roots of unity + (:func:`gr_dft_default_root`); the second takes a root of unity + *w* of the given *order*, which every component length must + divide, the component of length `m` using `w^{\mathrm{order}/m}`. + Any *flags* are passed on to the component plans. + +.. function:: void gr_dft_prod_precomp_clear(gr_dft_prod_pre_t P) + +.. function:: int gr_dft_prod_precomp(gr_ptr res, gr_srcptr vec, const gr_dft_prod_pre_t P, gr_ctx_t ctx) + int gr_dft_prod_inverse_precomp(gr_ptr res, gr_srcptr vec, const gr_dft_prod_pre_t P, gr_ctx_t ctx) + int gr_dft_prod(gr_ptr res, gr_srcptr vec, const ulong * cyc, slong num, gr_ctx_t ctx) + int gr_dft_prod_inverse(gr_ptr res, gr_srcptr vec, const ulong * cyc, slong num, gr_ctx_t ctx) + + Product DFT and its inverse (the latter including the `1/n` + normalization, applied once at the end rather than per axis). + Aliasing of *res* and *vec* is allowed. + +.. function:: void gr_dft_prod_precomp_nfixed_bound(double * peak, double * err_ulps, double in_mag, double in_err, const gr_dft_prod_pre_t P) + + Composes the fixed-point error bounds of the component + transforms, axis by axis. + +.. function:: void gr_dft_acb_prod(acb_ptr w, acb_srcptr v, const ulong * cyc, slong num, slong prec) + void gr_dft_acb_prod_inverse(acb_ptr w, acb_srcptr v, const ulong * cyc, slong num, slong prec) + int gr_dft_acb_prod_precomp_init(gr_dft_acb_prod_pre_t Q, const ulong * cyc, slong num, slong prec) + void gr_dft_acb_prod_precomp_clear(gr_dft_acb_prod_pre_t Q) + void gr_dft_acb_prod_precomp(acb_ptr w, acb_srcptr v, const gr_dft_acb_prod_pre_t Q, slong prec) + void gr_dft_acb_prod_inverse_precomp(acb_ptr w, acb_srcptr v, const gr_dft_acb_prod_pre_t Q, slong prec) + + Product DFT with complex ball input and output, using fixed-point + arithmetic internally with rigorous error bounds (the same + scaling and error analysis as :func:`gr_dft_acb`, with the + composed product bound), falling back to ball arithmetic when + fixed point does not apply. + +DFT on Dirichlet groups +------------------------------------------------------------------------------- + +Counterparts of :func:`acb_dirichlet_dft` and +:func:`acb_dirichlet_dft_index`: the group of Dirichlet characters +mod `q` is a product of cyclic groups, so its DFT is a product DFT +over the Conrey component sizes. + +.. function:: int gr_dft_dirichlet_index(gr_ptr w, gr_srcptr v, const dirichlet_group_t G, gr_ctx_t ctx) + int gr_dft_dirichlet(gr_ptr w, gr_srcptr v, const dirichlet_group_t G, gr_ctx_t ctx) + + DFT of *v* over the Dirichlet group *G*, in lexicographic Conrey + indexing (array size ``G->phi_q``) and number indexing (array + size ``G->q``) respectively, over an arbitrary ring with the + required roots of unity. + +.. function:: void gr_dft_acb_dirichlet_index(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec) + void gr_dft_acb_dirichlet(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec) + + The same transforms for complex ball input and output, routed + through the fixed-point product transform. diff --git a/src/acb_dft.h b/src/acb_dft.h index b987591735..02d5860ba2 100644 --- a/src/acb_dft.h +++ b/src/acb_dft.h @@ -1,5 +1,6 @@ /* Copyright (C) 2016 Pascal Molin + Copyright (C) 2026 Fredrik Johansson This file is part of FLINT. @@ -12,300 +13,30 @@ #ifndef ACB_DFT_H #define ACB_DFT_H -#include "limb_types.h" -#include "acb.h" +#include "acb_types.h" +#include "gr_dft.h" #ifdef __cplusplus extern "C" { #endif -void _acb_dft_naive(acb_ptr w, acb_srcptr v, slong dv, acb_srcptr z, slong dz, slong len, slong prec); -void acb_dft_naive(acb_ptr w, acb_srcptr v, slong len, slong prec); -void acb_dft_crt(acb_ptr w, acb_srcptr v, slong len, slong prec); -void acb_dft_cyc(acb_ptr w, acb_srcptr v, slong len, slong prec); -void acb_dft_inverse_cyc(acb_ptr w, acb_srcptr v, slong len, slong prec); -void acb_dft_rad2_inplace(acb_ptr v, int e, slong prec); -void acb_dft_rad2(acb_ptr w, acb_srcptr v, int e, slong prec); -void acb_dft_bluestein(acb_ptr w, acb_srcptr v, slong len, slong prec); -void acb_dft_prod(acb_ptr w, acb_srcptr v, slong * cyc, slong num, slong prec); - -void acb_dft_rad2_inplace_threaded(acb_ptr v, int e, slong prec); - -void acb_dft_convol_naive(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec); -void acb_dft_convol_dft(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec); -void acb_dft_convol_rad2(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec); -void acb_dft_convol_mullow(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec); -void acb_dft_convol(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec); - -#define CRT_MAX 15 -typedef struct -{ - int num; - nmod_t n; - slong m[CRT_MAX]; - ulong M[CRT_MAX]; - ulong vM[CRT_MAX]; -} -crt_struct; - -typedef crt_struct crt_t[1]; - -void crt_init(crt_t c, ulong n); - -void crt_decomp(acb_ptr y, acb_srcptr x, slong dx, const crt_t c, ulong len); -void crt_recomp(acb_ptr y, acb_srcptr x, const crt_t c, ulong len); - -typedef struct acb_dft_step_struct acb_dft_step_struct; -typedef acb_dft_step_struct * acb_dft_step_ptr; - -typedef struct -{ - slong n; - acb_ptr z; - int zclear; - slong num; - acb_dft_step_ptr cyc; -} -acb_dft_cyc_struct; - -typedef acb_dft_cyc_struct acb_dft_cyc_t[1]; - -typedef struct -{ - int e; - slong n; /* = 1 << e */ - slong dv; - slong nz; /* = n but could be bigger */ - acb_ptr z; -} -acb_dft_rad2_struct; - -typedef acb_dft_rad2_struct acb_dft_rad2_t[1]; - -typedef struct -{ - slong n; - slong dv; - acb_ptr z; /* z[k] = e(k^2/2n) */ - acb_ptr g; /* g[k] = dft( z ) */ - acb_dft_rad2_t rad2; -} -acb_dft_bluestein_struct; - -typedef acb_dft_bluestein_struct acb_dft_bluestein_t[1]; - -typedef struct -{ - slong n; - slong num; - acb_dft_step_ptr cyc; -} -acb_dft_prod_struct; - -typedef acb_dft_prod_struct acb_dft_prod_t[1]; - -typedef struct -{ - slong n; - crt_t c; - slong dv; - /* then a product */ - acb_dft_step_ptr cyc; -} -acb_dft_crt_struct; - -typedef acb_dft_crt_struct acb_dft_crt_t[1]; - -typedef struct -{ - slong n; - slong dv; - int zclear; - acb_ptr z; - slong dz; -} -acb_dft_naive_struct; - -typedef acb_dft_naive_struct acb_dft_naive_t[1]; - -typedef struct -{ - slong n; - int type; - union - { - acb_dft_rad2_t rad2; - acb_dft_cyc_t cyc; - acb_dft_prod_t prod; - acb_dft_crt_t crt; - acb_dft_naive_t naive; - acb_dft_bluestein_t bluestein; - } t; -} -acb_dft_pre_struct; +/* This module is a thin wrapper around gr_dft, which computes DFTs + over generic rings and provides the complex-ball functionality + through fixed-point arithmetic with rigorous error bounds + (gr_dft_acb), falling back to ball arithmetic when fixed point + does not apply. Product DFTs, convolutions and access to specific + algorithms are available in gr_dft. */ +typedef gr_dft_acb_pre_struct acb_dft_pre_struct; typedef acb_dft_pre_struct acb_dft_pre_t[1]; -/* covers both product and cyclic case */ -struct -acb_dft_step_struct -{ - /* [G:H] */ - slong m; - /* card H */ - slong M; - slong dv; /* = M for prod, also = M if cyc is reordered */ - /* pointer on some roots of unity, if needed */ - acb_srcptr z; - /* index of mM in z */ - slong dz; - /* dft to call on H */ - acb_dft_pre_t pre; - /* dft to call on G/H ? */ -}; - -/*typedef acb_dft_pre_struct acb_dft_pre_t[1];*/ - -#define DFT_VERB 0 - -enum -{ - DFT_NAIVE, DFT_CYC, DFT_PROD, DFT_CRT , DFT_RAD2 , DFT_CONV -}; - -void acb_dft_step(acb_ptr w, acb_srcptr v, acb_dft_step_ptr cyc, slong num, slong prec); - -void acb_dft_precomp(acb_ptr w, acb_srcptr v, const acb_dft_pre_t pre, slong prec); -void acb_dft_inverse_precomp(acb_ptr w, acb_srcptr v, const acb_dft_pre_t pre, slong prec); -void acb_dft_naive_precomp(acb_ptr w, acb_srcptr v, const acb_dft_naive_t pol, slong prec); -void acb_dft_cyc_precomp(acb_ptr w, acb_srcptr v, const acb_dft_cyc_t cyc, slong prec); - -void acb_dft_rad2_precomp_inplace(acb_ptr v, const acb_dft_rad2_t rad2, slong prec); -void acb_dft_rad2_precomp(acb_ptr w, acb_srcptr v, const acb_dft_rad2_t rad2, slong prec); -void acb_dft_crt_precomp(acb_ptr w, acb_srcptr v, const acb_dft_crt_t crt, slong prec); -void acb_dft_prod_precomp(acb_ptr w, acb_srcptr v, const acb_dft_prod_t prod, slong prec); -void acb_dft_bluestein_precomp(acb_ptr w, acb_srcptr v, const acb_dft_bluestein_t t, slong prec); - -void acb_dft_rad2_precomp_inplace_threaded(acb_ptr v, const acb_dft_rad2_t rad2, slong prec); - -void acb_dft_inverse_rad2_precomp_inplace(acb_ptr v, const acb_dft_rad2_t rad2, slong prec); -void acb_dft_inverse_rad2_precomp(acb_ptr w, acb_srcptr v, const acb_dft_rad2_t rad2, slong prec); -void acb_dft_convol_rad2_precomp(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, const acb_dft_rad2_t rad2, slong prec); - -void _acb_dft_precomp_init(acb_dft_pre_t pre, slong dv, acb_ptr z, slong dz, slong len, slong prec); -void acb_dft_precomp_init(acb_dft_pre_t pre, slong len, slong prec); -void acb_dft_precomp_clear(acb_dft_pre_t pre); - void acb_dft(acb_ptr w, acb_srcptr v, slong len, slong prec); void acb_dft_inverse(acb_ptr w, acb_srcptr v, slong len, slong prec); -acb_dft_step_ptr _acb_dft_steps_prod(slong * m, slong num, slong prec); - -FLINT_FORCE_INLINE void -acb_dft_prod_init(acb_dft_prod_t t, slong * cyc, slong num, slong prec) -{ - t->num = num; - t->cyc = _acb_dft_steps_prod(cyc, num, prec); -} - -void acb_dft_prod_clear(acb_dft_prod_t t); - -void _acb_dft_cyc_init_z_fac(acb_dft_cyc_t t, n_factor_t fac, slong dv, acb_ptr z, slong dz, slong len, slong prec); -void _acb_dft_cyc_init(acb_dft_cyc_t t, slong dv, slong len, slong prec); - -FLINT_FORCE_INLINE void -acb_dft_cyc_init(acb_dft_cyc_t t, slong len, slong prec) -{ - _acb_dft_cyc_init(t, 1, len, prec); -} - -void acb_dft_cyc_clear(acb_dft_cyc_t t); - -void _acb_dft_naive_init(acb_dft_naive_t pol, slong dv, acb_ptr z, slong dz, slong len, slong prec); - -FLINT_FORCE_INLINE void -acb_dft_naive_init(acb_dft_naive_t pol, slong len, slong prec) -{ - _acb_dft_naive_init(pol, 1, NULL, 0, len, prec); -} - -FLINT_FORCE_INLINE void -acb_dft_naive_clear(acb_dft_naive_t pol) -{ - if (pol->zclear) - _acb_vec_clear(pol->z, pol->n); -} - -void _acb_dft_rad2_init(acb_dft_rad2_t t, slong dv, int e, slong prec); - -FLINT_FORCE_INLINE void -acb_dft_rad2_init(acb_dft_rad2_t t, int e, slong prec) -{ - _acb_dft_rad2_init(t, 1, e, prec); -} - -FLINT_FORCE_INLINE void -acb_dft_rad2_clear(acb_dft_rad2_t t) -{ - _acb_vec_clear(t->z, t->nz); -} - -void _acb_dft_bluestein_init(acb_dft_bluestein_t t, slong dv, slong n, slong prec); - -FLINT_FORCE_INLINE void -acb_dft_bluestein_init(acb_dft_bluestein_t t, slong n, slong prec) -{ - _acb_dft_bluestein_init(t, 1, n, prec); -} - -FLINT_FORCE_INLINE void -acb_dft_bluestein_clear(acb_dft_bluestein_t t) -{ - if (t->n != 0) - { - _acb_vec_clear(t->z, t->n); - _acb_vec_clear(t->g, t->rad2->n); - acb_dft_rad2_clear(t->rad2); - } -} - -void _acb_dft_crt_init(acb_dft_crt_t crt, slong dv, slong len, slong prec); -void acb_dft_crt_init(acb_dft_crt_t crt, slong len, slong prec); -void acb_dft_crt_clear(acb_dft_crt_t crt); - -/* utils, could be moved elsewhere */ - -FLINT_FORCE_INLINE void -acb_swap_ri(acb_t x) -{ - arb_swap(acb_realref(x), acb_imagref(x)); -} - -FLINT_FORCE_INLINE void -acb_vec_swap_ri(acb_ptr v, slong len) -{ - slong k; - for (k = 0; k < len; k++) - acb_swap_ri(v + k); -} - -FLINT_FORCE_INLINE void -_acb_vec_kronecker_mul(acb_ptr z, acb_srcptr x, acb_srcptr y, slong len, slong prec) -{ - slong k; - for (k = 0; k < len; k++) - acb_mul(z + k, x + k, y + k, prec); -} - -FLINT_FORCE_INLINE void -_acb_vec_kronecker_mul_step(acb_ptr z, acb_srcptr x, acb_srcptr y, slong step, slong len, slong prec) -{ - slong k; - for (k = 0; k < len; k++) - { - acb_mul(z + k, x + k, y + k * step, prec); - } -} +void acb_dft_precomp_init(acb_dft_pre_t pre, slong len, slong prec); +void acb_dft_precomp_clear(acb_dft_pre_t pre); +void acb_dft_precomp(acb_ptr w, acb_srcptr v, const acb_dft_pre_t pre, slong prec); +void acb_dft_inverse_precomp(acb_ptr w, acb_srcptr v, const acb_dft_pre_t pre, slong prec); #ifdef __cplusplus } diff --git a/src/acb_dft/bluestein.c b/src/acb_dft/bluestein.c deleted file mode 100644 index b8c5bddfbb..0000000000 --- a/src/acb_dft/bluestein.c +++ /dev/null @@ -1,154 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "nmod.h" -#include "acb_dft.h" -#include "acb_modular.h" - -/* z[k] = z^(k^2), z a 2n-th root of unity */ -static void -_acb_vec_bluestein_factors(acb_ptr z, slong n, slong prec) -{ - /* this function is used mostly with prime-power n - * so the set of squares has index 2 only. - * computing an addition sequence does not considerably improve things */ - if (n < 30) - { - slong k, k2; - acb_ptr z2n; - nmod_t n2; - - z2n = _acb_vec_init(2 * n); - _acb_vec_unit_roots(z2n, -2 * n, 2 * n, prec); - nmod_init(&n2, FLINT_MAX(2 * n, 1)); - - for (k = 0, k2 = 0; k < n; k++) - { - acb_set(z + k, z2n + k2); - k2 = nmod_add(k2, 2 * k + 1, n2); - } - - _acb_vec_clear(z2n, 2 * n); - } - else - { - nmod_t n2; - slong k, k2, dk, wp; - slong * v, * s; - acb_ptr t; - s = flint_malloc(n * sizeof(slong)); - v = flint_malloc((n + 1)* sizeof(slong)); - t = _acb_vec_init(n + 1); - nmod_init(&n2, 2 * n); - - for (k = 0; k < n; k++) - v[k] = 0; - for (k = 0, k2 = 0, dk = 1; k < n; k++) - { - s[k] = k2; - if (k2 < n) - v[k2] = -1; - else - v[2 * n - k2] = -1; - - k2 = nmod_add(k2, dk, n2); - dk = nmod_add(dk, 2, n2); - } - acb_modular_fill_addseq(v, n); - - wp = prec + 6 + 2 * FLINT_BIT_COUNT(2 * n); - - acb_one(t + 0); - acb_unit_root(t + 1, 2 * n, wp); - acb_conj(t + 1, t + 1); - acb_set_si(t + n, -1); - for (k = 2; k < n; k++) - if (v[k]) - acb_mul(t + k, t + v[k], t + k - v[k], wp); - for (k = 0; k < n; k++) - { - if (s[k] <= n) - acb_set(z + k, t + s[k]); - else - acb_conj(z + k, t + 2 * n - s[k]); - } - _acb_vec_set_round(z, z, n, prec); - _acb_vec_clear(t, n + 1); - flint_free(s); - flint_free(v); - } -} - -void -_acb_dft_bluestein_init(acb_dft_bluestein_t t, slong dv, slong n, slong prec) -{ - acb_ptr z, g; - slong k, n2; - int e; - - t->n = n; - t->dv = dv; - - if (n == 0) - return; - - e = FLINT_CLOG2(2 * n - 1); - - if (DFT_VERB) - flint_printf("dft_bluestein: init z[2^%i]\n", e); - - acb_dft_rad2_init(t->rad2, e, prec); - - t->z = z = _acb_vec_init(n); - _acb_vec_bluestein_factors(t->z, n, prec); - - n2 = t->rad2->n; - t->g = g = _acb_vec_init(n2); - acb_one(g + 0); - for (k = 1; k < n; k++) - { - acb_conj(g + k, z + k); - acb_set(g + n2 - k, g + k); - } - acb_dft_rad2_precomp_inplace(g, t->rad2, prec); -} - -void -acb_dft_bluestein_precomp(acb_ptr w, acb_srcptr v, const acb_dft_bluestein_t t, slong prec) -{ - slong n = t->n, dv = t->dv, np; - acb_ptr fp; - - if (n == 0) - return; - - np = t->rad2->n; - fp = _acb_vec_init(np); - _acb_vec_kronecker_mul_step(fp, t->z, v, dv, n, prec); - - acb_dft_rad2_precomp_inplace(fp, t->rad2, prec); - _acb_vec_kronecker_mul(fp, t->g, fp, np, prec); - - acb_dft_inverse_rad2_precomp_inplace(fp, t->rad2, prec); - - _acb_vec_kronecker_mul(w, t->z, fp, n, prec); - - _acb_vec_clear(fp, np); -} - -void -acb_dft_bluestein(acb_ptr w, acb_srcptr v, slong len, slong prec) -{ - acb_dft_bluestein_t t; - acb_dft_bluestein_init(t, len, prec); - acb_dft_bluestein_precomp(w, v, t, prec); - acb_dft_bluestein_clear(t); -} diff --git a/src/acb_dft/convol.c b/src/acb_dft/convol.c deleted file mode 100644 index db6183b9fd..0000000000 --- a/src/acb_dft/convol.c +++ /dev/null @@ -1,35 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "acb_dft.h" - -static int use_dft(slong len, slong prec) -{ - slong l2 = len; - while (l2 >= 16) l2 >>= 1; - if (l2 < 11) - { - while (!(len & 1)) len >>= 1; - while (len % 3 == 0) len /= 3; - while (len % 5 == 0) len /= 5; - while (len % 7 == 0) len /= 7; - return (len == 1); - } - return 0; -} - -void acb_dft_convol(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec) -{ - if (use_dft(len, prec)) - acb_dft_convol_dft(w, f, g, len, prec); - else - acb_dft_convol_rad2(w, f, g, len, prec); -} diff --git a/src/acb_dft/convol_dft.c b/src/acb_dft/convol_dft.c deleted file mode 100644 index 4c20e0279f..0000000000 --- a/src/acb_dft/convol_dft.c +++ /dev/null @@ -1,40 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "acb_dft.h" - -static void -acb_dft_convol_dft_precomp(acb_ptr w, acb_srcptr f, acb_srcptr g, const acb_dft_pre_t pre, slong prec) -{ - acb_ptr fp, gp; - - fp = _acb_vec_init(pre->n); - gp = _acb_vec_init(pre->n); - - acb_dft_precomp(fp, f, pre, prec); - acb_dft_precomp(gp, g, pre, prec); - - _acb_vec_kronecker_mul(gp, gp, fp, pre->n, prec); - - acb_dft_inverse_precomp(w, gp, pre, prec); - - _acb_vec_clear(fp, pre->n); - _acb_vec_clear(gp, pre->n); -} - -void -acb_dft_convol_dft(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec) -{ - acb_dft_pre_t pre; - acb_dft_precomp_init(pre, len, prec); - acb_dft_convol_dft_precomp(w, f, g, pre, prec); - acb_dft_precomp_clear(pre); -} diff --git a/src/acb_dft/convol_mullow.c b/src/acb_dft/convol_mullow.c deleted file mode 100644 index e5e03a53f9..0000000000 --- a/src/acb_dft/convol_mullow.c +++ /dev/null @@ -1,31 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "acb_dft.h" -#include "acb_poly.h" - -void -acb_dft_convol_mullow(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec) -{ - /* TODO: should probably use (acb_struct *) arrays */ - acb_ptr gg, ww; - if (len == 0) - return; - gg = _acb_vec_init(2 * len - 1); - ww = _acb_vec_init(2 * len - 1); - _acb_vec_set(gg, g, len); - _acb_vec_set(gg + len, g, len - 1); - _acb_poly_mullow(ww, f, len, gg, 2 * len - 1, 2 * len - 1, prec); - _acb_vec_set(w, ww + len, len - 1); - acb_set(w + len - 1, ww + len - 1); - _acb_vec_clear(gg, 2 * len - 1); - _acb_vec_clear(ww, 2 * len - 1); -} diff --git a/src/acb_dft/convol_naive.c b/src/acb_dft/convol_naive.c deleted file mode 100644 index e1061b1976..0000000000 --- a/src/acb_dft/convol_naive.c +++ /dev/null @@ -1,30 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "acb_dft.h" - -void -acb_dft_convol_naive(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec) -{ - slong x, y; - for (x = 0; x < len; x ++) - { - acb_ptr wx; - acb_srcptr fx; - wx = w + x; - fx = f + x; - acb_zero(wx); - for (y = 0; y <= x; y++) - acb_addmul(wx, fx - y, g + y, prec); - for (; y < len; y++) - acb_addmul(wx, fx + (len - y), g + y, prec); - } -} diff --git a/src/acb_dft/convol_rad2.c b/src/acb_dft/convol_rad2.c deleted file mode 100644 index bc31f9ecc1..0000000000 --- a/src/acb_dft/convol_rad2.c +++ /dev/null @@ -1,87 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "acb_dft.h" - -/* assume np >= 2 * n - 1 */ -static void -acb_dft_convol_pad(acb_ptr fp, acb_ptr gp, acb_srcptr f, acb_srcptr g, slong n, slong np) -{ - slong k; - - if (np < 2 * n - 1) - { - flint_throw(FLINT_ERROR, "dft_convol_pad: overlapping padding %wd < 2*%wd-1\n", np, n); - } - - for (k = 0; k < n; k++) - acb_set(gp + k, g + k); - for (; k < np; k++) - acb_zero(gp + k); - - for (k = 0; k < n; k++) - acb_set(fp + k, f + k); - for (k = 1; k < n; k++) - acb_set(fp + np - k, f + n - k); - for (k = n; k <= np - n; k++) - acb_zero(fp + k); - -} - -void -acb_dft_convol_rad2_precomp(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, const acb_dft_rad2_t rad2, slong prec) -{ - slong np; - acb_ptr fp, gp; - np = rad2->n; - - if (len <= 0) - return; - - fp = _acb_vec_init(np); - gp = _acb_vec_init(np); - - if (len == np) - { - _acb_vec_set(fp, f, len); - _acb_vec_set(gp, g, len); - } - else - acb_dft_convol_pad(fp, gp, f, g, len, np); - - acb_dft_rad2_precomp_inplace(fp, rad2, prec); - acb_dft_rad2_precomp_inplace(gp, rad2, prec); - - _acb_vec_kronecker_mul(gp, gp, fp, np, prec); - - acb_dft_inverse_rad2_precomp_inplace(gp, rad2, prec); - - _acb_vec_set(w, gp, len); - _acb_vec_clear(fp, np); - _acb_vec_clear(gp, np); -} - -void -acb_dft_convol_rad2(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec) -{ - int e; - acb_dft_rad2_t dft; - /* catch power of 2 */ - if (len <= 0) - return; - else if ((len & (len - 1)) == 0) - e = FLINT_CLOG2(len); - else - e = FLINT_CLOG2(2 * len - 1); - acb_dft_rad2_init(dft, e, prec); - acb_dft_convol_rad2_precomp(w, f, g, len, dft, prec); - acb_dft_rad2_clear(dft); -} diff --git a/src/acb_dft/crt.c b/src/acb_dft/crt.c deleted file mode 100644 index 9d231a5cdd..0000000000 --- a/src/acb_dft/crt.c +++ /dev/null @@ -1,210 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "nmod.h" -#include "acb_dft.h" - -void -crt_init(crt_t c, ulong n) -{ - int k; - n_factor_t fac; - - n_factor_init(&fac); - if (n) - n_factor(&fac, n, 1); - else - fac.num = 0; - - nmod_init(&c->n, FLINT_MAX(n, 1)); - - c->num = fac.num; - for (k = 0; k < fac.num; k++) - { - c->m[k] = n_pow(fac.p[k], fac.exp[k]); - c->M[k] = n / c->m[k]; - c->vM[k] = c->M[k] * n_invmod(c->M[k] % c->m[k], c->m[k]); - /* - flint_printf("m[%ld]=%ld, M[%ld]=%wu, vM[%ld]=%wu\n", - k, c->m[k], k, c->M[k], k, c->vM[k]); - */ - } -} - -#if 0 -/* lexicographic index of crt elt j */ -static ulong -index_crt(crt_t c, ulong j) -{ - int k; - ulong res = 0; - for (k = 0; k < c->num; k ++) - res = res * c->m[k] + (j % c->m[k]); - return res; -} - -/* crt elt of lexicographic index i */ -static ulong -crt_index(crt_t c, ulong i) -{ - int k; - ulong j, res = 0; - for (k = 0; k < c->num; k ++) - { - j = i % c->m[k]; - i = i / c->m[k]; - res = nmod_add(res, j * c->vM[k], c->n); - } - return res; -} -/* for all elements can do fast conrey-like loop just adding vM[k] */ -static acb_ptr -reorder_to_crt(acb_srcptr v, crt_t c, ulong len) -{ - ulong k; - acb_ptr res; - res = flint_malloc(len * sizeof(acb_struct)); - for (k = 0; k < len; k++) - res[index_crt(c, k)] = v[k]; - return res; -} -static acb_ptr -reorder_from_crt(acb_srcptr v, crt_t c, ulong len) -{ - ulong k; - acb_ptr res; - res = flint_malloc(len * sizeof(acb_struct)); - for (k = 0; k < len; k++) - res[k] = v[index_crt(c, k)]; - return res; -} -#endif - -void -crt_decomp(acb_ptr y, acb_srcptr x, slong dv, const crt_t c, ulong len) -{ - int j, e[CRT_MAX]; - ulong k, l; - - for (j = 0; j < c->num; j++) - e[j] = 0; - - l = 0; - for(k = 0; k < len; k++) - { - acb_set(y + k, x + l * dv); - for (j = c->num - 1; j >= 0; e[j] = 0, j--) - { - e[j]++; l = nmod_add(l, c->vM[j], c->n); - if (e[j] < c->m[j]) - break; - } - } -} - -void -crt_recomp(acb_ptr y, acb_srcptr x, const crt_t c, ulong len) -{ - int j, e[CRT_MAX]; - ulong k, l; - - for (j = 0; j < c->num; j++) - e[j] = 0; - - l = 0; - for(k = 0; k < len; k++) - { - acb_set(y + l, x + k); - for (j = c->num - 1; j >= 0; e[j] = 0, j--) - { - e[j]++; l = nmod_add(l, c->M[j], c->n); - if (e[j] < c->m[j]) - break; - } - } -} - -void -_acb_dft_crt_init(acb_dft_crt_t crt, slong dv, slong len, slong prec) -{ - crt->n = len; - crt_init(crt->c, len); - crt->dv = dv; - crt->cyc = _acb_dft_steps_prod(crt->c->m, crt->c->num, prec); -} - -void -acb_dft_crt_init(acb_dft_crt_t crt, slong len, slong prec) -{ - crt->n = len; - crt_init(crt->c, len); - crt->dv = 1; - crt->cyc = _acb_dft_steps_prod(crt->c->m, crt->c->num, prec); -} - -void -acb_dft_crt_clear(acb_dft_crt_t crt) -{ - slong i; - for (i = 0; i < crt->c->num; i++) - acb_dft_precomp_clear(crt->cyc[i].pre); - flint_free(crt->cyc); -} - -void -acb_dft_crt_precomp(acb_ptr w, acb_srcptr v, const acb_dft_crt_t crt, slong prec) -{ - if (crt->n <= 1) - { - if (crt->n == 1) - acb_set(w, v); - } - else - { - acb_ptr t; - t = _acb_vec_init(crt->n); - if (w == v) - { - _acb_vec_set(t, v, crt->n); - v = t; - } - crt_decomp(w, v, crt->dv, crt->c, crt->n); - acb_dft_step(t, w, crt->cyc, crt->c->num, prec); - crt_recomp(w, t, crt->c, crt->n); - _acb_vec_clear(t, crt->n); - } -} - -void -acb_dft_crt(acb_ptr w, acb_srcptr v, slong len, slong prec) -{ - if (len <= 1) - { - if (len == 1) - acb_set(w, v); - } - else - { - crt_t c; - acb_ptr t; - t = _acb_vec_init(len); - if (w == v) - { - _acb_vec_set(t, v, len); - v = t; - } - crt_init(c, len); - crt_decomp(w, v, 1, c, len); - acb_dft_prod(t, w, c->m, c->num, prec); - crt_recomp(w, t, c, len); - _acb_vec_clear(t, len); - } -} diff --git a/src/acb_dft/cyc.c b/src/acb_dft/cyc.c deleted file mode 100644 index 55a321fa4a..0000000000 --- a/src/acb_dft/cyc.c +++ /dev/null @@ -1,121 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "ulong_extras.h" -#include "acb_dft.h" - -void -_acb_dft_cyc_init_z_fac(acb_dft_cyc_t t, n_factor_t fac, slong dv, acb_ptr z, slong dz, slong len, slong prec) -{ - slong i, j, num; - t->n = len; - num = 0; - for (i = 0; i < fac.num; i++) - num += fac.exp[i]; - t->num = num; - t->cyc = flint_malloc(num * sizeof(acb_dft_step_struct)); - - if (z == NULL) - { - if (DFT_VERB) - flint_printf("dft_cyc: init z[%wu]\n", t->n); - z = _acb_vec_init(t->n); - _acb_vec_unit_roots(z, -t->n, t->n, prec); - dz = 1; - t->zclear = 1; - } - else - { - if (DFT_VERB) - flint_printf("dft_cyc: roots of order %wu already computed\n", t->n); - t->zclear = 0; - } - - t->z = z; - - num = 0; - for (i = 0; i < fac.num; i++) - { - for (j = 0; j < fac.exp[i]; j++) - { - slong m, M; - m = fac.p[i]; - M = (len /= m); - t->cyc[num].m = m; - t->cyc[num].M = M; - t->cyc[num].dv = dv; - t->cyc[num].z = z; - t->cyc[num].dz = dz; - /* TODO: ugly, reorder should solve this */ - if (num == t->num - 1) - _acb_dft_precomp_init(t->cyc[num].pre, dv, z, dz, m, prec); - else - _acb_dft_precomp_init(t->cyc[num].pre, M, z, dz * M, m, prec); - dv *= m; - dz *= m; - num++; - } - } -} - -void -_acb_dft_cyc_init(acb_dft_cyc_t t, slong dv, slong len, slong prec) -{ - n_factor_t fac; - n_factor_init(&fac); - if (len) - n_factor(&fac, len, 1); - else - fac.num = 0; - _acb_dft_cyc_init_z_fac(t, fac, dv, NULL, 0, len, prec); -} - -void -acb_dft_cyc_clear(acb_dft_cyc_t t) -{ - slong i; - for (i = 0; i < t->num; i++) - acb_dft_precomp_clear(t->cyc[i].pre); - if (t->zclear) - _acb_vec_clear(t->z, t->n); - flint_free(t->cyc); -} - -void -acb_dft_cyc_precomp(acb_ptr w, acb_srcptr v, const acb_dft_cyc_t cyc, slong prec) -{ - if (cyc->num == 0) - { - if (cyc->n == 1) - acb_set(w, v); - } - else - acb_dft_step(w, v, cyc->cyc, cyc->num, prec); -} - -void -acb_dft_cyc(acb_ptr w, acb_srcptr v, slong len, slong prec) -{ - acb_dft_cyc_t cyc; - acb_dft_cyc_init(cyc, len, prec); - acb_dft_cyc_precomp(w, v, cyc, prec); - acb_dft_cyc_clear(cyc); -} - -void -acb_dft_inverse_cyc(acb_ptr w, acb_srcptr v, slong len, slong prec) -{ - /* divide before to keep v const */ - _acb_vec_scalar_div_ui(w, v, len, len, prec); - acb_vec_swap_ri(w, len); - acb_dft_cyc(w, w, len, prec); - acb_vec_swap_ri(w, len); -} diff --git a/src/acb_dft/impl.h b/src/acb_dft/impl.h deleted file mode 100644 index 85a4f44fc3..0000000000 --- a/src/acb_dft/impl.h +++ /dev/null @@ -1,19 +0,0 @@ -/* - Copyright (C) 2025 Albin Ahlbäck - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#ifndef ACB_DFT_IMPL_H -#define ACB_DFT_IMPL_H - -#include "acb_types.h" - -void acb_dft_rad2_reorder(acb_ptr v, slong n); - -#endif diff --git a/src/acb_dft/naive.c b/src/acb_dft/naive.c deleted file mode 100644 index c21bffff1a..0000000000 --- a/src/acb_dft/naive.c +++ /dev/null @@ -1,95 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "acb_dft.h" - -/* all roots are already computed, w != v */ -void -_acb_dft_naive(acb_ptr w, acb_srcptr v, slong dv, acb_srcptr z, slong dz, slong len, slong prec) -{ - slong i, j; - acb_ptr wi, v1 = NULL; - acb_srcptr vj; - - if (w == v) - { - flint_throw(FLINT_ERROR, "\n_acb_dft_naive: does not accept aliasing\n"); - } - - for (i = 0, wi = w; i < len; i++, wi++) - { - acb_zero(wi); - for (j = 0, vj = v; j < len; j++, vj += dv) - acb_addmul(wi, vj, z + dz * (i * j % len), prec); - } - - if (v1) - _acb_vec_clear(v1, len); -} - -void -acb_dft_naive_precomp(acb_ptr w, acb_srcptr v, const acb_dft_naive_t pol, slong prec) -{ - if (v == w) - { - acb_ptr v1 = _acb_vec_init(pol->n); - _acb_vec_set(v1, v, pol->n); - _acb_dft_naive(w, v1, pol->dv, pol->z, pol->dz, pol->n, prec); - _acb_vec_clear(v1, pol->n); - } - else - _acb_dft_naive(w, v, pol->dv, pol->z, pol->dz, pol->n, prec); -} - -void -acb_dft_naive(acb_ptr w, acb_srcptr v, slong len, slong prec) -{ - acb_ptr z, v1 = NULL; - - z = _acb_vec_init(len); - _acb_vec_unit_roots(z, -len, len, prec); - if (w == v) - { - v1 = _acb_vec_init(len); - _acb_vec_set(v1, v, len); - v = v1; - } - - _acb_dft_naive(w, v, 1, z, 1, len, prec); - - if (v1) - _acb_vec_clear(v1, len); - - _acb_vec_clear(z, len); -} - -void -_acb_dft_naive_init(acb_dft_naive_t pol, slong dv, acb_ptr z, slong dz, slong len, slong prec) -{ - pol->n = len; - pol->dv = dv; - - if (z == NULL) - { - if (DFT_VERB) - flint_printf("dft_naive: init z[%wd]\n",len); - pol->z = _acb_vec_init(len); - _acb_vec_unit_roots(pol->z, -len, len, prec); - pol->dz = 1; - pol->zclear = 1; - } - else - { - pol->z = z; - pol->dz = dz; - pol->zclear = 0; - } -} diff --git a/src/acb_dft/precomp.c b/src/acb_dft/precomp.c index 2b0575ad9d..3e94687276 100644 --- a/src/acb_dft/precomp.c +++ b/src/acb_dft/precomp.c @@ -1,5 +1,6 @@ /* Copyright (C) 2016 Pascal Molin + Copyright (C) 2026 Fredrik Johansson This file is part of FLINT. @@ -9,148 +10,70 @@ (at your option) any later version. See . */ -#include "ulong_extras.h" +#include "acb.h" #include "acb_dft.h" +/* Thin wrappers around gr_dft (see gr_dft_acb in gr_dft/acb.c), which + use fixed-point arithmetic with rigorous error bounds by default + and fall back to ball arithmetic when fixed point does not + apply. */ + void -_acb_dft_precomp_init(acb_dft_pre_t pre, slong dv, acb_ptr z, slong dz, slong len, slong prec) +acb_dft(acb_ptr w, acb_srcptr v, slong len, slong prec) { - pre->n = len; - if (len <= 1) - { - pre->type = DFT_NAIVE; - _acb_dft_naive_init(pre->t.naive, dv, z, dz, len, prec); - } - else if (n_is_prime(len)) - { - if (len < 100) - { - pre->type = DFT_NAIVE; - _acb_dft_naive_init(pre->t.naive, dv, z, dz, len, prec); - } - else - { - pre->type = DFT_CONV; - /* FIXME: do not recompute the same bluestein - * scheme if needed several times */ - _acb_dft_bluestein_init(pre->t.bluestein, dv, len, prec); - } - } - else - { - n_factor_t fac; + if (len <= 0) + return; - n_factor_init(&fac); - n_factor(&fac, len, 1); - - if (fac.num == 1) - { - /* TODO: could be p^e, or 2^e, but with dv shift */ - if (fac.p[0] == 2) - { - pre->type = DFT_RAD2; - _acb_dft_rad2_init(pre->t.rad2, dv, fac.exp[0], prec); - } - else - { - pre->type = DFT_CYC; - _acb_dft_cyc_init_z_fac(pre->t.cyc, fac, dv, z, dz, len, prec); - } - } - else - { - pre->type = DFT_CRT; - _acb_dft_crt_init(pre->t.crt, dv, len, prec); - } - } + gr_dft_acb(w, v, len, prec); } void -acb_dft_precomp_init(acb_dft_pre_t pre, slong len, slong prec) +acb_dft_inverse(acb_ptr w, acb_srcptr v, slong len, slong prec) { - _acb_dft_precomp_init(pre, 1, NULL, 0, len, prec); -} + if (len <= 0) + return; -void -acb_dft_precomp_clear(acb_dft_pre_t pre) -{ - switch (pre->type) - { - case DFT_NAIVE: - acb_dft_naive_clear(pre->t.naive); - break; - case DFT_CYC: - acb_dft_cyc_clear(pre->t.cyc); - break; - case DFT_PROD: - acb_dft_prod_clear(pre->t.prod); - break; - case DFT_CRT: - acb_dft_crt_clear(pre->t.crt); - break; - case DFT_RAD2: - acb_dft_rad2_clear(pre->t.rad2); - break; - case DFT_CONV: - acb_dft_bluestein_clear(pre->t.bluestein); - break; - default: - flint_throw(FLINT_ERROR, "acb_dft_clear: unknown strategy code %i\n", pre->type); - } + gr_dft_acb_inverse(w, v, len, prec); } void -acb_dft_precomp(acb_ptr w, acb_srcptr v, const acb_dft_pre_t pre, slong prec) +acb_dft_precomp_init(acb_dft_pre_t pre, slong len, slong prec) { - switch (pre->type) + if (len <= 0) { - case DFT_NAIVE: - acb_dft_naive_precomp(w, v, pre->t.naive, prec); - break; - case DFT_CYC: - acb_dft_cyc_precomp(w, v, pre->t.cyc, prec); - break; - case DFT_PROD: - acb_dft_prod_precomp(w, v, pre->t.prod, prec); - break; - case DFT_CRT: - acb_dft_crt_precomp(w, v, pre->t.crt, prec); - break; - case DFT_RAD2: - acb_dft_rad2_precomp(w, v, pre->t.rad2, prec); - break; - case DFT_CONV: - acb_dft_bluestein_precomp(w, v, pre->t.bluestein, prec); - break; - default: - flint_throw(FLINT_ERROR, "acb_dft_precomp: unknown strategy code %i\n", pre->type); + /* empty plan: transforms of length 0 are no-ops, and clearing + is safe */ + pre->n = 0; + pre->which = 0; + return; } + + if (gr_dft_acb_precomp_init(pre, len, prec) != GR_SUCCESS) + flint_throw(FLINT_ERROR, + "acb_dft_precomp_init: len = %wd, prec = %wd\n", len, prec); } void -acb_dft(acb_ptr w, acb_srcptr v, slong len, slong prec) +acb_dft_precomp_clear(acb_dft_pre_t pre) { - acb_dft_pre_t t; - acb_dft_precomp_init(t, len, prec); - acb_dft_precomp(w, v, t, prec); - acb_dft_precomp_clear(t); + gr_dft_acb_precomp_clear(pre); } void -acb_dft_inverse_precomp(acb_ptr w, acb_srcptr v, const acb_dft_pre_t pre, slong prec) +acb_dft_precomp(acb_ptr w, acb_srcptr v, const acb_dft_pre_t pre, slong prec) { - slong k; - acb_dft_precomp(w, v, pre, prec); - for (k = 1; 2 * k < pre->n; k++) - acb_swap(w + k, w + pre->n - k); - _acb_vec_scalar_div_ui(w, w, pre->n, pre->n, prec); + if (pre->n <= 0) + return; + + gr_dft_acb_precomp(w, v, pre, prec); } + void -acb_dft_inverse(acb_ptr w, acb_srcptr v, slong len, slong prec) +acb_dft_inverse_precomp(acb_ptr w, acb_srcptr v, const acb_dft_pre_t pre, + slong prec) { - slong k; - acb_dft(w, v, len, prec); - for (k = 1; 2 * k < len; k++) - acb_swap(w + k, w + len - k); - _acb_vec_scalar_div_ui(w, w, len, len, prec); + if (pre->n <= 0) + return; + + gr_dft_acb_inverse_precomp(w, v, pre, prec); } diff --git a/src/acb_dft/prod.c b/src/acb_dft/prod.c deleted file mode 100644 index 8acbc1b37a..0000000000 --- a/src/acb_dft/prod.c +++ /dev/null @@ -1,65 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "acb_dft.h" - -acb_dft_step_ptr -_acb_dft_steps_prod(slong * cyc, slong num, slong prec) -{ - slong i, len; - acb_dft_step_ptr s; - s = flint_malloc(num * sizeof(acb_dft_step_struct)); - - len = 1; - for (i = 0; i < num; i++) - len *= cyc[i]; - - for (i = 0; i < num; i++) - { - slong m, M; - m = cyc[i]; - M = (len /= m); - s[i].m = m; - s[i].M = M; - s[i].dv = M; - s[i].dz = 0; - s[i].z = NULL; - _acb_dft_precomp_init(s[i].pre, M, NULL, 0, m, prec); - } - - return s; -} - -void -acb_dft_prod_clear(acb_dft_prod_t t) -{ - slong i; - for (i = 0; i < t->num; i++) - acb_dft_precomp_clear(t->cyc[i].pre); - flint_free(t->cyc); -} - - -void -acb_dft_prod_precomp(acb_ptr w, acb_srcptr v, const acb_dft_prod_t prod, slong prec) -{ - if (prod->num >= 1) - acb_dft_step(w, v, prod->cyc, prod->num, prec); -} - -void -acb_dft_prod(acb_ptr w, acb_srcptr v, slong * cyc, slong num, slong prec) -{ - acb_dft_prod_t t; - acb_dft_prod_init(t, cyc, num, prec); - acb_dft_prod_precomp(w, v, t, prec); - acb_dft_prod_clear(t); -} diff --git a/src/acb_dft/profile/p-convol.c b/src/acb_dft/profile/p-convol.c deleted file mode 100644 index 2000a28dd0..0000000000 --- a/src/acb_dft/profile/p-convol.c +++ /dev/null @@ -1,136 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include -#include "acb_dft.h" -#include "profiler.h" - -#define LOG 0 -#define CSV 1 -#define JSON 2 - -typedef void (*do_f) (acb_ptr z, acb_srcptr x, acb_srcptr y, slong len, slong prec); - -int main(int argc, char *argv[]) -{ - int out; - flint_rand_t state; - - int l, nf = 5; - do_f func[5] = { - acb_dft_convol_naive, - acb_dft_convol_dft, - acb_dft_convol_rad2, - acb_dft_convol_mullow, - acb_dft_convol, - }; - char * name[5] = { - "naive", - "dft", - "rad2", - "mullow", - "default" - }; - - int i, ni = 14; - ulong len[14] = { 60, 5*5*5, 131, 2*3*5*7, 256, - 3*3*3*3*3, 2*2*2*43, 503, 2*2*2*101, 3001, - 8192, 13*11*61, 10000, 101*103 }; - char * factor[14] = { "2^2*3*5", "5^3", "131", "2*3*5*7", "2^8", - "3^5", "2^3*43", "503", "2^3*101", "3001", - "2^13", "11*13*61", "2^4*5^4", "101*103" }; - slong ml[14] = { 20, 10, 10, 5, 5, - 5, 3, 3, 2, 2, - 1, 1, 1, 1}; - - int j, nj = 4; - slong prec[4] = { 64, 128, 512, 1024 }; - slong mp[4] = { 40, 20, 5, 2 }; - - if (argc < 2) - out = LOG; - else if (!strcmp(argv[1], "json")) - out = JSON; - else if (!strcmp(argv[1], "csv")) - out = CSV; - else if (!strcmp(argv[1], "log")) - out = LOG; - else - { - printf("usage: %s [log|csv|json]\n", argv[0]); - flint_abort(); - } - - if (out == CSV) - flint_printf("# %-12s, %7s, %7s, %7s\n","name", "prec", "len", "factor", "time"); - - flint_rand_init(state); - - for (j = 0; j < nj; j++) - { - - for (i = 0; i < ni; i++) - { - slong k, nm; - acb_ptr u, v, w; - - nm = mp[j] * ml[i]; - - u = _acb_vec_init(len[i]); - v = _acb_vec_init(len[i]); - w = _acb_vec_init(len[i]); - - for (k = 0; k < len[i]; k++) - { - acb_randtest_precise(u + k, state, prec[j], 4); - acb_randtest_precise(v + k, state, prec[j], 4); - } - - if (out == LOG) - flint_printf("%wu * CONV(%wu = %s), prec %wu....\n", nm, len[i], factor[i], prec[j]); - - for (l = 0; l < nf; l++) - { - slong m; - - if (len[i] > 1000 && l == 0) - continue; - - if (out == LOG) - flint_printf("%-14s ... ", name[l]); - else if (out == CSV) - flint_printf("%-12s, %7d, %8s, %7d, ", name[l], prec[j], len[i], factor[i]); - else if (out == JSON) - flint_printf("{ \"name\": \"%s\", \"prec\": %d, \"len\": %d, \"factor\": %s, \"time\": ", - name[l], prec[j], len[i], factor[i]); - - TIMEIT_ONCE_START; - - for (m = 0; m < nm; m++) - func[l](w, u, v, len[i], prec[j]); - - TIMEIT_ONCE_STOP; - - if (out == JSON) - flint_printf("}\n"); - else - flint_printf("\n"); - } - - _acb_vec_clear(u, len[i]); - _acb_vec_clear(v, len[i]); - _acb_vec_clear(w, len[i]); - } - - } - flint_cleanup(); - return 0; -} diff --git a/src/acb_dft/profile/p-dft.c b/src/acb_dft/profile/p-dft.c deleted file mode 100644 index 81852c96cb..0000000000 --- a/src/acb_dft/profile/p-dft.c +++ /dev/null @@ -1,145 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include -#include "acb_dft.h" -#include "profiler.h" - -#define LOG 0 -#define CSV 1 -#define JSON 2 - -typedef void (*do_f) (acb_ptr y, acb_srcptr x, slong len, slong prec); - -int main(int argc, char *argv[]) -{ - int out; - flint_rand_t state; - - int l, nf = 5; - do_f func[5] = { - acb_dft_naive, - acb_dft_crt, - acb_dft_cyc, - acb_dft_bluestein, - acb_dft - }; - char * name[6] = { - "naive", - "crt", - "cyc", - "bluestein", - "default", - "precomp" - }; - - int i, ni = 14; - ulong len[14] = { 60, 5*5*5, 131, 2*3*5*7, 256, - 3*3*3*3*3, 2*2*2*43, 503, 2*2*2*101, 3001, - 8192, 13*11*61, 10000, 101*103 }; - char * factor[14] = { "2^2*3*5", "5^3", "131", "2*3*5*7", "2^8", - "3^5", "2^3*43", "503", "2^3*101", "3001", - "2^13", "11*13*61", "2^4*5^4", "101*103" }; - slong ml[14] = { 20, 10, 20, 20, 10, - 10, 10, 5, 5, 5, - 3, 3, 3, 3}; - - int j, nj = 4; - slong prec[4] = { 64, 128, 512, 1024 }; - slong mp[4] = { 40, 20, 5, 2 }; - - if (argc < 2) - out = LOG; - else if (!strcmp(argv[1], "json")) - out = JSON; - else if (!strcmp(argv[1], "csv")) - out = CSV; - else if (!strcmp(argv[1], "log")) - out = LOG; - else - { - printf("usage: %s [log|csv|json]\n", argv[0]); - flint_abort(); - } - - if (out == CSV) - flint_printf("# %-12s, %7s, %7s, %7s\n","name", "prec", "len", "factor", "time"); - - flint_rand_init(state); - - for (j = 0; j < nj; j++) - { - - for (i = 0; i < ni; i++) - { - slong k, nm; - acb_ptr u, v; - - nm = mp[j] * ml[i]; - - u = _acb_vec_init(len[i]); - v = _acb_vec_init(len[i]); - - for (k = 0; k < len[i]; k++) - acb_randtest_precise(u + k, state, prec[j], 4); - - if (out == LOG) - flint_printf("%wu * DFT(%wu = %s), prec %wu....\n", nm, len[i], factor[i], prec[j]); - - for (l = 0; l <= nf; l++) - { - slong m; - - if (len[i] > 300 && l == 0) - continue; - - if (out == LOG) - flint_printf("%-14s ... ", name[l]); - else if (out == CSV) - flint_printf("%-12s, %7d, %8s, %7d, ", name[l], prec[j], len[i], factor[i]); - else if (out == JSON) - flint_printf("{ \"name\": \"%s\", \"prec\": %d, \"len\": %d, \"factor\": %s, \"time\": ", - name[l], prec[j], len[i], factor[i]); - - TIMEIT_ONCE_START; - - if (l < nf) - { - for (m = 0; m < nm; m++) - func[l](v, u, len[i], prec[j]); - } - else - { - acb_dft_pre_t t; - acb_dft_precomp_init(t, len[i], prec[j]); - - for (m = 0; m < nm; m++) - acb_dft_precomp(v, u, t, prec[j]); - - acb_dft_precomp_clear(t); - } - - TIMEIT_ONCE_STOP; - - if (out == JSON) - flint_printf("}\n"); - else - flint_printf("\n"); - } - - _acb_vec_clear(u, len[i]); - _acb_vec_clear(v, len[i]); - } - - } - flint_cleanup(); - return 0; -} diff --git a/src/acb_dft/rad2.c b/src/acb_dft/rad2.c deleted file mode 100644 index e018b69682..0000000000 --- a/src/acb_dft/rad2.c +++ /dev/null @@ -1,121 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "acb_dft.h" -#include "acb_dft/impl.h" - -/* swap each element with one with bit-reversed index */ -void acb_dft_rad2_reorder(acb_ptr v, slong n) -{ - slong i, j, k, l; - - for (i = 0, l = n>>1; i < l; i++) - { - /* j = bit reversal of i */ - for (k = 1, j = 0; k < n; k <<= 1) - { - j <<= 1; - if (i & k) - j |= 1; - } - if (i < j) - acb_swap(v + i, v + j); - else if (i > j) - acb_swap(v + n - 1 - i, v + n - 1 - j); - i++, j |= l; - acb_swap(v + i, v + j); - } - -} - -void -_acb_dft_rad2_init(acb_dft_rad2_t t, slong dv, int e, slong prec) -{ - if (e < 0) - { - flint_throw(FLINT_ERROR, "acb_dft_rad2_init: need e >= 0"); - } - t->e = e; - t->n = WORD(1) << e; - t->dv = dv; - t->nz = t->n >> 1; - t->z = _acb_vec_init(t->nz); - /* set n/2 roots of order n */ - _acb_vec_unit_roots(t->z, -t->n, t->nz, prec); -} - -/* remark: can use same rad2 with smaller power of 2 */ -void -acb_dft_rad2_precomp_inplace(acb_ptr v, const acb_dft_rad2_t rad2, slong prec) -{ - if (flint_get_num_threads() > 1 && rad2-> e > 9) - { - acb_dft_rad2_precomp_inplace_threaded(v, rad2, prec); - } - else - { - slong j, k, l; - slong n = rad2->n, nz = rad2->nz; - acb_ptr p, vend = v + n, w = rad2->z; - acb_t tmp; - acb_init(tmp); - - acb_dft_rad2_reorder(v, n); - - for (k = 1, l = nz; k < n; k <<= 1, l >>= 1) - for (p = v; p < vend; p += k) - for (j = 0; j < nz; j += l, p++) - { - acb_mul(tmp, p + k, w + j, prec); - acb_sub(p + k, p + 0, tmp, prec); - acb_add(p + 0, p + 0, tmp, prec); - } - - acb_clear(tmp); - } -} - -void -acb_dft_inverse_rad2_precomp_inplace(acb_ptr v, const acb_dft_rad2_t rad2, slong prec) -{ - slong k, n = rad2->n; - acb_dft_rad2_precomp_inplace(v, rad2, prec); - _acb_vec_scalar_mul_2exp_si(v, v, n, - rad2->e); - for (k = 1; k < n / 2; k++) - acb_swap(v + k, v + n - k); -} - -void -acb_dft_rad2_inplace(acb_ptr v, int e, slong prec) -{ - acb_dft_rad2_t rad2; - acb_dft_rad2_init(rad2, e, prec); - acb_dft_rad2_precomp_inplace(v, rad2, prec); - acb_dft_rad2_clear(rad2); -} - -void -acb_dft_rad2_precomp(acb_ptr w, acb_srcptr v, const acb_dft_rad2_t rad2, slong prec) -{ - slong k; - for (k = 0; k < rad2->n; k++, v += rad2->dv) - acb_set(w + k, v + 0); - acb_dft_rad2_precomp_inplace(w, rad2, prec); -} - -void -acb_dft_rad2(acb_ptr w, acb_srcptr v, int e, slong prec) -{ - acb_dft_rad2_t rad2; - acb_dft_rad2_init(rad2, e, prec); - acb_dft_rad2_precomp(w, v, rad2, prec); - acb_dft_rad2_clear(rad2); -} diff --git a/src/acb_dft/rad2_threaded.c b/src/acb_dft/rad2_threaded.c deleted file mode 100644 index 5d7544a4db..0000000000 --- a/src/acb_dft/rad2_threaded.c +++ /dev/null @@ -1,129 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - Copyright (C) 2020 D.H.J. Polymath - Copyright (C) 2023 Albin Ahlbäck - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "thread_pool.h" -#include "thread_support.h" -#include "acb_dft.h" -#include "acb_dft/impl.h" - -typedef struct -{ - acb_ptr v; - acb_ptr vend; - slong k; - slong l; - slong jstart; - slong jend; - acb_srcptr w; - slong prec; -} -_worker_arg; - -static void -_acb_dft_rad2_thread(void * arg_ptr) -{ - _worker_arg arg = *((_worker_arg *) arg_ptr); - slong j, rstart, pstep; - acb_ptr p, r; - acb_t tmp; - - acb_init(tmp); - rstart = arg.jstart / arg.l; - pstep = 2 * arg.k; - - for (p = arg.v; p < arg.vend; p += pstep) - { - for (r = p + rstart, j = arg.jstart; j < arg.jend; j += arg.l, r++) - { - acb_mul(tmp, r + arg.k, arg.w + j, arg.prec); - acb_sub(r + arg.k, r + 0, tmp, arg.prec); - acb_add(r + 0, r + 0, tmp, arg.prec); - } - } - - acb_clear(tmp); - flint_cleanup(); - return; -} - -void -acb_dft_rad2_precomp_inplace_threaded(acb_ptr v, const acb_dft_rad2_t rad2, slong prec) -{ - slong num_threads, num_workers; - thread_pool_handle * handles; - _worker_arg * args; - - slong t, logt, logk, logl; - slong n = rad2->n; - slong nz = rad2->nz; /* always n/2 ? */ - slong logn = rad2->e; - - num_workers = flint_request_threads(&handles, nz); - num_threads = num_workers + 1; - - for (logt = 0; WORD(1) << (logt + 1) <= num_threads; logt++); - t = WORD(1) << logt; - - args = FLINT_ARRAY_ALLOC(t, _worker_arg); - - acb_dft_rad2_reorder(v, n); - - for (logk = 0, logl = logn - 1; logk < logn; logk += 1, logl -= 1) - { - slong i, j, p; - slong logpstep = logk + 1 + FLINT_MAX(0, logl - logt); - slong logjstep = logl + FLINT_MIN(logk, logn - 1 - logt); - slong pstep = WORD(1) << logpstep; - slong jstep = WORD(1) << logjstep; - i = 0; - for (p = 0; p < n; p += pstep) - { - for (j = 0; j < nz ; j += jstep) - { - args[i].v = v + p; - args[i].vend = v + p + pstep; - args[i].jstart = j; - args[i].jend = j + jstep; - args[i].k = WORD(1) << logk; - args[i].l = WORD(1) << logl; - args[i].w = rad2->z; - args[i].prec = prec; - - if (i != num_workers) - thread_pool_wake(global_thread_pool, handles[i], 0, _acb_dft_rad2_thread, &args[i]); - else - _acb_dft_rad2_thread(&args[i]); - - i++; - } - } - - if (i != t) - flint_throw(FLINT_ERROR, "unequal i=%wd, t=%wd in %s\n", i, t, __func__); - - for (i = 0; i < num_workers; i++) - thread_pool_wait(global_thread_pool, handles[i]); - } - - flint_give_back_threads(handles, num_workers); - flint_free(args); -} - -void -acb_dft_rad2_inplace_threaded(acb_ptr v, int e, slong prec) -{ - acb_dft_rad2_t rad2; - acb_dft_rad2_init(rad2, e, prec); - acb_dft_rad2_precomp_inplace_threaded(v, rad2, prec); - acb_dft_rad2_clear(rad2); -} diff --git a/src/acb_dft/step.c b/src/acb_dft/step.c deleted file mode 100644 index 1145733446..0000000000 --- a/src/acb_dft/step.c +++ /dev/null @@ -1,84 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "acb_dft.h" - -#define REORDER 0 - -void -acb_dft_step(acb_ptr w, acb_srcptr v, acb_dft_step_ptr cyc, slong num, slong prec) -{ - acb_dft_step_struct c; - if (num == 0) - return; - c = cyc[0]; - if (num == 1) - { - acb_dft_precomp(w, v, c.pre, prec); - /*_acb_dft_base(w, v, c.dv, c.z, c.dz, c.m, prec);*/ - } - else - { - slong i, j; - slong m = c.m, M = c.M, dv = c.dv, dz = c.dz; - acb_srcptr z = c.z; - acb_ptr t; -#if REORDER - acb_ptr w2; -#endif - - t = _acb_vec_init(m * M); - - if (w == v) - { - _acb_vec_set(t, v, m * M); - v = t; - } - - /* m DFT of size M */ - for (i = 0; i < m; i++) - acb_dft_step(w + i * M, v + i * dv, cyc + 1, num - 1, prec); - - /* twiddle if non trivial product */ - if (c.z != NULL) - { - acb_ptr wi; - for (wi = w + M, i = 1; i < m; i++, wi += M) - for (j = 1; j < M; j++) - { - if (DFT_VERB) - flint_printf("z[%wu*%wu]",dz,i*j); - acb_mul(wi + j, wi + j, z + dz * i * j, prec); - } - if (DFT_VERB) - flint_printf("\n"); - } - -#if REORDER - /* reorder w to avoid dv shifts in next DFT */ - w2 = flint_malloc(m * M * sizeof(acb_struct)); - for (j = 0; j < M; j++) - for (i = 0; i < m; i++) - w2[j + M * i] = w[i + m * j]; -#endif - - /* M DFT of size m */ - for (j = 0; j < M; j++) - acb_dft_precomp(t + m * j, w + j, c.pre, prec); - - /* reorder */ - for (i = 0; i < m; i++) - for (j = 0; j < M; j++) - acb_set(w + j + M * i, t + i + m * j); - - _acb_vec_clear(t, m * M); - } -} diff --git a/src/acb_dft/test/main.c b/src/acb_dft/test/main.c index abb80b3894..38905333c2 100644 --- a/src/acb_dft/test/main.c +++ b/src/acb_dft/test/main.c @@ -12,7 +12,6 @@ /* Include functions *********************************************************/ #include "t-accuracy.c" -#include "t-convol.c" #include "t-dft.c" /* Array of test functions ***************************************************/ @@ -20,7 +19,6 @@ test_struct tests[] = { TEST_FUNCTION(acb_dft_accuracy), - TEST_FUNCTION(acb_dft_convol), TEST_FUNCTION(acb_dft) }; diff --git a/src/acb_dft/test/t-accuracy.c b/src/acb_dft/test/t-accuracy.c index 0874a5575d..70247daa5e 100644 --- a/src/acb_dft/test/t-accuracy.c +++ b/src/acb_dft/test/t-accuracy.c @@ -1,5 +1,6 @@ /* Copyright (C) 2026 Edgar Costa + Copyright (C) 2026 Fredrik Johansson This file is part of FLINT. @@ -13,28 +14,20 @@ #include "acb.h" #include "acb_dft.h" -/* Regression guard for the acb_dft plan root/factor tables (issue #2709). - - Every acb_dft plan builds its roots of unity by raising a base root to - powers. When the base root was built at the call precision, the powering - amplified its error roughly in proportion to the root index, so the stored - roots near the top of the table carried about (index) * 2^-prec of error and - the transform inherited it: the output radius grew about linearly in n - (about 2490x looser than achievable at n = 2^16, prec = 128 for rad2, and - more for the bluestein factor table used at prime lengths). - - With exact integer input (zero input radius) the output radius is purely the - transform error, so we isolate the plan table by holding the transform - arithmetic at prec while building the plan once at prec and once at prec+64, - and comparing the output radii. A loose table makes rad(plan=prec) grow with - n relative to rad(plan=prec+64); a tight table keeps the ratio a small - constant. Comparing acb_dft at prec vs prec+64 would instead vary the - arithmetic and is dominated by the working-precision gap, so we always vary - only the plan precision. acb_dft_naive is intentionally not checked: it has - its own O(n) summation widening, independent of the root table. */ - -#define RATIO_MAX 8.0 /* fixed: rad2 ~1.5x, bluestein ~2-4x; broken: >= 24x */ -#define GROWTH_MAX 4.0 /* fixed: ~1x; broken: ~54x over 2^8..2^14 */ +/* Accuracy guard for the acb_dft plans (historically issue #2709, + where loosely built root tables made the output radius grow about + linearly in n). The module is now a wrapper around gr_dft, whose + fixed-point path guarantees an absolute output error below + 2^-(prec+2) relative to the input scale; with exact integer input + (zero input radius) the output radius is purely the transform + error, so we assert directly that the radius stays below + 2^-(prec-MARGIN) times the largest output at every length. This + guards the entire limb-selection and error-bound machinery: any + regression that lets the error grow with n shows up as a shrinking + margin. A plan built at prec+64 must also stay consistent (overlap) + with the plan built at prec. */ + +#define MARGIN 16 static void fill_pattern(acb_ptr v, slong n) @@ -49,30 +42,6 @@ fill_pattern(acb_ptr v, slong n) } } -static double -max_rad(acb_srcptr v, slong n) -{ - arb_t m, r; - double d; - slong i; - - arb_init(m); - arb_init(r); - for (i = 0; i < n; i++) - { - arb_get_rad_arb(r, acb_realref(v + i)); - if (arb_gt(r, m)) - arb_set(m, r); - arb_get_rad_arb(r, acb_imagref(v + i)); - if (arb_gt(r, m)) - arb_set(m, r); - } - d = arf_get_d(arb_midref(m), ARF_RND_NEAR); - arb_clear(m); - arb_clear(r); - return d; -} - static void check_overlap(acb_srcptr a, acb_srcptr b, slong n, const char * what) { @@ -88,166 +57,103 @@ check_overlap(acb_srcptr a, acb_srcptr b, slong n, const char * what) } } +/* max radius <= 2^-(prec-MARGIN) max |midpoint|, in mag arithmetic + (robust at high precision, where the radii underflow doubles) */ static void -check_ratio(double ratio, double bound, const char * what) -{ - if (ratio > bound) - { - flint_printf("FAIL (looseness): %s : ratio %g > %g (see issue #2709)\n", - what, ratio, bound); - flint_abort(); - } -} - -/* rad2 plan built at prec vs prec+64, both transforms run at prec. */ -static double -rad2_plan_ratio(int e, slong prec) +check_radius(acb_srcptr w, slong n, slong prec, const char * what) { - slong n = WORD(1) << e; - acb_ptr in = _acb_vec_init(n); - acb_ptr lo = _acb_vec_init(n); - acb_ptr hi = _acb_vec_init(n); - acb_dft_rad2_t plo, phi; - double rlo, rhi; - - fill_pattern(in, n); - acb_dft_rad2_init(plo, e, prec); - acb_dft_rad2_init(phi, e, prec + 64); - acb_dft_rad2_precomp(lo, in, plo, prec); - acb_dft_rad2_precomp(hi, in, phi, prec); - acb_dft_rad2_clear(plo); - acb_dft_rad2_clear(phi); - - check_overlap(lo, hi, n, "rad2 plan prec vs prec+64"); - rlo = max_rad(lo, n); - rhi = max_rad(hi, n); - - _acb_vec_clear(in, n); - _acb_vec_clear(lo, n); - _acb_vec_clear(hi, n); - return (rhi > 0.0) ? rlo / rhi : 1.0; -} + mag_t rad, mid, t; + slong i; -/* bluestein plan (any length >= 30) built at prec vs prec+64, run at prec. */ -static double -bluestein_plan_ratio(slong n, slong prec) -{ - acb_ptr in = _acb_vec_init(n); - acb_ptr lo = _acb_vec_init(n); - acb_ptr hi = _acb_vec_init(n); - acb_dft_bluestein_t plo, phi; - double rlo, rhi; + mag_init(rad); + mag_init(mid); + mag_init(t); - fill_pattern(in, n); - acb_dft_bluestein_init(plo, n, prec); - acb_dft_bluestein_init(phi, n, prec + 64); - acb_dft_bluestein_precomp(lo, in, plo, prec); - acb_dft_bluestein_precomp(hi, in, phi, prec); - acb_dft_bluestein_clear(plo); - acb_dft_bluestein_clear(phi); + for (i = 0; i < n; i++) + { + mag_max(rad, rad, arb_radref(acb_realref(w + i))); + mag_max(rad, rad, arb_radref(acb_imagref(w + i))); + arf_get_mag(t, arb_midref(acb_realref(w + i))); + mag_max(mid, mid, t); + arf_get_mag(t, arb_midref(acb_imagref(w + i))); + mag_max(mid, mid, t); + } - check_overlap(lo, hi, n, "bluestein plan prec vs prec+64"); - rlo = max_rad(lo, n); - rhi = max_rad(hi, n); + mag_mul_2exp_si(mid, mid, -(prec - MARGIN)); + if (mag_cmp(rad, mid) > 0) + { + flint_printf("FAIL (radius): %s, n = %wd, prec = %wd\n", + what, n, prec); + flint_printf("max rad = "); mag_printd(rad, 5); + flint_printf(", bound = "); mag_printd(mid, 5); + flint_printf("\n"); + flint_abort(); + } - _acb_vec_clear(in, n); - _acb_vec_clear(lo, n); - _acb_vec_clear(hi, n); - return (rhi > 0.0) ? rlo / rhi : 1.0; + mag_clear(rad); + mag_clear(mid); + mag_clear(t); } -/* the public acb_dft dispatcher at prec vs a tight reference plan (same kind, - built at prec+64, run at prec). kind 0 = rad2 (power of two), 1 = bluestein. */ -static double -dispatcher_ratio(slong n, int e, slong prec, int kind) +/* one length: the one-shot and the plans at prec and prec + 64 all + agree, and every output radius meets the precision bound */ +static void +check_length(slong n, slong prec) { acb_ptr in = _acb_vec_init(n); - acb_ptr wd = _acb_vec_init(n); - acb_ptr ref = _acb_vec_init(n); - double rd, rr; + acb_ptr w = _acb_vec_init(n); + acb_ptr wp = _acb_vec_init(n); + acb_dft_pre_t plo, phi; fill_pattern(in, n); - acb_dft(wd, in, n, prec); - if (kind == 0) - { - acb_dft_rad2_t p; - acb_dft_rad2_init(p, e, prec + 64); - acb_dft_rad2_precomp(ref, in, p, prec); - acb_dft_rad2_clear(p); - } - else - { - acb_dft_bluestein_t p; - acb_dft_bluestein_init(p, n, prec + 64); - acb_dft_bluestein_precomp(ref, in, p, prec); - acb_dft_bluestein_clear(p); - } + acb_dft(w, in, n, prec); + check_radius(w, n, prec, "acb_dft"); - check_overlap(wd, ref, n, "acb_dft dispatcher vs tight reference"); - rd = max_rad(wd, n); - rr = max_rad(ref, n); + acb_dft_precomp_init(plo, n, prec); + acb_dft_precomp_init(phi, n, prec + 64); + acb_dft_precomp(wp, in, plo, prec); + check_overlap(w, wp, n, "one-shot vs plan"); + check_radius(wp, n, prec, "plan at prec"); + acb_dft_precomp(wp, in, phi, prec); + check_overlap(w, wp, n, "plan prec vs prec+64"); + acb_dft_precomp_clear(plo); + acb_dft_precomp_clear(phi); _acb_vec_clear(in, n); - _acb_vec_clear(wd, n); - _acb_vec_clear(ref, n); - return (rr > 0.0) ? rd / rr : 1.0; + _acb_vec_clear(w, n); + _acb_vec_clear(wp, n); } TEST_FUNCTION_START(acb_dft_accuracy, state) { slong prec = 128; int es[5] = { 8, 11, 12, 14, 16 }; - int ne = 4; /* default tops out at 2^14 (~0.1s) */ - int ehi, e_disp, i; + int ne = 4; /* default tops out at 2^14 */ + int i; slong nb; - double r_lo, r_hi = 0.0; if (flint_test_multiplier() > 1) - ne = 5; /* exercise the headline n = 2^16 too */ - ehi = es[ne - 1]; + ne = 5; /* exercise n = 2^16 too */ - /* (A) rad2 plan tightness, isolated, plus (B) the ratio must not grow - with n (a ratio of ratios, so per-platform rounding cancels). */ - r_lo = rad2_plan_ratio(es[0], prec); + /* powers of two across the size range: the radius bound holding + uniformly in n is the growth guard */ for (i = 0; i < ne; i++) - { - double ratio = rad2_plan_ratio(es[i], prec); + check_length(WORD(1) << es[i], prec); - check_ratio(ratio, RATIO_MAX, "rad2 plan"); - if (es[i] == ehi) - r_hi = ratio; - } - if (r_hi > GROWTH_MAX * r_lo) - { - flint_printf("FAIL (growth): rad2 ratio(2^%d)=%g vs ratio(2^%d)=%g " - "(factor %g > %g) (see issue #2709)\n", - ehi, r_hi, es[0], r_lo, - r_lo > 0 ? r_hi / r_lo : 0.0, GROWTH_MAX); - flint_abort(); - } - - /* tiny n and small precision exercise the worst guard margin */ + /* tiny n */ for (i = 1; i <= 3; i++) - check_ratio(rad2_plan_ratio(i, prec), RATIO_MAX, "rad2 plan tiny n"); - check_ratio(rad2_plan_ratio(10, 8), RATIO_MAX, "rad2 plan small prec"); + check_length(WORD(1) << i, prec); + + /* small and large working precision */ + check_length(1024, 8); + check_length(1009, 256); - /* (.5) bluestein factor table: a prime >= 100 and a prime in [30, 100) */ + /* prime lengths */ nb = (flint_test_multiplier() > 1) ? 16411 : 1009; - check_ratio(bluestein_plan_ratio(nb, prec), RATIO_MAX, "bluestein large prime"); - check_ratio(bluestein_plan_ratio(97, prec), RATIO_MAX, "bluestein n=97"); - check_ratio(bluestein_plan_ratio(29, prec), RATIO_MAX, "bluestein n=29"); - - /* prec > 128 forces heap-allocated mantissas: this case additionally lets - `make valgrind` catch a regression of the np-length scratch clear in - acb_dft_bluestein_precomp (which only leaks above the inline limb cap). */ - check_ratio(bluestein_plan_ratio(1009, 256), RATIO_MAX, "bluestein prec 256"); - - /* the public acb_dft dispatcher is tight on both routed paths */ - e_disp = (flint_test_multiplier() > 1) ? 16 : 14; - check_ratio(dispatcher_ratio(WORD(1) << e_disp, e_disp, prec, 0), RATIO_MAX, - "acb_dft power of two"); - check_ratio(dispatcher_ratio(nb, 0, prec, 1), RATIO_MAX, "acb_dft prime"); + check_length(nb, prec); + check_length(97, prec); + check_length(29, prec); TEST_FUNCTION_END(state); } diff --git a/src/acb_dft/test/t-convol.c b/src/acb_dft/test/t-convol.c deleted file mode 100644 index c3122d4f48..0000000000 --- a/src/acb_dft/test/t-convol.c +++ /dev/null @@ -1,135 +0,0 @@ -/* - Copyright (C) 2016 Pascal Molin - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "test_helpers.h" -#include "acb_dft.h" - -/* Defined in t-convol.c and t-dft.c. */ -#ifndef acb_vec_printd_index -# define acb_vec_printd_index acb_vec_printd_index -static void acb_vec_printd_index(acb_srcptr vec, slong len, slong digits) -{ - slong i; - for (i = 0; i < len; i++) - { - flint_printf("[%wd] ", i); - acb_printd(vec + i, digits); - flint_printf("\n"); - } -} -#endif - -/* Both of these are defined in t-convol.c and t-dft.c. However, do_f is not the - * same. */ -#define do_f hogusbogusname -typedef void (*do_f) (acb_ptr z, acb_srcptr x, acb_srcptr y, slong len, slong prec); - -#ifndef check_vec_eq_prec -#define check_vec_eq_prec check_vec_eq_prec -void -check_vec_eq_prec(acb_srcptr w1, acb_srcptr w2, slong len, slong prec, slong digits, ulong q, char d[], char f1[], char f2[]) -{ - slong i; - - for (i = 0; i < len; i++) - { - if (!acb_overlaps(w1 + i, w2 + i)) - { - flint_printf("FAIL\n\n"); - flint_printf("q = %wu, size = %wd\n", q, len); - flint_printf("\nDFT %sdiffer from index %wd / %wd \n", d, i, len); - flint_printf("\n%s =\n", f1); - acb_vec_printd_index(w1, len, digits); - flint_printf("\n%s =\n", f2); - acb_vec_printd_index(w2, len, digits); - flint_printf("\n\n"); - flint_abort(); - } - else if (!acb_is_zero(w1+i) && (acb_rel_accuracy_bits(w1 + i) < 30 - || acb_rel_accuracy_bits(w2 + i) < 30)) - { - flint_printf("FAIL\n\n"); - flint_printf("q = %wu\n", q); - flint_printf("\nDFT inaccurate from index %wd / %wd \n", i, len); - flint_printf("\n%s =\n", f1); - acb_printd(w1 + i, digits); - flint_printf("\n%s =\n", f2); - acb_printd(w2 + i, digits); - flint_printf("\nerrors %wd & %wd [prec = %wd]\n", - acb_rel_accuracy_bits(w1 + i), - acb_rel_accuracy_bits(w2 + i), prec); - flint_abort(); - } - } -} -#endif - -TEST_FUNCTION_START(acb_dft_convol, state) -{ - slong k; - slong prec = 100, digits = 30; - slong nq = 17; - ulong q[17] = { 0, 1, 2, 3, 4, 5, 6, 23, 10, 15, 30, 59, 256, 308, 335, 344, 961}; - ulong nr = 5; - - slong f, nf = 4; - do_f func[4] = { acb_dft_convol_naive, acb_dft_convol_rad2, acb_dft_convol_dft, acb_dft_convol_mullow }; - char * name[4] = { "naive", "rad2", "dft", "mullow" }; - - for (k = 0; k < nq + nr; k++) - { - slong i, len; - acb_ptr z1, z2, x, y; - - if (k < nq) - len = q[k]; - else - len = n_randint(state, 2000); - - z1 = _acb_vec_init(len); - z2 = _acb_vec_init(len); - x = _acb_vec_init(len); - y = _acb_vec_init(len); - - for (i = 0; i < len; i++) - { -#if 1 - acb_set_si(x + i, n_randint(state, 4 * len)); - acb_set_si(y + i, n_randint(state, 4 * len)); -#else - acb_set_si_si(x + i, n_randint(state, 4 * len), n_randint(state, 4 * len)); - acb_set_si_si(y + i, n_randint(state, 4 * len), n_randint(state, 4 * len)); -#endif - } - - for (f = 0; f < nf; f++) - { - - acb_ptr z = (f == 0) ? z1 : z2; - - func[f](z, x, y, len, prec); - - if (f == 0) - continue; - - check_vec_eq_prec(z1, z2, len, prec, digits, len, "", name[0], name[f]); - - } - - _acb_vec_clear(x, len); - _acb_vec_clear(y, len); - _acb_vec_clear(z1, len); - _acb_vec_clear(z2, len); - } - - TEST_FUNCTION_END(state); -} -#undef do_f diff --git a/src/acb_dft/test/t-dft.c b/src/acb_dft/test/t-dft.c index 0a5eac2b9e..6689442256 100644 --- a/src/acb_dft/test/t-dft.c +++ b/src/acb_dft/test/t-dft.c @@ -1,5 +1,6 @@ /* Copyright (C) 2016 Pascal Molin + Copyright (C) 2026 Fredrik Johansson This file is part of FLINT. @@ -10,185 +11,108 @@ */ #include "test_helpers.h" +#include "arb.h" +#include "acb.h" #include "acb_dft.h" -/* Defined in t-convol.c and t-dft.c. */ -#ifndef acb_vec_printd_index -# define acb_vec_printd_index acb_vec_printd_index -static void acb_vec_printd_index(acb_srcptr vec, slong len, slong digits) +/* The module is a thin wrapper around gr_dft, whose test suite + stress-tests the algorithms; here we verify that every public + function behaves as documented (definition, inverse, plan reuse, + aliasing, trivial lengths). */ + +/* naive X_k = sum_j v_j e(-jk/n) directly from the definition */ +static void +_dft_naive(acb_ptr w, acb_srcptr v, slong n, slong prec) { - slong i; - for (i = 0; i < len; i++) + slong j, k; + acb_ptr roots; + acb_t t; + + roots = _acb_vec_init(n); + acb_init(t); + for (j = 0; j < n; j++) { - flint_printf("[%wd] ", i); - acb_printd(vec + i, digits); - flint_printf("\n"); + arb_set_si(acb_realref(t), -2 * j); + arb_div_ui(acb_realref(t), acb_realref(t), n, prec + 10); + arb_zero(acb_imagref(t)); + acb_exp_pi_i(roots + j, t, prec + 10); } -} -#endif - -/* Both of these are defined in t-convol.c and t-dft.c. However, do_f is not the - * same. */ -#define do_f rajtantajtan -typedef void (*do_f) (acb_ptr w, acb_srcptr v, slong len, slong prec); - -#ifndef check_vec_eq_prec -#define check_vec_eq_prec check_vec_eq_prec -void -check_vec_eq_prec(acb_srcptr w1, acb_srcptr w2, slong len, slong prec, slong digits, ulong q, char d[], char f1[], char f2[]) -{ - slong i; - for (i = 0; i < len; i++) + for (k = 0; k < n; k++) { - if (!acb_overlaps(w1 + i, w2 + i)) - { - flint_printf("FAIL\n\n"); - flint_printf("q = %wu, size = %wd\n", q, len); - flint_printf("\nDFT %sdiffer from index %wd / %wd \n", d, i, len); - flint_printf("\n%s =\n", f1); - acb_vec_printd_index(w1, len, digits); - flint_printf("\n%s =\n", f2); - acb_vec_printd_index(w2, len, digits); - flint_printf("\n\n"); - flint_abort(); - } - else if (!acb_is_zero(w1+i) && (acb_rel_accuracy_bits(w1 + i) < 30 - || acb_rel_accuracy_bits(w2 + i) < 30)) - { - flint_printf("FAIL\n\n"); - flint_printf("q = %wu\n", q); - flint_printf("\nDFT inaccurate from index %wd / %wd \n", i, len); - flint_printf("\n%s =\n", f1); - acb_printd(w1 + i, digits); - flint_printf("\n%s =\n", f2); - acb_printd(w2 + i, digits); - flint_printf("\nerrors %wd & %wd [prec = %wd]\n", - acb_rel_accuracy_bits(w1 + i), - acb_rel_accuracy_bits(w2 + i), prec); - flint_abort(); - } + acb_zero(w + k); + for (j = 0; j < n; j++) + acb_addmul(w + k, roots + (j * k) % n, v + j, prec); } + + _acb_vec_clear(roots, n); + acb_clear(t); } -#endif TEST_FUNCTION_START(acb_dft, state) { slong k; - slong prec = 100, digits = 30; - slong nq = 19; - ulong q[19] = { 0, 1, 2, 3, 4, 5, 6, 23, 10, 15, 16, 30, 59, 125, 308, 335, 525, 961, 1225}; - slong nr = 5; + slong lens[9] = { 1, 2, 3, 4, 6, 8, 12, 30, 59 }; - slong f, nf = 5; - do_f func[5] = { acb_dft_naive, acb_dft_cyc, acb_dft_crt, acb_dft_bluestein, acb_dft }; - char * name[5] = { "naive", "cyc", "crt", "bluestein", "default" }; - - /* cyclic dft */ - for (k = 0; k < nq + nr; k++) + for (k = 0; k < 9; k++) { - slong i, len, f0; - acb_ptr v, w1, w2, w3; - - if (k < nq) - len = q[k]; - else - len = n_randint(state, 1000); + slong len = lens[k], i; + slong prec = 64 + n_randint(state, 128); + acb_ptr v, w1, w2; + acb_dft_pre_t pre; v = _acb_vec_init(len); w1 = _acb_vec_init(len); w2 = _acb_vec_init(len); - w3 = _acb_vec_init(len); for (i = 0; i < len; i++) - acb_set_si_si(v + i, i, 3 - i); - - /* avoid naive for long transforms */ - f0 = (len > 50); - - for (f = f0; f < nf; f++) - { - - acb_ptr w = (f == f0) ? w1 : w2; - - if (DFT_VERB) - flint_printf("\n%s %wu\n", name[f], len); + acb_randtest_precise(v + i, state, prec, 0); - /* compute DFT */ - func[f](w, v, len, prec); - - if (len < 500) - { - /* check aliasing */ - _acb_vec_set(w3, v, len); - func[f](w3, w3, len, prec); - - check_vec_eq_prec(w1, w3, len, prec, digits, len, "alias ", name[0], name[f]); - } + /* the one-shot transform computes the DFT */ + _dft_naive(w1, v, len, prec); + acb_dft(w2, v, len, prec); + for (i = 0; i < len; i++) + if (!acb_overlaps(w1 + i, w2 + i)) + TEST_FUNCTION_FAIL("dft vs definition, len = %wd i = %wd\n", + len, i); + + /* inverse roundtrip, aliased */ + _acb_vec_set(w2, v, len); + acb_dft(w2, w2, len, prec); + acb_dft_inverse(w2, w2, len, prec); + for (i = 0; i < len; i++) + if (!acb_overlaps(w2 + i, v + i)) + TEST_FUNCTION_FAIL("roundtrip, len = %wd i = %wd\n", len, i); - if (f > f0) - { - /* check non aliased */ - check_vec_eq_prec(w1, w2, len, prec, digits, len, "no alias ", name[0], name[f]); - } - else - { - /* check inverse */ - acb_dft_inverse(w2, w1, len, prec); - check_vec_eq_prec(v, w2, len, prec, digits, len, "inverse ", "original", "inverse"); - } - } + /* the plan interface agrees with the one-shots and is + reusable */ + acb_dft_precomp_init(pre, len, prec); + acb_dft_precomp(w2, v, pre, prec); + for (i = 0; i < len; i++) + if (!acb_overlaps(w1 + i, w2 + i)) + TEST_FUNCTION_FAIL("precomp, len = %wd i = %wd\n", len, i); + acb_dft_inverse_precomp(w2, w2, pre, prec); + for (i = 0; i < len; i++) + if (!acb_overlaps(w2 + i, v + i)) + TEST_FUNCTION_FAIL("precomp roundtrip, len = %wd i = %wd\n", + len, i); + acb_dft_precomp_clear(pre); _acb_vec_clear(v, len); _acb_vec_clear(w1, len); _acb_vec_clear(w2, len); - _acb_vec_clear(w3, len); } - /* radix2 dft */ - for (k = 0; k < 11; k++) + /* len = 0 is a no-op for every entry point */ { - slong n = 1 << k, j; - acb_ptr v, w1, w2; - v = w2 = _acb_vec_init(n); - w1 = _acb_vec_init(n); - - for (j = 0; j < n; j++) - acb_set_si_si(v + j, j, j + 2); - - acb_dft_cyc(w1, v, n, prec); - acb_dft_rad2_inplace(w2, k, prec); - - check_vec_eq_prec(w1, w2, n, prec, digits, n, "rad2 ", "cyc", "rad2"); - - _acb_vec_clear(v, n); - _acb_vec_clear(w1, n); - - } - - /* multi-threaded radix2 dft */ - for (k = 0; k < 11; k++) - { - slong n = 1 << k, j; - acb_ptr v, w1, w2; - v = w2 = _acb_vec_init(n); - w1 = _acb_vec_init(n); - - flint_set_num_threads(k % 5 + 1); - - for (j = 0; j < n; j++) - acb_set_si_si(v + j, j, j + 2); - - acb_dft_cyc(w1, v, n, prec); - acb_dft_rad2_inplace_threaded(w2, k, prec); - - check_vec_eq_prec(w1, w2, n, prec, digits, n, "rad2 ", "cyc", "rad2"); - - _acb_vec_clear(v, n); - _acb_vec_clear(w1, n); - + acb_dft_pre_t pre; + acb_dft(NULL, NULL, 0, 64); + acb_dft_inverse(NULL, NULL, 0, 64); + acb_dft_precomp_init(pre, 0, 64); + acb_dft_precomp(NULL, NULL, pre, 64); + acb_dft_inverse_precomp(NULL, NULL, pre, 64); + acb_dft_precomp_clear(pre); } TEST_FUNCTION_END(state); } -#undef do_f diff --git a/src/acb_dirichlet/dft.c b/src/acb_dirichlet/dft.c index c7eb8542cc..010cc9345d 100644 --- a/src/acb_dirichlet/dft.c +++ b/src/acb_dirichlet/dft.c @@ -1,5 +1,6 @@ /* Copyright (C) 2016 Pascal Molin + Copyright (C) 2026 Fredrik Johansson This file is part of FLINT. @@ -9,61 +10,25 @@ (at your option) any later version. See . */ -#include "acb_dft.h" #include "acb_dirichlet.h" +#include "gr_dft.h" + +/* Thin wrappers around the gr_dft module (see gr_dft/dirichlet.c), + which computes the product DFT over the Conrey components in + fixed-point arithmetic with rigorous error bounds when possible. */ /* dft, lexicographic conrey indexing, array size G->phi_q */ void -acb_dirichlet_dft_index(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec) +acb_dirichlet_dft_index(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, + slong prec) { - if (G->phi_q == 1) - { - acb_set(w, v); - } - else - { - slong k, l, * cyc; - - - cyc = flint_malloc(G->num * sizeof(slong)); - for (k = 0, l = G->num - 1; l >= 0; k++, l--) - cyc[k] = G->P[k].phi.n; - - acb_dft_prod(w, v, cyc, G->num, prec); - flint_free(cyc); - } + gr_dft_acb_dirichlet_index(w, v, G, prec); } /* dft, number indexing, array size G->q */ void -acb_dirichlet_dft(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec) +acb_dirichlet_dft(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, + slong prec) { - ulong i, len; - acb_ptr t1, t2; - dirichlet_char_t x; - - len = G->phi_q; - t1 = flint_malloc(len * sizeof(acb_struct)); - - dirichlet_char_init(x, G); - dirichlet_char_one(x, G); - for (i = 0; i < len; i++) - { - t1[i] = v[x->n]; - dirichlet_char_next(x, G); - }; - - t2 = _acb_vec_init(len); - acb_dirichlet_dft_index(t2, t1, G, prec); - - dirichlet_char_one(x, G); - for (i = 0; i < len; i++) - { - acb_set(w + x->n, t2 + i); - dirichlet_char_next(x, G); - }; - - _acb_vec_clear(t2, len); - dirichlet_char_clear(x); - flint_free(t1); + gr_dft_acb_dirichlet(w, v, G, prec); } diff --git a/src/acb_dirichlet/platt_multieval.c b/src/acb_dirichlet/platt_multieval.c index 1ef69699b5..2fb1d337dc 100644 --- a/src/acb_dirichlet/platt_multieval.c +++ b/src/acb_dirichlet/platt_multieval.c @@ -10,11 +10,20 @@ */ #include "fmpz_vec.h" +#include "acb.h" #include "acb_dirichlet.h" #include "acb_dirichlet/impl.h" #include "arb_hypgeom.h" #include "acb_dft.h" +FLINT_FORCE_INLINE void +_acb_vec_kronecker_mul(acb_ptr z, acb_srcptr x, acb_srcptr y, slong len, slong prec) +{ + slong k; + for (k = 0; k < len; k++) + acb_mul(z + k, x + k, y + k, prec); +} + static void _acb_dot_arb(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, arb_srcptr y, slong ystep, diff --git a/src/acb_dirichlet/test/t-dft.c b/src/acb_dirichlet/test/t-dft.c index 0d9e82c3aa..05cfe86566 100644 --- a/src/acb_dirichlet/test/t-dft.c +++ b/src/acb_dirichlet/test/t-dft.c @@ -1,5 +1,6 @@ /* Copyright (C) 2016 Pascal Molin + Copyright (C) 2026 Fredrik Johansson This file is part of FLINT. @@ -10,109 +11,112 @@ */ #include "test_helpers.h" +#include "arb.h" +#include "acb.h" #include "acb_dirichlet.h" -typedef void (*do_f) (acb_ptr w, acb_srcptr v, slong len, slong prec); - -void -check_vec_eq_prec(acb_srcptr w1, acb_srcptr w2, slong len, slong prec, slong digits, ulong q, char f1[], char f2[]) -{ - slong i; - - for (i = 0; i < len; i++) - { - if (!acb_overlaps(w1 + i, w2 + i)) - { - flint_printf("FAIL\n\n"); - flint_printf("q = %wu, size = %wd\n", q, len); - flint_printf("\nDFT differ from index %wd / %wd \n", i, len); - flint_printf("\n%s =\n", f1); - _acb_vec_printd(w1, len, digits); - flint_printf("\n%s =\n", f2); - _acb_vec_printd(w2, len, digits); - flint_printf("\n\n"); - flint_abort(); - } - else if (acb_rel_accuracy_bits(w1 + i) < 30 - || acb_rel_accuracy_bits(w2 + i) < 30) - { - flint_printf("FAIL\n\n"); - flint_printf("q = %wu\n", q); - flint_printf("\nDFT inaccurate from index %wd / %wd \n", i, len); - flint_printf("\nnaive =\n"); - acb_printd(w1 + i, digits); - flint_printf("\nfast =\n"); - acb_printd(w2 + i, digits); - flint_printf("\nerrors %wd & %wd [prec = %wd]\n", - acb_rel_accuracy_bits(w1 + i), - acb_rel_accuracy_bits(w2 + i), prec); - flint_abort(); - } - } -} +/* acb_dirichlet_dft and acb_dirichlet_dft_index are thin wrappers + around gr_dft, whose test suite validates the transforms against + independent references; here we verify that the wrappers behave as + documented: dft_index computes the pairing sum on a small group, + and dft is its conjugation by the Conrey gather/scatter. */ TEST_FUNCTION_START(acb_dirichlet_dft, state) { slong k; - slong prec = 100, digits = 30; - slong nq = 13; - ulong q[13] = { 2, 3, 4, 5, 6, 23, 10, 15, 30, 59, 308, 335, 961}; + ulong q[6] = { 1, 2, 5, 8, 12, 30 }; - /* Dirichlet group DFT */ - for (k = 0; k < nq; k++) + for (k = 0; k < 6; k++) { - slong i, j, len; + slong prec = 100; dirichlet_group_t G; dirichlet_char_t x, y; - acb_dirichlet_roots_t roots; - acb_t chiy; - acb_ptr v, w1, w2; + acb_ptr v, vq, w1, w2, t1; + acb_t chiy, u; + ulong i, j, len; dirichlet_group_init(G, q[k]); - len = G->phi_q; + v = _acb_vec_init(len); + vq = _acb_vec_init(q[k]); w1 = _acb_vec_init(len); - w2 = _acb_vec_init(len); + w2 = _acb_vec_init(FLINT_MAX(q[k], len)); + t1 = _acb_vec_init(len); acb_init(chiy); - acb_dirichlet_roots_init(roots, G->expo, len, prec); - + acb_init(u); dirichlet_char_init(x, G); dirichlet_char_init(y, G); for (i = 0; i < len; i++) acb_randtest_precise(v + i, state, prec, 0); + for (i = 0; i < q[k]; i++) + acb_randtest_precise(vq + i, state, prec, 0); - /* naive */ + /* dft_index against the pairing sum + w[x] = sum_y conj(chi_x(y)) v[y] */ dirichlet_char_one(x, G); - for (i = 0; i < len; i++) + for (i = 0; i < len; i++) { acb_zero(w1 + i); dirichlet_char_one(y, G); for (j = 0; j < len; j++) { - acb_dirichlet_root(chiy, roots, dirichlet_pairing_char(G, x, y), prec); - acb_conj(chiy, chiy); + arb_set_si(acb_realref(u), + -2 * (slong) dirichlet_pairing_char(G, x, y)); + arb_div_ui(acb_realref(u), acb_realref(u), G->expo, + prec + 10); + arb_zero(acb_imagref(u)); + acb_exp_pi_i(chiy, u, prec + 10); acb_addmul(w1 + i, chiy, v + j, prec); dirichlet_char_next(y, G); } dirichlet_char_next(x, G); } - /* dft */ acb_dirichlet_dft_index(w2, v, G, prec); + for (i = 0; i < len; i++) + if (!acb_overlaps(w1 + i, w2 + i)) + TEST_FUNCTION_FAIL("dft_index, q = %wu i = %wu\n", q[k], i); - check_vec_eq_prec(w1, w2, len, prec, digits, q[k], "naive", "group"); + /* dft is dft_index conjugated by the Conrey gather/scatter */ + if (q[k] > 1) + { + dirichlet_char_one(x, G); + for (i = 0; i < len; i++) + { + acb_set(t1 + i, vq + x->n); + dirichlet_char_next(x, G); + } + acb_dirichlet_dft_index(w1, t1, G, prec); + + acb_dirichlet_dft(w2, vq, G, prec); + dirichlet_char_one(x, G); + for (i = 0; i < len; i++) + { + if (!acb_overlaps(w1 + i, w2 + x->n)) + TEST_FUNCTION_FAIL("dft, q = %wu i = %wu\n", q[k], i); + dirichlet_char_next(x, G); + } + } + else + { + /* q = 1: identity on the single entry */ + acb_dirichlet_dft(w2, vq, G, prec); + if (!acb_equal(w2, vq)) + TEST_FUNCTION_FAIL("dft, q = 1\n"); + } _acb_vec_clear(v, len); + _acb_vec_clear(vq, q[k]); _acb_vec_clear(w1, len); - _acb_vec_clear(w2, len); - acb_dirichlet_roots_clear(roots); - + _acb_vec_clear(w2, FLINT_MAX(q[k], len)); + _acb_vec_clear(t1, len); + acb_clear(chiy); + acb_clear(u); dirichlet_char_clear(x); dirichlet_char_clear(y); dirichlet_group_clear(G); - acb_clear(chiy); } TEST_FUNCTION_END(state); diff --git a/src/acb_theta/ql_jet_fd.c b/src/acb_theta/ql_jet_fd.c index cd8a496116..823e2e99ce 100644 --- a/src/acb_theta/ql_jet_fd.c +++ b/src/acb_theta/ql_jet_fd.c @@ -14,7 +14,7 @@ #include "acb.h" #include "arb_mat.h" #include "acb_mat.h" -#include "acb_dft.h" +#include "gr_dft.h" #include "acb_theta.h" static void @@ -126,7 +126,7 @@ acb_theta_jet_finite_diff(acb_ptr dth, const arf_t eps, const arf_t err, acb_ptr aux; arb_t t, e; slong * tups; - slong * cyc; + ulong * cyc; slong j, i, l; slong k; @@ -134,13 +134,13 @@ acb_theta_jet_finite_diff(acb_ptr dth, const arf_t eps, const arf_t err, arb_init(t); arb_init(e); tups = flint_malloc(g * nbjet * sizeof(slong)); - cyc = flint_malloc(g * sizeof(slong)); + cyc = flint_malloc(g * sizeof(ulong)); for (j = 0; j < g; j++) { cyc[j] = ord + 1; } - acb_dft_prod(aux, val, cyc, g, prec); + gr_dft_acb_prod(aux, val, cyc, g, prec); arb_set_si(t, nbaux); _acb_vec_scalar_div_arb(aux, aux, nbaux, t, prec); diff --git a/src/gr_dft.h b/src/gr_dft.h index 130e04f81a..767d783b6c 100644 --- a/src/gr_dft.h +++ b/src/gr_dft.h @@ -17,6 +17,7 @@ #include "gr.h" #include "arf_types.h" #include "acb_types.h" +#include "dirichlet.h" #ifdef __cplusplus extern "C" { @@ -227,6 +228,79 @@ void gr_dft_acb_precomp(acb_ptr w, acb_srcptr v, const gr_dft_acb_pre_t Q, slong void gr_dft_acb_inverse_precomp(acb_ptr w, acb_srcptr v, const gr_dft_acb_pre_t Q, slong prec); int _gr_dft_acb_precomp(acb_ptr w, acb_srcptr v, int inverse, const gr_dft_acb_pre_t Q, slong prec); + +/* DFT on products of cyclic groups ******************************************/ + +/* Multidimensional DFT over Z/cyc[0] x ... x Z/cyc[num-1], row-major + indexing (the first component varies slowest). Realized as one + cyclic transform along each axis; plans for repeated component + lengths are shared. */ +typedef struct +{ + ulong n; /* product of the component lengths */ + slong num; + ulong * cyc; + slong num_plans; /* distinct component lengths */ + slong * plan_of; /* axis -> plan index */ + gr_dft_pre_struct * plans; + int flags; + slong serial_block; +} +gr_dft_prod_pre_struct; + +typedef gr_dft_prod_pre_struct gr_dft_prod_pre_t[1]; + +WARN_UNUSED_RESULT int gr_dft_prod_precomp_init(gr_dft_prod_pre_t P, const ulong * cyc, slong num, int flags, gr_ctx_t ctx); +WARN_UNUSED_RESULT int gr_dft_prod_precomp_init_root(gr_dft_prod_pre_t P, gr_srcptr w, ulong order, const ulong * cyc, slong num, int flags, gr_ctx_t ctx); +void gr_dft_prod_precomp_clear(gr_dft_prod_pre_t P); +WARN_UNUSED_RESULT int gr_dft_prod_precomp(gr_ptr res, gr_srcptr vec, const gr_dft_prod_pre_t P, gr_ctx_t ctx); +WARN_UNUSED_RESULT int gr_dft_prod_inverse_precomp(gr_ptr res, gr_srcptr vec, const gr_dft_prod_pre_t P, gr_ctx_t ctx); +WARN_UNUSED_RESULT int gr_dft_prod(gr_ptr res, gr_srcptr vec, const ulong * cyc, slong num, gr_ctx_t ctx); +WARN_UNUSED_RESULT int gr_dft_prod_inverse(gr_ptr res, gr_srcptr vec, const ulong * cyc, slong num, gr_ctx_t ctx); +void gr_dft_prod_precomp_set_serial_block(gr_dft_prod_pre_t P, slong serial_block); +void gr_dft_prod_precomp_nfixed_bound(double * peak, double * err_ulps, double in_mag, double in_err, const gr_dft_prod_pre_t P); +WARN_UNUSED_RESULT int _gr_dft_prod_precomp_raw(gr_ptr res, gr_srcptr vec, int inverse, const gr_dft_prod_pre_t P, gr_ctx_t ctx); +int _gr_dft_prod_precomp_init_layout(gr_dft_prod_pre_t P, const ulong * cyc, slong num, int flags, int complex_mode); +int _gr_dft_prod_precomp_realize(gr_dft_prod_pre_t P, gr_ctx_struct * real_ctx, gr_ctx_t ctx); + +/* acb wrapper for the product DFT, using fixed-point arithmetic + internally with rigorous error bounds, falling back to ball + arithmetic when fixed point does not apply */ +typedef struct +{ + slong n; + slong num; + ulong * cyc; + slong prec; + int which; + slong nl; + slong p1e; + double in_mag; + double errulps; + gr_ctx_struct rctx[1]; + gr_ctx_struct cctx[1]; + gr_dft_prod_pre_t P; +} +gr_dft_acb_prod_pre_struct; + +typedef gr_dft_acb_prod_pre_struct gr_dft_acb_prod_pre_t[1]; + +void gr_dft_acb_prod(acb_ptr w, acb_srcptr v, const ulong * cyc, slong num, slong prec); +void gr_dft_acb_prod_inverse(acb_ptr w, acb_srcptr v, const ulong * cyc, slong num, slong prec); +int _gr_dft_acb_prod(acb_ptr w, acb_srcptr v, const ulong * cyc, slong num, int inverse, int which, slong prec); +WARN_UNUSED_RESULT int gr_dft_acb_prod_precomp_init(gr_dft_acb_prod_pre_t Q, const ulong * cyc, slong num, slong prec); +void gr_dft_acb_prod_precomp_clear(gr_dft_acb_prod_pre_t Q); +void gr_dft_acb_prod_precomp(acb_ptr w, acb_srcptr v, const gr_dft_acb_prod_pre_t Q, slong prec); +void gr_dft_acb_prod_inverse_precomp(acb_ptr w, acb_srcptr v, const gr_dft_acb_prod_pre_t Q, slong prec); +int _gr_dft_acb_prod_precomp(acb_ptr w, acb_srcptr v, int inverse, const gr_dft_acb_prod_pre_t Q, slong prec); + +/* DFT on Dirichlet groups (counterparts of acb_dirichlet_dft and + acb_dirichlet_dft_index) */ +WARN_UNUSED_RESULT int gr_dft_dirichlet_index(gr_ptr w, gr_srcptr v, const dirichlet_group_t G, gr_ctx_t ctx); +WARN_UNUSED_RESULT int gr_dft_dirichlet(gr_ptr w, gr_srcptr v, const dirichlet_group_t G, gr_ctx_t ctx); +void gr_dft_acb_dirichlet_index(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec); +void gr_dft_acb_dirichlet(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec); + /* Internal functions *******************************************************/ WARN_UNUSED_RESULT int _gr_dft_mul_root(gr_ptr res, gr_srcptr x, ulong e, int inverse, gr_ptr rtmp, const gr_dft_pre_t P); diff --git a/src/gr_dft/acb.c b/src/gr_dft/acb.c index e7785ab98b..bf94c3fd32 100644 --- a/src/gr_dft/acb.c +++ b/src/gr_dft/acb.c @@ -86,22 +86,23 @@ _gr_dft_acb_ball(acb_ptr w, acb_srcptr v, slong n, int inverse, slong prec) return status; } +/* Scaling and limb-count selection shared by the cyclic and product + fixed-point wrappers, parameterized by the plan's error-bound + evaluator. On failure the caller owns cleaning up its plan. */ +typedef void (*_gr_dft_nfixed_bound_fn)(double * peak, double * err_ulps, + double in_mag, double in_err, const void * plan); + static int -_gr_dft_acb_nfixed_select(gr_dft_acb_pre_struct * Q, slong n, slong prec) +_gr_dft_acb_select_core(slong * nl_out, slong * p1e_out, double * in_mag_out, + double * errulps_out, _gr_dft_nfixed_bound_fn bound, + const void * plan, slong prec) { - int status = GR_SUCCESS; - double peak1, err0, peak; + double peak1, err0, peak, in_mag, errulps; slong nl, p1e; - prec = FLINT_MAX(prec, 2); - - status = _gr_dft_precomp_init_layout(Q->P, n, GR_DFT_ALG_AUTO, 0, 1); - if (status != GR_SUCCESS) - return status; - { int p1exp; - gr_dft_precomp_nfixed_bound(&peak1, &err0, 1.0, 0.0, Q->P); + bound(&peak1, &err0, 1.0, 0.0, plan); /* the bound of a unit input is at least 1 and, for sane transform lengths, far below 2^500; refuse defensively if @@ -113,44 +114,75 @@ _gr_dft_acb_nfixed_select(gr_dft_acb_pre_struct * Q, slong n, slong prec) range: an underflow to zero there would silently understate the bound. */ if (!(peak1 >= 1.0 && peak1 <= ldexp(1.0, 500))) - { - gr_dft_precomp_clear(Q->P); return GR_UNABLE; - } frexp(peak1, &p1exp); /* peak1 < 2^p1exp */ p1e = p1exp; } - Q->p1e = p1e; /* p1exp lies in [1, 1024] here, so this neither under- nor overflows */ - Q->in_mag = ldexp(1.0, -p1e - 1); + in_mag = ldexp(1.0, -p1e - 1); - gr_dft_precomp_nfixed_bound(&peak, &Q->errulps, Q->in_mag, 1.0, Q->P); - if (!(peak < 1.0) || - !(Q->errulps >= 1.0 && Q->errulps <= DBL_MAX)) - { - gr_dft_precomp_clear(Q->P); + bound(&peak, &errulps, in_mag, 1.0, plan); + if (!(peak < 1.0) || !(errulps >= 1.0 && errulps <= DBL_MAX)) return GR_UNABLE; - } /* guard bits: absolute roundoff errulps 2^(-64 nl) should stay below 2^-(prec+2) in_mag */ nl = (prec + FLINT_BITS - 1) / FLINT_BITS; { - double need = prec + 2.0 + log2(Q->errulps) - log2(Q->in_mag); + double need = prec + 2.0 + log2(errulps) - log2(in_mag); slong nl2 = (slong) (need / FLINT_BITS) + 1; nl = FLINT_MAX(nl, nl2); } if (nl > GR_DFT_NFIXED_MAX_NLIMBS) + return GR_UNABLE; + + *nl_out = nl; + *p1e_out = p1e; + *in_mag_out = in_mag; + *errulps_out = errulps; + return GR_SUCCESS; +} + +static void +_gr_dft_nfixed_bound_cyc(double * peak, double * err_ulps, double in_mag, + double in_err, const void * plan) +{ + gr_dft_precomp_nfixed_bound(peak, err_ulps, in_mag, in_err, + (const gr_dft_pre_struct *) plan); +} + +static void +_gr_dft_nfixed_bound_prod(double * peak, double * err_ulps, double in_mag, + double in_err, const void * plan) +{ + gr_dft_prod_precomp_nfixed_bound(peak, err_ulps, in_mag, in_err, + (const gr_dft_prod_pre_struct *) plan); +} + +static int +_gr_dft_acb_nfixed_select(gr_dft_acb_pre_struct * Q, slong n, slong prec) +{ + int status = GR_SUCCESS; + double peak, errulps2; + + prec = FLINT_MAX(prec, 2); + + status = _gr_dft_precomp_init_layout(Q->P, n, GR_DFT_ALG_AUTO, 0, 1); + if (status != GR_SUCCESS) + return status; + + status = _gr_dft_acb_select_core(&Q->nl, &Q->p1e, &Q->in_mag, + &Q->errulps, _gr_dft_nfixed_bound_cyc, Q->P, prec); + if (status != GR_SUCCESS) { gr_dft_precomp_clear(Q->P); - return GR_UNABLE; + return status; } - Q->nl = nl; - status |= gr_dft_ctx_init_nfixed(Q->rctx, nl); - status |= gr_dft_ctx_init_nfixed_complex(Q->cctx, nl); + status |= gr_dft_ctx_init_nfixed(Q->rctx, Q->nl); + status |= gr_dft_ctx_init_nfixed_complex(Q->cctx, Q->nl); status |= _gr_dft_precomp_realize(Q->P, Q->rctx, Q->cctx); if (status != GR_SUCCESS) { @@ -163,7 +195,8 @@ _gr_dft_acb_nfixed_select(gr_dft_acb_pre_struct * Q, slong n, slong prec) /* re-evaluate the error bound with the realized plan (actual complex-multiplication constant and root-table error), for tighter output radii */ - gr_dft_precomp_nfixed_bound(&peak, &Q->errulps, Q->in_mag, 1.0, Q->P); + gr_dft_precomp_nfixed_bound(&peak, &errulps2, Q->in_mag, 1.0, Q->P); + Q->errulps = errulps2; if (!(peak < 1.0) || !(Q->errulps >= 1.0 && Q->errulps <= DBL_MAX)) { @@ -483,13 +516,13 @@ _acb_dft_scan(mag_t R, arf_t t, acb_srcptr v, slong n) } static int -_gr_dft_acb_nfixed_apply(acb_ptr w, acb_srcptr v, int inverse, - const gr_dft_acb_pre_struct * Q, slong prec) +_gr_dft_acb_nfixed_run(acb_ptr w, acb_srcptr v, int inverse, + slong n, slong nl, slong p1e, double errulps, + gr_ctx_struct * rctx, gr_ctx_struct * cctx, + int (*transform)(gr_ptr, gr_srcptr, int, const void *, gr_ctx_t), + const void * plan, slong prec) { int status = GR_SUCCESS; - gr_ctx_struct * rctx = (gr_ctx_struct *) Q->rctx; - gr_ctx_struct * cctx = (gr_ctx_struct *) Q->cctx; - slong n = Q->n; gr_ptr x, y; mag_t R, errmag; arf_t t; @@ -543,7 +576,7 @@ _gr_dft_acb_nfixed_apply(acb_ptr w, acb_srcptr v, int inverse, } /* input scale 2^-e with (m 2^-e) peak(1) <= 1/2 */ - e = EM + Q->p1e + 1; + e = EM + p1e + 1; rsz = rctx->sizeof_elem; @@ -556,17 +589,14 @@ _gr_dft_acb_nfixed_apply(acb_ptr w, acb_srcptr v, int inverse, /* scale (exactly) and truncate the midpoints, writing the fraction limbs directly */ - _acb_dft_conv_run(0, v, (nn_ptr) x, NULL, NULL, 2 * n, e, Q->nl, + _acb_dft_conv_run(0, v, (nn_ptr) x, NULL, NULL, 2 * n, e, nl, rsz, 0, prec, NULL); - if (inverse) - status |= _gr_dft_precomp_raw(y, x, 1, Q->P, cctx); - else - status |= gr_dft_precomp(y, x, Q->P, cctx); + status |= transform(y, x, inverse, plan, cctx); /* absolute roundoff error, back at the input scale */ - mag_set_d(errmag, Q->errulps); - mag_mul_2exp_si(errmag, errmag, -Q->nl * FLINT_BITS + e); + mag_set_d(errmag, errulps); + mag_mul_2exp_si(errmag, errmag, -nl * FLINT_BITS + e); mag_add(errmag, errmag, R); /* the raw inverse omits the 1/n normalization: for power-of-two @@ -594,7 +624,7 @@ _gr_dft_acb_nfixed_apply(acb_ptr w, acb_srcptr v, int inverse, } _acb_dft_conv_run(1, NULL, NULL, w, (nn_srcptr) y, 2 * n, eout, - Q->nl, rsz, den, prec, errmag); + nl, rsz, den, prec, errmag); } flint_free(x); @@ -608,6 +638,37 @@ _gr_dft_acb_nfixed_apply(acb_ptr w, acb_srcptr v, int inverse, return status; } + +static int +_gr_dft_acb_transform_cyc(gr_ptr y, gr_srcptr x, int inverse, + const void * plan, gr_ctx_t cctx) +{ + const gr_dft_pre_struct * P = plan; + + if (inverse) + return _gr_dft_precomp_raw(y, x, 1, P, cctx); + else + return gr_dft_precomp(y, x, P, cctx); +} + +static int +_gr_dft_acb_transform_prod(gr_ptr y, gr_srcptr x, int inverse, + const void * plan, gr_ctx_t cctx) +{ + const gr_dft_prod_pre_struct * P = plan; + + return _gr_dft_prod_precomp_raw(y, x, inverse, P, cctx); +} + +static int +_gr_dft_acb_nfixed_apply(acb_ptr w, acb_srcptr v, int inverse, + const gr_dft_acb_pre_struct * Q, slong prec) +{ + return _gr_dft_acb_nfixed_run(w, v, inverse, Q->n, Q->nl, Q->p1e, + Q->errulps, (gr_ctx_struct *) Q->rctx, (gr_ctx_struct *) Q->cctx, + _gr_dft_acb_transform_cyc, Q->P, prec); +} + static int _gr_dft_acb_nfixed(acb_ptr w, acb_srcptr v, slong n, int inverse, slong prec) { @@ -787,3 +848,273 @@ gr_dft_acb_inverse(acb_ptr w, acb_srcptr v, slong n, slong prec) if (n > 0) GR_IGNORE(_gr_dft_acb(w, v, n, 1, 0, prec)); } + +/* Product-group variant *******************************************/ + +static int +_gr_dft_acb_prod_ball(acb_ptr w, acb_srcptr v, const ulong * cyc, slong num, + int inverse, slong prec) +{ + int status; + gr_ctx_t ctx; + + gr_ctx_init_complex_acb(ctx, prec); + if (inverse) + status = gr_dft_prod_inverse((gr_ptr) w, (gr_srcptr) v, cyc, num, ctx); + else + status = gr_dft_prod((gr_ptr) w, (gr_srcptr) v, cyc, num, ctx); + gr_ctx_clear(ctx); + + return status; +} + +static int +_gr_dft_acb_prod_nfixed_select(gr_dft_acb_prod_pre_struct * Q, + const ulong * cyc, slong num, slong prec) +{ + int status = GR_SUCCESS; + double peak, errulps2; + + prec = FLINT_MAX(prec, 2); + + status = _gr_dft_prod_precomp_init_layout(Q->P, cyc, num, 0, 1); + if (status != GR_SUCCESS) + { + gr_dft_prod_precomp_clear(Q->P); + return status; + } + + status = _gr_dft_acb_select_core(&Q->nl, &Q->p1e, &Q->in_mag, + &Q->errulps, _gr_dft_nfixed_bound_prod, Q->P, prec); + if (status != GR_SUCCESS) + { + gr_dft_prod_precomp_clear(Q->P); + return status; + } + + status |= gr_dft_ctx_init_nfixed(Q->rctx, Q->nl); + status |= gr_dft_ctx_init_nfixed_complex(Q->cctx, Q->nl); + status |= _gr_dft_prod_precomp_realize(Q->P, Q->rctx, Q->cctx); + if (status != GR_SUCCESS) + { + gr_dft_prod_precomp_clear(Q->P); + gr_ctx_clear(Q->rctx); + gr_ctx_clear(Q->cctx); + return status; + } + + /* re-evaluate with the realized component plans */ + gr_dft_prod_precomp_nfixed_bound(&peak, &errulps2, Q->in_mag, 1.0, Q->P); + Q->errulps = errulps2; + if (!(peak < 1.0) || + !(Q->errulps >= 1.0 && Q->errulps <= DBL_MAX)) + { + gr_dft_prod_precomp_clear(Q->P); + gr_ctx_clear(Q->rctx); + gr_ctx_clear(Q->cctx); + return GR_UNABLE; + } + + return GR_SUCCESS; +} + +static int +_gr_dft_acb_prod_nfixed_apply(acb_ptr w, acb_srcptr v, int inverse, + const gr_dft_acb_prod_pre_struct * Q, slong prec) +{ + return _gr_dft_acb_nfixed_run(w, v, inverse, Q->n, Q->nl, Q->p1e, + Q->errulps, (gr_ctx_struct *) Q->rctx, (gr_ctx_struct *) Q->cctx, + _gr_dft_acb_transform_prod, Q->P, prec); +} + +int +gr_dft_acb_prod_precomp_init(gr_dft_acb_prod_pre_t Q, const ulong * cyc, + slong num, slong prec) +{ + int status = GR_SUCCESS; + slong a, n = 1; + + for (a = 0; a < num; a++) + { + if (cyc[a] == 0 || (ulong) n > (WORD_MAX / 8) / cyc[a]) + return GR_DOMAIN; + n *= (slong) cyc[a]; + } + + prec = FLINT_MAX(prec, 2); + + Q->n = n; + Q->num = num; + Q->prec = prec; + Q->cyc = flint_malloc(FLINT_MAX(num, 1) * sizeof(ulong)); + for (a = 0; a < num; a++) + Q->cyc[a] = cyc[a]; + + status = _gr_dft_acb_prod_nfixed_select(Q, cyc, num, prec); + if (status == GR_SUCCESS) + { + Q->which = 2; + return GR_SUCCESS; + } + + /* ball arithmetic */ + gr_ctx_init_real_arb(Q->rctx, prec); + gr_ctx_init_complex_acb(Q->cctx, prec); + status = gr_dft_prod_precomp_init(Q->P, cyc, num, 0, + (gr_ctx_struct *) Q->cctx); + if (status != GR_SUCCESS) + { + gr_dft_prod_precomp_clear(Q->P); + gr_ctx_clear(Q->rctx); + gr_ctx_clear(Q->cctx); + flint_free(Q->cyc); + Q->cyc = NULL; + Q->which = 0; + return status; + } + Q->which = 1; + + return GR_SUCCESS; +} + +void +gr_dft_acb_prod_precomp_clear(gr_dft_acb_prod_pre_t Q) +{ + if (Q->which != 0) + { + gr_dft_prod_precomp_clear(Q->P); + gr_ctx_clear(Q->rctx); + gr_ctx_clear(Q->cctx); + Q->which = 0; + } + flint_free(Q->cyc); + Q->cyc = NULL; +} + +int +_gr_dft_acb_prod_precomp(acb_ptr w, acb_srcptr v, int inverse, + const gr_dft_acb_prod_pre_t Q, slong prec) +{ + int status; + slong n = Q->n; + + prec = FLINT_MAX(prec, 2); + + if (Q->which == 2) + { + status = _gr_dft_acb_prod_nfixed_apply(w, v, inverse, Q, prec); + + /* inputs the fixed-point path cannot handle: fall back to a + one-shot ball transform */ + if (status != GR_SUCCESS) + status = _gr_dft_acb_prod_ball(w, v, Q->cyc, Q->num, inverse, + prec); + } + else if (Q->which == 1) + { + if (inverse) + status = gr_dft_prod_inverse_precomp((gr_ptr) w, (gr_srcptr) v, + Q->P, (gr_ctx_struct *) Q->cctx); + else + status = gr_dft_prod_precomp((gr_ptr) w, (gr_srcptr) v, + Q->P, (gr_ctx_struct *) Q->cctx); + } + else + { + status = GR_UNABLE; + } + + if (status != GR_SUCCESS) + { + slong j; + for (j = 0; j < n; j++) + acb_indeterminate(w + j); + } + + return status; +} + +void +gr_dft_acb_prod_precomp(acb_ptr w, acb_srcptr v, + const gr_dft_acb_prod_pre_t Q, slong prec) +{ + GR_IGNORE(_gr_dft_acb_prod_precomp(w, v, 0, Q, prec)); +} + +void +gr_dft_acb_prod_inverse_precomp(acb_ptr w, acb_srcptr v, + const gr_dft_acb_prod_pre_t Q, slong prec) +{ + GR_IGNORE(_gr_dft_acb_prod_precomp(w, v, 1, Q, prec)); +} + +int +_gr_dft_acb_prod(acb_ptr w, acb_srcptr v, const ulong * cyc, slong num, + int inverse, int which, slong prec) +{ + int status = GR_SUCCESS; + gr_dft_acb_prod_pre_struct Q; + slong a, n = 1; + + for (a = 0; a < num; a++) + { + if (cyc[a] == 0 || (ulong) n > (WORD_MAX / 8) / cyc[a]) + return GR_DOMAIN; + n *= (slong) cyc[a]; + } + + prec = FLINT_MAX(prec, 2); + + if (which == 0 || which == 2) + { + Q.n = n; + Q.num = num; + Q.prec = prec; + Q.which = 2; + + status = _gr_dft_acb_prod_nfixed_select(&Q, cyc, num, prec); + if (status == GR_SUCCESS) + { + status = _gr_dft_acb_prod_nfixed_apply(w, v, inverse, &Q, prec); + + gr_dft_prod_precomp_clear(Q.P); + gr_ctx_clear(Q.rctx); + gr_ctx_clear(Q.cctx); + } + + if (status == GR_SUCCESS || which == 2) + { + if (status != GR_SUCCESS) + { + slong j; + for (j = 0; j < n; j++) + acb_indeterminate(w + j); + } + return status; + } + } + + status = _gr_dft_acb_prod_ball(w, v, cyc, num, inverse, prec); + if (status != GR_SUCCESS) + { + slong j; + for (j = 0; j < n; j++) + acb_indeterminate(w + j); + } + + return status; +} + +void +gr_dft_acb_prod(acb_ptr w, acb_srcptr v, const ulong * cyc, slong num, + slong prec) +{ + GR_IGNORE(_gr_dft_acb_prod(w, v, cyc, num, 0, 0, prec)); +} + +void +gr_dft_acb_prod_inverse(acb_ptr w, acb_srcptr v, const ulong * cyc, + slong num, slong prec) +{ + GR_IGNORE(_gr_dft_acb_prod(w, v, cyc, num, 1, 0, prec)); +} diff --git a/src/gr_dft/dirichlet.c b/src/gr_dft/dirichlet.c new file mode 100644 index 0000000000..11fefb183e --- /dev/null +++ b/src/gr_dft/dirichlet.c @@ -0,0 +1,149 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include "dirichlet.h" +#include "acb.h" +#include "gr.h" +#include "gr_vec.h" +#include "gr_dft.h" + +/* Generic-ring counterparts of acb_dirichlet_dft_index and + acb_dirichlet_dft: the group of Dirichlet characters mod q is the + product of the cyclic component groups given by the Conrey + decomposition, so its DFT is a product DFT over the component + sizes, with a gather/scatter between number indexing and + lexicographic Conrey indexing for the non-index variant. */ + +static void +_gr_dft_dirichlet_cyc(ulong * cyc, const dirichlet_group_t G) +{ + slong k; + + for (k = 0; k < G->num; k++) + cyc[k] = G->P[k].phi.n; +} + +/* dft, lexicographic conrey indexing, array size G->phi_q */ +int +gr_dft_dirichlet_index(gr_ptr w, gr_srcptr v, const dirichlet_group_t G, + gr_ctx_t ctx) +{ + int status; + ulong * cyc; + + if (G->phi_q == 1) + return gr_set(w, v, ctx); + + cyc = flint_malloc(G->num * sizeof(ulong)); + _gr_dft_dirichlet_cyc(cyc, G); + status = gr_dft_prod(w, v, cyc, G->num, ctx); + flint_free(cyc); + + return status; +} + +/* dft, number indexing, array size G->q */ +int +gr_dft_dirichlet(gr_ptr w, gr_srcptr v, const dirichlet_group_t G, + gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + slong sz = ctx->sizeof_elem; + ulong i, len = G->phi_q; + gr_ptr t; + dirichlet_char_t x; + + /* the trivial group: the Conrey number of the unique character + is 1, outside an array of size q = 1, so the generic + gather/scatter below must not run */ + if (G->q == 1) + return gr_set(w, v, ctx); + + t = gr_heap_init_vec(len, ctx); + dirichlet_char_init(x, G); + + dirichlet_char_one(x, G); + for (i = 0; i < len; i++) + { + status |= gr_set(GR_ENTRY(t, i, sz), GR_ENTRY(v, x->n, sz), ctx); + dirichlet_char_next(x, G); + } + + status |= gr_dft_dirichlet_index(t, t, G, ctx); + + dirichlet_char_one(x, G); + for (i = 0; i < len; i++) + { + status |= gr_set(GR_ENTRY(w, x->n, sz), GR_ENTRY(t, i, sz), ctx); + dirichlet_char_next(x, G); + } + + dirichlet_char_clear(x); + gr_heap_clear_vec(t, len, ctx); + + return status; +} + +/* acb variants routed through the fixed-point product transform */ +void +gr_dft_acb_dirichlet_index(acb_ptr w, acb_srcptr v, + const dirichlet_group_t G, slong prec) +{ + ulong * cyc; + + if (G->phi_q == 1) + { + acb_set(w, v); + return; + } + + cyc = flint_malloc(G->num * sizeof(ulong)); + _gr_dft_dirichlet_cyc(cyc, G); + gr_dft_acb_prod(w, v, cyc, G->num, prec); + flint_free(cyc); +} + +void +gr_dft_acb_dirichlet(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, + slong prec) +{ + ulong i, len = G->phi_q; + acb_ptr t; + dirichlet_char_t x; + + if (G->q == 1) + { + acb_set(w, v); + return; + } + + t = _acb_vec_init(len); + dirichlet_char_init(x, G); + + dirichlet_char_one(x, G); + for (i = 0; i < len; i++) + { + acb_set(t + i, v + x->n); + dirichlet_char_next(x, G); + } + + gr_dft_acb_dirichlet_index(t, t, G, prec); + + dirichlet_char_one(x, G); + for (i = 0; i < len; i++) + { + acb_set(w + x->n, t + i); + dirichlet_char_next(x, G); + } + + dirichlet_char_clear(x); + _acb_vec_clear(t, len); +} diff --git a/src/gr_dft/prod.c b/src/gr_dft/prod.c new file mode 100644 index 0000000000..5acf03bbeb --- /dev/null +++ b/src/gr_dft/prod.c @@ -0,0 +1,466 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include "thread_pool.h" +#include "thread_support.h" +#include "gr.h" +#include "gr_vec.h" +#include "gr_dft.h" + +/* DFT on a product of cyclic groups Z/cyc[0] x ... x Z/cyc[num-1], + with row-major indexing (the first component varies slowest), as in + acb_dft_prod: one cyclic DFT along each axis, with no twiddle + factors in between. Along axis a, the lines have stride equal to + the suffix product of the later component lengths; strided lines + are gathered into a contiguous buffer, transformed with the + component plan, and scattered back, while the lines of the last + axis (stride 1) are transformed in place. The lines of one axis are + independent, which is what the threaded version parallelizes over. + + Component plans are shared between axes of equal length (Dirichlet + groups routinely repeat component sizes). */ + +static void +_gr_dft_prod_zero(gr_dft_prod_pre_t P) +{ + P->n = 1; + P->num = 0; + P->cyc = NULL; + P->num_plans = 0; + P->plan_of = NULL; + P->plans = NULL; + P->flags = 0; + P->serial_block = 0; +} + +int +_gr_dft_prod_precomp_init_layout(gr_dft_prod_pre_t P, const ulong * cyc, + slong num, int flags, int complex_mode) +{ + int status = GR_SUCCESS; + slong a, b; + + _gr_dft_prod_zero(P); + + P->num = num; + P->flags = flags; + P->cyc = flint_malloc(FLINT_MAX(num, 1) * sizeof(ulong)); + P->plan_of = flint_malloc(FLINT_MAX(num, 1) * sizeof(slong)); + P->plans = flint_malloc(FLINT_MAX(num, 1) * sizeof(gr_dft_pre_struct)); + + for (a = 0; a < num; a++) + { + P->cyc[a] = cyc[a]; + P->n *= cyc[a]; + + P->plan_of[a] = -1; + for (b = 0; b < a; b++) + { + if (cyc[b] == cyc[a]) + { + P->plan_of[a] = P->plan_of[b]; + break; + } + } + + if (P->plan_of[a] == -1) + { + P->plan_of[a] = P->num_plans; + status |= _gr_dft_precomp_init_layout(P->plans + P->num_plans, + cyc[a], GR_DFT_ALG_AUTO, flags, complex_mode); + P->num_plans++; + } + } + + return status; +} + +int +_gr_dft_prod_precomp_realize(gr_dft_prod_pre_t P, gr_ctx_struct * real_ctx, + gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + slong k; + + for (k = 0; k < P->num_plans && status == GR_SUCCESS; k++) + status |= _gr_dft_precomp_realize(P->plans + k, real_ctx, ctx); + + return status; +} + +int +gr_dft_prod_precomp_init(gr_dft_prod_pre_t P, const ulong * cyc, slong num, + int flags, gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + slong a, b; + + _gr_dft_prod_zero(P); + + P->num = num; + P->flags = flags; + P->cyc = flint_malloc(FLINT_MAX(num, 1) * sizeof(ulong)); + P->plan_of = flint_malloc(FLINT_MAX(num, 1) * sizeof(slong)); + P->plans = flint_malloc(FLINT_MAX(num, 1) * sizeof(gr_dft_pre_struct)); + + for (a = 0; a < num && status == GR_SUCCESS; a++) + { + P->cyc[a] = cyc[a]; + P->n *= cyc[a]; + + P->plan_of[a] = -1; + for (b = 0; b < a; b++) + { + if (cyc[b] == cyc[a]) + { + P->plan_of[a] = P->plan_of[b]; + break; + } + } + + if (P->plan_of[a] == -1) + { + P->plan_of[a] = P->num_plans; + status |= gr_dft_precomp_init(P->plans + P->num_plans, + cyc[a], GR_DFT_ALG_AUTO, flags, ctx); + P->num_plans++; + } + } + + return status; +} + +/* Variant for rings without canonical roots of unity: w is a root of + unity of the given order, which every component length must divide; + the component of length m uses w^(order/m). */ +int +gr_dft_prod_precomp_init_root(gr_dft_prod_pre_t P, gr_srcptr w, ulong order, + const ulong * cyc, slong num, int flags, gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + slong a, b; + gr_ptr wa; + + _gr_dft_prod_zero(P); + + for (a = 0; a < num; a++) + if (cyc[a] == 0 || order % cyc[a] != 0) + return GR_DOMAIN; + + P->num = num; + P->flags = flags; + P->cyc = flint_malloc(FLINT_MAX(num, 1) * sizeof(ulong)); + P->plan_of = flint_malloc(FLINT_MAX(num, 1) * sizeof(slong)); + P->plans = flint_malloc(FLINT_MAX(num, 1) * sizeof(gr_dft_pre_struct)); + + GR_TMP_INIT(wa, ctx); + + for (a = 0; a < num && status == GR_SUCCESS; a++) + { + P->cyc[a] = cyc[a]; + P->n *= cyc[a]; + + P->plan_of[a] = -1; + for (b = 0; b < a; b++) + { + if (cyc[b] == cyc[a]) + { + P->plan_of[a] = P->plan_of[b]; + break; + } + } + + if (P->plan_of[a] == -1) + { + P->plan_of[a] = P->num_plans; + status |= gr_pow_ui(wa, w, order / cyc[a], ctx); + if (status == GR_SUCCESS) + status |= gr_dft_precomp_init_root(P->plans + P->num_plans, + wa, cyc[a], GR_DFT_ALG_AUTO, flags, ctx); + P->num_plans++; + } + } + + GR_TMP_CLEAR(wa, ctx); + return status; +} + +void +gr_dft_prod_precomp_clear(gr_dft_prod_pre_t P) +{ + slong k; + + for (k = 0; k < P->num_plans; k++) + gr_dft_precomp_clear(P->plans + k); + + flint_free(P->cyc); + flint_free(P->plan_of); + flint_free(P->plans); + _gr_dft_prod_zero(P); +} + +void +gr_dft_prod_precomp_set_serial_block(gr_dft_prod_pre_t P, slong serial_block) +{ + P->serial_block = serial_block; +} + +/* Compose the fixed-point error bounds of the component transforms: + the output of one axis, exact copies aside, is the input of the + next. */ +void +gr_dft_prod_precomp_nfixed_bound(double * peak, double * err_ulps, + double in_mag, double in_err, const gr_dft_prod_pre_t P) +{ + double mag = in_mag, err = in_err, mag2, err2; + slong a; + + for (a = 0; a < P->num; a++) + { + if (P->cyc[a] == 1) + continue; + + gr_dft_precomp_nfixed_bound(&mag2, &err2, mag, err, + P->plans + P->plan_of[a]); + mag = mag2; + err = err2; + } + + *peak = mag; + *err_ulps = err; +} + +/* transform the contiguous lines [llo, lhi) of one axis; tmp is a + line buffer of m elements for the strided (inner > 1) case */ +static int +_gr_dft_prod_lines(gr_ptr x, int inverse, ulong m, ulong inner, + ulong llo, ulong lhi, const gr_dft_pre_struct * plan, gr_ptr tmp, + gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + slong sz = ctx->sizeof_elem; + ulong L, t; + + if (inner == 1) + { + for (L = llo; L < lhi; L++) + { + gr_ptr line = GR_ENTRY(x, (slong) (L * m), sz); + status |= _gr_dft_precomp_raw(line, line, inverse, plan, ctx); + } + } + else + { + for (L = llo; L < lhi; L++) + { + ulong o = L / inner, i = L % inner; + gr_ptr base = GR_ENTRY(x, (slong) (o * m * inner + i), sz); + + for (t = 0; t < m; t++) + status |= gr_set(GR_ENTRY(tmp, t, sz), + GR_ENTRY(base, (slong) (t * inner), sz), ctx); + + status |= _gr_dft_precomp_raw(tmp, tmp, inverse, plan, ctx); + + for (t = 0; t < m; t++) + status |= gr_set(GR_ENTRY(base, (slong) (t * inner), sz), + GR_ENTRY(tmp, t, sz), ctx); + } + } + + return status; +} + +typedef struct +{ + gr_ptr x; + int inverse; + ulong m; + ulong inner; + ulong llo; + ulong lhi; + const gr_dft_pre_struct * plan; + gr_ctx_struct * ctx; + int status; +} +_gr_dft_prod_work_t; + +static void +_gr_dft_prod_worker(void * arg) +{ + _gr_dft_prod_work_t * a = arg; + gr_ctx_struct * ctx = a->ctx; + gr_ptr tmp = NULL; + + if (a->inner > 1) + tmp = gr_heap_init_vec(a->m, ctx); + + a->status = _gr_dft_prod_lines(a->x, a->inverse, a->m, a->inner, + a->llo, a->lhi, a->plan, tmp, ctx); + + if (tmp != NULL) + gr_heap_clear_vec(tmp, a->m, ctx); +} + +int +_gr_dft_prod_precomp_raw(gr_ptr res, gr_srcptr vec, int inverse, + const gr_dft_prod_pre_t P, gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + ulong n = P->n; + slong a; + ulong inner; + thread_pool_handle * handles = NULL; + slong num_workers = 0; + slong blk, max_chunks; + + if (res != vec) + status |= _gr_vec_set(res, vec, n, ctx); + + if (P->num == 0 || n <= 1) + return status; + + blk = (P->serial_block > 0) ? P->serial_block + : GR_DFT_SERIAL_BLOCK_DEFAULT; + max_chunks = (slong) (n / (ulong) blk); + + if (max_chunks >= 2 && flint_get_num_threads() > 1 && + gr_ctx_is_threadsafe(ctx) == T_TRUE) + num_workers = flint_request_threads(&handles, max_chunks); + + /* suffix products give the line strides */ + inner = n; + for (a = 0; a < P->num && status == GR_SUCCESS; a++) + { + ulong m = P->cyc[a]; + ulong nlines; + const gr_dft_pre_struct * plan = P->plans + P->plan_of[a]; + + inner /= m; + + if (m == 1) + continue; + + nlines = n / m; + + if (num_workers < 1 || nlines < 2) + { + gr_ptr tmp = NULL; + + if (inner > 1) + tmp = gr_heap_init_vec(m, ctx); + status |= _gr_dft_prod_lines(res, inverse, m, inner, + 0, nlines, plan, tmp, ctx); + if (tmp != NULL) + gr_heap_clear_vec(tmp, m, ctx); + } + else + { + slong i, nchunks = FLINT_MIN(num_workers + 1, (slong) nlines); + _gr_dft_prod_work_t * args; + + args = flint_malloc(nchunks * sizeof(_gr_dft_prod_work_t)); + + for (i = 0; i < nchunks; i++) + { + args[i].x = res; + args[i].inverse = inverse; + args[i].m = m; + args[i].inner = inner; + args[i].llo = (nlines * (ulong) i) / (ulong) nchunks; + args[i].lhi = (nlines * (ulong) (i + 1)) / (ulong) nchunks; + args[i].plan = plan; + args[i].ctx = ctx; + args[i].status = GR_SUCCESS; + + if (i < nchunks - 1) + thread_pool_wake(global_thread_pool, handles[i], 0, + _gr_dft_prod_worker, &args[i]); + } + + _gr_dft_prod_worker(&args[nchunks - 1]); + + for (i = 0; i < nchunks - 1; i++) + thread_pool_wait(global_thread_pool, handles[i]); + + for (i = 0; i < nchunks; i++) + status |= args[i].status; + + flint_free(args); + } + } + + if (handles != NULL) + flint_give_back_threads(handles, num_workers); + + return status; +} + +int +gr_dft_prod_precomp(gr_ptr res, gr_srcptr vec, const gr_dft_prod_pre_t P, + gr_ctx_t ctx) +{ + return _gr_dft_prod_precomp_raw(res, vec, 0, P, ctx); +} + +/* the raw inverse omits the 1/n normalization of each component; + apply the full 1/n once */ +int +gr_dft_prod_inverse_precomp(gr_ptr res, gr_srcptr vec, + const gr_dft_prod_pre_t P, gr_ctx_t ctx) +{ + int status; + gr_ptr c; + + status = _gr_dft_prod_precomp_raw(res, vec, 1, P, ctx); + + if (status == GR_SUCCESS && P->n > 1) + { + GR_TMP_INIT(c, ctx); + status |= gr_set_ui(c, P->n, ctx); + status |= gr_inv(c, c, ctx); + if (status == GR_SUCCESS) + status |= _gr_vec_mul_scalar(res, res, (slong) P->n, c, ctx); + GR_TMP_CLEAR(c, ctx); + } + + return status; +} + +int +gr_dft_prod(gr_ptr res, gr_srcptr vec, const ulong * cyc, slong num, + gr_ctx_t ctx) +{ + int status; + gr_dft_prod_pre_t P; + + status = gr_dft_prod_precomp_init(P, cyc, num, 0, ctx); + if (status == GR_SUCCESS) + status = gr_dft_prod_precomp(res, vec, P, ctx); + gr_dft_prod_precomp_clear(P); + + return status; +} + +int +gr_dft_prod_inverse(gr_ptr res, gr_srcptr vec, const ulong * cyc, slong num, + gr_ctx_t ctx) +{ + int status; + gr_dft_prod_pre_t P; + + status = gr_dft_prod_precomp_init(P, cyc, num, 0, ctx); + if (status == GR_SUCCESS) + status = gr_dft_prod_inverse_precomp(res, vec, P, ctx); + gr_dft_prod_precomp_clear(P); + + return status; +} diff --git a/src/gr_dft/profile/p-gr_dft_acb.c b/src/gr_dft/profile/p-gr_dft_acb.c index 6ef3af6eaa..f42916a497 100644 --- a/src/gr_dft/profile/p-gr_dft_acb.c +++ b/src/gr_dft/profile/p-gr_dft_acb.c @@ -13,13 +13,16 @@ #include #include "profiler.h" #include "acb.h" -#include "acb_dft.h" #include "gr.h" #include "gr_dft.h" -/* Compares acb_dft against the drop-in gr_dft_acb, reporting the - precomputation and per-transform times separately. Optional - arguments: max_n max_prec. */ +/* Compares the two internal paths of the gr_dft_acb drop-in: ball + arithmetic (a Karatsuba plan over arb/acb contexts, as used by the + fallback) against the default fixed-point path, reporting the + precomputation and per-transform times separately. The path column + shows which arithmetic the automatic plan selected at this length + and precision (it falls back to ball when the fixed-point limb + count would be excessive). Optional arguments: max_n max_prec. */ static double time_ms(void (*f)(void *), void * arg, slong min_reps) @@ -47,26 +50,39 @@ typedef struct acb_srcptr v; slong n; slong prec; - acb_dft_pre_struct * apre; + gr_ctx_struct * rctx; + gr_ctx_struct * cctx; + gr_dft_pre_struct * bpre; gr_dft_acb_pre_struct * qpre; } args_t; -static void f_acb_pre(void * p) +/* ball path: contexts and Karatsuba plan, as in the which = 1 branch + of gr_dft_acb_precomp_init */ +static void f_ball_pre(void * p) { args_t * a = p; - acb_dft_pre_t pre; - acb_dft_precomp_init(pre, a->n, a->prec); - acb_dft_precomp_clear(pre); + gr_ctx_t rctx, cctx; + gr_dft_pre_t pre; + + gr_ctx_init_real_arb(rctx, a->prec); + gr_ctx_init_complex_acb(cctx, a->prec); + GR_IGNORE(gr_dft_precomp_init_karatsuba(pre, a->n, GR_DFT_ALG_AUTO, 0, + rctx, cctx)); + gr_dft_precomp_clear(pre); + gr_ctx_clear(rctx); + gr_ctx_clear(cctx); } -static void f_acb_dft(void * p) +static void f_ball_dft(void * p) { args_t * a = p; - acb_dft_precomp(a->w, a->v, a->apre, a->prec); + GR_IGNORE(gr_dft_precomp((gr_ptr) a->w, (gr_srcptr) a->v, a->bpre, + a->cctx)); } -static void f_gr_pre(void * p) +/* fixed-point path: the automatic drop-in plan */ +static void f_nfixed_pre(void * p) { args_t * a = p; gr_dft_acb_pre_t Q; @@ -74,7 +90,7 @@ static void f_gr_pre(void * p) gr_dft_acb_precomp_clear(Q); } -static void f_gr_dft(void * p) +static void f_nfixed_dft(void * p) { args_t * a = p; gr_dft_acb_precomp(a->w, a->v, (const gr_dft_acb_pre_struct *) a->qpre, @@ -109,18 +125,20 @@ main(int argc, char * argv[]) continue; flint_printf("prec = %wd\n", prec); - flint_printf("%8s %12s %12s %12s %12s %10s\n", "n", - "acb pre", "acb dft", "gr pre", "gr dft", "dft ratio"); + flint_printf("%8s %12s %12s %12s %12s %10s %7s\n", "n", + "ball pre", "ball dft", "auto pre", "auto dft", "dft ratio", + "path"); for (ni = 0; ns[ni] != 0; ni++) { ulong n = ns[ni]; slong j; acb_ptr v, w; - acb_dft_pre_t apre; + gr_ctx_t rctx, cctx; + gr_dft_pre_t bpre; gr_dft_acb_pre_t qpre; args_t a; - double t_acb_pre, t_acb_dft, t_gr_pre, t_gr_dft; + double t_ball_pre, t_ball_dft, t_nf_pre, t_nf_dft; if ((slong) n > max_n) continue; @@ -130,26 +148,34 @@ main(int argc, char * argv[]) for (j = 0; j < (slong) n; j++) acb_urandom(v + j, state, prec); - acb_dft_precomp_init(apre, n, prec); + gr_ctx_init_real_arb(rctx, prec); + gr_ctx_init_complex_acb(cctx, prec); + GR_IGNORE(gr_dft_precomp_init_karatsuba(bpre, n, + GR_DFT_ALG_AUTO, 0, rctx, cctx)); GR_IGNORE(gr_dft_acb_precomp_init(qpre, n, prec)); a.w = w; a.v = v; a.n = n; a.prec = prec; - a.apre = apre; + a.rctx = rctx; + a.cctx = cctx; + a.bpre = bpre; a.qpre = qpre; - t_acb_pre = time_ms(f_acb_pre, &a, 1); - t_acb_dft = time_ms(f_acb_dft, &a, 1); - t_gr_pre = time_ms(f_gr_pre, &a, 1); - t_gr_dft = time_ms(f_gr_dft, &a, 1); + t_ball_pre = time_ms(f_ball_pre, &a, 1); + t_ball_dft = time_ms(f_ball_dft, &a, 1); + t_nf_pre = time_ms(f_nfixed_pre, &a, 1); + t_nf_dft = time_ms(f_nfixed_dft, &a, 1); - flint_printf("%8wu %12.4g %12.4g %12.4g %12.4g %9.2fx\n", - n, t_acb_pre, t_acb_dft, t_gr_pre, t_gr_dft, - t_acb_dft / t_gr_dft); + flint_printf("%8wu %12.4g %12.4g %12.4g %12.4g %9.2fx %7s\n", + n, t_ball_pre, t_ball_dft, t_nf_pre, t_nf_dft, + t_ball_dft / t_nf_dft, + (qpre->which == 2) ? "nfixed" : "ball"); - acb_dft_precomp_clear(apre); + gr_dft_precomp_clear(bpre); + gr_ctx_clear(rctx); + gr_ctx_clear(cctx); gr_dft_acb_precomp_clear(qpre); _acb_vec_clear(v, n); _acb_vec_clear(w, n); diff --git a/src/gr_dft/test/main.c b/src/gr_dft/test/main.c index 6f5574267e..6fd1c51968 100644 --- a/src/gr_dft/test/main.c +++ b/src/gr_dft/test/main.c @@ -15,7 +15,9 @@ #include "t-dft.c" #include "t-dft_inverse.c" #include "t-acb_dft.c" +#include "t-dirichlet_dft.c" #include "t-nfixed.c" +#include "t-prod.c" /* Array of test functions ***************************************************/ @@ -24,7 +26,9 @@ test_struct tests[] = TEST_FUNCTION(gr_dft), TEST_FUNCTION(gr_dft_inverse), TEST_FUNCTION(gr_dft_acb), + TEST_FUNCTION(gr_dft_dirichlet), TEST_FUNCTION(gr_dft_nfixed), + TEST_FUNCTION(gr_dft_prod), }; /* main function *************************************************************/ diff --git a/src/gr_dft/test/t-acb_dft.c b/src/gr_dft/test/t-acb_dft.c index 9c55cfaa52..3c6cf05a0e 100644 --- a/src/gr_dft/test/t-acb_dft.c +++ b/src/gr_dft/test/t-acb_dft.c @@ -13,7 +13,6 @@ #include "test_helpers.h" #include "ulong_extras.h" #include "acb.h" -#include "acb_dft.h" #include "gr.h" #include "gr_dft.h" @@ -46,13 +45,13 @@ TEST_FUNCTION_START(gr_dft_acb, state) acb_get_mid(v + j, v + j); } - /* reference: acb_dft at elevated precision on the midpoints */ + /* reference: this module's ball path (plain gr_dft over an + acb context, itself tested against naive references in + t-dft) at elevated precision on the midpoints */ for (j = 0; j < n; j++) acb_get_mid(ref + j, v + j); - if (inverse) - acb_dft_inverse(ref, ref, n, refprec); - else - acb_dft(ref, ref, n, refprec); + if (_gr_dft_acb(ref, ref, n, inverse, 1, refprec) != GR_SUCCESS) + TEST_FUNCTION_FAIL("reference failed, n = %wd\n", n); /* both internal paths must contain the reference and agree */ status = _gr_dft_acb(w1, v, n, inverse, 2, prec); diff --git a/src/gr_dft/test/t-dirichlet_dft.c b/src/gr_dft/test/t-dirichlet_dft.c new file mode 100644 index 0000000000..50969c4130 --- /dev/null +++ b/src/gr_dft/test/t-dirichlet_dft.c @@ -0,0 +1,224 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include "test_helpers.h" +#include "arb.h" +#include "acb.h" +#include "dirichlet.h" +#include "gr.h" +#include "gr_vec.h" +#include "gr_dft.h" + +/* Independent naive reference (adapted from acb_dirichlet's t-dft.c, + but not using acb_dft or acb_dirichlet, which are to become + wrappers around this module): + + w[x] = sum_y conj(chi_x(y)) v[y] + + over the characters in lexicographic Conrey order, with the + pairing values chi_x(y) = e(dirichlet_pairing_char(G, x, y) / + G->expo) evaluated through a directly computed table of roots of + unity. */ +static void +_dirichlet_dft_index_naive(acb_ptr w, acb_srcptr v, + const dirichlet_group_t G, slong prec) +{ + slong i, j, len = G->phi_q; + ulong expo = G->expo, k; + acb_ptr roots; + acb_t t; + dirichlet_char_t x, y; + + /* roots[k] = e(-k / expo), the conjugate pairing values */ + roots = _acb_vec_init(expo); + acb_init(t); + for (k = 0; k < expo; k++) + { + arb_set_si(acb_realref(t), -2 * (slong) k); + arb_div_ui(acb_realref(t), acb_realref(t), expo, prec + 10); + arb_zero(acb_imagref(t)); + acb_exp_pi_i(roots + k, t, prec + 10); + } + + dirichlet_char_init(x, G); + dirichlet_char_init(y, G); + + dirichlet_char_one(x, G); + for (i = 0; i < len; i++) + { + acb_zero(w + i); + dirichlet_char_one(y, G); + for (j = 0; j < len; j++) + { + acb_addmul(w + i, + roots + dirichlet_pairing_char(G, x, y), v + j, prec); + dirichlet_char_next(y, G); + } + dirichlet_char_next(x, G); + } + + dirichlet_char_clear(x); + dirichlet_char_clear(y); + _acb_vec_clear(roots, expo); + acb_clear(t); +} + +TEST_FUNCTION_START(gr_dft_dirichlet, state) +{ + /* moduli covering odd, even (q_even = 2, 4 and >= 8, the last + giving two even Conrey components), prime, prime power and + composite structure */ + slong k; + slong nq = 16; + ulong q[16] = { 2, 3, 4, 5, 6, 8, 10, 15, 16, 23, 24, 30, 59, 308, + 335, 961 }; + + for (k = 0; k < nq; k++) + { + slong prec = 100; + dirichlet_group_t G; + acb_ptr v, vq, w1, w2; + ulong i, len; + gr_ctx_t ctx; + dirichlet_char_t x; + int status = GR_SUCCESS; + + dirichlet_group_init(G, q[k]); + len = G->phi_q; + + v = _acb_vec_init(len); + vq = _acb_vec_init(q[k]); + w1 = _acb_vec_init(len); + w2 = _acb_vec_init(FLINT_MAX(q[k], len)); + + for (i = 0; i < len; i++) + acb_randtest_precise(v + i, state, prec, 0); + for (i = 0; i < q[k]; i++) + acb_randtest_precise(vq + i, state, prec, 0); + + _dirichlet_dft_index_naive(w1, v, G, prec); + + /* conrey-indexed transform: the generic gr entry point on an + acb context, and the fixed-point wrapper */ + gr_ctx_init_complex_acb(ctx, prec); + status |= gr_dft_dirichlet_index((gr_ptr) w2, (gr_srcptr) v, G, ctx); + if (status != GR_SUCCESS) + TEST_FUNCTION_FAIL("index status, q = %wu\n", q[k]); + for (i = 0; i < len; i++) + if (!acb_overlaps(w1 + i, w2 + i)) + TEST_FUNCTION_FAIL("generic index vs naive, q = %wu " + "i = %wu\n", q[k], i); + gr_ctx_clear(ctx); + + gr_dft_acb_dirichlet_index(w2, v, G, prec); + for (i = 0; i < len; i++) + { + if (!acb_overlaps(w1 + i, w2 + i)) + TEST_FUNCTION_FAIL("acb index vs naive, q = %wu i = %wu\n", + q[k], i); + if (acb_rel_accuracy_bits(w2 + i) < 30 && len > 1) + TEST_FUNCTION_FAIL("index accuracy, q = %wu i = %wu " + "(%wd bits)\n", q[k], i, + acb_rel_accuracy_bits(w2 + i)); + } + + /* number-indexed transform against its definition: gather the + values in Conrey order, transform in index order, scatter */ + { + acb_ptr t1; + + t1 = _acb_vec_init(len); + dirichlet_char_init(x, G); + + dirichlet_char_one(x, G); + for (i = 0; i < len; i++) + { + acb_set(t1 + i, vq + x->n); + dirichlet_char_next(x, G); + } + _dirichlet_dft_index_naive(w1, t1, G, prec); + + gr_ctx_init_complex_acb(ctx, prec); + status |= gr_dft_dirichlet((gr_ptr) w2, (gr_srcptr) vq, G, ctx); + if (status != GR_SUCCESS) + TEST_FUNCTION_FAIL("status, q = %wu\n", q[k]); + gr_ctx_clear(ctx); + + dirichlet_char_one(x, G); + for (i = 0; i < len; i++) + { + if (!acb_overlaps(w1 + i, w2 + x->n)) + TEST_FUNCTION_FAIL("generic number-indexed, q = %wu " + "i = %wu\n", q[k], i); + dirichlet_char_next(x, G); + } + + gr_dft_acb_dirichlet(w2, vq, G, prec); + dirichlet_char_one(x, G); + for (i = 0; i < len; i++) + { + if (!acb_overlaps(w1 + i, w2 + x->n)) + TEST_FUNCTION_FAIL("acb number-indexed, q = %wu " + "i = %wu\n", q[k], i); + dirichlet_char_next(x, G); + } + + dirichlet_char_clear(x); + _acb_vec_clear(t1, len); + } + + _acb_vec_clear(v, len); + _acb_vec_clear(vq, q[k]); + _acb_vec_clear(w1, len); + _acb_vec_clear(w2, FLINT_MAX(q[k], len)); + dirichlet_group_clear(G); + } + + /* q = 1: the transforms reduce to the identity on the single + entry (the Conrey number of the trivial character is 1, so the + generic gather/scatter must not run) */ + { + dirichlet_group_t G; + acb_t a, b; + gr_ctx_t ctx; + int status = GR_SUCCESS; + + dirichlet_group_init(G, 1); + acb_init(a); + acb_init(b); + acb_set_d_d(a, 1.25, -0.5); + + gr_ctx_init_complex_acb(ctx, 64); + status |= gr_dft_dirichlet((gr_ptr) b, (gr_srcptr) a, G, ctx); + if (status != GR_SUCCESS || !acb_equal(a, b)) + TEST_FUNCTION_FAIL("q = 1 generic\n"); + status |= gr_dft_dirichlet_index((gr_ptr) b, (gr_srcptr) a, G, ctx); + if (status != GR_SUCCESS || !acb_equal(a, b)) + TEST_FUNCTION_FAIL("q = 1 generic index\n"); + gr_ctx_clear(ctx); + + acb_zero(b); + gr_dft_acb_dirichlet(b, a, G, 64); + if (!acb_equal(a, b)) + TEST_FUNCTION_FAIL("q = 1 acb\n"); + + acb_zero(b); + gr_dft_acb_dirichlet_index(b, a, G, 64); + if (!acb_equal(a, b)) + TEST_FUNCTION_FAIL("q = 1 index\n"); + + acb_clear(a); + acb_clear(b); + dirichlet_group_clear(G); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/gr_dft/test/t-prod.c b/src/gr_dft/test/t-prod.c new file mode 100644 index 0000000000..c7fa801b94 --- /dev/null +++ b/src/gr_dft/test/t-prod.c @@ -0,0 +1,244 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include "test_helpers.h" +#include "acb.h" +#include "gr.h" +#include "gr_vec.h" +#include "gr_dft.h" +#include "ulong_extras.h" + +/* reference: X[k] = sum_j v[j] prod_a w_a^(j_a k_a) with w_a the + canonical root of order cyc[a], directly from the definition */ +static int +_prod_dft_naive(gr_ptr res, gr_srcptr vec, const ulong * cyc, slong num, + gr_srcptr w, ulong order, gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + slong sz = ctx->sizeof_elem; + ulong n = 1, j, k; + slong a; + gr_ptr roots, t; + + for (a = 0; a < num; a++) + n *= cyc[a]; + + /* roots[a] = canonical root of order cyc[a] */ + roots = gr_heap_init_vec(FLINT_MAX(num, 1), ctx); + for (a = 0; a < num; a++) + { + if (w == NULL) + status |= gr_dft_default_root(GR_ENTRY(roots, a, sz), + cyc[a], ctx); + else + status |= gr_pow_ui(GR_ENTRY(roots, a, sz), w, + order / cyc[a], ctx); + } + + GR_TMP_INIT(t, ctx); + + for (k = 0; k < n; k++) + { + gr_ptr acc = GR_ENTRY(res, k, sz); + status |= gr_zero(acc, ctx); + + for (j = 0; j < n; j++) + { + ulong jj = j, kk = k, inner = n; + + status |= gr_set(t, GR_ENTRY(vec, j, sz), ctx); + for (a = 0; a < num; a++) + { + ulong m = cyc[a], ja, ka; + gr_ptr p; + + inner /= m; + ja = (jj / inner) % m; + ka = (kk / inner) % m; + + GR_TMP_INIT(p, ctx); + status |= gr_pow_ui(p, GR_ENTRY(roots, a, sz), + (ja * ka) % m, ctx); + status |= gr_mul(t, t, p, ctx); + GR_TMP_CLEAR(p, ctx); + } + status |= gr_add(acc, acc, t, ctx); + } + } + + GR_TMP_CLEAR(t, ctx); + gr_heap_clear_vec(roots, FLINT_MAX(num, 1), ctx); + return status; +} + +TEST_FUNCTION_START(gr_dft_prod, state) +{ + slong iter; + + /* deterministic sweep over ring, shape, and threading */ + { + static const ulong shapes[][5] = { + { 2, 0 }, { 4, 0 }, { 2, 3, 0 }, { 3, 5, 0 }, { 2, 2, 2, 2, 0 }, + { 4, 4, 0 }, { 6, 5, 0 }, { 1, 3, 1, 0 }, { 7, 0 }, { 8, 3, 0 }, + }; + slong si, ring, threads; + + for (ring = 0; ring < 2; ring++) + for (si = 0; si < (slong) (sizeof(shapes) / sizeof(shapes[0])); si++) + for (threads = 1; threads <= 3; threads += 2) + { + gr_ctx_t ctx; + ulong cyc[5], n = 1, j; + slong num = 0, sz; + gr_ptr x, y, z; + gr_dft_prod_pre_t P; + int status = GR_SUCCESS; + + while (shapes[si][num] != 0) + { + cyc[num] = shapes[si][num]; + n *= cyc[num]; + num++; + } + + flint_set_num_threads(threads); + + if (ring == 0) + gr_ctx_init_complex_acb(ctx, 128); + else + { + /* prime q = 1 mod lcm: q = 1 mod 840 covers all shapes */ + gr_ctx_init_nmod(ctx, 68041); /* 68041 = 81 * 840 + 1 */ + } + sz = ctx->sizeof_elem; + + x = gr_heap_init_vec(n, ctx); + y = gr_heap_init_vec(n, ctx); + z = gr_heap_init_vec(n, ctx); + for (j = 0; j < n; j++) + status |= gr_set_ui(GR_ENTRY(x, j, sz), + 1 + (j * j * 17 + j) % 97, ctx); + + if (ring == 0) + { + status |= _prod_dft_naive(y, x, cyc, num, NULL, 0, ctx); + status |= gr_dft_prod_precomp_init(P, cyc, num, 0, ctx); + } + else + { + /* w = generator^((q-1)/840): a root of order 840, + divisible by every component length in the sweep */ + gr_ptr w; + GR_TMP_INIT(w, ctx); + status |= gr_set_ui(w, n_primitive_root_prime(68041), ctx); + status |= gr_pow_ui(w, w, (68041 - 1) / 840, ctx); + status |= _prod_dft_naive(y, x, cyc, num, w, 840, ctx); + status |= gr_dft_prod_precomp_init_root(P, w, 840, cyc, + num, 0, ctx); + GR_TMP_CLEAR(w, ctx); + } + gr_dft_prod_precomp_set_serial_block(P, 1); + status |= gr_dft_prod_precomp(z, x, P, ctx); + + if (status != GR_SUCCESS) + TEST_FUNCTION_FAIL("status, ring %wd shape %wd\n", ring, si); + + for (j = 0; j < n; j++) + if (gr_equal(GR_ENTRY(z, j, sz), GR_ENTRY(y, j, sz), ctx) + == T_FALSE) + TEST_FUNCTION_FAIL("forward, ring %wd shape %wd j %wu\n", + ring, si, j); + + /* aliased forward, then inverse roundtrip */ + status |= _gr_vec_set(z, x, n, ctx); + status |= gr_dft_prod_precomp(z, z, P, ctx); + status |= gr_dft_prod_inverse_precomp(z, z, P, ctx); + if (status != GR_SUCCESS) + TEST_FUNCTION_FAIL("roundtrip status, shape %wd\n", si); + for (j = 0; j < n; j++) + if (gr_equal(GR_ENTRY(z, j, sz), GR_ENTRY(x, j, sz), ctx) + == T_FALSE) + TEST_FUNCTION_FAIL("roundtrip, ring %wd shape %wd j %wu\n", + ring, si, j); + + gr_dft_prod_precomp_clear(P); + gr_heap_clear_vec(x, n, ctx); + gr_heap_clear_vec(y, n, ctx); + gr_heap_clear_vec(z, n, ctx); + gr_ctx_clear(ctx); + flint_set_num_threads(1); + } + } + + /* acb wrapper: fixed-point vs ball containment (both enclose the + exact transform, so they must overlap), plus roundtrip */ + for (iter = 0; iter < 30 * flint_test_multiplier(); iter++) + { + ulong cyc[4], n = 1, j; + slong num = 1 + n_randint(state, 3), a; + slong prec = 32 + n_randint(state, 300); + acb_ptr v, w1, w2; + + for (a = 0; a < num; a++) + { + cyc[a] = 1 + n_randint(state, 8); + n *= cyc[a]; + } + + v = _acb_vec_init(n); + w1 = _acb_vec_init(n); + w2 = _acb_vec_init(n); + + for (j = 0; j < n; j++) + { + acb_randtest(v + j, state, prec, 4); + if (n_randint(state, 8) == 0) + acb_mul_2exp_si(v + j, v + j, + (slong) n_randint(state, 200) - 100); + } + + if (_gr_dft_acb_prod(w1, v, cyc, num, 0, 2, prec) == GR_SUCCESS) + { + if (_gr_dft_acb_prod(w2, v, cyc, num, 0, 1, prec) != GR_SUCCESS) + TEST_FUNCTION_FAIL("ball path failed\n"); + for (j = 0; j < n; j++) + if (!acb_overlaps(w1 + j, w2 + j)) + TEST_FUNCTION_FAIL("containment, n %wu j %wu\n", n, j); + } + + /* roundtrip through the public entries */ + gr_dft_acb_prod(w1, v, cyc, num, prec); + gr_dft_acb_prod_inverse(w2, w1, cyc, num, prec); + for (j = 0; j < n; j++) + if (!acb_contains(w2 + j, v + j) && !acb_overlaps(w2 + j, v + j)) + TEST_FUNCTION_FAIL("acb roundtrip, n %wu j %wu\n", n, j); + + /* plan interface must agree with the one-shot */ + { + gr_dft_acb_prod_pre_t Q; + if (gr_dft_acb_prod_precomp_init(Q, cyc, num, prec) + == GR_SUCCESS) + { + gr_dft_acb_prod_precomp(w2, v, Q, prec); + for (j = 0; j < n; j++) + if (!acb_overlaps(w1 + j, w2 + j)) + TEST_FUNCTION_FAIL("plan vs one-shot, j %wu\n", j); + gr_dft_acb_prod_precomp_clear(Q); + } + } + + _acb_vec_clear(v, n); + _acb_vec_clear(w1, n); + _acb_vec_clear(w2, n); + } + + TEST_FUNCTION_END(state); +} From ea153f8d1d8308fb98c90de2f5c4ae114479b096 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Thu, 9 Jul 2026 18:56:22 +0200 Subject: [PATCH 2/2] Fixes --- CMakeLists.txt | 2 +- doc/source/gr_dft.rst | 55 ++++++++++++++++++++++++++----------------- src/gr_dft/prod.c | 16 ++++--------- 3 files changed, 39 insertions(+), 34 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e7a29faf54..30e8b55c10 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 diff --git a/doc/source/gr_dft.rst b/doc/source/gr_dft.rst index f4a81c72d5..1f9b4326be 100644 --- a/doc/source/gr_dft.rst +++ b/doc/source/gr_dft.rst @@ -48,11 +48,12 @@ The *alg* parameter of the plan constructors selects the algorithm: decomposes a transform of length `n = n_1 n_2` with `n_1, n_2 \approx \sqrt{n}` into column transforms of length `n_2`, pointwise twiddle multiplications, row transforms of length `n_1`, - and a final transposition. This is cache-friendlier than plain - Cooley-Tukey for transforms too large to fit in cache. The row and + and a final transposition. The row and column transforms currently use radix-2 Cooley-Tukey sub-plans. - This is the algorithm supporting multithreading (see below), and - ``GR_DFT_ALG_AUTO`` selects it for large power-of-two lengths. + Only available on explicit request: ``GR_DFT_ALG_AUTO`` never + selects it, since in measurements up to `n = 2^{20}` its cache + blocking did not outperform the plain transforms, serially or + threaded (see the multithreading section below). * ``GR_DFT_ALG_SPLIT`` -- recursive split-radix, performing about `\tfrac{1}{3} n \log_2 n` multiplications by nontrivial roots of unity when multiplication by `w^{n/4}` (playing the role of `-i`) @@ -95,10 +96,12 @@ The *alg* parameter of the plan constructors selects the algorithm: transforms plus a pointwise multiplication by the precomputed transformed kernel). All chirp factors are entries of the root table. The power-of-two sub-plan requires a root of unity for the - convolution length, obtained with :func:`gr_dft_default_root`; - over rings where this fails (e.g. finite fields without roots of - the required order), Bluestein is unavailable and plans fall back - to direct prime kernels. + convolution length, obtained with :func:`gr_dft_default_root` (or, + over the fixed-point contexts, built directly by the canonical + fixed-point table constructor); over rings where no such root is + available (e.g. finite fields without roots of the required order), + Bluestein is unavailable and plans fall back to direct prime + kernels. The *flags* parameter is a bitwise combination of: @@ -302,15 +305,18 @@ The values `1`, `-1`, `i`, `-i` are not representable and are rounded to magnitude `1 - \operatorname{ulp}`; the transforms never multiply by these values explicitly. Consequently the contexts are not rings in the strict sense: ``is_one`` and ``is_neg_one`` return ``T_FALSE``, -and integers of absolute value greater than 1 cannot be assigned. In -particular :func:`gr_dft_inverse_precomp`, which scales by `1/n`, -returns ``GR_DOMAIN`` over these contexts; use the unscaled inverse -:func:`_gr_dft_precomp_raw` and incorporate the scaling into the -caller's own normalization. For the same reason the Bluestein -algorithm is not available over fixed-point contexts (its kernel -normalization and convolution intermediates do not fit the -representation); prime and prime-power lengths automatically fall -back to the direct kernels of the mixed-radix algorithm. +and integers of absolute value greater than 1 cannot be assigned. +Division by an unsigned integer is supported (a truncating +``mpn`` division per component, with error below 1 ulp), so the +scaled inverse :func:`gr_dft_inverse_precomp` works over these +contexts; the unscaled :func:`_gr_dft_precomp_raw` remains available +when the caller prefers to fold the `1/n` into its own +normalization, as the ``acb`` drop-in transforms do. The Bluestein +algorithm is likewise supported: the `1/\mathrm{conv\_len}` scaling +of its inverse sub-transform is folded into the precomputed +transformed kernel, together with an extra factor `1/2` keeping the +kernel entries inside the representable range `|t| < 1`, which is +undone with an exact doubling during the pointwise multiplication. The root table of a fixed-point plan is built from a primitive root computed via ``arb`` at elevated internal precision, truncated to @@ -578,8 +584,7 @@ Plan construction are decidable in the ring), returning ``GR_DOMAIN`` on a definite failure. - A sanity check `w^{n/2} = -1` is performed; if this can be decided - and is false, ``GR_DOMAIN`` is returned. If the return status is + If the return status is not ``GR_SUCCESS``, the plan is left in a cleared state and must not be used (calling :func:`gr_dft_precomp_clear` is harmless). @@ -628,8 +633,13 @@ Transforms Sets *res* to the forward respectively inverse DFT of *vec*, both of length `n`, using the plan *P*. Aliasing of *res* and *vec* is - allowed. The inverse transform includes the scaling by `1/n` and - returns ``GR_DOMAIN`` if `n` is not invertible in the ring. + allowed. The inverse transform includes the scaling by `1/n`, + performed as a division by `n` on the raw result (so it also + succeeds in rings where `n` has no inverse whenever the ring can + carry out the division, and over the fixed-point contexts, where + it is a truncating division per component); the division's error + status, ``GR_DOMAIN`` in a ring where it fails, is passed + through. .. function:: int gr_dft(gr_ptr res, gr_srcptr vec, ulong n, gr_ctx_t ctx) int gr_dft_inverse(gr_ptr res, gr_srcptr vec, ulong n, gr_ctx_t ctx) @@ -735,7 +745,8 @@ and the ring is marked thread-safe. int gr_dft_prod_inverse(gr_ptr res, gr_srcptr vec, const ulong * cyc, slong num, gr_ctx_t ctx) Product DFT and its inverse (the latter including the `1/n` - normalization, applied once at the end rather than per axis). + normalization, performed as a single division by `n` at the end + rather than per axis, as in :func:`gr_dft_inverse_precomp`). Aliasing of *res* and *vec* is allowed. .. function:: void gr_dft_prod_precomp_nfixed_bound(double * peak, double * err_ulps, double in_mag, double in_err, const gr_dft_prod_pre_t P) diff --git a/src/gr_dft/prod.c b/src/gr_dft/prod.c index 5acf03bbeb..039393e6a1 100644 --- a/src/gr_dft/prod.c +++ b/src/gr_dft/prod.c @@ -412,25 +412,19 @@ gr_dft_prod_precomp(gr_ptr res, gr_srcptr vec, const gr_dft_prod_pre_t P, } /* the raw inverse omits the 1/n normalization of each component; - apply the full 1/n once */ + scale by the full 1/n once at the end (as in + gr_dft_inverse_precomp, using scalar division so that exact + division succeeds even in rings without an inverse of n) */ int gr_dft_prod_inverse_precomp(gr_ptr res, gr_srcptr vec, const gr_dft_prod_pre_t P, gr_ctx_t ctx) { int status; - gr_ptr c; status = _gr_dft_prod_precomp_raw(res, vec, 1, P, ctx); - if (status == GR_SUCCESS && P->n > 1) - { - GR_TMP_INIT(c, ctx); - status |= gr_set_ui(c, P->n, ctx); - status |= gr_inv(c, c, ctx); - if (status == GR_SUCCESS) - status |= _gr_vec_mul_scalar(res, res, (slong) P->n, c, ctx); - GR_TMP_CLEAR(c, ctx); - } + if (P->n >= 2) + status |= _gr_vec_div_scalar_ui(res, res, (slong) P->n, P->n, ctx); return status; }