gr_dft.h

Sets w to the default n-th root of unity: 1 for n \le 1, -1 for n = 2, and e^{-2\pi i/n} otherwise. Returns GR_DOMAIN or GR_UNABLE if the ring cannot provide e^{-2\pi i/n}, for example over \mathbb{Z}/m\mathbb{Z}, finite fields, real numbers and acf.

Initialise P with a plan for transforms of length n \ge 1 over ctx with the algorithm alg and the flags flags. The first function uses the principal n-th root of unity w, the other two the default root e^{-2\pi i/n}. The third also enables complex mode when real_ctx is not NULL: ctx must then be a complex ring whose elements are stored as a real and an imaginary part from real_ctx, such as acb over arb. On failure the plan is left cleared.

Frees the memory used by P, including its sub-plans. It is harmless to call it after a failed initialisation.

Writes to perm (an array of n entries) the permutation by which entry i of the forward output holds X_{\mathit{perm}[i]}; the inverse transform expects X_{\mathit{perm}[i]} in entry i of its input. This is the identity unless the plan has GR_DFT_SCRAMBLED set.

Print a description of the plan to out or to standard output: the algorithm and decomposition of each node, with one indented block per sub-plan, the tables present, and the threading settings if they are not the defaults.

Set res to the forward or inverse transform of vec with the plan P. Both vectors have length n and may be the same. ctx must be the ring of the plan. The inverse transform includes the division by n, and returns GR_DOMAIN if the ring cannot perform it.

Set res to the forward or inverse transform of vec with respect to the default root, using a temporary plan from gr_dft_precomp_init with GR_DFT_ALG_AUTO. Return GR_DOMAIN or GR_UNABLE if the ring has no default root of order n.

Attaches num_threads worker threads to P and its sub-plans; every later transform uses them instead of requesting threads from the pool. The handles are borrowed: the caller obtains them, for example with flint_request_threads, and gives them back with flint_give_back_threads after the plan is cleared. Passing NULL and 0 detaches them.

Sets the threading granularity of P and its sub-plans: each pass or phase is divided into at most n / \mathit{serial\_block} chunks, so a transform with n < 2 \cdot \mathit{serial\_block} runs serially. Zero restores the default GR_DFT_SERIAL_BLOCK_DEFAULT.

Initialise P for the product DFT with the num component lengths cyc. The first function uses the default roots w_a = e^{-2\pi i/c_a}. The second takes a root w of the given order, which every component length must divide, and uses w_a = w^{\mathit{order}/c_a}. The flags are passed to the component plans. After any return value, P must be cleared with gr_dft_prod_precomp_clear.

Frees the memory used by P.

Set res to the forward or inverse product DFT of vec, both of length n = c_0 \cdots c_{r-1}, which may be the same. The last two use a temporary plan with default roots. The inverse divides by n once at the end.

Sets the threading granularity of the product transform: the lines of one axis are distributed over at most n / \mathit{serial\_block} chunks. It does not change the component plans. Zero restores the default GR_DFT_SERIAL_BLOCK_DEFAULT.

Set w to the DFT of v over the characters of G. In gr_dft_dirichlet_index the vectors have G->phi_q entries in the lexicographic order of Conrey logarithms (that of dirichlet_char_next). In gr_dft_dirichlet they have G->q entries indexed by Conrey numbers; entries at indices not coprime to q are ignored in v and not written in w. The vectors may be the same. Both use default roots, so they need a ring with e^{-2\pi i/\varphi_k} such as acb, and return GR_UNABLE or GR_DOMAIN otherwise.

Set w to the forward or inverse transform (including the division by n) of v, both of length n, with working precision prec. The vectors may be the same, and nothing is done for n \le 0. _gr_dft_acb forces the method with which (0 automatic, 1 ball arithmetic, 2 fixed point) and returns GR_UNABLE if the forced fixed-point path cannot handle the input, in which case w is set to indeterminate values.

