hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-06 09:51:39 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

fix dft on group

Browse source
This commit is contained in:
Pascal 2016-08-04 16:32:31 +02:00
parent bf8d58cbf2
commit 9441c452ec
2 changed files with 87 additions and 4 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

16
acb_dirichlet/dft_fast.c
View file

@ -202,7 +202,10 @@ _acb_dft_prod(acb_ptr w, acb_srcptr v, dft_prod_step * cyc, slong num, slong pre
{ {
dft_prod_step c; dft_prod_step c;
if (num == 0) if (num == 0)
{
flint_printf("error: reached num = 0 in dft_prod\n");
abort(); /* or just copy v to w */ abort(); /* or just copy v to w */
}
c = cyc[0]; c = cyc[0];
if (num == 1) if (num == 1)
{ {
@ -221,7 +224,7 @@ _acb_dft_prod(acb_ptr w, acb_srcptr v, dft_prod_step * cyc, slong num, slong pre
{ {
_acb_dft_prod(wi, vi, cyc + 1, num - 1, prec); _acb_dft_prod(wi, vi, cyc + 1, num - 1, prec);
wi += M; wi += M;
vi += dv; vi += dv; /* here = M */
} }
/* after first pass */ /* after first pass */
for (j = 0; j < M; j++) for (j = 0; j < M; j++)
@ -237,7 +240,14 @@ _acb_dft_prod(acb_ptr w, acb_srcptr v, dft_prod_step * cyc, slong num, slong pre
void void
acb_dirichlet_dft_prod_precomp(acb_ptr w, acb_srcptr v, acb_dft_prod_t t, slong prec) acb_dirichlet_dft_prod_precomp(acb_ptr w, acb_srcptr v, acb_dft_prod_t t, slong prec)
{ {
_acb_dft_prod(w, v, t->cyc, t->num, prec); if (t->num == 0)
{
acb_set(w + 0, v + 0);
}
else
{
_acb_dft_prod(w, v, t->cyc, t->num, prec);
}
} }
void void
@ -258,7 +268,7 @@ acb_dirichlet_dft_prod_init(acb_dft_prod_t t, slong * cyc, slong num, slong prec
len /= m; len /= m;
t->cyc[i].m = m; t->cyc[i].m = m;
t->cyc[i].M = len; t->cyc[i].M = len;
t->cyc[i].dv = dv; t->cyc[i].dv = len;
_acb_dirichlet_dft_cyc_init(t->cyc[i].pre, m, len, prec); _acb_dirichlet_dft_cyc_init(t->cyc[i].pre, m, len, prec);
dv *= m; dv *= m;
} }

75
acb_dirichlet/test/t-dft.c
View file

@ -36,6 +36,7 @@ int main()
flint_printf("dft...."); flint_printf("dft....");
fflush(stdout); fflush(stdout);
/* cyclic dft */
for (k = 0; k < nq; k++) for (k = 0; k < nq; k++)
{ {
slong i; slong i;
@ -64,11 +65,83 @@ int main()
abort(); abort();
} }
} }
_acb_vec_clear(v, q[k]); _acb_vec_clear(v, q[k]);
_acb_vec_clear(w1, q[k]); _acb_vec_clear(w1, q[k]);
_acb_vec_clear(w2, q[k]); _acb_vec_clear(w2, q[k]);
} }
/* Dirichlet group DFT */
for (k = 0; k < nq; k++)
{
slong i, j, len;
acb_dirichlet_group_t G;
acb_dirichlet_conrey_t x, y;
acb_t chiy;
acb_ptr v, w1, w2;
acb_dirichlet_group_init(G, q[k]);
len = G->phi_q;
v = _acb_vec_init(len);
w1 = _acb_vec_init(len);
w2 = _acb_vec_init(len);
acb_dirichlet_conrey_init(x, G);
acb_dirichlet_conrey_one(x, G);
for (i = 0; i < len; i++)
{
acb_set_si(v + i, x->n);
acb_dirichlet_conrey_next(x, G);
}
/* naive */
acb_init(chiy);
acb_dirichlet_conrey_init(y, G);
acb_dirichlet_conrey_one(x, G);
for (i = 0; i < len; i++)
{
acb_zero(w1 + i);
acb_dirichlet_conrey_one(y, G);
for (j = 0; j < len; j++)
{
acb_dirichlet_pairing_conrey(chiy, G, x, y, prec);
acb_addmul(w1 + i, chiy, v + j, prec);
acb_dirichlet_conrey_next(y, G);
}
acb_dirichlet_conrey_next(x, G);
}
acb_clear(chiy);
acb_dirichlet_conrey_clear(y);
acb_dirichlet_conrey_clear(x);
/* dft */
acb_dirichlet_dft_conrey(w2, v, G, prec);
for (i = 0; i < len; i++)
{
if (!acb_overlaps(w1 + i, w2 + i))
{
flint_printf("FAIL\n\n");
flint_printf("q = %wu\n", q[k]);
flint_printf("v [size %wu]\n", len);
acb_vec_printd(v, len, 10);
flint_printf("\nDFT differ from index %ld / %ld \n", i, len);
flint_printf("\nnaive =\n");
acb_vec_printd(w1, len, 10);
flint_printf("\nfast =\n");
acb_vec_printd(w2, len, 10);
flint_printf("\n\n");
abort();
}
}
_acb_vec_clear(v, len);
_acb_vec_clear(w1, len);
_acb_vec_clear(w2, len);
acb_dirichlet_group_clear(G);
}
flint_printf("PASS\n"); flint_printf("PASS\n");
} }