implement Hardy Z-function for Dirichlet L-functions

This commit is contained in:
Fredrik Johansson 2016-11-16 22:49:45 +01:00
parent 776acc62ab
commit 66741901ad
7 changed files with 415 additions and 23 deletions

View file

@ -126,6 +126,15 @@ void acb_dirichlet_l(acb_t res, const acb_t s, const dirichlet_group_t G, const
void acb_dirichlet_l_jet(acb_ptr res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, int deflate, slong len, slong prec);
void acb_dirichlet_hardy_theta(acb_ptr res, const acb_t t,
const dirichlet_group_t G, const dirichlet_char_t chi,
slong len, slong prec);
void acb_dirichlet_hardy_z(acb_t res, const acb_t t,
const dirichlet_group_t G, const dirichlet_char_t chi,
slong len, slong prec);
/* utils */
ACB_DIRICHLET_INLINE void

126
acb_dirichlet/hardy_theta.c Normal file
View file

@ -0,0 +1,126 @@
/*
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"
#include "acb_poly.h"
void
acb_dirichlet_hardy_theta(acb_ptr res, const acb_t t,
const dirichlet_group_t G, const dirichlet_char_t chi,
slong len, slong prec)
{
acb_struct y[2];
arb_t c;
slong k;
int parity;
ulong q;
if (len <= 0)
return;
if (t == res)
{
acb_init(y);
acb_set(y, t);
acb_dirichlet_hardy_theta(res, y, G, chi, len, prec);
acb_clear(y);
return;
}
if (G == NULL)
{
parity = 0;
q = 1;
}
else
{
parity = dirichlet_parity_char(G, chi);
q = G->q;
if (q != dirichlet_conductor_char(G, chi))
{
flint_printf("hardy theta: need primitive character\n");
abort();
}
}
if (!acb_is_finite(t))
{
_acb_vec_indeterminate(res, len);
return;
}
acb_init(y + 0);
acb_init(y + 1);
arb_init(c);
/* res = log gamma((s+parity)/2), s = 0.5+i(t+x) */
acb_mul_onei(y, t);
arb_set_d(c, 0.5 + parity);
arb_add(acb_realref(y), acb_realref(y), c, prec);
acb_mul_2exp_si(y, y, -1);
acb_onei(y + 1);
acb_mul_2exp_si(y + 1, y + 1, -1);
_acb_poly_lgamma_series(res, y, FLINT_MIN(len, 2), len, prec);
if (acb_is_real(t))
{
for (k = 0; k < len; k++)
{
arb_set(acb_realref(res + k), acb_imagref(res + k));
arb_zero(acb_imagref(res + k));
}
}
else
{
acb_ptr v = _acb_vec_init(len);
/* v = log gamma(((1-s)+parity)/2), s = 0.5+i(t+x) */
acb_div_onei(y, t);
arb_set_d(c, 0.5 + parity);
arb_add(acb_realref(y), acb_realref(y), c, prec);
acb_mul_2exp_si(y, y, -1);
acb_neg(y + 1, y + 1);
_acb_poly_lgamma_series(v, y, FLINT_MIN(len, 2), len, prec);
_acb_vec_sub(res, res, v, len, prec);
for (k = 0; k < len; k++)
{
acb_div_onei(res + k, res + k);
acb_mul_2exp_si(res + k, res + k, -1);
}
_acb_vec_clear(v, len);
}
/* (t+x) [-(1/2) log(pi/q)] */
arb_const_pi(c, prec);
arb_div_ui(c, c, q, prec);
arb_log(c, c, prec);
arb_mul_2exp_si(c, c, -1);
acb_submul_arb(res, t, c, prec);
if (len > 1)
acb_sub_arb(res + 1, res + 1, c, prec);
/* i log(eps) / 2 */
if (q > 1)
{
acb_dirichlet_root_number(y, G, chi, prec);
acb_arg(c, y, prec);
arb_mul_2exp_si(c, c, -1);
arb_sub(acb_realref(res), acb_realref(res), c, prec);
}
acb_clear(y + 0);
acb_clear(y + 1);
arb_clear(c);
}

62
acb_dirichlet/hardy_z.c Normal file
View file