Precomputed versions of gr_dft_acb and gr_dft_acb_inverse for repeated transforms of length n \ge 1. The precision given to gr_dft_acb_precomp_init determines the accuracy; the prec of the transforms only sets the rounding of the outputs. The initialisation returns GR_DOMAIN for n \le 0; Q may be cleared after any return status.

The product DFT of Products of cyclic groups for acb vectors, computed like gr_dft_acb. _gr_dft_acb_prod returns GR_DOMAIN if a component length is zero or the total length is too large.

Precomputed versions of gr_dft_acb_prod and gr_dft_acb_prod_inverse. The initialisation returns GR_DOMAIN if a component length is zero or the total length exceeds about WORD_MAX/8; Q may be cleared after any return status.

The transforms of gr_dft_dirichlet_index and gr_dft_dirichlet for acb vectors, computed with gr_dft_acb_prod. They implement acb_dirichlet_dft_index and acb_dirichlet_dft (Discrete Fourier transforms on characters).

Initialise ctx as the real or complex fixed-point context with nlimbs fraction limbs. Return GR_UNABLE unless 1 \le \mathit{nlimbs} \le GR_DFT_NFIXED_MAX_NLIMBS. The contexts are freed with gr_ctx_clear.

Convert between arf values and elements of a real fixed-point context. Setting truncates toward zero, with error below 1 ulp, and saturates values with |x| \ge 1 to \pm(1 - \mathrm{ulp}); it returns GR_DOMAIN only for non-finite x. Getting is exact. The parts of a complex element are converted by treating them as two consecutive real elements.

For a forward transform or an inverse transform without the division by n, executed in fixed-point arithmetic with the plan P, given inputs of complex modulus at most in_mag with errors of at most in_err_ulps ulp in each part: set peak to a bound for the modulus of every input, intermediate and output value, and err_ulps to a bound for the error of each part of the output, in ulp. There is no overflow if peak is less than 1.

Initialise P as a layout for length n with the given algorithm and flags (complex_mode tells the automatic choice whether the plan will be realised in complex mode), and realise a layout over ctx with default roots, in complex mode if real_ctx is not NULL. A layout can be passed to gr_dft_precomp_nfixed_bound, gr_dft_precomp_print and gr_dft_precomp_clear. On failure, both functions leave P cleared.

The same for product plans. P must be cleared with gr_dft_prod_precomp_clear after any return value of either function.

Apply the plan, forward or inverse, without the division by n of the inverse transform. res and vec may be the same.

The kernels of the algorithms, each for a plan of the matching algorithm, without the division by n. The radix-2 and four-step kernels work in place on x (the radix-2 kernels on the entries x[i * stride]); the others write to res, which must not alias vec. The _threaded versions use threads and may only be called on a top-level plan; the others never use threads.

The first sets res to w^e x (or w^{-e} x if inverse is set), for 0 \le e < n and a plan that stores the full root table. Radix-2 plans, and split-radix plans in complex mode, of length at least 4 and without Karatsuba tables store per-stage tables instead. The second sets res to w x for a plan constant w, using the complex Karatsuba entry wtab3 (three real elements c, d - c, d + c for w = c + di) if it is not NULL. res must not alias x. In complex mode, rtmp is an initialised element of the real ring used as scratch space.

Permutes the entries x[i * stride], 0 \le i < 2^{\mathit{depth}}, by bit reversal of the index, in place.

Return 1 if ctx is a real or complex fixed-point context, and 0 otherwise.

The root table of a fixed-point plan and its error bounds. The first function sets roots to approximations of e^{-2\pi i j/n}, 0 \le j < n, over the complex fixed-point context ctx and writes a bound for their errors, in ulp, to err_ulps. The second bounds the error of the entries of a table of length n built with products whose rounding error is at most kmul ulp in modulus. The third returns the error bound per part of one complex multiplication in ctx: 2, 4 or 6 ulp.

Multiply two elements of a complex fixed-point context with four or with three real multiplications, regardless of the limb count. They exist for tuning the crossover between the two, 22 limbs, with the profile program p-gr_dft_nfixed.