mirror of
https://github.com/vale981/arb
synced 2025-03-04 17:01:40 -05:00
temporary remove acb_dft and hide theta internals
This commit is contained in:
Pascal
parent
0428ebf806
commit
30525ec26a
23 changed files with 1 additions and 1988 deletions
|
|
@ -13,7 +13,7 @@ QUIET_AR = @echo ' ' AR ' ' $@;
|
|||
AT=@
|
||||
|
||||
BUILD_DIRS = fmpr arf mag arb arb_mat arb_poly arb_calc acb acb_mat acb_poly \
|
||||
acb_calc acb_hypgeom acb_modular dirichlet acb_dft acb_dirichlet arb_hypgeom bernoulli hypgeom \
|
||||
acb_calc acb_hypgeom acb_modular dirichlet acb_dirichlet arb_hypgeom bernoulli hypgeom \
|
||||
fmpz_extras bool_mat partitions dlog \
|
||||
$(EXTRA_BUILD_DIRS)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,67 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dft.h"
|
||||
|
||||
void
|
||||
acb_dft_bluestein_init(acb_dft_bluestein_t t, slong n, slong prec)
|
||||
{
|
||||
|
||||
nmod_t n2;
|
||||
slong k, k2;
|
||||
acb_ptr z2n;
|
||||
int e = n_clog(2 * n - 1, 2);
|
||||
acb_dft_rad2_init(t->rad2, e, prec);
|
||||
z2n = _acb_vec_init(2 * n);
|
||||
_acb_vec_nth_roots(z2n, 2 * n, prec);
|
||||
nmod_init(&n2, 2 * n);
|
||||
t->n = n;
|
||||
t->z = _acb_vec_init(n);
|
||||
for (k = 0, k2 = 0; k < n; k++)
|
||||
{
|
||||
acb_conj(t->z + k, z2n + k2);
|
||||
k2 = nmod_add(k2, 2 * k + 1, n2);
|
||||
}
|
||||
_acb_vec_clear(z2n, 2 * n);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_bluestein_precomp(acb_ptr w, acb_srcptr v, const acb_dft_bluestein_t t, slong prec)
|
||||
{
|
||||
slong n = t->n;
|
||||
acb_ptr vz, wz, z;
|
||||
z = t->z;
|
||||
/* TODO: allocate directly length 2^e and pad */
|
||||
flint_printf("\n\n====================\n\nv\n");
|
||||
acb_vec_printd_index(v, n, 10);
|
||||
vz = _acb_vec_init(n);
|
||||
acb_vec_kronecker_mul_conj(vz, z, v, n, prec);
|
||||
flint_printf("\nvz\n");
|
||||
acb_vec_printd_index(vz, n, 10);
|
||||
wz = _acb_vec_init(n);
|
||||
acb_dft_convol_rad2_precomp(wz, vz, z, n, t->rad2, prec);
|
||||
flint_printf("\nwz\n");
|
||||
acb_vec_printd_index(wz, n, 10);
|
||||
acb_vec_kronecker_mul_conj(w, z, wz, n, prec);
|
||||
flint_printf("\nw\n");
|
||||
acb_vec_printd_index(w, n, 10);
|
||||
_acb_vec_clear(wz, n);
|
||||
_acb_vec_clear(vz, n);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_bluestein(acb_ptr w, acb_srcptr v, slong len, slong prec)
|
||||
{
|
||||
acb_dft_bluestein_t t;
|
||||
acb_dft_bluestein_init(t, len, prec);
|
||||
acb_dft_bluestein_precomp(w, v, t, prec);
|
||||
acb_dft_bluestein_clear(t);
|
||||
}
|
||||
|
|
@ -1,105 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dft.h"
|
||||
|
||||
/* assume np >= 2 * n - 1 */
|
||||
void
|
||||
acb_dft_convol_pad(acb_ptr fp, acb_ptr gp, acb_srcptr f, acb_srcptr g, slong n, slong np)
|
||||
{
|
||||
slong k;
|
||||
|
||||
if (np < 2 * n - 1)
|
||||
{
|
||||
flint_printf("dft_convol_pad: overlapping padding %ld < 2*%ld-1\n", np, n);
|
||||
abort();
|
||||
}
|
||||
|
||||
for (k = 0; k < n; k++)
|
||||
acb_set(gp + k, g + k);
|
||||
for (; k < np; k++)
|
||||
acb_zero(gp + k);
|
||||
|
||||
for (k = 0; k < n; k++)
|
||||
acb_set(fp + k, f + k);
|
||||
for (k = 1; k < n; k++)
|
||||
acb_set(fp + np - k, f + n - k);
|
||||
for (k = n; k <= np - n; k++)
|
||||
acb_zero(fp + k);
|
||||
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_inverse_cyc(acb_ptr w, acb_srcptr v, slong len, slong prec)
|
||||
{
|
||||
/* divide before to keep v const */
|
||||
_acb_vec_scalar_div_ui(w, v, len, len, prec);
|
||||
acb_vec_swap_ri(w, len);
|
||||
acb_dft_cyc(w, w, len, prec);
|
||||
acb_vec_swap_ri(w, len);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_convol_rad2_precomp(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, const acb_dft_rad2_t rad2, slong prec)
|
||||
{
|
||||
slong np;
|
||||
acb_ptr fp, gp;
|
||||
np = rad2->n;
|
||||
|
||||
flint_printf("\nf\n");
|
||||
acb_vec_printd_index(f, len, 10);
|
||||
flint_printf("\ng\n");
|
||||
acb_vec_printd_index(g, len, 10);
|
||||
|
||||
fp = _acb_vec_init(np);
|
||||
gp = _acb_vec_init(np);
|
||||
acb_dft_convol_pad(fp, gp, f, g, len, np);
|
||||
|
||||
flint_printf("\nF\n");
|
||||
acb_vec_printd_index(fp, np, 10);
|
||||
flint_printf("\nG\n");
|
||||
acb_vec_printd_index(gp, np, 10);
|
||||
|
||||
acb_dft_rad2_precomp(fp, rad2, prec);
|
||||
|
||||
flint_printf("\nDFT F\n");
|
||||
acb_vec_printd_index(fp, np, 10);
|
||||
|
||||
acb_dft_rad2_precomp(gp, rad2, prec);
|
||||
|
||||
flint_printf("\nDFT G\n");
|
||||
acb_vec_printd_index(gp, np, 10);
|
||||
|
||||
_acb_vec_kronecker_mul(gp, gp, fp, np, prec);
|
||||
|
||||
flint_printf("\n(DFT F)(DFT G)=DFT(F*G)\n");
|
||||
acb_vec_printd_index(gp, np, 10);
|
||||
|
||||
acb_dft_inverse_rad2_precomp(gp, rad2, prec);
|
||||
|
||||
flint_printf("\nF*G\n");
|
||||
acb_vec_printd_index(gp, np, 10);
|
||||
|
||||
_acb_vec_set(w, gp, len);
|
||||
_acb_vec_clear(fp, np);
|
||||
_acb_vec_clear(gp, np);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_convol_rad2(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec)
|
||||
{
|
||||
int e;
|
||||
acb_dft_rad2_t dft;
|
||||
e = n_clog(2 * len - 1, 2);
|
||||
acb_dft_rad2_init(dft, e, prec);
|
||||
acb_dft_convol_rad2_precomp(w, f, g, len, dft, prec);
|
||||
acb_dft_rad2_clear(dft);
|
||||
}
|
||||
|
|
@ -1,31 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dft.h"
|
||||
|
||||
void
|
||||
acb_dft_convol_naive(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec)
|
||||
{
|
||||
slong x, y;
|
||||
for (x = 0; x < len; x ++)
|
||||
{
|
||||
acb_ptr wx;
|
||||
acb_srcptr fx, gy;
|
||||
wx = w + x;
|
||||
fx = f + x;
|
||||
gy = g;
|
||||
acb_zero(wx);
|
||||
for (y = 0; y <= x; y++)
|
||||
acb_addmul(wx, fx - y, g + y, prec);
|
||||
for (; y < len; y++)
|
||||
acb_addmul(wx, fx + (len - y), g + y, prec);
|
||||
}
|
||||
}
|
||||
|
|
@ -1,176 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dft.h"
|
||||
|
||||
void
|
||||
crt_init(crt_t c, ulong n)
|
||||
{
|
||||
int k;
|
||||
n_factor_t fac;
|
||||
|
||||
n_factor_init(&fac);
|
||||
n_factor(&fac, n, 1);
|
||||
|
||||
nmod_init(&c->n, n);
|
||||
|
||||
c->num = fac.num;
|
||||
for (k = 0; k < fac.num; k++)
|
||||
{
|
||||
c->m[k] = n_pow(fac.p[k], fac.exp[k]);
|
||||
c->M[k] = n / c->m[k];
|
||||
c->vM[k] = c->M[k] * n_invmod(c->M[k] % c->m[k], c->m[k]);
|
||||
/*
|
||||
flint_printf("m[%ld]=%ld, M[%ld]=%wu, vM[%ld]=%wu\n",
|
||||
k, c->m[k], k, c->M[k], k, c->vM[k]);
|
||||
*/
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
crt_print(const crt_t c)
|
||||
{
|
||||
slong k;
|
||||
if (c->num == 0)
|
||||
{
|
||||
flint_printf("trivial group\n");
|
||||
abort();
|
||||
}
|
||||
for (k = 0; k < c->num; k++)
|
||||
flint_printf("Z/%wuZ ", c->m[k]);
|
||||
flint_printf("\n");
|
||||
}
|
||||
|
||||
#if 0
|
||||
/* lexicographic index of crt elt j */
|
||||
static ulong
|
||||
index_crt(crt_t c, ulong j)
|
||||
{
|
||||
int k;
|
||||
ulong res = 0;
|
||||
for (k = 0; k < c->num; k ++)
|
||||
res = res * c->m[k] + (j % c->m[k]);
|
||||
return res;
|
||||
}
|
||||
|
||||
/* crt elt of lexicographic index i */
|
||||
static ulong
|
||||
crt_index(crt_t c, ulong i)
|
||||
{
|
||||
int k;
|
||||
ulong j, res = 0;
|
||||
for (k = 0; k < c->num; k ++)
|
||||
{
|
||||
j = i % c->m[k];
|
||||
i = i / c->m[k];
|
||||
res = nmod_add(res, j * c->vM[k], c->n);
|
||||
}
|
||||
return res;
|
||||
}
|
||||
/* for all elements can do fast conrey-like loop just adding vM[k] */
|
||||
static acb_ptr
|
||||
reorder_to_crt(acb_srcptr v, crt_t c, ulong len)
|
||||
{
|
||||
ulong k;
|
||||
acb_ptr res;
|
||||
res = flint_malloc(len * sizeof(acb_struct));
|
||||
for (k = 0; k < len; k++)
|
||||
res[index_crt(c, k)] = v[k];
|
||||
return res;
|
||||
}
|
||||
static acb_ptr
|
||||
reorder_from_crt(acb_srcptr v, crt_t c, ulong len)
|
||||
{
|
||||
ulong k;
|
||||
acb_ptr res;
|
||||
res = flint_malloc(len * sizeof(acb_struct));
|
||||
for (k = 0; k < len; k++)
|
||||
res[k] = v[index_crt(c, k)];
|
||||
return res;
|
||||
}
|
||||
#endif
|
||||
|
||||
void
|
||||
crt_decomp(acb_ptr y, acb_srcptr x, slong dv, const crt_t c, ulong len)
|
||||
{
|
||||
int j, e[CRT_MAX];
|
||||
ulong k, l;
|
||||
|
||||
for (j = 0; j < c->num; j++)
|
||||
e[j] = 0;
|
||||
|
||||
l = 0;
|
||||
for(k = 0; k < len; k++)
|
||||
{
|
||||
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);
|
||||
if (e[j] < c->m[j])
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
crt_recomp(acb_ptr y, acb_srcptr x, const crt_t c, ulong len)
|
||||
{
|
||||
int j, e[CRT_MAX];
|
||||
ulong k, l;
|
||||
|
||||
for (j = 0; j < c->num; j++)
|
||||
e[j] = 0;
|
||||
|
||||
l = 0;
|
||||
for(k = 0; k < len; k++)
|
||||
{
|
||||
acb_set(y + l, x + k);
|
||||
for (j = c->num - 1; j >= 0; e[j] = 0, j--)
|
||||
{
|
||||
e[j]++; l = nmod_add(l, c->M[j], c->n);
|
||||
if (e[j] < c->m[j])
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
_acb_dft_crt_init(acb_dft_crt_t crt, slong dv, slong len, slong prec)
|
||||
{
|
||||
crt->n = len;
|
||||
crt_init(crt->c, len);
|
||||
crt->dv = dv;
|
||||
crt->cyc = _acb_dft_steps_prod(crt->c->m, crt->c->num, prec);
|
||||
}
|
||||
|
||||
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);
|
||||
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);
|
||||
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);
|
||||
}
|
||||
|
|
@ -1,93 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dft.h"
|
||||
|
||||
void
|
||||
_acb_dft_cyc_init_z_fac(acb_dft_cyc_t t, n_factor_t fac, slong dv, acb_ptr z, slong dz, slong len, slong prec)
|
||||
{
|
||||
slong i, j, num;
|
||||
t->n = len;
|
||||
num = 0;
|
||||
for (i = 0; i < fac.num; i++)
|
||||
num += fac.exp[i];
|
||||
t->num = num;
|
||||
t->cyc = flint_malloc(num * sizeof(acb_dft_step_struct));
|
||||
|
||||
if (z == NULL)
|
||||
{
|
||||
z = _acb_vec_init(t->n);
|
||||
_acb_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;
|
||||
|
||||
num = 0;
|
||||
for (i = 0; i < fac.num; i++)
|
||||
{
|
||||
for (j = 0; j < fac.exp[i]; j++)
|
||||
{
|
||||
slong m, M;
|
||||
m = fac.p[i];
|
||||
M = (len /= m);
|
||||
t->cyc[num].m = m;
|
||||
t->cyc[num].M = M;
|
||||
t->cyc[num].dv = dv;
|
||||
t->cyc[num].z = z;
|
||||
t->cyc[num].dz = dz;
|
||||
/* TODO: ugly, reorder should solve this */
|
||||
if (num == t->num - 1)
|
||||
_acb_dft_precomp_init(t->cyc[num].pre, dv, z, dz, m, prec);
|
||||
else
|
||||
_acb_dft_precomp_init(t->cyc[num].pre, M, z, dz * M, m, prec);
|
||||
dv *= m;
|
||||
dz *= m;
|
||||
num++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
_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);
|
||||
_acb_dft_cyc_init_z_fac(t, fac, dv, NULL, 0, len, prec);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_cyc_clear(acb_dft_cyc_t t)
|
||||
{
|
||||
slong i;
|
||||
for (i = 0; i < t->num; i++)
|
||||
acb_dft_precomp_clear(t->cyc[i].pre);
|
||||
if (t->zclear)
|
||||
_acb_vec_clear(t->z, t->n);
|
||||
flint_free(t->cyc);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_cyc(acb_ptr w, acb_srcptr v, slong len, slong prec)
|
||||
{
|
||||
acb_dft_cyc_t cyc;
|
||||
acb_dft_cyc_init(cyc, len, prec);
|
||||
acb_dft_cyc_precomp(w, v, cyc, prec);
|
||||
acb_dft_cyc_clear(cyc);
|
||||
}
|
||||
|
|
@ -1,278 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dirichlet.h"
|
||||
|
||||
/* shallow extract */
|
||||
static acb_ptr
|
||||
vec_extract(acb_srcptr v, slong step, slong len)
|
||||
{
|
||||
slong k, l;
|
||||
acb_ptr res;
|
||||
res = flint_malloc(len * sizeof(acb_struct));
|
||||
for (k = 0, l = 0; k < len; k++, l+=step)
|
||||
res[k] = v[l];
|
||||
return res;
|
||||
}
|
||||
|
||||
void
|
||||
_acb_dft_base(acb_ptr w, acb_srcptr v, slong dv, acb_srcptr z, slong dz, slong n, slong prec)
|
||||
{
|
||||
acb_ptr v1, z1;
|
||||
v1 = vec_extract(v, dv, n);
|
||||
z1 = vec_extract(z, dz, n);
|
||||
_acb_dirichlet_dft_pol(w, v1, z1, n, prec);
|
||||
flint_free(v1);
|
||||
flint_free(z1);
|
||||
}
|
||||
|
||||
typedef struct
|
||||
{
|
||||
slong m;
|
||||
slong M;
|
||||
slong dv;
|
||||
acb_srcptr z;
|
||||
slong dz;
|
||||
}
|
||||
dft_cyc_step;
|
||||
|
||||
typedef struct
|
||||
{
|
||||
slong n;
|
||||
acb_ptr z;
|
||||
slong num;
|
||||
dft_cyc_step * cyc;
|
||||
}
|
||||
acb_dft_cyc_struct;
|
||||
|
||||
typedef acb_dft_cyc_struct acb_dft_cyc_t[1];
|
||||
|
||||
typedef struct
|
||||
{
|
||||
slong m;
|
||||
slong M;
|
||||
slong dv;
|
||||
acb_dft_cyc_t pre;
|
||||
}
|
||||
dft_prod_step;
|
||||
|
||||
typedef struct
|
||||
{
|
||||
slong n;
|
||||
slong num;
|
||||
dft_prod_step * cyc;
|
||||
}
|
||||
acb_dft_prod_struct;
|
||||
|
||||
typedef acb_dft_prod_struct acb_dft_prod_t[1];
|
||||
|
||||
void
|
||||
_acb_dirichlet_dft_cyc_init(acb_dft_cyc_t t, slong len, slong dv, slong prec)
|
||||
{
|
||||
slong i, j, num, dz;
|
||||
n_factor_t fac;
|
||||
acb_ptr z;
|
||||
|
||||
t->n = len;
|
||||
|
||||
n_factor_init(&fac);
|
||||
n_factor(&fac, len, 0);
|
||||
num = 0;
|
||||
for (i = 0; i < fac.num; i++)
|
||||
num += fac.exp[i];
|
||||
t->num = num;
|
||||
t->cyc = flint_malloc(num * sizeof(dft_cyc_step));
|
||||
t->z = z = _acb_vec_init(t->n);
|
||||
acb_dirichlet_vec_nth_roots(z, t->n, prec);
|
||||
|
||||
num = 0;
|
||||
dz = 1;
|
||||
for (i = 0; i < fac.num; i++)
|
||||
{
|
||||
for (j = 0; j < fac.exp[i]; j++)
|
||||
{
|
||||
slong m = fac.p[i];
|
||||
t->cyc[num].m = m;
|
||||
t->cyc[num].M = (len /= m);
|
||||
t->cyc[num].dv = dv;
|
||||
t->cyc[num].z = z;
|
||||
t->cyc[num].dz = dz;
|
||||
dv *= m;
|
||||
dz *= m;
|
||||
num++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
static void
|
||||
acb_dirichlet_dft_cyc_init(acb_dft_cyc_t t, slong len, slong prec)
|
||||
{
|
||||
_acb_dirichlet_dft_cyc_init(t, len, 1, prec);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dirichlet_dft_cyc_clear(acb_dft_cyc_t t)
|
||||
{
|
||||
_acb_vec_clear(t->z, t->n);
|
||||
flint_free(t->cyc);
|
||||
}
|
||||
|
||||
void
|
||||
_acb_dft_cyc(acb_ptr w, acb_srcptr v, dft_cyc_step * cyc, slong num, slong prec)
|
||||
{
|
||||
dft_cyc_step c;
|
||||
if (num == 0)
|
||||
abort(); /* or just copy v to w */
|
||||
c = cyc[0];
|
||||
if (num == 1)
|
||||
{
|
||||
_acb_dft_base(w, v, c.dv, c.z, c.dz, c.m, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
slong i, j;
|
||||
slong m = c.m, M = c.M, dv = c.dv, dz = c.dz;
|
||||
acb_srcptr z = c.z, vi;
|
||||
acb_ptr t, wi;
|
||||
|
||||
t = _acb_vec_init(m * M);
|
||||
wi = w; vi = v;
|
||||
for (i = 0; i < m; i++)
|
||||
{
|
||||
_acb_dft_cyc(wi, vi, cyc + 1, num - 1, prec);
|
||||
if (i)
|
||||
{
|
||||
for (j = 1; j < M; j++)
|
||||
acb_mul(wi + j, wi + j, z + dz * i * j, prec);
|
||||
}
|
||||
wi += M;
|
||||
vi += dv;
|
||||
}
|
||||
/* after first pass */
|
||||
for (j = 0; j < M; j++)
|
||||
_acb_dft_base(t + m * j, w + j, M, z, dz * M, m, prec);
|
||||
/* reorder */
|
||||
for (i = 0; i < m; i++)
|
||||
for (j = 0; j < M; j++)
|
||||
acb_set(w + j + M * i, t + i + m * j);
|
||||
_acb_vec_clear(t, m * M);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_dirichlet_dft_cyc_precomp(acb_ptr w, acb_srcptr v, acb_dft_cyc_t t, slong prec)
|
||||
{
|
||||
_acb_dft_cyc(w, v, t->cyc, t->num, prec);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dirichlet_dft_fast(acb_ptr w, acb_srcptr v, slong len, slong prec)
|
||||
{
|
||||
acb_dft_cyc_t t;
|
||||
acb_dirichlet_dft_cyc_init(t, len, prec);
|
||||
acb_dirichlet_dft_cyc_precomp(w, v, t, prec);
|
||||
acb_dirichlet_dft_cyc_clear(t);
|
||||
}
|
||||
|
||||
void
|
||||
_acb_dft_prod(acb_ptr w, acb_srcptr v, dft_prod_step * cyc, slong num, slong prec)
|
||||
{
|
||||
dft_prod_step c;
|
||||
if (num == 0)
|
||||
{
|
||||
flint_printf("error: reached num = 0 in dft_prod\n");
|
||||
abort(); /* or just copy v to w */
|
||||
}
|
||||
c = cyc[0];
|
||||
if (num == 1)
|
||||
{
|
||||
acb_dirichlet_dft_cyc_precomp(w, v, c.pre, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
slong i, j;
|
||||
slong m = c.m, M = c.M, dv = c.dv;
|
||||
acb_srcptr vi;
|
||||
acb_ptr t, wi;
|
||||
|
||||
t = _acb_vec_init(m * M);
|
||||
wi = w; vi = v;
|
||||
for (i = 0; i < m; i++)
|
||||
{
|
||||
_acb_dft_prod(wi, vi, cyc + 1, num - 1, prec);
|
||||
wi += M;
|
||||
vi += dv; /* here = M */
|
||||
}
|
||||
/* after first pass */
|
||||
for (j = 0; j < M; j++)
|
||||
acb_dirichlet_dft_cyc_precomp(t + m * j, w + j, c.pre, prec);
|
||||
/* reorder */
|
||||
for (i = 0; i < m; i++)
|
||||
for (j = 0; j < M; j++)
|
||||
acb_set(w + j + M * i, t + i + m * j);
|
||||
_acb_vec_clear(t, m * M);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_dirichlet_dft_prod_precomp(acb_ptr w, acb_srcptr v, acb_dft_prod_t t, slong prec)
|
||||
{
|
||||
if (t->num == 0)
|
||||
{
|
||||
acb_set(w + 0, v + 0);
|
||||
}
|
||||
else
|
||||
{
|
||||
_acb_dft_prod(w, v, t->cyc, t->num, prec);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_dirichlet_dft_prod_init(acb_dft_prod_t t, slong * cyc, slong num, slong prec)
|
||||
{
|
||||
slong i, len, dv;
|
||||
|
||||
t->num = num;
|
||||
t->cyc = flint_malloc(num * sizeof(dft_prod_step));
|
||||
|
||||
len = 1; dv = 1;
|
||||
for (i = 0; i < num; i++)
|
||||
len *= cyc[i];
|
||||
|
||||
for (i = 0; i < num; i++)
|
||||
{
|
||||
slong m = cyc[i];
|
||||
len /= m;
|
||||
t->cyc[i].m = m;
|
||||
t->cyc[i].M = len;
|
||||
t->cyc[i].dv = len;
|
||||
_acb_dirichlet_dft_cyc_init(t->cyc[i].pre, m, len, prec);
|
||||
dv *= m;
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_dirichlet_dft_prod_clear(acb_dft_prod_t t)
|
||||
{
|
||||
slong i;
|
||||
for (i = 0; i < t->num; i++)
|
||||
acb_dirichlet_dft_cyc_clear(t->cyc[i].pre);
|
||||
flint_free(t->cyc);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dirichlet_dft_prod(acb_ptr w, acb_srcptr v, slong * cyc, slong num, slong prec)
|
||||
{
|
||||
acb_dft_prod_t t;
|
||||
acb_dirichlet_dft_prod_init(t, cyc, num, prec);
|
||||
acb_dirichlet_dft_prod_precomp(w, v, t, prec);
|
||||
acb_dirichlet_dft_prod_clear(t);
|
||||
}
|
||||
|
|
@ -1,62 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_poly.h"
|
||||
#include "acb_dft.h"
|
||||
|
||||
/* all roots are already computed */
|
||||
void
|
||||
_acb_dft_pol(acb_ptr w, acb_srcptr v, slong dv, acb_srcptr z, slong dz, slong len, slong prec)
|
||||
{
|
||||
slong i, j;
|
||||
acb_ptr wi;
|
||||
acb_srcptr vj;
|
||||
|
||||
for (i = 0, wi = w; i < len; i++, wi++)
|
||||
{
|
||||
acb_zero(wi);
|
||||
for (j = 0, vj = v; j < len; j++, vj += dv)
|
||||
acb_addmul(wi, vj, z + dz * (i * j % len), prec);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_pol(acb_ptr w, acb_srcptr v, slong len, slong prec)
|
||||
{
|
||||
acb_ptr z;
|
||||
z = _acb_vec_init(len);
|
||||
_acb_vec_nth_roots(z, len, prec);
|
||||
_acb_dft_pol(w, v, 1, z, 1, len, prec);
|
||||
_acb_vec_clear(z, len);
|
||||
}
|
||||
|
||||
void
|
||||
_acb_dft_pol_init(acb_dft_pol_t pol, slong dv, acb_ptr z, slong dz, slong len, slong prec)
|
||||
{
|
||||
pol->n = len;
|
||||
pol->dv = dv;
|
||||
|
||||
if (z == NULL)
|
||||
{
|
||||
if (DFT_VERB)
|
||||
flint_printf("warning: init z[%ld] in dft_pol, should be avoided\n",len);
|
||||
pol->z = _acb_vec_init(len);
|
||||
_acb_vec_nth_roots(pol->z, len, prec);
|
||||
pol->dz = 1;
|
||||
pol->zclear = 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
pol->z = z;
|
||||
pol->dz = dz;
|
||||
pol->zclear = 0;
|
||||
}
|
||||
}
|
||||
|
|
@ -1,101 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dft.h"
|
||||
|
||||
void
|
||||
_acb_dft_precomp_init(acb_dft_pre_t pre, slong dv, acb_ptr z, slong dz, slong len, slong prec)
|
||||
{
|
||||
if (len <= 1)
|
||||
{
|
||||
flint_printf("precomp init: trivial group. abort.\n");
|
||||
abort();
|
||||
}
|
||||
if (n_is_prime(len))
|
||||
{
|
||||
/* TODO: need convolution if len is large */
|
||||
pre->type = DFT_POL;
|
||||
_acb_dft_pol_init(pre->t.pol, dv, z, dz, len, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
n_factor_t fac;
|
||||
|
||||
n_factor_init(&fac);
|
||||
n_factor(&fac, len, 1);
|
||||
|
||||
if (fac.num == 1)
|
||||
{
|
||||
/* TODO: could be p^e, or 2^e */
|
||||
pre->type = DFT_CYC;
|
||||
_acb_dft_cyc_init_z_fac(pre->t.cyc, fac, dv, z, dz, len, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
pre->type = DFT_CRT;
|
||||
_acb_dft_crt_init(pre->t.crt, dv, len, prec);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_precomp_init(acb_dft_pre_t pre, slong len, slong prec)
|
||||
{
|
||||
_acb_dft_precomp_init(pre, 1, NULL, 0, len, prec);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_precomp_clear(acb_dft_pre_t pre)
|
||||
{
|
||||
switch (pre->type)
|
||||
{
|
||||
case DFT_POL:
|
||||
acb_dft_pol_clear(pre->t.pol);
|
||||
break;
|
||||
case DFT_CYC:
|
||||
acb_dft_cyc_clear(pre->t.cyc);
|
||||
break;
|
||||
case DFT_PROD:
|
||||
acb_dft_prod_clear(pre->t.prod);
|
||||
break;
|
||||
case DFT_CRT:
|
||||
acb_dft_crt_clear(pre->t.crt);
|
||||
break;
|
||||
default:
|
||||
flint_printf("acb_dft_clear: unknown strategy code %i\n", pre->type);
|
||||
abort();
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_precomp(acb_ptr w, acb_srcptr v, const acb_dft_pre_t pre, slong prec)
|
||||
{
|
||||
switch (pre->type)
|
||||
{
|
||||
case DFT_POL:
|
||||
acb_dft_pol_precomp(w, v, pre->t.pol, prec);
|
||||
break;
|
||||
case DFT_CYC:
|
||||
acb_dft_cyc_precomp(w, v, pre->t.cyc, prec);
|
||||
break;
|
||||
case DFT_PROD:
|
||||
acb_dft_prod_precomp(w, v, pre->t.prod, prec);
|
||||
break;
|
||||
case DFT_CRT:
|
||||
acb_dft_crt_precomp(w, v, pre->t.crt, prec);
|
||||
break;
|
||||
default:
|
||||
flint_printf("acb_dft_precomp: unknown strategy code %i\n", pre->type);
|
||||
abort();
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
|
@ -1,67 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dft.h"
|
||||
|
||||
acb_dft_step_ptr
|
||||
_acb_dft_steps_prod(slong * cyc, slong num, slong prec)
|
||||
{
|
||||
slong i, len;
|
||||
acb_dft_step_ptr s;
|
||||
s = flint_malloc(num * sizeof(acb_dft_step_struct));
|
||||
|
||||
len = 1;
|
||||
for (i = 0; i < num; i++)
|
||||
len *= cyc[i];
|
||||
|
||||
for (i = 0; i < num; i++)
|
||||
{
|
||||
slong m, M;
|
||||
m = cyc[i];
|
||||
M = (len /= m);
|
||||
s[i].m = m;
|
||||
s[i].M = M;
|
||||
s[i].dv = M;
|
||||
s[i].dz = 0;
|
||||
s[i].z = NULL;
|
||||
_acb_dft_precomp_init(s[i].pre, M, NULL, 0, m, prec);
|
||||
}
|
||||
|
||||
return s;
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_prod_clear(acb_dft_prod_t t)
|
||||
{
|
||||
slong i;
|
||||
for (i = 0; i < t->num; i++)
|
||||
acb_dft_precomp_clear(t->cyc[i].pre);
|
||||
flint_free(t->cyc);
|
||||
}
|
||||
|
||||
|
||||
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
|
||||
acb_dft_step(w, v, prod->cyc, prod->num, prec);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_prod(acb_ptr w, acb_srcptr v, slong * cyc, slong num, slong prec)
|
||||
{
|
||||
acb_dft_prod_t t;
|
||||
acb_dft_prod_init(t, cyc, num, prec);
|
||||
acb_dft_prod_precomp(w, v, t, prec);
|
||||
acb_dft_prod_clear(t);
|
||||
}
|
||||
|
|
@ -1,80 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dft.h"
|
||||
|
||||
/* swap each element with one with bit-reversed index */
|
||||
void
|
||||
acb_dft_rad2_reorder(acb_ptr v, slong n)
|
||||
{
|
||||
slong i, j, k, l;
|
||||
|
||||
for (i = 0, l = n>>1; i < l; i++)
|
||||
{
|
||||
/* j = bit reversal of i */
|
||||
for (k = 1, j = 0; k < n; k <<= 1)
|
||||
{
|
||||
j <<= 1;
|
||||
if (i & k)
|
||||
j |= 1;
|
||||
}
|
||||
if (i < j)
|
||||
acb_swap(v + i, v + j);
|
||||
else if (i > j)
|
||||
acb_swap(v + n - 1 - i, v + n - 1 - j);
|
||||
i++, j |= l;
|
||||
acb_swap(v + i, v + j);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* remark: can use same rad2 with smaller power of 2 */
|
||||
void
|
||||
acb_dft_rad2_precomp(acb_ptr v, const acb_dft_rad2_t rad2, slong prec)
|
||||
{
|
||||
slong j, k, l;
|
||||
slong n = rad2->n, nz = rad2->nz;
|
||||
acb_ptr p, vend = v + n, w = rad2->z;
|
||||
acb_t tmp;
|
||||
acb_init(tmp);
|
||||
|
||||
acb_dft_rad2_reorder(v, n);
|
||||
|
||||
for (k = 1, l = nz; k < n; k <<= 1, l >>= 1)
|
||||
for (p = v; p < vend; p += k)
|
||||
for (j = 0; j < nz; j += l, p++)
|
||||
{
|
||||
acb_mul(tmp, p + k, w + j, prec);
|
||||
acb_sub(p + k, p + 0, tmp, prec);
|
||||
acb_add(p + 0, p + 0, tmp, prec);
|
||||
}
|
||||
|
||||
acb_clear(tmp);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_inverse_rad2_precomp(acb_ptr v, const acb_dft_rad2_t rad2, slong prec)
|
||||
{
|
||||
slong k, n = rad2->n;
|
||||
acb_dft_rad2_precomp(v, rad2, prec);
|
||||
_acb_vec_scalar_mul_2exp_si(v, v, n, - rad2->e);
|
||||
for (k = 1; k < n / 2; k++)
|
||||
acb_swap(v + k, v + n - k);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dft_rad2(acb_ptr v, int e, slong prec)
|
||||
{
|
||||
acb_dft_rad2_t rad2;
|
||||
acb_dft_rad2_init(rad2, e, prec);
|
||||
acb_dft_rad2_precomp(v, rad2, prec);
|
||||
acb_dft_rad2_clear(rad2);
|
||||
}
|
||||
|
|
@ -1,75 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dft.h"
|
||||
|
||||
#define REORDER 0
|
||||
|
||||
void
|
||||
acb_dft_step(acb_ptr w, acb_srcptr v, acb_dft_step_ptr cyc, slong num, slong prec)
|
||||
{
|
||||
acb_dft_step_struct c;
|
||||
if (num == 0)
|
||||
{
|
||||
flint_printf("error: reached num = 0 in acb_dft_step\n");
|
||||
abort(); /* or just copy v to w */
|
||||
}
|
||||
c = cyc[0];
|
||||
if (num == 1)
|
||||
{
|
||||
acb_dft_precomp(w, v, c.pre, prec);
|
||||
/*_acb_dft_base(w, v, c.dv, c.z, c.dz, c.m, prec);*/
|
||||
}
|
||||
else
|
||||
{
|
||||
slong i, j;
|
||||
slong m = c.m, M = c.M, dv = c.dv, dz = c.dz;
|
||||
acb_srcptr z = c.z;
|
||||
acb_ptr t;
|
||||
#if REORDER
|
||||
acb_ptr w2;
|
||||
#endif
|
||||
|
||||
t = _acb_vec_init(m * M);
|
||||
|
||||
/* m DFT of size M */
|
||||
for (i = 0; i < m; i++)
|
||||
acb_dft_step(w + i * M, v + i * dv, cyc + 1, num - 1, prec);
|
||||
|
||||
/* twiddle if non trivial product */
|
||||
if (c.z != NULL)
|
||||
{
|
||||
acb_ptr wi;
|
||||
for (wi = w + M, i = 1; i < m; i++, wi += M)
|
||||
for (j = 1; j < M; j++)
|
||||
acb_mul(wi + j, wi + j, z + dz * i * j, prec);
|
||||
}
|
||||
|
||||
#if REORDER
|
||||
/* 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++)
|
||||
acb_dft_precomp(t + m * j, w + j, c.pre, prec);
|
||||
|
||||
/* reorder */
|
||||
for (i = 0; i < m; i++)
|
||||
for (j = 0; j < M; j++)
|
||||
acb_set(w + j + M * i, t + i + m * j);
|
||||
|
||||
_acb_vec_clear(t, m * M);
|
||||
}
|
||||
}
|
||||
|
|
@ -1,118 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dft.h"
|
||||
|
||||
typedef void (*do_f) (acb_ptr z, acb_srcptr x, acb_srcptr y, slong len, slong prec);
|
||||
|
||||
void
|
||||
check_vec_eq_prec(acb_srcptr w1, acb_srcptr w2, slong len, slong prec, slong digits, ulong q, char f1[], char f2[])
|
||||
{
|
||||
slong i;
|
||||
|
||||
for (i = 0; i < len; i++)
|
||||
{
|
||||
if (!acb_overlaps(w1 + i, w2 + i))
|
||||
{
|
||||
flint_printf("FAIL\n\n");
|
||||
flint_printf("q = %wu, size = %wu\n", q, len);
|
||||
flint_printf("\nDFT differ from index %ld / %ld \n", i, len);
|
||||
flint_printf("\n%s =\n", f1);
|
||||
acb_vec_printd_index(w1, len, digits);
|
||||
flint_printf("\n%s =\n", f2);
|
||||
acb_vec_printd_index(w2, len, digits);
|
||||
flint_printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
else if (!acb_is_zero(w1+i) && (acb_rel_accuracy_bits(w1 + i) < 30
|
||||
|| acb_rel_accuracy_bits(w2 + i) < 30))
|
||||
{
|
||||
flint_printf("FAIL\n\n");
|
||||
flint_printf("q = %wu\n", q);
|
||||
flint_printf("\nDFT inaccurate from index %ld / %ld \n", i, len);
|
||||
flint_printf("\nnaive =\n");
|
||||
acb_printd(w1 + i, digits);
|
||||
flint_printf("\nfast =\n");
|
||||
acb_printd(w2 + i, digits);
|
||||
flint_printf("\nerrors %ld & %ld [prec = %wu]\n",
|
||||
acb_rel_accuracy_bits(w1 + i),
|
||||
acb_rel_accuracy_bits(w2 + i), prec);
|
||||
abort();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
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};
|
||||
flint_rand_t state;
|
||||
|
||||
slong f, nf = 2;
|
||||
do_f func[2] = { acb_dft_convol_naive, acb_dft_convol_rad2 };
|
||||
char * name[4] = { "naive", "rad2" };
|
||||
|
||||
flint_printf("convol....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
for (k = 0; k < nq; k++)
|
||||
{
|
||||
slong i;
|
||||
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]);
|
||||
|
||||
for (i = 0; i < q[k]; 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]));
|
||||
*/
|
||||
}
|
||||
|
||||
for (f = 0; f < nf; f++)
|
||||
{
|
||||
|
||||
acb_ptr z = (f == 0) ? z1 : z2;
|
||||
|
||||
func[f](z, x, y, q[k], prec);
|
||||
|
||||
if (f == 0)
|
||||
continue;
|
||||
|
||||
check_vec_eq_prec(z1, z2, q[k], prec, digits, q[k], 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]);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -1,127 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dft.h"
|
||||
|
||||
typedef void (*do_f) (acb_ptr w, acb_srcptr v, slong len, slong prec);
|
||||
|
||||
void
|
||||
check_vec_eq_prec(acb_srcptr w1, acb_srcptr w2, slong len, slong prec, slong digits, ulong q, char f1[], char f2[])
|
||||
{
|
||||
slong i;
|
||||
|
||||
for (i = 0; i < len; i++)
|
||||
{
|
||||
if (!acb_overlaps(w1 + i, w2 + i))
|
||||
{
|
||||
flint_printf("FAIL\n\n");
|
||||
flint_printf("q = %wu, size = %wu\n", q, len);
|
||||
flint_printf("\nDFT differ from index %ld / %ld \n", i, len);
|
||||
flint_printf("\n%s =\n", f1);
|
||||
acb_vec_printd_index(w1, len, digits);
|
||||
flint_printf("\n%s =\n", f2);
|
||||
acb_vec_printd_index(w2, len, digits);
|
||||
flint_printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
else if (acb_rel_accuracy_bits(w1 + i) < 30
|
||||
|| acb_rel_accuracy_bits(w2 + i) < 30)
|
||||
{
|
||||
flint_printf("FAIL\n\n");
|
||||
flint_printf("q = %wu\n", q);
|
||||
flint_printf("\nDFT inaccurate from index %ld / %ld \n", i, len);
|
||||
flint_printf("\nnaive =\n");
|
||||
acb_printd(w1 + i, digits);
|
||||
flint_printf("\nfast =\n");
|
||||
acb_printd(w2 + i, digits);
|
||||
flint_printf("\nerrors %ld & %ld [prec = %wu]\n",
|
||||
acb_rel_accuracy_bits(w1 + i),
|
||||
acb_rel_accuracy_bits(w2 + i), prec);
|
||||
abort();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
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};
|
||||
flint_rand_t state;
|
||||
|
||||
slong f, nf = 4;
|
||||
do_f func[4] = { acb_dft_pol, acb_dft_cyc, acb_dft_crt , acb_dft_bluestein };
|
||||
char * name[4] = { "pol", "cyc", "crt", "bluestein" };
|
||||
|
||||
flint_printf("dft....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
/* cyclic dft */
|
||||
for (k = 0; k < nq; k++)
|
||||
{
|
||||
slong i;
|
||||
acb_ptr v, w1, w2;
|
||||
|
||||
v = _acb_vec_init(q[k]);
|
||||
w1 = _acb_vec_init(q[k]);
|
||||
w2 = _acb_vec_init(q[k]);
|
||||
|
||||
for (i = 0; i < q[k]; i++)
|
||||
acb_set_si(v + i, i);
|
||||
|
||||
for (f = 0; f < nf; f++)
|
||||
{
|
||||
|
||||
acb_ptr w = (f == 0) ? w1 : w2;
|
||||
|
||||
func[f](w, v, q[k], prec);
|
||||
|
||||
if (f == 0)
|
||||
continue;
|
||||
|
||||
check_vec_eq_prec(w1, w2, q[k], prec, digits, q[k], name[0], name[f]);
|
||||
|
||||
}
|
||||
|
||||
_acb_vec_clear(v, q[k]);
|
||||
_acb_vec_clear(w1, q[k]);
|
||||
_acb_vec_clear(w2, q[k]);
|
||||
}
|
||||
|
||||
/* radix2 dft */
|
||||
for (k = 1; k < 12; k++)
|
||||
{
|
||||
slong n = 1 << k, j;
|
||||
acb_ptr v, w1, w2;
|
||||
v = w2 = _acb_vec_init(n);
|
||||
w1 = _acb_vec_init(n);
|
||||
|
||||
for (j = 0; j < n; j++)
|
||||
acb_set_si(v + k, k);
|
||||
|
||||
acb_dft_pol(w1, v, n, prec);
|
||||
acb_dft_rad2(v, k, prec);
|
||||
|
||||
check_vec_eq_prec(w1, w2, n, prec, digits, n, "pol", "rad2");
|
||||
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
|
|
@ -88,12 +88,6 @@ void acb_dirichlet_jacobi_sum(acb_t res, const dirichlet_group_t G, const dirich
|
|||
void acb_dirichlet_jacobi_sum_ui(acb_t res, const dirichlet_group_t G, ulong a, ulong b, slong prec);
|
||||
|
||||
void acb_dirichlet_l_hurwitz(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
|
||||
void acb_dirichlet_l_vec_hurwitz(acb_ptr res, const acb_t s, const dirichlet_group_t G, slong prec);
|
||||
|
||||
/* Discrete Fourier Transform */
|
||||
|
||||
void acb_dirichlet_dft_conrey(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec);
|
||||
void acb_dirichlet_dft(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec);
|
||||
|
||||
/* utils */
|
||||
|
||||
|
|
|
|||
|
|
@ -1,60 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dft.h"
|
||||
#include "acb_dirichlet.h"
|
||||
|
||||
/* dft, lexicographic conrey indexing, array size G->phi_q */
|
||||
void
|
||||
acb_dirichlet_dft_conrey(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec)
|
||||
{
|
||||
slong k, l, * cyc;
|
||||
cyc = flint_malloc(G->num * sizeof(slong));
|
||||
for (k = 0, l = G->num - 1; l >= 0; k++, l--)
|
||||
cyc[k] = G->P[k].phi;
|
||||
|
||||
acb_dft_prod(w, v, cyc, G->num, prec);
|
||||
flint_free(cyc);
|
||||
}
|
||||
|
||||
/* dft, number indexing, array size G->q */
|
||||
void
|
||||
acb_dirichlet_dft(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec)
|
||||
{
|
||||
ulong i, len;
|
||||
acb_ptr t1, t2;
|
||||
dirichlet_conrey_t x;
|
||||
|
||||
len = G->phi_q;
|
||||
t1 = flint_malloc(len * sizeof(acb_struct));
|
||||
|
||||
dirichlet_conrey_init(x, G);
|
||||
dirichlet_conrey_one(x, G);
|
||||
for (i = 0; i < len; i++)
|
||||
{
|
||||
t1[i] = v[x->n];
|
||||
dirichlet_conrey_next(x, G);
|
||||
};
|
||||
|
||||
t2 = _acb_vec_init(len);
|
||||
acb_dirichlet_dft_conrey(t2, t1, G, prec);
|
||||
|
||||
dirichlet_conrey_one(x, G);
|
||||
for (i = 0; i < len; i++)
|
||||
{
|
||||
acb_set(w + x->n, t2 + i);
|
||||
dirichlet_conrey_next(x, G);
|
||||
};
|
||||
|
||||
_acb_vec_clear(t2, len);
|
||||
dirichlet_conrey_clear(x);
|
||||
flint_free(t1);
|
||||
}
|
||||
|
|
@ -1,65 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dirichlet.h"
|
||||
#include "acb_poly.h"
|
||||
|
||||
void
|
||||
acb_dirichlet_l_vec_hurwitz(acb_ptr res, const acb_t s,
|
||||
const dirichlet_group_t G, slong prec)
|
||||
{
|
||||
acb_t a, qs;
|
||||
acb_ptr zeta, z;
|
||||
dirichlet_conrey_t cn;
|
||||
int deflate;
|
||||
|
||||
/* remove pole in Hurwitz zeta at s = 1 */
|
||||
deflate = acb_is_one(s);
|
||||
|
||||
dirichlet_conrey_init(cn, G);
|
||||
acb_init(qs);
|
||||
acb_init(a);
|
||||
|
||||
prec += n_clog(G->phi_q, 2);
|
||||
|
||||
acb_set_ui(qs, G->q);
|
||||
acb_neg(a, s);
|
||||
acb_pow(qs, qs, a, prec);
|
||||
|
||||
zeta = z = _acb_vec_init(G->phi_q);
|
||||
dirichlet_conrey_one(cn, G);
|
||||
do {
|
||||
|
||||
acb_set_ui(a, cn->n);
|
||||
acb_div_ui(a, a, G->q, prec);
|
||||
|
||||
if (!deflate)
|
||||
acb_hurwitz_zeta(z, s, a, prec);
|
||||
else
|
||||
_acb_poly_zeta_cpx_series(z, s, a, 1, 1, prec);
|
||||
|
||||
acb_mul(z, z, qs, prec);
|
||||
|
||||
z++;
|
||||
} while (dirichlet_conrey_next(cn, G) >= 0);
|
||||
|
||||
acb_dirichlet_dft_conrey(res, zeta, G, prec);
|
||||
|
||||
/* restore pole for the principal character */
|
||||
if (deflate)
|
||||
acb_indeterminate(res);
|
||||
|
||||
dirichlet_conrey_clear(cn);
|
||||
_acb_vec_clear(zeta, G->phi_q);
|
||||
acb_clear(qs);
|
||||
acb_clear(a);
|
||||
}
|
||||
|
||||
|
|
@ -1,125 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dirichlet.h"
|
||||
|
||||
typedef void (*do_f) (acb_ptr w, acb_srcptr v, slong len, slong prec);
|
||||
|
||||
void
|
||||
check_vec_eq_prec(acb_srcptr w1, acb_srcptr w2, slong len, slong prec, slong digits, ulong q, char f1[], char f2[])
|
||||
{
|
||||
slong i;
|
||||
|
||||
for (i = 0; i < len; i++)
|
||||
{
|
||||
if (!acb_overlaps(w1 + i, w2 + i))
|
||||
{
|
||||
flint_printf("FAIL\n\n");
|
||||
flint_printf("q = %wu, size = %wu\n", q, len);
|
||||
flint_printf("\nDFT differ from index %ld / %ld \n", i, len);
|
||||
flint_printf("\n%s =\n", f1);
|
||||
acb_vec_printd(w1, len, digits);
|
||||
flint_printf("\n%s =\n", f2);
|
||||
acb_vec_printd(w2, len, digits);
|
||||
flint_printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
else if (acb_rel_accuracy_bits(w1 + i) < 30
|
||||
|| acb_rel_accuracy_bits(w2 + i) < 30)
|
||||
{
|
||||
flint_printf("FAIL\n\n");
|
||||
flint_printf("q = %wu\n", q);
|
||||
flint_printf("\nDFT inaccurate from index %ld / %ld \n", i, len);
|
||||
flint_printf("\nnaive =\n");
|
||||
acb_printd(w1 + i, digits);
|
||||
flint_printf("\nfast =\n");
|
||||
acb_printd(w2 + i, digits);
|
||||
flint_printf("\nerrors %ld & %ld [prec = %wu]\n",
|
||||
acb_rel_accuracy_bits(w1 + i),
|
||||
acb_rel_accuracy_bits(w2 + i), prec);
|
||||
abort();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
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};
|
||||
flint_rand_t state;
|
||||
|
||||
flint_printf("dft....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
/* Dirichlet group DFT */
|
||||
for (k = 0; k < nq; k++)
|
||||
{
|
||||
slong i, j, len;
|
||||
dirichlet_group_t G;
|
||||
dirichlet_conrey_t x, y;
|
||||
acb_t chiy;
|
||||
acb_ptr v, w1, w2;
|
||||
|
||||
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);
|
||||
|
||||
dirichlet_conrey_init(x, G);
|
||||
dirichlet_conrey_one(x, G);
|
||||
for (i = 0; i < len; i++)
|
||||
acb_randtest_precise(v + i, state, prec, 0);
|
||||
|
||||
/* naive */
|
||||
acb_init(chiy);
|
||||
dirichlet_conrey_init(y, G);
|
||||
dirichlet_conrey_one(x, G);
|
||||
for (i = 0; i < len; i++)
|
||||
{
|
||||
acb_zero(w1 + i);
|
||||
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);
|
||||
dirichlet_conrey_next(y, G);
|
||||
}
|
||||
dirichlet_conrey_next(x, G);
|
||||
}
|
||||
acb_clear(chiy);
|
||||
dirichlet_conrey_clear(y);
|
||||
dirichlet_conrey_clear(x);
|
||||
|
||||
/* dft */
|
||||
acb_dirichlet_dft_conrey(w2, v, G, prec);
|
||||
|
||||
check_vec_eq_prec(w1, w2, len, prec, digits, q[k], "naive", "group");
|
||||
|
||||
_acb_vec_clear(v, len);
|
||||
_acb_vec_clear(w1, len);
|
||||
_acb_vec_clear(w2, len);
|
||||
|
||||
dirichlet_group_clear(G);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
|
@ -14,77 +14,6 @@
|
|||
#define nq 5
|
||||
#define nx 3
|
||||
|
||||
void
|
||||
test_dft(ulong q)
|
||||
{
|
||||
ulong i;
|
||||
slong prec = 100;
|
||||
dirichlet_group_t G;
|
||||
dirichlet_conrey_t x;
|
||||
dirichlet_char_t chi;
|
||||
acb_t s, z;
|
||||
acb_ptr v;
|
||||
|
||||
dirichlet_group_init(G, q);
|
||||
dirichlet_conrey_init(x, G);
|
||||
dirichlet_char_init(chi, G);
|
||||
|
||||
acb_init(s);
|
||||
acb_one(s);
|
||||
acb_div_si(s, s, 2, prec);
|
||||
|
||||
v = _acb_vec_init(G->phi_q);
|
||||
|
||||
/* all at once */
|
||||
acb_dirichlet_l_vec_hurwitz(v, s, G, prec);
|
||||
|
||||
/* check with complete loop */
|
||||
|
||||
i = 0;
|
||||
acb_init(z);
|
||||
dirichlet_conrey_one(x, G);
|
||||
do {
|
||||
|
||||
dirichlet_char_conrey(chi, G, x);
|
||||
acb_dirichlet_l_hurwitz(z, s, G, chi, prec);
|
||||
|
||||
if (!acb_overlaps(z, v + i))
|
||||
{
|
||||
flint_printf("\n L value differ");
|
||||
flint_printf("\nL(1/2, %wu) single = ", x->n);
|
||||
acb_printd(z, 20);
|
||||
flint_printf("\nL(1/2, %wu) multi = ", x->n);
|
||||
acb_printd(v + i, 20);
|
||||
flint_printf("\n\n");
|
||||
acb_vec_printd(v, G->phi_q, 10);
|
||||
flint_printf("\n\n");
|
||||
}
|
||||
else if (acb_rel_accuracy_bits(z) < prec - 8
|
||||
|| acb_rel_accuracy_bits(v + i) < prec - 8)
|
||||
{
|
||||
flint_printf("FAIL\n\n");
|
||||
flint_printf("q = %wu\n", q);
|
||||
flint_printf("\nL(1/2,chi_%wu(%wu,)) inaccurate\n", q, x->n);
|
||||
flint_printf("\nsingle =\n");
|
||||
acb_printd(z, 30);
|
||||
flint_printf("\ndft =\n");
|
||||
acb_printd(v + i, 30);
|
||||
flint_printf("\nerrors %ld & %ld [prec = %wu]\n",
|
||||
acb_rel_accuracy_bits(z),
|
||||
acb_rel_accuracy_bits(v + i), prec);
|
||||
abort();
|
||||
}
|
||||
|
||||
i++;
|
||||
} while (dirichlet_conrey_next(x, G) >= 0);
|
||||
|
||||
acb_clear(s);
|
||||
_acb_vec_clear(v, G->phi_q);
|
||||
dirichlet_char_clear(chi);
|
||||
dirichlet_conrey_clear(x);
|
||||
dirichlet_group_clear(G);
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
|
||||
|
|
@ -204,9 +133,6 @@ int main()
|
|||
dirichlet_char_clear(chi);
|
||||
dirichlet_group_clear(G);
|
||||
|
||||
/* test using dft */
|
||||
test_dft(q[i]);
|
||||
|
||||
}
|
||||
acb_clear(ref);
|
||||
acb_clear(res);
|
||||
|
|
|
|||
|
|
@ -1,56 +0,0 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dirichlet.h"
|
||||
|
||||
void
|
||||
acb_dirichlet_vec_mellin_arb(acb_ptr res, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, const arb_t t, slong n, slong prec)
|
||||
{
|
||||
slong k;
|
||||
arb_t tk, xt, stk, st;
|
||||
acb_ptr a;
|
||||
mag_t e;
|
||||
a = _acb_vec_init(len);
|
||||
acb_dirichlet_chi_vec(a, G, chi, len, prec);
|
||||
if (chi->parity)
|
||||
{
|
||||
for (k = 2; k < len; k++)
|
||||
acb_mul_si(a + k, a + k, k, prec);
|
||||
}
|
||||
arb_init(tk);
|
||||
arb_init(xt);
|
||||
arb_init(st);
|
||||
arb_init(stk);
|
||||
mag_init(e);
|
||||
|
||||
arb_sqrt(st, t, prec);
|
||||
arb_one(tk);
|
||||
arb_one(stk);
|
||||
for (k = 0; k < n; k++)
|
||||
{
|
||||
_acb_dirichlet_theta_argument_at_arb(xt, G->q, tk, prec);
|
||||
mag_tail_kexpk2_arb(e, xt, len);
|
||||
arb_neg(xt, xt);
|
||||
arb_exp(xt, xt, prec);
|
||||
/* TODO: reduce len */
|
||||
acb_dirichlet_qseries_arb(res + k, a, xt, len, prec);
|
||||
acb_add_error_mag(res + k, e);
|
||||
acb_mul_arb(res + k, res + k, stk, prec);
|
||||
arb_mul(tk, tk, t, prec);
|
||||
arb_mul(stk, stk, st, prec);
|
||||
}
|
||||
mag_clear(e);
|
||||
arb_clear(xt);
|
||||
arb_clear(tk);
|
||||
arb_clear(stk);
|
||||
arb_clear(st);
|
||||
_acb_vec_clear(a, len);
|
||||
}
|
||||
|
|
@ -1,110 +0,0 @@
|
|||
.. _acb-dft:
|
||||
|
||||
**acb_dft.h** -- Discrete Fourier Transform on finite abelian groups
|
||||
===================================================================================
|
||||
|
||||
*Warning: the interfaces in this module are experimental and may change
|
||||
without notice.*
|
||||
|
||||
Let *G* be a finite abelian group, and `\chi` a character of *G*.
|
||||
For any map `f:G\to\mathbb C`, the discrete fourier transform
|
||||
`\hat f:\hat G\to \mathbb C` is defined by
|
||||
|
||||
.. math::
|
||||
|
||||
\hat f(\chi) = \sum_{x\in G}\overline{\chi(x)}f(x)
|
||||
|
||||
|
||||
Fast Fourier Transform techniques allow to compute efficiently
|
||||
all values `\hat f(\chi)` by reusing common computations.
|
||||
|
||||
Specifically, if `H\triangleleft G` is a subgroup of size `M` and index
|
||||
`[G:H]=m`, then writing `f_x(h)=f(xh)` the translate of `f` by representatives
|
||||
`x` of `G/H`, one has a decomposition
|
||||
|
||||
.. math::
|
||||
|
||||
\hat f(\chi) = \sum_{x\in G/H} \overline{\chi(x)} \hat{f_x}(\chi_{H})
|
||||
|
||||
so that the DFT on `G` can be computed using `m` DFT on `H` (of
|
||||
appropriate translates of `f`), then `M` DFT on `G/H`, one for
|
||||
each restriction `\chi_{H}`.
|
||||
|
||||
This decomposition can be done recursively.
|
||||
|
||||
Note that by inversion formula
|
||||
|
||||
.. math::
|
||||
|
||||
\widehat{\hat f}(\chi) = \#G\times f(\chi^{-1})
|
||||
|
||||
it is straightforward to recover `f` from its DFT `\hat f`.
|
||||
|
||||
DFT on Z/nZ
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
If `G=\mathbb Z/n\mathbb Z`, we compute the DFT according to the usual convention
|
||||
|
||||
.. math::
|
||||
|
||||
w_x = \sum_{y\bmod n} v_y e^{-\frac{2iπ}nxy}
|
||||
|
||||
.. function:: void acb_dirichlet_dft_pol(acb_ptr w, acb_srcptr v, slong n, slong prec)
|
||||
|
||||
.. function:: void acb_dirichlet_dft_crt(acb_ptr w, acb_srcptr v, slong n, slong prec)
|
||||
|
||||
.. function:: void acb_dirichlet_dft_cyc(acb_ptr w, acb_srcptr v, slong n, slong prec)
|
||||
|
||||
Set *w* to the DFT of *v* of length *len*.
|
||||
|
||||
The first variant uses the naive `O(n^2)` algorithm.
|
||||
|
||||
The second one uses CRT to express `Z/nZ` as a product of cyclic groups.
|
||||
|
||||
The *cyc* version uses each prime factor of `m` of `n` to decompose with
|
||||
the subgroup `H=m\mathbb Z/n\mathbb Z`.
|
||||
|
||||
|
||||
DFT on products
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
A finite abelian group is isomorphic to a product of cyclic components
|
||||
|
||||
.. math::
|
||||
|
||||
G = \bigoplus_{i=1}^r \mathbb Z/n_i\mathbb Z
|
||||
|
||||
then a character is a product of characters of all components and the DFT reads
|
||||
|
||||
.. math::
|
||||
|
||||
\hat f(x_1,\dots x_r) = \sum_{y_1\dots y_r} f(y_1,\dots y_r)
|
||||
e^{-2iπ\sum\frac{x_i y_i}{n_i}}
|
||||
|
||||
We assume that `f` is given by a vector of length `\prod n_i` corresponding
|
||||
to a lexicographic ordering of the values `y_1,\dots y_r`, and the computation
|
||||
returns the same indexing for values of `\hat f`.
|
||||
|
||||
.. function:: void acb_dirichlet_dft_prod(acb_ptr w, acb_srcptr v, slong * cyc, slong num, slong prec);
|
||||
|
||||
Computes the DFT on the group product of *num* cyclic components of sizes *cyc*.
|
||||
|
||||
Precomputations
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
Convolution
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
For functions `f` and `g` on `G` we consider the convolution
|
||||
|
||||
.. math::
|
||||
|
||||
(f \star g)(x) = \sum_{y\in G} f(x-y)g(y)
|
||||
|
||||
which satisfies
|
||||
|
||||
.. math::
|
||||
|
||||
\widehat{f \star g}(\chi) = \hat f(\chi)\hat g(\chi)
|
||||
|
||||
|
||||
|
|
@ -34,38 +34,6 @@ Character evaluation
|
|||
|
||||
Compute the *nv* first Dirichlet values.
|
||||
|
||||
Roots of unity
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
.. type:: acb_dirichlet_powers_struct
|
||||
|
||||
.. type:: acb_dirichlet_powers_t
|
||||
|
||||
This structure allows to compute *n* powers of a fixed root of unity of order *m*
|
||||
using precomputations. Extremal cases are
|
||||
|
||||
- all powers are stored: `O(m)` initialization + storage, `O(n)` evaluations
|
||||
|
||||
- nothing stored: `O(1)` initialization + storage, `O(\log(m)n)` evaluations
|
||||
|
||||
- `k` step decomposition: `O(k m^{\frac1k})` init + storage, `O((k-1)n)` evaluations.
|
||||
|
||||
Currently, only baby-step giant-step decomposition (i.e. `k=2`)
|
||||
is implemented, allowing to obtain each power using one multiplication.
|
||||
|
||||
.. function:: void acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong num, slong prec)
|
||||
|
||||
Initialize the powers structure for *num* evaluations of powers of the root of unity
|
||||
of order *order*.
|
||||
|
||||
.. function:: void acb_dirichlet_powers_clear(acb_dirichlet_powers_t t)
|
||||
|
||||
Clears *t*.
|
||||
|
||||
.. function:: void acb_dirichlet_power(acb_t z, const acb_dirichlet_powers_t t, ulong n, slong prec)
|
||||
|
||||
Sets *z* to `x^n` where *t* contains precomputed powers of `x`.
|
||||
|
||||
Gauss and Jacobi sums
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
|
|
@ -163,77 +131,6 @@ For `\Re(t)>0` we write `x(t)=\exp(-\frac{\pi}{N}t^2)` and define
|
|||
should be used, which is not done automatically (to avoid recomputing the
|
||||
Gauss sum).
|
||||
|
||||
.. function:: ulong acb_dirichlet_theta_length(ulong q, const arb_t t, slong prec)
|
||||
|
||||
Compute the number of terms to be summed in the theta series of argument *t*
|
||||
so that the tail is less than `2^{-\mathrm{prec}}`.
|
||||
|
||||
.. function:: void acb_dirichlet_qseries_powers_naive(acb_t res, const arb_t x, int p, const ulong * a, const acb_dirichlet_powers_t z, slong len, slong prec)
|
||||
|
||||
.. function:: void acb_dirichlet_qseries_powers_smallorder(acb_t res, const arb_t x, int p, const ulong * a, const acb_dirichlet_powers_t z, slong len, slong prec)
|
||||
|
||||
Compute the series `\sum n^p z^{a_n} x^{n^2}` for exponent list *a*,
|
||||
precomputed powers *z* and parity *p* (being 0 or 1).
|
||||
|
||||
The *naive* version sums the series as defined, while the *smallorder*
|
||||
variant evaluates the series on the quotient ring by a cyclotomic polynomial
|
||||
before evaluating at the root of unity, ignoring its argument *z*.
|
||||
|
||||
Discrete Fourier Transforms (DFT)
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
Let *G* be a finite abelian group, and `\chi` a character of *G*.
|
||||
For any map `f:G\to\mathbb C`, the discrete fourier transform
|
||||
`\hat f:\hat G\to \mathbb C` is defined by
|
||||
|
||||
.. math::
|
||||
|
||||
\hat f(\chi) = \sum_{x\in G}\chi(x)f(x)
|
||||
|
||||
Fast Fourier Transform techniques allow to compute efficiently
|
||||
all values `\hat f(\chi)`.
|
||||
|
||||
For a Dirichlet group `G` modulo `q`, we take advantage
|
||||
of the Conrey isomorphism `G \to \hat G` to consider the
|
||||
the Fourier transform on Conrey labels as
|
||||
|
||||
.. math::
|
||||
|
||||
g(a) = \sum_{b\bmod q}\chi_q(a,b)f(b)
|
||||
|
||||
|
||||
.. function:: void acb_dirichlet_dft_conrey(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec)
|
||||
|
||||
Compute the DFT of *v* using Conrey indices.
|
||||
This function assumes *v* and *w* are vectors
|
||||
of size *G->phi_q*, whose values correspond to a lexicographic ordering
|
||||
of Conrey logs (as obtained using :func:`dirichlet_conrey_next` or
|
||||
by :func:`dirichlet_index_conrey`).
|
||||
|
||||
For example, if `q=15`, the Conrey elements are stored in following
|
||||
order
|
||||
|
||||
======= ============= =====================
|
||||
index log = [e,f] number = 7^e 11^f
|
||||
======= ============= =====================
|
||||
0 [0, 0] 1
|
||||
1 [0, 1] 7
|
||||
2 [0, 2] 4
|
||||
3 [0, 3] 13
|
||||
4 [0, 4] 1
|
||||
5 [1, 0] 11
|
||||
6 [1, 1] 2
|
||||
7 [1, 2] 14
|
||||
8 [1, 3] 8
|
||||
9 [1, 4] 11
|
||||
======= ============= =====================
|
||||
|
||||
.. function:: void acb_dirichlet_dft(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec)
|
||||
|
||||
Compute the DFT of *v* using Conrey numbers.
|
||||
This function assumes *v* and *w* are vectors of size *G->q*.
|
||||
All values at index not coprime to *G->q* are ignored.
|
||||
|
||||
Euler products
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
|
|
@ -301,10 +198,3 @@ L-functions
|
|||
is used to avoid poles.
|
||||
|
||||
This formula is slow for large *q*.
|
||||
|
||||
.. function:: void acb_dirichlet_l_vec_hurwitz(acb_ptr res, const acb_t s, const dirichlet_group_t G, slong prec)
|
||||
|
||||
Compute all values `L(s,\chi)` for `\chi` mod `q`, by Hurwitz formula and
|
||||
discrete Fourier transform.
|
||||
*res* is assumed to have length *G->phi_q* and values are stored by lexicographically ordered
|
||||
Conrey logs. See :func:`acb_dirichlet_dft_conrey`.
|
||||
|
|
|
|||
|
|
@ -123,7 +123,6 @@ modules.
|
|||
bernoulli.rst
|
||||
hypgeom.rst
|
||||
partitions.rst
|
||||
acb_dft.rst
|
||||
|
||||
Calculus
|
||||
::::::::::::::::::::::::::::::::::::
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue