diff --git a/acb_dirichlet.h b/acb_dirichlet.h index d03c59af..28704fab 100644 --- a/acb_dirichlet.h +++ b/acb_dirichlet.h @@ -316,7 +316,7 @@ typedef crt_struct crt_t[1]; void crt_init(crt_t c, ulong n); -void crt_decomp(acb_ptr y, acb_srcptr x, const crt_t c, ulong len); +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_dirichlet_dft_step_struct acb_dirichlet_dft_step_struct; @@ -348,6 +348,7 @@ typedef struct { slong n; crt_t c; + slong dv; /* then a product */ acb_dirichlet_dft_step_ptr cyc; } @@ -402,6 +403,8 @@ acb_dirichlet_dft_step_struct /*typedef acb_dirichlet_dft_pre_struct acb_dirichlet_dft_pre_t[1];*/ +#define DFT_VERB 0 + enum { DFT_POL, DFT_CYC, DFT_PROD, DFT_CRT /*, DFT_2E, DFT_CONV */ @@ -462,18 +465,17 @@ ACB_DIRICHLET_INLINE void acb_dirichlet_dft_pol_clear(acb_dirichlet_dft_pol_t pol) { if (pol->zclear) - { - flint_printf(" ## clearing pol [len=%ld]....", pol->n); _acb_vec_clear(pol->z, pol->n); - flint_printf("done\n"); - } } +void _acb_dirichlet_dft_crt_init(acb_dirichlet_dft_crt_t crt, slong dv, slong len, slong prec); + ACB_DIRICHLET_INLINE void acb_dirichlet_dft_crt_init(acb_dirichlet_dft_crt_t crt, slong len, slong prec) { crt->n = len; crt_init(crt->c, len); + crt->dv = 1; crt->cyc = _acb_dirichlet_dft_steps_prod(crt->c->m, crt->c->num, prec); } diff --git a/acb_dirichlet/dft_crt.c b/acb_dirichlet/dft_crt.c index e210056f..0f8eaea5 100644 --- a/acb_dirichlet/dft_crt.c +++ b/acb_dirichlet/dft_crt.c @@ -41,10 +41,9 @@ crt_print(const crt_t c) slong k; if (c->num == 0) { - flint_printf("trivial group, absurd\n"); + flint_printf("trivial group\n"); abort(); } - flint_printf("crt decomp "); for (k = 0; k < c->num; k++) flint_printf("Z/%wuZ ", c->m[k]); flint_printf("\n"); @@ -100,7 +99,7 @@ reorder_from_crt(acb_srcptr v, crt_t c, ulong len) #endif void -crt_decomp(acb_ptr y, acb_srcptr x, const crt_t c, ulong len) +crt_decomp(acb_ptr y, acb_srcptr x, slong dv, const crt_t c, ulong len) { int j, e[CRT_MAX]; ulong k, l; @@ -111,8 +110,7 @@ crt_decomp(acb_ptr y, acb_srcptr x, const crt_t c, ulong len) l = 0; for(k = 0; k < len; k++) { - /*flint_printf("set y[%wu] = x[%wu]\n", k, l);*/ - acb_set(y + k, x + l); + 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); @@ -145,17 +143,12 @@ crt_recomp(acb_ptr y, acb_srcptr x, const crt_t c, ulong len) } void -acb_dirichlet_dft_crt(acb_ptr w, acb_srcptr v, slong len, slong prec) +_acb_dirichlet_dft_crt_init(acb_dirichlet_dft_crt_t crt, slong dv, slong len, slong prec) { - crt_t c; - acb_ptr t; - t = _acb_vec_init(len); - crt_init(c, len); - crt_print(c); - crt_decomp(w, v, c, len); - acb_dirichlet_dft_prod(t, w, c->m, c->num, prec); - crt_recomp(w, t, c, len); - _acb_vec_clear(t, len); + crt->n = len; + crt_init(crt->c, len); + crt->dv = dv; + crt->cyc = _acb_dirichlet_dft_steps_prod(crt->c->m, crt->c->num, prec); } void @@ -163,9 +156,21 @@ acb_dirichlet_dft_crt_precomp(acb_ptr w, acb_srcptr v, const acb_dirichlet_dft_c { acb_ptr t; t = _acb_vec_init(crt->n); - crt_print(crt->c); - crt_decomp(w, v, crt->c, crt->n); + crt_decomp(w, v, crt->dv, crt->c, crt->n); acb_dirichlet_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_dirichlet_dft_crt(acb_ptr w, acb_srcptr v, slong len, slong prec) +{ + crt_t c; + acb_ptr t; + t = _acb_vec_init(len); + crt_init(c, len); + crt_decomp(w, v, 1, c, len); + acb_dirichlet_dft_prod(t, w, c->m, c->num, prec); + crt_recomp(w, t, c, len); + _acb_vec_clear(t, len); +} diff --git a/acb_dirichlet/dft_cyc.c b/acb_dirichlet/dft_cyc.c index 83ed005c..65601651 100644 --- a/acb_dirichlet/dft_cyc.c +++ b/acb_dirichlet/dft_cyc.c @@ -24,14 +24,17 @@ _acb_dirichlet_dft_cyc_init_z_fac(acb_dirichlet_dft_cyc_t t, n_factor_t fac, slo if (z == NULL) { - flint_printf("dft_cyc: compute roots of order %wu\n", t->n); z = _acb_vec_init(t->n); acb_dirichlet_vec_nth_roots(z, 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; diff --git a/acb_dirichlet/dft_pol.c b/acb_dirichlet/dft_pol.c index 2efefcfb..eebf657c 100644 --- a/acb_dirichlet/dft_pol.c +++ b/acb_dirichlet/dft_pol.c @@ -44,10 +44,10 @@ _acb_dirichlet_dft_pol_init(acb_dirichlet_dft_pol_t pol, slong dv, acb_ptr z, sl pol->n = len; pol->dv = dv; - /* warning: if set here, must be cleared somewhere */ if (z == NULL) { - flint_printf("warning: init z in dft_pol, should be avoided\n"); + if (DFT_VERB) + flint_printf("warning: init z[%ld] in dft_pol, should be avoided\n",len); pol->z = _acb_vec_init(len); acb_dirichlet_vec_nth_roots(pol->z, len, prec); pol->dz = 1; @@ -59,5 +59,4 @@ _acb_dirichlet_dft_pol_init(acb_dirichlet_dft_pol_t pol, slong dv, acb_ptr z, sl pol->dz = dz; pol->zclear = 0; } - flint_printf("init dft_pol len = %ld, dz = %ld, dv = %ld\n", pol->n, pol->dz, pol->dv); } diff --git a/acb_dirichlet/dft_precomp.c b/acb_dirichlet/dft_precomp.c index c1171e12..70e46e34 100644 --- a/acb_dirichlet/dft_precomp.c +++ b/acb_dirichlet/dft_precomp.c @@ -41,7 +41,7 @@ _acb_dirichlet_dft_precomp_init(acb_dirichlet_dft_pre_t pre, slong dv, acb_ptr z else { pre->type = DFT_CRT; - acb_dirichlet_dft_crt_init(pre->t.crt, len, prec); + _acb_dirichlet_dft_crt_init(pre->t.crt, dv, len, prec); } } } diff --git a/acb_dirichlet/dft_step.c b/acb_dirichlet/dft_step.c index 775b2344..290f9115 100644 --- a/acb_dirichlet/dft_step.c +++ b/acb_dirichlet/dft_step.c @@ -54,13 +54,12 @@ acb_dirichlet_dft_step(acb_ptr w, acb_srcptr v, acb_dirichlet_dft_step_ptr cyc, } #if REORDER - /* reorder w to avoid shifting in next DFT */ + /* 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++) diff --git a/acb_dirichlet/test/t-dft.c b/acb_dirichlet/test/t-dft.c index 282e56e4..b258ab5c 100644 --- a/acb_dirichlet/test/t-dft.c +++ b/acb_dirichlet/test/t-dft.c @@ -19,8 +19,7 @@ int main() 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[13] = { 20, 3, 4, 5, 6, 23, 10, 15, 30, 59, 308, 335, 961}; + ulong q[13] = { 2, 3, 4, 5, 6, 23, 10, 15, 30, 59, 308, 335, 961}; flint_rand_t state; slong f, nf = 3; @@ -50,7 +49,6 @@ int main() acb_ptr w = (f == 0) ? w1 : w2; - flint_printf("\ndo dft %s on Z/%wuZ\n", name[f], q[k]); func[f](w, v, q[k], prec); if (f == 0)