From c2fb3c8221d0b52af726c7b661d5645764125c26 Mon Sep 17 00:00:00 2001 From: Pascal Date: Sun, 22 Oct 2017 22:44:24 +0200 Subject: [PATCH] dft: fix len=0,1 + test --- acb_dft.h | 7 ++--- acb_dft/crt.c | 61 +++++++++++++++++++++++++++--------------- acb_dft/cyc.c | 17 +++++++++++- acb_dft/precomp.c | 6 ++--- acb_dft/prod.c | 4 +-- acb_dft/step.c | 3 ++- acb_dft/test/t-dft.c | 10 +++---- doc/source/acb_dft.rst | 2 +- 8 files changed, 69 insertions(+), 41 deletions(-) diff --git a/acb_dft.h b/acb_dft.h index 76957127..cc681778 100644 --- a/acb_dft.h +++ b/acb_dft.h @@ -179,11 +179,8 @@ void acb_dft_step(acb_ptr w, acb_srcptr v, acb_dft_step_ptr cyc, slong num, slon 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); -ACB_DFT_INLINE void -acb_dft_cyc_precomp(acb_ptr w, acb_srcptr v, const acb_dft_cyc_t cyc, slong prec) -{ - acb_dft_step(w, v, cyc->cyc, cyc->num, 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); diff --git a/acb_dft/crt.c b/acb_dft/crt.c index 98c4dae1..07283fe3 100644 --- a/acb_dft/crt.c +++ b/acb_dft/crt.c @@ -18,7 +18,10 @@ crt_init(crt_t c, ulong n) n_factor_t fac; n_factor_init(&fac); - n_factor(&fac, n, 1); + if (n) + n_factor(&fac, n, 1); + else + fac.num = 0; nmod_init(&c->n, n); @@ -172,33 +175,49 @@ acb_dft_crt_clear(acb_dft_crt_t crt) void acb_dft_crt_precomp(acb_ptr w, acb_srcptr v, const acb_dft_crt_t crt, slong prec) { - acb_ptr t; - t = _acb_vec_init(crt->n); - if (w == v) + if (crt->n <= 1) { - _acb_vec_set(t, v, crt->n); - v = t; + 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); } - 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) { - crt_t c; - acb_ptr t; - t = _acb_vec_init(len); - if (w == v) + if (len <= 1) { - _acb_vec_set(t, v, len); - v = t; + 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); } - 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/acb_dft/cyc.c b/acb_dft/cyc.c index 2b55fe9c..2e922fa3 100644 --- a/acb_dft/cyc.c +++ b/acb_dft/cyc.c @@ -70,7 +70,10 @@ _acb_dft_cyc_init(acb_dft_cyc_t t, slong dv, slong len, slong prec) { n_factor_t fac; n_factor_init(&fac); - n_factor(&fac, len, 0); + if (len) + n_factor(&fac, len, 1); + else + fac.num = 0; _acb_dft_cyc_init_z_fac(t, fac, dv, NULL, 0, len, prec); } @@ -85,6 +88,18 @@ acb_dft_cyc_clear(acb_dft_cyc_t t) 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) { diff --git a/acb_dft/precomp.c b/acb_dft/precomp.c index f72d9f8a..300d9a39 100644 --- a/acb_dft/precomp.c +++ b/acb_dft/precomp.c @@ -17,10 +17,10 @@ _acb_dft_precomp_init(acb_dft_pre_t pre, slong dv, acb_ptr z, slong dz, slong le pre->n = len; if (len <= 1) { - flint_printf("precomp init: trivial group. abort.\n"); - abort(); + pre->type = DFT_NAIVE; + _acb_dft_naive_init(pre->t.naive, dv, z, dz, len, prec); } - if (n_is_prime(len)) + else if (n_is_prime(len)) { if (len < 100) { diff --git a/acb_dft/prod.c b/acb_dft/prod.c index 0a4f2356..8941c5fc 100644 --- a/acb_dft/prod.c +++ b/acb_dft/prod.c @@ -51,9 +51,7 @@ acb_dft_prod_clear(acb_dft_prod_t t) void acb_dft_prod_precomp(acb_ptr w, acb_srcptr v, const acb_dft_prod_t prod, slong prec) { - if (prod->num == 0) - acb_set(w + 0, v + 0); - else + if (prod->num >= 1) acb_dft_step(w, v, prod->cyc, prod->num, prec); } diff --git a/acb_dft/step.c b/acb_dft/step.c index 3da6768b..f66dcfc7 100644 --- a/acb_dft/step.c +++ b/acb_dft/step.c @@ -19,8 +19,9 @@ acb_dft_step(acb_ptr w, acb_srcptr v, acb_dft_step_ptr cyc, slong num, slong pre acb_dft_step_struct c; if (num == 0) { + /* only possible if len = 0 */ flint_printf("error: reached num = 0 in acb_dft_step\n"); - abort(); /* or just copy v to w */ + return; } c = cyc[0]; if (num == 1) diff --git a/acb_dft/test/t-dft.c b/acb_dft/test/t-dft.c index d32758c9..506d80e4 100644 --- a/acb_dft/test/t-dft.c +++ b/acb_dft/test/t-dft.c @@ -54,8 +54,8 @@ int main() { slong k; slong prec = 100, digits = 30; - slong nq = 17; - ulong q[17] = { 2, 3, 4, 5, 6, 23, 10, 15, 16, 30, 59, 125, 308, 335, 525, 961, 1225}; + 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}; flint_rand_t state; slong f, nf = 5; @@ -103,10 +103,8 @@ int main() check_vec_eq_prec(w1, w2, q[k], prec, digits, q[k], "no alias", name[0], name[f]); } - /* acb_dft_inverse(w2, w1, q[k], prec); - check_vec_eq_prec(v, w2, q[k], prec, digits, q[k], "original", "inverse"); - */ + check_vec_eq_prec(v, w2, q[k], prec, digits, q[k], "inverse", "original", "inverse"); _acb_vec_clear(v, q[k]); _acb_vec_clear(w1, q[k]); @@ -115,7 +113,7 @@ int main() } /* radix2 dft */ - for (k = 1; k < 12; k++) + for (k = 0; k < 12; k++) { slong n = 1 << k, j; acb_ptr v, w1, w2; diff --git a/doc/source/acb_dft.rst b/doc/source/acb_dft.rst index c55ec97e..02e195b1 100644 --- a/doc/source/acb_dft.rst +++ b/doc/source/acb_dft.rst @@ -225,7 +225,7 @@ CRT decomposition Sets *w* to the DFT of *v* of size *t->n*, using the CRT decomposition scheme *t*. -Cooley-Tuckey decomposition +Cooley-Tukey decomposition ............................................................................... .. function:: void acb_dft_cyc(acb_ptr w, acb_srcptr v, slong n, slong prec)