@ -0,0 +1,62 @@
/*
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"
#include "acb_poly.h"
void
acb_dirichlet_hardy_z(acb_t res, const acb_t t,
const dirichlet_group_t G, const dirichlet_char_t chi,
slong len, slong prec)
{
acb_ptr v, w;
slong k;
int is_real;
if (len <= 0)
return;
v = _acb_vec_init(len);
w = _acb_vec_init(len);
is_real = acb_is_real(t);
/* v = exp(i theta(t+x)) */
acb_dirichlet_hardy_theta(v, t, G, chi, len, prec);
_acb_vec_scalar_mul_onei(v, v, len);
_acb_poly_exp_series(v, v, len, len, prec);
/* w = L(1/2 + i (t+x)) */
acb_one(w);
acb_mul_2exp_si(w, w, -1);
arb_sub(acb_realref(w), acb_realref(w), acb_imagref(t), prec);
arb_set(acb_imagref(w), acb_realref(t));
acb_dirichlet_l_jet(w, w, G, chi, 0, len, prec);
for (k = 0; k < len; k++)
{
if (k % 4 == 1)
acb_mul_onei(w + k, w + k);
else if (k % 4 == 2)
acb_neg(w + k, w + k);
else if (k % 4 == 3)
acb_div_onei(w + k, w + k);
}
_acb_poly_mullow(res, v, len, w, len, len, prec);
if (is_real)
for (k = 0; k < len; k++)
arb_zero(acb_imagref(res + k));
_acb_vec_clear(v, len);
_acb_vec_clear(w, len);
}

View file

