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

acb_dirichlet roots of unity: several algorithm improvements; rename, test, document

Browse source
This commit is contained in:
Fredrik Johansson 2016-11-22 21:48:29 +01:00
parent 4291b7e331
commit 3383033cea
11 changed files with 343 additions and 115 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

19
acb_dirichlet.h
View file

@ -67,31 +67,32 @@ void acb_dirichlet_eta(acb_t res, const acb_t s, slong prec);
void acb_dirichlet_pairing(acb_t res, const dirichlet_group_t G, ulong m, ulong n, slong prec);
void acb_dirichlet_pairing_char(acb_t res, const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b, slong prec);
/* precompute powers of a root of unity */
/* precompute roots of unity */
typedef struct
{
ulong order;
ulong reduced_order;
acb_t z;
slong size;
slong depth;
acb_ptr * Z;
int use_pow;
}
acb_dirichlet_powers_struct;
acb_dirichlet_roots_struct;
typedef acb_dirichlet_powers_struct acb_dirichlet_powers_t[1];
typedef acb_dirichlet_roots_struct acb_dirichlet_roots_t[1];
void _acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong size, slong depth, slong prec);
void acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong num, slong prec);
void acb_dirichlet_powers_clear(acb_dirichlet_powers_t t);
void acb_dirichlet_power(acb_t z, const acb_dirichlet_powers_t t, ulong n, slong prec);
void acb_dirichlet_roots_init(acb_dirichlet_roots_t t, ulong order, slong num, slong prec);
void acb_dirichlet_roots_clear(acb_dirichlet_roots_t t);
void acb_dirichlet_root(acb_t z, const acb_dirichlet_roots_t t, ulong n, slong prec);
void acb_dirichlet_chi(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, ulong n, slong prec);
void acb_dirichlet_chi_vec(acb_ptr v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv, slong prec);
void acb_dirichlet_arb_quadratic_powers(arb_ptr v, slong nv, const arb_t x, slong prec);
void acb_dirichlet_qseries_arb(acb_t res, acb_srcptr a, const arb_t x, slong len, slong prec);
void acb_dirichlet_qseries_arb_powers_naive(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec);
void acb_dirichlet_qseries_arb_powers_smallorder(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec);
void acb_dirichlet_qseries_arb_powers_naive(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_roots_t z, slong len, slong prec);
void acb_dirichlet_qseries_arb_powers_smallorder(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_roots_t z, slong len, slong prec);
ulong acb_dirichlet_theta_length_d(ulong q, double x, slong prec);
ulong acb_dirichlet_theta_length(ulong q, const arb_t x, slong prec);

8
acb_dirichlet/chi_vec.c
View file

@ -16,23 +16,23 @@ acb_dirichlet_chi_vec(acb_ptr v, const dirichlet_group_t G, const dirichlet_char
{
slong k;
ulong * a, order;
acb_dirichlet_powers_t t;
acb_dirichlet_roots_t t;
a = flint_malloc(nv * sizeof(ulong));
order = dirichlet_order_char(G, chi);
dirichlet_chi_vec_order(a, G, chi, order, nv);
acb_dirichlet_powers_init(t, order, nv, prec);
acb_dirichlet_roots_init(t, order, nv, prec);
acb_zero(v + 0);
for (k = 0; k < nv; k++)
{
if (a[k] != DIRICHLET_CHI_NULL)
acb_dirichlet_power(v + k, t, a[k], prec);
acb_dirichlet_root(v + k, t, a[k], prec);
else
acb_zero(v + k);
}
acb_dirichlet_powers_clear(t);
acb_dirichlet_roots_clear(t);
flint_free(a);
}

36
acb_dirichlet/power.c
View file

@ -1,36 +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_power(acb_t z, const acb_dirichlet_powers_t t, ulong n, slong prec)
{
if (!t->depth)
{
acb_pow_ui(z, t->z, n, prec);
}
else
{
slong k;
ulong r;
r = n % t->size;
n = n / t->size;
acb_set(z, t->Z[0] + r);
for (k = 1; k < t->depth && n; k++)
{
r = n % t->size;
n = n / t->size;
acb_mul(z, z, t->Z[k] + r, prec);
}
}
}

51
acb_dirichlet/powers_init.c
View file

@ -1,51 +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_powers_init(acb_dirichlet_powers_t t, ulong order, slong size, slong depth, slong prec)
{
slong k;
t->order = order;
acb_init(t->z);
acb_nth_root(t->z, order, prec);
t->size = size;
t->depth = depth;
if (depth)
{
acb_struct * z;
z = t->z + 0;
t->Z = flint_malloc(depth * sizeof(acb_ptr));
for (k = 0; k < depth; k++)
{
t->Z[k] = _acb_vec_init(size);
_acb_vec_set_powers(t->Z[k], z, size, prec);
z = t->Z[k] + 1;
}
}
else
{
t->Z = NULL;
}
}
void
acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong num, slong prec)
{
slong size, depth;
depth = (num > order) ? 1 : n_flog(order, num);
size = n_root(order, depth) + 1;
_acb_dirichlet_powers_init(t, order, size, depth, prec);
}

