dft: also test random length

This commit is contained in:
Pascal 2017-10-23 09:45:01 +02:00
parent c2fb3c8221
commit 0610898c74
3 changed files with 48 additions and 46 deletions

View file

@ -18,11 +18,7 @@ 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");
return;
}
c = cyc[0];
if (num == 1)
{

View file

@ -50,14 +50,14 @@ check_vec_eq_prec(acb_srcptr w1, acb_srcptr w2, slong len, slong prec, slong dig
}
}
int main()
{
slong k;
slong prec = 100, digits = 30;
slong nq = 15;
ulong q[15] = { 2, 3, 4, 5, 6, 23, 10, 15, 30, 59, 256, 308, 335, 344, 961};
slong nq = 17;
ulong q[17] = { 0, 1, 2, 3, 4, 5, 6, 23, 10, 15, 30, 59, 256, 308, 335, 344, 961};
ulong nr = 5;
flint_rand_t state;
slong f, nf = 3;
@ -69,24 +69,25 @@ int main()
flint_randinit(state);
for (k = 0; k < nq; k++)
for (k = 0; k < nq + nr; k++)
{
slong i;
slong i, len;
acb_ptr z1, z2, x, y;
z1 = _acb_vec_init(q[k]);
z2 = _acb_vec_init(q[k]);
x = _acb_vec_init(q[k]);
y = _acb_vec_init(q[k]);
if (k < nq)
len = q[k];
else
len = n_randint(state, 2000);
for (i = 0; i < q[k]; i++)
z1 = _acb_vec_init(len);
z2 = _acb_vec_init(len);
x = _acb_vec_init(len);
y = _acb_vec_init(len);
for (i = 0; i < len; i++)
{
acb_set_si(x + i, q[k] - i);
acb_set_si(y + i, i * i);
/*
acb_set_si(x + i, n_randint(state, q[k]));
acb_set_si(y + i, n_randint(state, q[k]));
*/
acb_set_si(x + i, n_randint(state, len));
acb_set_si(y + i, n_randint(state, len));
}
for (f = 0; f < nf; f++)
@ -94,19 +95,19 @@ int main()
acb_ptr z = (f == 0) ? z1 : z2;
func[f](z, x, y, q[k], prec);
func[f](z, x, y, len, prec);
if (f == 0)
continue;
check_vec_eq_prec(z1, z2, q[k], prec, digits, q[k], name[0], name[f]);
check_vec_eq_prec(z1, z2, len, prec, digits, len, name[0], name[f]);
}
_acb_vec_clear(x, q[k]);
_acb_vec_clear(y, q[k]);
_acb_vec_clear(z1, q[k]);
_acb_vec_clear(z2, q[k]);
_acb_vec_clear(x, len);
_acb_vec_clear(y, len);
_acb_vec_clear(z1, len);
_acb_vec_clear(z2, len);
}
flint_randclear(state);

View file

@ -56,6 +56,7 @@ int main()
slong prec = 100, digits = 30;
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};
slong nr = 10;
flint_rand_t state;
slong f, nf = 5;
@ -68,17 +69,22 @@ int main()
flint_randinit(state);
/* cyclic dft */
for (k = 0; k < nq; k++)
for (k = 0; k < nq + nr; k++)
{
slong i;
slong i, len;
acb_ptr v, w1, w2, w3;
v = _acb_vec_init(q[k]);
w1 = _acb_vec_init(q[k]);
w2 = _acb_vec_init(q[k]);
w3 = _acb_vec_init(q[k]);
if (k < nq)
len = q[k];
else
len = n_randint(state, 2000);
for (i = 0; i < q[k]; i++)
v = _acb_vec_init(len);
w1 = _acb_vec_init(len);
w2 = _acb_vec_init(len);
w3 = _acb_vec_init(len);
for (i = 0; i < len; i++)
acb_set_si_si(v + i, i, 3 - i);
for (f = 0; f < nf; f++)
@ -87,29 +93,29 @@ int main()
acb_ptr w = (f == 0) ? w1 : w2;
if (DFT_VERB)
flint_printf("\n%s %wu\n", name[f], q[k]);
flint_printf("\n%s %wu\n", name[f], len);
func[f](w, v, q[k], prec);
func[f](w, v, len, prec);
/* check aliasing */
_acb_vec_set(w3, v, q[k]);
func[f](w3, w3, q[k], prec);
_acb_vec_set(w3, v, len);
func[f](w3, w3, len, prec);
check_vec_eq_prec(w1, w3, q[k], prec, digits, q[k], "alias", name[0], name[f]);
check_vec_eq_prec(w1, w3, len, prec, digits, len, "alias", name[0], name[f]);
if (f == 0)
continue;
check_vec_eq_prec(w1, w2, q[k], prec, digits, q[k], "no alias", name[0], name[f]);
check_vec_eq_prec(w1, w2, len, prec, digits, len, "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], "inverse", "original", "inverse");
acb_dft_inverse(w2, w1, len, prec);
check_vec_eq_prec(v, w2, len, prec, digits, len, "inverse", "original", "inverse");
_acb_vec_clear(v, q[k]);
_acb_vec_clear(w1, q[k]);
_acb_vec_clear(w2, q[k]);
_acb_vec_clear(w3, q[k]);
_acb_vec_clear(v, len);
_acb_vec_clear(w1, len);
_acb_vec_clear(w2, len);
_acb_vec_clear(w3, len);
}
/* radix2 dft */
@ -138,4 +144,3 @@ int main()
flint_printf("PASS\n");
return EXIT_SUCCESS;
}