diff --git a/acb_dirichlet/dft_convol_fft.c b/acb_dirichlet/dft_convol_fft.c new file mode 100644 index 00000000..23a3b72d --- /dev/null +++ b/acb_dirichlet/dft_convol_fft.c @@ -0,0 +1,58 @@ +/* + Copyright (C) 2016 Pascal Molin + + This file is part of Arb. + + Arb 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 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "acb_dirichlet.h" + +/* assume np >= 2 * n - 1 */ +void +acb_dirichlet_dft_convol_pad(acb_ptr fp, acb_ptr gp, acb_srcptr f, acb_srcptr g, slong n, slong np) +{ + slong k; + for (k = 0; k < n; k++) + acb_set(gp + k, g + 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_dirichlet_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_dirichlet_dft_cyc(w, w, len, prec); + acb_vec_swap_ri(w, len); +} + +void +acb_dirichlet_dft_convol_fft(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec) +{ + int e; + slong k, np; + acb_ptr fp, gp; + acb_dirichlet_dft_rad2_t dft; + e = n_clog(2 * len + 1, 2); + acb_dirichlet_dft_rad2_init(dft, e, prec); + np = dft->n; + fp = _acb_vec_init(np); + gp = _acb_vec_init(np); + acb_dirichlet_dft_convol_pad(fp, gp, f, g, len, np); + acb_dirichlet_dft_rad2_precomp(fp, dft, prec); + acb_dirichlet_dft_rad2_precomp(gp, dft, prec); + _acb_vec_kronecker_mul(gp, gp, fp, np, prec); + acb_dirichlet_dft_inverse_rad2_precomp(gp, dft, prec); + for (k = 0; k < len; k++) + acb_set(w + k, gp + k); +} diff --git a/acb_dirichlet/dft_convol_naive.c b/acb_dirichlet/dft_convol_naive.c new file mode 100644 index 00000000..be7f2c8e --- /dev/null +++ b/acb_dirichlet/dft_convol_naive.c @@ -0,0 +1,31 @@ +/* + Copyright (C) 2016 Pascal Molin + + This file is part of Arb. + + Arb 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 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "acb_dirichlet.h" + +void +acb_dirichlet_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, gy; + wx = w + x; + fx = f + x; + gy = g; + acb_zero(wx); + for (y = 0; y <= x; y++) + acb_addmul(wx, fx--, gy++, prec); + for (; y < len; y++) + acb_addmul(wx, f + x - y, g + y, prec); + } +}