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

fix dv in dft_crt

Browse source
This commit is contained in:
Pascal 2016-09-30 21:22:06 +02:00
parent 97fcc27e6a
commit a7a4d316a9
7 changed files with 38 additions and 32 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

12
acb_dirichlet.h
View file

@ -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);
}

39
acb_dirichlet/dft_crt.c
View file

@ -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);
}

5
acb_dirichlet/dft_cyc.c
View file

@ -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;

5
acb_dirichlet/dft_pol.c
View file

@ -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);
}

2
acb_dirichlet/dft_precomp.c
View file

@ -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);
}
}
}

3
acb_dirichlet/dft_step.c
View file

@ -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++)

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

@ -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)