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

refactor and add errors to theta series

Browse source
This commit is contained in:
Pascal 2016-09-08 18:01:09 +02:00
parent f626fe53b5
commit 245b3f9442
10 changed files with 235 additions and 190 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

10
acb_dirichlet.h
View file

@ -234,14 +234,16 @@ void acb_dirichlet_chi(acb_t res, const acb_dirichlet_group_t G, const acb_diric
void acb_dirichlet_chi_vec(acb_ptr v, const acb_dirichlet_group_t G, const acb_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_eval_arb(acb_t res, acb_srcptr a, const arb_t x, slong len, 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(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_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);
void acb_dirichlet_arb_theta_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_arb_theta_smallorder(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec);
void mag_tail_kexpk2_arb(mag_t res, const arb_t a, ulong n);
void acb_dirichlet_chi_theta_arb(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t t, slong prec);
void _acb_dirichlet_theta_argument_at_arb(arb_t xt, ulong q, const arb_t t, slong prec);
void acb_dirichlet_theta_arb(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t t, slong prec);
void acb_dirichlet_ui_theta_arb(acb_t res, const acb_dirichlet_group_t G, ulong a, const arb_t t, slong prec);
void acb_dirichlet_gauss_sum_naive(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec);

50
acb_dirichlet/arb_theta_naive.c
View file

@ -1,50 +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_arb_theta_naive(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec)
{
slong k;
arb_t xk2, dx, x2;
acb_t zk;
arb_init(xk2);
arb_init(dx);
arb_init(x2);
acb_init(zk);
arb_set(dx, x);
arb_set(xk2, dx);
arb_mul(x2, dx, dx, prec);
acb_set_arb(res, xk2);
/* TODO: reduce prec */
for (k = 2; k < len; k++)
{
arb_mul(dx, dx, x2, prec);
arb_mul(xk2, xk2, dx, prec);
if (a[k] != ACB_DIRICHLET_CHI_NULL)
{
acb_dirichlet_power(zk, z, a[k], prec);
if (parity)
acb_mul_si(zk, zk, k, prec);
acb_addmul_arb(res, zk, xk2, prec);
}
}
arb_clear(xk2);
arb_clear(x2);
arb_clear(dx);
acb_clear(zk);
}

131
acb_dirichlet/chi_theta_arb.c
View file

@ -1,131 +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"
/* x(t) = exp(-Pi / q * t^2) */
static void
acb_dirichlet_arb_theta_xt(arb_t xt, ulong q, const arb_t t, slong prec)
{
arb_const_pi(xt, prec);
arb_div_ui(xt, xt, q, prec);
arb_mul(xt, xt, t, prec);
arb_mul(xt, xt, t, prec);
arb_neg(xt, xt);
arb_exp(xt, xt, prec);
}
void
acb_dirichlet_chi_theta_arb_smallorder(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t xt, slong len, slong prec)
{
ulong * a;
acb_dirichlet_powers_t z;
a = flint_malloc(len * sizeof(ulong));
acb_dirichlet_ui_chi_vec(a, G, chi, len);
_acb_dirichlet_powers_init(z, chi->order.n, 0, 0, prec);
acb_dirichlet_arb_theta_smallorder(res, xt, chi->parity, a, z, len, prec);
acb_dirichlet_powers_clear(z);
flint_free(a);
}
void
acb_dirichlet_chi_theta_arb_series(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t xt, slong len, slong prec)
{
acb_ptr a;
a = _acb_vec_init(len);
acb_dirichlet_chi_vec(a, G, chi, len, prec);
if (chi->parity)
{
slong k;
for (k = 2; k < len; k++)
acb_mul_si(a + k, a + k, k, prec);
}
acb_dirichlet_qseries_eval_arb(res, a, xt, len, prec);
_acb_vec_clear(a, len);
}
void
acb_dirichlet_vec_mellin_series(acb_ptr res, const acb_dirichlet_group_t G, const acb_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;
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);
arb_sqrt(st, t, prec);
arb_one(tk);
arb_one(stk);
for (k = 0; k < n; k++)
{
acb_dirichlet_arb_theta_xt(xt, G->q, tk, prec);
/* TODO: reduce len */
acb_dirichlet_qseries_eval_arb(res + k, a, xt, len, prec);
acb_mul_arb(res + k, res + k, stk, prec);
arb_mul(tk, tk, t, prec);
arb_mul(stk, stk, st, prec);
}
arb_clear(xt);
arb_clear(tk);
arb_clear(stk);
arb_clear(st);
_acb_vec_clear(a, len);
}
void
acb_dirichlet_chi_theta_arb_naive(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t xt, slong len, slong prec)
{
ulong * a;
acb_dirichlet_powers_t z;
a = flint_malloc(len * sizeof(ulong));
acb_dirichlet_ui_chi_vec(a, G, chi, len);
acb_dirichlet_powers_init(z, chi->order.n, len, prec);
acb_dirichlet_arb_theta_naive(res, xt, chi->parity, a, z, len, prec);
acb_dirichlet_powers_clear(z);
flint_free(a);
}
void
acb_dirichlet_chi_theta_arb(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t t, slong prec)
{
slong len;
arb_t xt;
len = acb_dirichlet_theta_length(G->q, t, prec);
arb_init(xt);
acb_dirichlet_arb_theta_xt(xt, G->q, t, prec);
/* TODO: tune this limit */
if (chi->order.n < 30)
acb_dirichlet_chi_theta_arb_smallorder(res, G, chi, xt, len, prec);
else
acb_dirichlet_chi_theta_arb_naive(res, G, chi, xt, len, prec);
arb_clear(xt);
}

2
acb_dirichlet/gauss_sum_theta.c
View file

@ -28,7 +28,7 @@ acb_dirichlet_gauss_sum_theta(acb_t res, const acb_dirichlet_group_t G, const ac
}
arb_one(x);
acb_dirichlet_chi_theta_arb(res, G, chi, x, prec);
acb_dirichlet_theta_arb(res, G, chi, x, prec);
acb_init(eps);
acb_conj(eps, res);
acb_div(eps, res, eps, prec);

4
acb_dirichlet/qseries_eval_arb.c → acb_dirichlet/qseries_arb.c
View file

@ -11,9 +11,9 @@
#include "acb_dirichlet.h"
/* assume a[0] = 0 */
/* sum a_k x^(k^2), assume a[0] = 0 */
void
acb_dirichlet_qseries_eval_arb(acb_t res, acb_srcptr a, const arb_t x, slong len, slong prec)
acb_dirichlet_qseries_arb(acb_t res, acb_srcptr a, const arb_t x, slong len, slong prec)
{
slong k;
arb_t xk2, dx, x2;

39
acb_dirichlet/arb_theta_smallorder.c → acb_dirichlet/qseries_arb_powers.c
View file

@ -12,9 +12,46 @@
#include "acb_dirichlet.h"
#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)
{
slong k;
arb_t xk2, dx, x2;
acb_t zk;
arb_init(xk2);
arb_init(dx);
arb_init(x2);
acb_init(zk);
arb_set(dx, x);
arb_set(xk2, dx);
arb_mul(x2, dx, dx, prec);
acb_set_arb(res, xk2);
/* TODO: reduce prec */
for (k = 2; k < len; k++)
{
arb_mul(dx, dx, x2, prec);
arb_mul(xk2, xk2, dx, prec);
if (a[k] != ACB_DIRICHLET_CHI_NULL)
{
acb_dirichlet_power(zk, z, a[k], prec);
if (parity)
acb_mul_si(zk, zk, k, prec);
acb_addmul_arb(res, zk, xk2, prec);
}
}
arb_clear(xk2);
arb_clear(x2);
arb_clear(dx);
acb_clear(zk);
}
/* small order, multiply by chi at the end */
void
acb_dirichlet_arb_theta_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(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec)
{
slong k;
ulong order = z->order;

93
acb_dirichlet/theta_arb.c Normal file
View file

@ -0,0 +1,93 @@
/*
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_theta_arb_smallorder(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t xt, slong len, slong prec)
{
ulong * a;
acb_dirichlet_powers_t z;
a = flint_malloc(len * sizeof(ulong));
acb_dirichlet_ui_chi_vec(a, G, chi, len);
_acb_dirichlet_powers_init(z, chi->order.n, 0, 0, prec);
acb_dirichlet_qseries_arb_powers(res, xt, chi->parity, a, z, len, prec);
acb_dirichlet_powers_clear(z);
flint_free(a);
}
void
_acb_dirichlet_theta_arb_series(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t xt, slong len, slong prec)
{
acb_ptr a;
a = _acb_vec_init(len);
acb_dirichlet_chi_vec(a, G, chi, len, prec);
if (chi->parity)
{
slong k;
for (k = 2; k < len; k++)
acb_mul_si(a + k, a + k, k, prec);
}
acb_dirichlet_qseries_arb(res, a, xt, len, prec);
_acb_vec_clear(a, len);
}
void
_acb_dirichlet_theta_arb_naive(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t xt, slong len, slong prec)
{
ulong * a;
acb_dirichlet_powers_t z;
a = flint_malloc(len * sizeof(ulong));
acb_dirichlet_ui_chi_vec(a, G, chi, len);
acb_dirichlet_powers_init(z, chi->order.n, len, prec);
acb_dirichlet_qseries_arb_powers_naive(res, xt, chi->parity, a, z, len, prec);
acb_dirichlet_powers_clear(z);
flint_free(a);
}
void
acb_dirichlet_theta_arb(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t t, slong prec)
{
slong len;
arb_t xt;
mag_t e;
len = acb_dirichlet_theta_length(G->q, t, prec);
arb_init(xt);
_acb_dirichlet_theta_argument_at_arb(xt, G->q, t, prec);
mag_init(e);
mag_tail_kexpk2_arb(e, xt, len);
arb_neg(xt, xt);
arb_exp(xt, xt, prec);
/* TODO: tune this limit */
if (chi->order.n < 30)
_acb_dirichlet_theta_arb_smallorder(res, G, chi, xt, len, prec);
else
_acb_dirichlet_theta_arb_naive(res, G, chi, xt, len, prec);
arb_add_error_mag(acb_realref(res), e);
arb_add_error_mag(acb_imagref(res), e);
mag_clear(e);
arb_clear(xt);
}

38
acb_dirichlet/theta_length.c
View file

@ -37,3 +37,41 @@ acb_dirichlet_theta_length(ulong q, const arb_t t, slong prec)
arf_clear(at);
return len;
}
/* bound for sum_{k>n} k*exp(-a k^2) */
void
mag_tail_kexpk2_arb(mag_t res, const arb_t a, ulong n)
{
mag_t m;
mag_init(m);
arb_get_mag_lower(m, a);
/* a < 1/4 ? */
if (mag_cmp_2exp_si(m, -2) <= 0)
{
mag_t c;
mag_init(c);
mag_mul_ui(c, m, 2);
mag_addmul(c, c, c);
mag_mul_ui(res, m, n*n-n+1);
mag_expinv(res, res);
mag_div(res, res, c);
mag_clear(c);
}
else
{
mag_mul_ui(res, m, n*n-n-1);
mag_expinv(res, res);
mag_mul_ui(res, res, 2);
}
mag_clear(m);
}
/* a(t) = Pi / q * t^2 */
void
_acb_dirichlet_theta_argument_at_arb(arb_t xt, ulong q, const arb_t t, slong prec)
{
arb_const_pi(xt, prec);
arb_div_ui(xt, xt, q, prec);
arb_mul(xt, xt, t, prec);
arb_mul(xt, xt, t, prec);
}

2
acb_dirichlet/ui_theta_arb.c
View file

@ -19,7 +19,7 @@ acb_dirichlet_ui_theta_arb(acb_t res, const acb_dirichlet_group_t G, ulong a, co
acb_dirichlet_char_init(chi, G);
acb_dirichlet_char(chi, G, a);
acb_dirichlet_chi_theta_arb(res, G, chi, t, prec);
acb_dirichlet_theta_arb(res, G, chi, t, prec);
acb_dirichlet_char_clear(chi);
}

56
acb_dirichlet/vec_mellin_arb.c Normal file
View file

@ -0,0 +1,56 @@
/*
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 acb_dirichlet_group_t G, const acb_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);
}