13
acb_dirichlet/qseries_arb_powers.c
View file

@ -13,7 +13,7 @@
#include "acb_poly.h"
void
acb_dirichlet_qseries_arb_powers_naive(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec)
acb_dirichlet_qseries_arb_powers_naive(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_roots_t roots, slong len, slong prec)
{
slong k;
arb_t xk2, dx, x2;
@ -36,7 +36,7 @@ acb_dirichlet_qseries_arb_powers_naive(acb_t res, const arb_t x, int parity, con
arb_mul(xk2, xk2, dx, prec);
if (a[k] != DIRICHLET_CHI_NULL)
{
acb_dirichlet_power(zk, z, a[k], prec);
acb_dirichlet_root(zk, roots, a[k], prec);
if (parity)
acb_mul_si(zk, zk, k, prec);
acb_addmul_arb(res, zk, xk2, prec);
@ -51,10 +51,10 @@ acb_dirichlet_qseries_arb_powers_naive(acb_t res, const arb_t x, int parity, con
/* small order, multiply by chi at the end */
void
acb_dirichlet_qseries_arb_powers_smallorder(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec)
acb_dirichlet_qseries_arb_powers_smallorder(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_roots_t roots, slong len, slong prec)
{
slong k;
ulong order = z->order;
ulong order = roots->order;
arb_t xk2, kxk2, dx, x2;
acb_ptr t;
arb_init(xk2);
@ -90,8 +90,9 @@ acb_dirichlet_qseries_arb_powers_smallorder(acb_t res, const arb_t x, int parity
}
}
/* now Hörner */
_acb_poly_evaluate(res, t, order, z->z, prec);
/* now Horner */
acb_dirichlet_root(res, roots, 1, prec);
_acb_poly_evaluate(res, t, order, res, prec);
_acb_vec_clear(t, order);
arb_clear(xk2);

103
acb_dirichlet/root.c Normal file
View file

@ -0,0 +1,103 @@
/*
Copyright (C) 2016 Pascal Molin
Copyright (C) 2016 Fredrik Johansson
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_root(acb_t z, const acb_dirichlet_roots_t t, ulong k, slong prec)
{
ulong n = t->order;
int swap, flip, conjugate;
slong wp;
if (k > n)
k %= n;
swap = flip = conjugate = 0;
if (k > n / 2)
{
conjugate = 1;
k = n - k;
}
if (n % 2 == 0 && k > n / 4)
{
flip = 1;
k = n / 2 - k;
}
if (n % 4 == 0 && k > n / 8)
{
swap = 1;
k = n / 4 - k;
}
wp = prec + 6 + 2 * FLINT_BIT_COUNT(t->reduced_order);
if (k == 0)
{
acb_one(z);
}
else if (t->depth == 0)
{
if (t->use_pow)
{
acb_pow_ui(z, t->z, k, wp);
acb_set_round(z, z, prec);
}
else
{
ulong r;
fmpq_t t;
fmpq_init(t);
r = n_gcd(n, 2 * k); /* no overflow since since k <= n / 2 */
fmpz_set_ui(fmpq_numref(t), (2 * k) / r);
fmpz_set_ui(fmpq_denref(t), n / r);
arb_sin_cos_pi_fmpq(acb_imagref(z), acb_realref(z), t, prec);
fmpq_clear(t);
}
}
else if (t->depth == 1)
{
acb_set_round(z, t->Z[0] + k, prec);
}
else
{
slong j;
ulong r;
r = k % t->size;
k = k / t->size;
acb_set(z, t->Z[0] + r);
for (j = 1; j < t->depth && k != 0; j++)
{
r = k % t->size;
k = k / t->size;
acb_mul(z, z, t->Z[j] + r, wp);
}
if (k != 0)
abort();
acb_set_round(z, z, prec);
}
if (swap)
arb_swap(acb_realref(z), acb_imagref(z));
if (flip)
arb_neg(acb_realref(z), acb_realref(z));
if (conjugate)
arb_neg(acb_imagref(z), acb_imagref(z));
}

5
acb_dirichlet/powers_clear.c → acb_dirichlet/roots_clear.c
View file

@ -12,13 +12,14 @@
#include "acb_dirichlet.h"
void
acb_dirichlet_powers_clear(acb_dirichlet_powers_t t)
acb_dirichlet_roots_clear(acb_dirichlet_roots_t t)
{
slong k;
for (k = 0; k < t->depth; k++)
_acb_vec_clear(t->Z[k], t->size);
_acb_vec_clear(t->Z[k], t->size + 1);
flint_free(t->Z);
acb_clear(t->z);
}

108
acb_dirichlet/roots_init.c Normal file
View file

@ -0,0 +1,108 @@
/*
Copyright (C) 2016 Pascal Molin
Copyright (C) 2016 Fredrik Johansson
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_roots_init(acb_dirichlet_roots_t t, ulong order, slong num, slong prec)
{
slong k, size, depth, best_depth, wp;
ulong reduced_order;
double cost, best_cost;
/* exploit 90 deg symmetries */
if (order % 4 == 0)
reduced_order = order / 8 + 1;
else if (order % 2 == 0)
reduced_order = order / 4 + 1;
else
reduced_order = order / 2 + 1;
wp = prec + 6 + 2 * FLINT_BIT_COUNT(reduced_order);
t->order = order;
t->reduced_order = reduced_order;
t->use_pow = 0;
if (reduced_order <= 2 || num <= 2)
{
depth = 0;
size = 0;
}
else
{
/* At depth = d and reduced_order = n, for k evaluations we need about
log_2(n) * k muls if d == 0
d * n^(1/d) + (d-1) * k muls if d > 0 */
best_cost = FLINT_BIT_COUNT(reduced_order) * (double) num;
best_depth = 0;
for (depth = 1; depth <= 4; depth++)
{
size = n_root(reduced_order, depth) + 1;
/* limit memory usage */
if (depth * _acb_vec_estimate_allocated_bytes(size, wp) > 1e9)
continue;
cost = depth * (double) size + (depth - 1) * (double) num;
if (cost < best_cost)
{
best_depth = depth;
best_cost = cost;
}
}
depth = best_depth;
size = n_root(reduced_order, depth) + 1;
}
t->size = size;
t->depth = depth;
acb_init(t->z);
if (depth != 0)
{
acb_struct * z;
acb_nth_root(t->z, order, wp);
z = t->z;
t->Z = flint_malloc(depth * sizeof(acb_ptr));
/* todo: at the last level, we could avoid computing
entries that will never be reached */
for (k = 0; k < depth; k++)
{
t->Z[k] = _acb_vec_init(size + 1);
_acb_vec_set_powers(t->Z[k], z, size + 1, wp);
z = t->Z[k] + size;
}
}
else
{
/* this tuning could be improved */
if (reduced_order < 30)
t->use_pow = 1;
else if (reduced_order < 100)
t->use_pow = (prec >= 512);
else if (reduced_order < 10000)
t->use_pow = (prec >= 4096);
else
t->use_pow = (prec >= 16384);
if (t->use_pow)
acb_nth_root(t->z, order, wp);
t->Z = NULL;
}
}