@ -0,0 +1,122 @@
/*
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("hardy_z....");
fflush(stdout);
flint_randinit(state);
/* test self-consistency */
for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
{
acb_t s, s2;
dirichlet_group_t G;
dirichlet_char_t chi;
acb_ptr vec1, vec2;
slong len1, len2;
slong prec1, prec2;
ulong q, k;
slong i;
len1 = n_randint(state, 6);
len2 = n_randint(state, 6);
prec1 = 2 + n_randint(state, 100);
prec2 = 2 + n_randint(state, 100);
do {
q = 1 + n_randint(state, 30);
} while (q % 4 == 2);
dirichlet_group_init(G, q);
dirichlet_char_init(chi, G);
do {
k = n_randint(state, n_euler_phi(q));
dirichlet_char_index(chi, G, k);
} while (dirichlet_conductor_char(G, chi) != q);
acb_init(s);
acb_init(s2);
vec1 = _acb_vec_init(len1);
vec2 = _acb_vec_init(len2);
acb_randtest(s, state, 2 + n_randint(state, 200), 2);
acb_randtest(s2, state, 2 + n_randint(state, 200), 2);
acb_sub(s2, s2, s2, 200);
acb_add(s2, s, s2, 200);
acb_dirichlet_hardy_z(vec1, s, G, chi, len1, prec1);
acb_dirichlet_hardy_z(vec2, s2, G, chi, len2, prec2);
for (i = 0; i < FLINT_MIN(len1, len2); i++)
{
if (!acb_overlaps(vec1 + i, vec2 + i))
{
flint_printf("FAIL: overlap\n\n");
flint_printf("iter = %wd q = %wu k = %wu i = %wd\n\n", iter, q, k, i);
flint_printf("s = "); acb_printn(s, 50, 0); flint_printf("\n\n");
flint_printf("r1 = "); acb_printn(vec1 + i, 50, 0); flint_printf("\n\n");
flint_printf("r2 = "); acb_printn(vec2 + i, 50, 0); flint_printf("\n\n");
abort();
}
}
if (arb_contains_zero(acb_imagref(s)))
{
for (i = 0; i < len1; i++)
{
if (!arb_contains_zero(acb_imagref(vec1 + i)))
{
flint_printf("FAIL: real 1\n\n");
flint_printf("iter = %wd q = %wu k = %wu i = %wd\n\n", iter, q, k, i);
flint_printf("s = "); acb_printn(s, 50, 0); flint_printf("\n\n");
flint_printf("r1 = "); acb_printn(vec1 + i, 50, 0); flint_printf("\n\n");
abort();
}
}
}
if (arb_contains_zero(acb_imagref(s2)))
{
for (i = 0; i < len2; i++)
{
if (!arb_contains_zero(acb_imagref(vec2 + i)))
{
flint_printf("FAIL: real 1\n\n");
flint_printf("iter = %wd q = %wu k = %wu i = %wd\n\n", iter, q, k, i);
flint_printf("s = "); acb_printn(s, 50, 0); flint_printf("\n\n");
flint_printf("r1 = "); acb_printn(vec2 + i, 50, 0); flint_printf("\n\n");
abort();
}
}
}
dirichlet_char_clear(chi);
dirichlet_group_clear(G);
acb_clear(s);
acb_clear(s2);
_acb_vec_clear(vec1, len1);
_acb_vec_clear(vec2, len2);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

View file

@ -89,7 +89,7 @@ The coefficients `C_k(p)` in the asymptotic part of the expansion
are expressed in terms of certain auxiliary coefficients `d_j^{(k)}`
and `F^{(j)}(p)`.
Because of artificial discontinuities, *s* should be exact inside
the evaluation (automatic reduction to the exact case is not yet implemented).
the evaluation.
.. function:: void acb_dirichlet_zeta_rs_f_coeffs(acb_ptr f, const arb_t p, slong n, slong prec)
@ -373,3 +373,26 @@ L-functions
gives the regular part of the Laurent expansion.
When *chi* is non-principal, *deflate* has no effect.
Hardy Z-functions
-------------------------------------------------------------------------------
.. function:: void acb_dirichlet_hardy_theta(acb_ptr res, const acb_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)
Computes the phase function used to construct the Z-function.
We have
.. math ::
\theta(t) = -\frac{t}{2} \log(\pi/q) - \frac{\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`
is the root number as computed by :func:`acb_dirichlet_root_number`.
The first *len* terms in the Taylor expansion are written to the output.
.. function:: void acb_dirichlet_hardy_z(acb_t res, const acb_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)
Computes the Hardy Z-function, also known as the Riemann-Siegel Z-function
`Z(t) = e^{i \theta(t)} L(1/2+it)`, which is real-valued for real *t*.
The first *len* terms in the Taylor expansion are written to the output.

View file

@ -243,7 +243,7 @@ encoding the zero value.
.. function:: ulong dirichlet_chi(const dirichlet_group_t G, const dirichlet_char_t chi, ulong n)
Compute that value `\chi(n)` as the exponent modulo *G->expo*.
Compute the value `\chi(n)` as the exponent modulo *G->expo*.
.. function:: void dirichlet_chi_vec(ulong * v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv)
@ -264,7 +264,7 @@ Character operations
.. function:: void dirichlet_char_pow(dirichlet_char_t c, const dirichlet_group_t G, const dirichlet_char_t a, ulong n)
Take the power of some character.
Take the power of a character.
.. function:: void dirichlet_char_lift(dirichlet_char_t chi_G, const dirichlet_group_t G, const dirichlet_char_t chi_H, const dirichlet_group_t H)
@ -278,3 +278,4 @@ Character operations
This requires `c(\chi_G)\mid q_H\mid q_G`, where `c(\chi_G)` is the
conductor of `\chi_G` and `q_G, q_H` are the moduli of G and H.

View file

@ -3,25 +3,51 @@
#include <string.h>
#include "arb_calc.h"
#include "acb_hypgeom.h"
#include "acb_dirichlet.h"
#include "flint/profiler.h"
slong eval_count = 0;
typedef struct
{
dirichlet_group_t * G;
dirichlet_char_t * chi;
}
z_param_struct;
int
z_function(arb_ptr out, const arb_t inp, void * params, slong order, slong prec)
{
arb_struct x[2];
z_param_struct * par = params;
arb_init(x);
arb_init(x + 1);
if (par->G == NULL)
{
arb_struct x[2];
arb_set(x, inp);
arb_one(x + 1);
arb_init(x);
arb_init(x + 1);
_arb_poly_riemann_siegel_z_series(out, x, FLINT_MIN(2, order), order, prec);
arb_set(x, inp);
arb_one(x + 1);
arb_clear(x);
arb_clear(x + 1);
_arb_poly_riemann_siegel_z_series(out, x, FLINT_MIN(2, order), order, prec);
arb_clear(x);
arb_clear(x + 1);
}
else
{
acb_ptr tmp;
slong k;
tmp = _acb_vec_init(order);
acb_set_arb(tmp, inp);
acb_dirichlet_hardy_z(tmp, tmp, *(par->G), *(par->chi), order, prec);
for (k = 0; k < order; k++)
arb_set(out + k, acb_realref(tmp + k));
_acb_vec_clear(tmp, order);
}
eval_count++;
return 0;
@ -152,7 +178,11 @@ int main(int argc, char *argv[])
arf_interval_ptr blocks;
arb_calc_func_t function;
int * info;
int params;
void * params;
int param1;
z_param_struct param2;
dirichlet_group_t G;
dirichlet_char_t chi;
slong digits, low_prec, high_prec, i, num, found_roots, found_unknown;
slong maxdepth, maxeval, maxfound;
int refine;
@ -166,7 +196,7 @@ int main(int argc, char *argv[])
flint_printf("real_roots function a b [-refine d] [-verbose] "
"[-maxdepth n] [-maxeval n] [-maxfound n] [-prec n]\n");
flint_printf("available functions:\n");
flint_printf(" 0 Z(x), Riemann-Siegel Z-function\n");
flint_printf(" 0 Z(x), Z-function (Riemann zeta or Dirichlet L-function)\n");
flint_printf(" 1 sin(x)\n");
flint_printf(" 2 sin(x^2)\n");
flint_printf(" 3 sin(1/x)\n");
@ -174,15 +204,20 @@ int main(int argc, char *argv[])
flint_printf(" 5 Ai'(x), Airy function\n");
flint_printf(" 6 Bi(x), Airy function\n");
flint_printf(" 7 Bi'(x), Airy function\n");
flint_printf("With 0, specify optional Dirichlet character with [-character q n]\n");
return 1;
}
params = 0;
param1 = 0;
param2.G = NULL;
param2.chi = NULL;
params = &param1;
switch (atoi(argv[1]))
{
case 0:
function = z_function;
params = &param2;
break;
case 1:
function = sin_x;
@ -195,19 +230,19 @@ int main(int argc, char *argv[])
break;
case 4:
function = airy;
params = 0;
param1 = 0;
break;
case 5:
function = airy;
params = 1;
param1 = 1;
break;
case 6:
function = airy;
params = 2;
param1 = 2;
break;
case 7:
function = airy;
params = 3;
param1 = 3;
break;
default:
flint_printf("require a function 0-7\n");
@ -257,6 +292,14 @@ int main(int argc, char *argv[])
{
low_prec = atol(argv[i+1]);
}
else if (!strcmp(argv[i], "-character"))
{
dirichlet_group_init(G, atol(argv[i+1]));
dirichlet_char_init(chi, G);
dirichlet_char_log(chi, G, atol(argv[i+2]));
param2.G = &G;
param2.chi = &chi;
}
}
high_prec = digits * 3.32192809488736 + 10;
@ -280,7 +323,7 @@ int main(int argc, char *argv[])
TIMEIT_ONCE_START
num = arb_calc_isolate_roots(&blocks, &info, function,
&params, interval, maxdepth, maxeval, maxfound, low_prec);
params, interval, maxdepth, maxeval, maxfound, low_prec);
for (i = 0; i < num; i++)
{
@ -302,7 +345,7 @@ int main(int argc, char *argv[])
continue;
if (arb_calc_refine_root_bisect(t,
function, &params, blocks + i, 5, low_prec)
function, params, blocks + i, 5, low_prec)
!= ARB_CALC_SUCCESS)
{
flint_printf("warning: some bisection steps failed!\n");
@ -316,7 +359,7 @@ int main(int argc, char *argv[])
}
if (arb_calc_refine_root_bisect(blocks + i,
function, &params, t, 5, low_prec)
function, params, t, 5, low_prec)
!= ARB_CALC_SUCCESS)
{
flint_printf("warning: some bisection steps failed!\n");
@ -330,10 +373,10 @@ int main(int argc, char *argv[])
}
arf_interval_get_arb(v, t, high_prec);
arb_calc_newton_conv_factor(C, function, &params, v, low_prec);
arb_calc_newton_conv_factor(C, function, params, v, low_prec);
arf_interval_get_arb(w, blocks + i, high_prec);
if (arb_calc_refine_root_newton(z, function, &params,
if (arb_calc_refine_root_newton(z, function, params,
w, v, C, 10, high_prec) != ARB_CALC_SUCCESS)
{
flint_printf("warning: some newton steps failed!\n");
@ -357,6 +400,12 @@ int main(int argc, char *argv[])
flint_free(blocks);
flint_free(info);
if (param2.G != NULL)
{
dirichlet_group_clear(G);
dirichlet_char_clear(chi);
}
arf_interval_clear(t);
arf_interval_clear(interval);
arf_clear(C);