69
acb_dirichlet/test/t-roots.c Normal file
View file

@ -0,0 +1,69 @@
/*
Copyright (C) 2016 Fredrik Johansson
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"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("roots....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
{
acb_t x, y, z;
acb_dirichlet_roots_t roots;
ulong n, k;
slong prec;
slong iter2;
n = 1 + n_randint(state, 1000);
prec = 2 + n_randint(state, 200);
acb_init(x);
acb_init(y);
acb_init(z);
acb_dirichlet_roots_init(roots, n, n_randtest(state), prec);
acb_nth_root(y, n, prec);
for (iter2 = 0; iter2 <= FLINT_MIN(n, 20); iter2++)
{
k = n_randint(state, 2 * n);
acb_dirichlet_root(x, roots, k, prec);
acb_pow_ui(z, y, k, prec);
if (!acb_overlaps(x, z))
{
flint_printf("FAIL: overlap\n\n");
flint_printf("iter = %wd n = %wu k = %wu prec = %wd\n\n", iter, n, k, prec);
flint_printf("x = "); acb_printn(x, 30, 0); flint_printf("\n\n");
flint_printf("z = "); acb_printn(z, 30, 0); flint_printf("\n\n");
abort();
}
}
acb_dirichlet_roots_clear(roots);
acb_clear(x);
acb_clear(y);
acb_clear(z);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

12
acb_dirichlet/theta_arb.c
View file

@ -18,16 +18,16 @@ _acb_dirichlet_theta_arb_smallorder(acb_t res, const dirichlet_group_t G, const
ulong order;
ulong * a;
int parity;
acb_dirichlet_powers_t z;
acb_dirichlet_roots_t z;
parity = dirichlet_parity_char(G, chi);
order = dirichlet_order_char(G, chi);
a = flint_malloc(len * sizeof(ulong));
dirichlet_chi_vec_order(a, G, chi, order, len);
_acb_dirichlet_powers_init(z, order, 0, 0, prec);
acb_dirichlet_roots_init(z, order, 0, prec);
acb_dirichlet_qseries_arb_powers_smallorder(res, xt, parity, a, z, len, prec);
acb_dirichlet_powers_clear(z);
acb_dirichlet_roots_clear(z);
flint_free(a);
}
@ -52,7 +52,7 @@ void
_acb_dirichlet_theta_arb_naive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, const arb_t xt, slong len, slong prec)
{
ulong order, * a;
acb_dirichlet_powers_t z;
acb_dirichlet_roots_t z;
int parity;
parity = dirichlet_parity_char(G, chi);
@ -60,11 +60,11 @@ _acb_dirichlet_theta_arb_naive(acb_t res, const dirichlet_group_t G, const diric
a = flint_malloc(len * sizeof(ulong));
dirichlet_chi_vec_order(a, G, chi, order, len);
acb_dirichlet_powers_init(z, order, len, prec);
acb_dirichlet_roots_init(z, order, len, prec);
acb_dirichlet_qseries_arb_powers_naive(res, xt, parity, a, z, len, prec);
acb_dirichlet_powers_clear(z);
acb_dirichlet_roots_clear(z);
flint_free(a);
}

34
doc/source/acb_dirichlet.rst
View file

@ -22,6 +22,38 @@ The code in other modules for computing the Riemann zeta function,
Hurwitz zeta function and polylogarithm will possibly be migrated to this
module in the future.
Roots of unity
-------------------------------------------------------------------------------
.. type:: acb_dirichlet_roots_struct
.. type:: acb_dirichlet_roots_t
.. function:: void acb_dirichlet_roots_init(acb_dirichlet_roots_t roots, ulong n, slong num, slong prec)
Initializes *roots* with precomputed data for fast evaluation of roots of
unity `e^{2\pi i k/n}` of a fixed order *n*. The precomputation is
optimized for *num* evaluations.
For very small *num*, only the single root `e^{2\pi i/n}` will be
precomputed, which can then be raised to a power. For small *prec*
and large *n*, this method might even skip precomputing this single root
if it estimates that evaluating roots of unity from scratch will be faster
than powering.
If *num* is large enough, the whole set of roots in the first quadrant
will be precomputed at once. However, this is automatically avoided for
large *n* if too much memory would be used. For intermediate *num*,
baby-step giant-step tables are computed.
.. function:: void acb_dirichlet_roots_clear(acb_dirichlet_roots_t roots)
Clears the structure.
.. function:: void acb_dirichlet_root(acb_t res, const acb_dirichlet_roots_t roots, ulong k, slong prec)
Computes `e^{2\pi i k/n}`.
Truncated L-series and power sums
-------------------------------------------------------------------------------
@ -383,7 +415,7 @@ Hardy Z-functions
.. math ::
\theta(t) = -\frac{t}{2} \log(\pi/q) - \frac{\epsilon}{2}
\theta(t) = -\frac{t}{2} \log(\pi/q) - \frac{i \log(\epsilon)}{2}
+ \frac{\log \Gamma((s+\delta)/2) - \log \Gamma((1-s+\delta)/2)}{2i}
where `s = 1/2+it`, `\delta` is the parity of *chi*, and `\epsilon`