From 986f7fab2faf0e3e75ab50370f7239eaf78bebb8 Mon Sep 17 00:00:00 2001 From: Pascal Date: Wed, 23 Mar 2016 12:53:08 +0100 Subject: [PATCH] fix bug on char_next + add log precomputations --- acb_dirichlet.h | 9 ++ acb_dirichlet/arb_quadratic_powers.c | 51 +++++++++ acb_dirichlet/char_conductor.c | 56 ++++++++++ acb_dirichlet/char_is_odd.c | 40 +++++++ acb_dirichlet/char_next.c | 3 +- acb_dirichlet/char_next_primitive.c | 9 +- acb_dirichlet/chi_vec_loop.c | 44 ++++---- acb_dirichlet/conrey_exp.c | 3 +- acb_dirichlet/conrey_first_primitive.c | 20 ++-- acb_dirichlet/conrey_log.c | 31 ++++-- acb_dirichlet/conrey_mul.c | 4 +- acb_dirichlet/conrey_next.c | 24 ++--- acb_dirichlet/conrey_next_primitive.c | 9 +- acb_dirichlet/conrey_one.c | 5 +- acb_dirichlet/conrey_print.c | 5 +- acb_dirichlet/group_dlog_precompute.c | 42 ++++++++ acb_dirichlet/group_init.c | 3 + acb_dirichlet/test/t-thetanull.c | 143 +++++++++++++++++++++++++ acb_dirichlet/zeta.c | 36 +++++++ dlog.h | 1 + 20 files changed, 468 insertions(+), 70 deletions(-) create mode 100644 acb_dirichlet/arb_quadratic_powers.c create mode 100644 acb_dirichlet/char_conductor.c create mode 100644 acb_dirichlet/char_is_odd.c create mode 100644 acb_dirichlet/group_dlog_precompute.c create mode 100644 acb_dirichlet/test/t-thetanull.c create mode 100644 acb_dirichlet/zeta.c diff --git a/acb_dirichlet.h b/acb_dirichlet.h index 993d240c..3279cd39 100644 --- a/acb_dirichlet.h +++ b/acb_dirichlet.h @@ -23,6 +23,7 @@ extern "C" { typedef struct { ulong q; /* modulus */ + nmod_t mod; /* modulus with precomputed inverse */ ulong q_even; /* even part of modulus */ ulong q_odd; /* odd part of modulus */ ulong phi_q; /* phi(q) = group size */ @@ -35,6 +36,7 @@ typedef struct ulong * generators; /* generator for each prime p[k] lifted mod q */ ulong * phi; /* phi(k) = phi(p[k]^e[k]) */ ulong * PHI; /* PHI(k) = expo / phi(k) */ + dlog_precomp_t * dlog; /* precomputed data for discrete log mod p^e */ } acb_dirichlet_group_struct; @@ -101,6 +103,9 @@ void acb_dirichlet_char_conrey(acb_dirichlet_char_t chi, const acb_dirichlet_gro void acb_dirichlet_char_normalize(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G); void acb_dirichlet_char_denormalize(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G); +ulong acb_dirichlet_char_conductor(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi); +int acb_dirichlet_char_is_odd(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi); + void acb_dirichlet_char_one(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G); void acb_dirichlet_char_first_primitive(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G); @@ -115,6 +120,10 @@ void acb_dirichlet_chi_vec(ulong *v, ulong nv, const acb_dirichlet_group_t G, co void acb_dirichlet_char_vec(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n, slong prec); + +void acb_dirichlet_zeta(acb_t res, ulong order, slong prec); +void acb_dirichlet_arb_quadratic_powers(arb_ptr v, slong nv, const arb_t x, slong prec); + #ifdef __cplusplus } #endif diff --git a/acb_dirichlet/arb_quadratic_powers.c b/acb_dirichlet/arb_quadratic_powers.c new file mode 100644 index 00000000..9e2f1bee --- /dev/null +++ b/acb_dirichlet/arb_quadratic_powers.c @@ -0,0 +1,51 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2016 Pascal Molin + +******************************************************************************/ + +#include "acb_dirichlet.h" + +void +acb_dirichlet_arb_quadratic_powers(arb_ptr v, slong nv, const arb_t x, slong prec) +{ + slong i; + arb_t dx, x2; + arb_init(dx); + arb_init(x2); + arb_set(dx, x); + arb_mul(x2, x, x, prec); + for (i = 0; i < nv; i++) + { + if (i == 0) + arb_one(v + i); + else if (i == 1) + arb_set_round(v + i, x, prec); + else + { + arb_mul(dx, dx, x2, prec); + arb_mul(v + i, v + i - 1, dx, prec); + } + } + arb_clear(dx); + arb_clear(x2); +} diff --git a/acb_dirichlet/char_conductor.c b/acb_dirichlet/char_conductor.c new file mode 100644 index 00000000..62e96f0a --- /dev/null +++ b/acb_dirichlet/char_conductor.c @@ -0,0 +1,56 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2016 Pascal Molin + +******************************************************************************/ + +#include "acb_dirichlet.h" + +ulong +acb_dirichlet_char_conductor(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi) +{ + int k, f; + ulong cond = 1; + if (G->neven >= 1 && chi->expo[0] == 1) + cond = 4; + if (G->neven == 2 && chi->expo[1]) + { + ulong l2 = chi->expo[1]; + f = n_remove(&l2, 2); + cond = G->primepowers[1] >> f; + } + for (k = G->neven; k < G->neven; k++) + { + if (chi->expo[k]) + { + ulong p, lp; + p = G->primes[k]; + lp = chi->expo[k]; + f = n_remove(&lp, p); + if (f) + cond *= n_pow(p, G->exponents[k] - f); + else + cond *= G->primepowers[k]; + } + } + return cond; +} diff --git a/acb_dirichlet/char_is_odd.c b/acb_dirichlet/char_is_odd.c new file mode 100644 index 00000000..cb5d0502 --- /dev/null +++ b/acb_dirichlet/char_is_odd.c @@ -0,0 +1,40 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2016 Pascal Molin + +******************************************************************************/ + +#include "acb_dirichlet.h" + +int +acb_dirichlet_char_is_odd(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi) +{ + slong k, odd = 0; + for (k = 0; k < G->num; k++) + { + if (k == 1 && G->neven == 2) + continue; + if (chi->expo[k] % 2 == 1) + odd = 1 - odd; + } + return odd; +} diff --git a/acb_dirichlet/char_next.c b/acb_dirichlet/char_next.c index 0ab819da..5b7efb16 100644 --- a/acb_dirichlet/char_next.c +++ b/acb_dirichlet/char_next.c @@ -35,8 +35,7 @@ acb_dirichlet_char_next(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G) /* update index */ for (k=0; k < G->num ; k++) { - /* chi->n = n_mulmod(chi->n, G->generators[k], G->q); */ - chi->n = chi->n * G->generators[k] % G->q; + chi->n = nmod_mul(chi->n, G->generators[k], G->mod); chi->expo[k] += G->PHI[k]; if (chi->expo[k] < G->expo) break; diff --git a/acb_dirichlet/char_next_primitive.c b/acb_dirichlet/char_next_primitive.c index eb3b5b9d..cd3bda76 100644 --- a/acb_dirichlet/char_next_primitive.c +++ b/acb_dirichlet/char_next_primitive.c @@ -36,8 +36,7 @@ acb_dirichlet_char_next_primitive(acb_dirichlet_char_t chi, const acb_dirichlet_ k = 0; if (G->neven == 2) { - /* chi->n = n_mulmod(chi->n, G->generators[0], G->q); */ - chi->n = chi->n * G->generators[0] % G->q; + chi->n = nmod_mul(chi->n, G->generators[0], G->mod); if (++chi->expo[0] < G->expo) return 0; chi->expo[0] = 0; @@ -45,13 +44,11 @@ acb_dirichlet_char_next_primitive(acb_dirichlet_char_t chi, const acb_dirichlet_ } for (; k < G->num ; k++) { - /* chi->n = n_mulmod(chi->n, G->generators[k], G->q); */ - chi->n = chi->n * G->generators[k] % G->q; + chi->n = nmod_mul(chi->n, G->generators[k], G->mod); chi->expo[k] += G->PHI[k]; if (chi->expo[k] % G->primes[k] == 0) { - /* chi->n = n_mulmod(chi->n, G->generators[k], G->q); */ - chi->n = chi->n * G->generators[k] % G->q; + chi->n = nmod_mul(chi->n, G->generators[k], G->mod); chi->expo[k] += G->PHI[k]; } if (chi->expo[k] < G->expo) diff --git a/acb_dirichlet/chi_vec_loop.c b/acb_dirichlet/chi_vec_loop.c index 40d1f950..e2dec1db 100644 --- a/acb_dirichlet/chi_vec_loop.c +++ b/acb_dirichlet/chi_vec_loop.c @@ -29,24 +29,28 @@ void acb_dirichlet_chi_vec_loop(ulong *v, ulong nv, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi) { - int j; - ulong t, k; - acb_dirichlet_conrey_t x; - acb_dirichlet_conrey_init(x, G); - acb_dirichlet_conrey_one(x, G); - t = v[1] = 0; - while ( (j = acb_dirichlet_conrey_next(x, G)) < G->num ) - { - /* exponents were modified up to j */ - for (k = 0; k < j; k++) - t = (t + chi->expo[k] * x->log[k]) % chi->order; - if (x->n < nv) - v[x->n] = t; - } - /* fix result outside primes */ - acb_dirichlet_vec_set_null(v, nv, G); - /* copy outside modulus */ - for (k = G->q + 1; k < nv ; k++ ) - v[k] = v[k - G->q]; - acb_dirichlet_conrey_clear(x); + int j; + ulong t, k; + acb_dirichlet_conrey_t x; + acb_dirichlet_conrey_init(x, G); + acb_dirichlet_conrey_one(x, G); + + for (k = 0; k < nv; k++) + v[k] = CHI_NULL; + + t = v[1] = 0; + while ( (j = acb_dirichlet_conrey_next(x, G)) < G->num ) + { + /* exponents were modified up to j */ + for (k = 0; k <= j; k++) + t = (t + chi->expo[k]) % chi->order; + if (x->n < nv) + v[x->n] = t; + } + /* fix result outside primes */ + /*acb_dirichlet_vec_set_null(v, nv, G);*/ + /* copy outside modulus */ + for (k = G->q; k < nv ; k++ ) + v[k] = v[k - G->q]; + acb_dirichlet_conrey_clear(x); } diff --git a/acb_dirichlet/conrey_exp.c b/acb_dirichlet/conrey_exp.c index 3dd6c7b1..af268814 100644 --- a/acb_dirichlet/conrey_exp.c +++ b/acb_dirichlet/conrey_exp.c @@ -30,8 +30,7 @@ acb_dirichlet_conrey_exp(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G { ulong k, n = 1; for (k = G->neven; k < G->num; k++) - /* n = n_mulmod(n, n_powmod(G->generators[k], x->log[k], G->q), G->q); */ - n = n * n_powmod(G->generators[k], x->log[k], G->q) % G->q; + n = nmod_mul(n, nmod_pow_ui(G->generators[k], x->log[k], G->mod), G->mod); x->n = n; return n; } diff --git a/acb_dirichlet/conrey_first_primitive.c b/acb_dirichlet/conrey_first_primitive.c index 488a6186..50913511 100644 --- a/acb_dirichlet/conrey_first_primitive.c +++ b/acb_dirichlet/conrey_first_primitive.c @@ -26,15 +26,23 @@ #include "acb_dirichlet.h" void -acb_dirichlet_conrey_first_primitive(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G) { +acb_dirichlet_conrey_first_primitive(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G) +{ ulong k; if (G->q % 4 == 2) { - flint_printf("Exception (acb_dirichlet_conrey_first_primitive). no primitive element mod %wu.\n",G->q); - abort(); + flint_printf("Exception (acb_dirichlet_conrey_first_primitive). No primitive element mod %wu.\n",G->q); + abort(); } + x->n = 1; for (k = 0; k < G->num ; k++) - x->log[k] = 1; - if (G->neven == 2) - x->log[0] = 0; + { + if (k == 0 && G->neven == 2) + x->log[k] = 0; + else + { + x->n = nmod_mul(x->n, G->generators[k], G->mod); + x->log[k] = 1; + } + } } diff --git a/acb_dirichlet/conrey_log.c b/acb_dirichlet/conrey_log.c index bc41fb02..a7fd1401 100644 --- a/acb_dirichlet/conrey_log.c +++ b/acb_dirichlet/conrey_log.c @@ -35,20 +35,31 @@ acb_dirichlet_conrey_log(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G ulong k, pk, gk; /* even part */ if (G->neven >= 1) - x->log[0] = (m % 4 == 3); - if (G->neven == 2) { - ulong q_even = G->q_even; - ulong g2 = 5; - ulong m2 = (m % 4 == 3) ? -m % q_even : m % q_even; - x->log[1] = n_discrete_log_bsgs(m2, g2, q_even); + x->log[0] = (m % 4 == 3); + if (G->neven == 2) + { + ulong m2 = (x->log[0]) ? -m % G->q_even : m % G->q_even; + if (G->dlog == NULL) + x->log[1] = n_discrete_log_bsgs(m2, 5, G->q_even); + else + x->log[1] = dlog_precomp(G->dlog[1], m2); + } } /* odd part */ - for (k = G->neven; k < G->num; k++) + if (G->dlog == NULL) { - pk = G->primepowers[k]; - gk = G->generators[k] % pk; - x->log[k] = n_discrete_log_bsgs(m % pk, gk, pk); + for (k = G->neven; k < G->num; k++) + { + pk = G->primepowers[k]; + gk = G->generators[k] % pk; + x->log[k] = n_discrete_log_bsgs(m % pk, gk, pk); + } + } + else + { + for (k = G->neven; k < G->num; k++) + x->log[k] = dlog_precomp(G->dlog[k], m % G->primepowers[k]); } /* keep value m */ x->n = m; diff --git a/acb_dirichlet/conrey_mul.c b/acb_dirichlet/conrey_mul.c index c0c7d95c..cdac5acd 100644 --- a/acb_dirichlet/conrey_mul.c +++ b/acb_dirichlet/conrey_mul.c @@ -30,6 +30,6 @@ acb_dirichlet_conrey_mul(acb_dirichlet_conrey_t c, const acb_dirichlet_group_t G { ulong k; for (k = 0; k < G->num ; k++) - c->log[k] = a->log[k] + b->log[k] % G->phi[k]; - c->n = a->n * b->n % G->q; + c->log[k] = (a->log[k] + b->log[k]) % G->phi[k]; + c->n = nmod_mul(a->n, b->n, G->mod); } diff --git a/acb_dirichlet/conrey_next.c b/acb_dirichlet/conrey_next.c index ab5b6ecf..5b4be26c 100644 --- a/acb_dirichlet/conrey_next.c +++ b/acb_dirichlet/conrey_next.c @@ -28,16 +28,16 @@ int acb_dirichlet_conrey_next(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G) { - /* update index */ - ulong k; - for (k=0; k < G->num ; k++) - { - /* x->n = n_mulmod(x->n, G->generators[k], G->q); */ - x->n = x->n * G->generators[k] % G->q; - if (++x->log[k] < G->phi[k]) - break; - x->log[k] = 0; - } - /* return last index modified */ - return k; + /* update index */ + int k; + for (k=0; k < G->num; k++) + { + x->n = nmod_mul(x->n, G->generators[k], G->mod); + x->log[k] += 1; + if (x->log[k] < G->phi[k]) + break; + x->log[k] = 0; + } + /* return last index modified */ + return k; } diff --git a/acb_dirichlet/conrey_next_primitive.c b/acb_dirichlet/conrey_next_primitive.c index acae8a76..90671bd6 100644 --- a/acb_dirichlet/conrey_next_primitive.c +++ b/acb_dirichlet/conrey_next_primitive.c @@ -33,8 +33,7 @@ acb_dirichlet_conrey_next_primitive(acb_dirichlet_conrey_t x, const acb_dirichle ulong k = 0; if (G->neven == 2) { - /* x->n = n_mulmod(x->n, G->generators[0], G->q); */ - x->n = x->n * G->generators[0] % G->q; + x->n = nmod_mul(x->n, G->generators[0], G->mod); if (++x->log[0] == 1) return 0; x->log[0] = 0; @@ -42,12 +41,10 @@ acb_dirichlet_conrey_next_primitive(acb_dirichlet_conrey_t x, const acb_dirichle } for (; k < G->num ; k++) { - /* x->n = n_mulmod(x->n, G->generators[k], G->q); */ - x->n = x->n * G->generators[k] % G->q; + x->n = nmod_mul(x->n, G->generators[k], G->mod); if (++x->log[k] % G->primes[k] == 0) { - /* x->n = n_mulmod(x->n, G->generators[k], G->q); */ - x->n = x->n * G->generators[k] % G->q; + x->n = nmod_mul(x->n, G->generators[k], G->mod); ++x->log[k]; } if (x->log[k] < G->phi[k]) diff --git a/acb_dirichlet/conrey_one.c b/acb_dirichlet/conrey_one.c index 36e23b58..d5eac769 100644 --- a/acb_dirichlet/conrey_one.c +++ b/acb_dirichlet/conrey_one.c @@ -26,9 +26,10 @@ #include "acb_dirichlet.h" void -acb_dirichlet_conrey_one(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G) { +acb_dirichlet_conrey_one(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G) +{ ulong k; for (k = 0; k < G->num ; k++) - x->log[k] = 0; + x->log[k] = 0; x->n = 1; } diff --git a/acb_dirichlet/conrey_print.c b/acb_dirichlet/conrey_print.c index 266bbf23..92aeee2d 100644 --- a/acb_dirichlet/conrey_print.c +++ b/acb_dirichlet/conrey_print.c @@ -29,8 +29,9 @@ void acb_dirichlet_conrey_print(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x) { int k; - flint_printf("[%wu",x->log[0]); + if (G->num) + flint_printf("[%wu", x->log[0]); for (k = 1; k < G->num; k++) - flint_printf(", %wu",x->log[k]); + flint_printf(", %wu", x->log[k]); flint_printf("]\n"); } diff --git a/acb_dirichlet/group_dlog_precompute.c b/acb_dirichlet/group_dlog_precompute.c new file mode 100644 index 00000000..f16bd567 --- /dev/null +++ b/acb_dirichlet/group_dlog_precompute.c @@ -0,0 +1,42 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2016 Pascal Molin + +******************************************************************************/ + +#include "acb_dirichlet.h" + +void +acb_dirichlet_group_dlog_precompute(acb_dirichlet_group_t G, ulong num) +{ + slong k; + G->dlog = flint_malloc(G->num * sizeof(dlog_precomp_t)); + for (k = G->neven; k < G->num; k++) + { + ulong p, e, pe, a; + p = G->primes[k]; + e = G->exponents[k]; + pe = G->primepowers[k]; + a = G->generators[k] % pe; + dlog_precomp_modpe_init(G->dlog[k], a, p, e, pe, num); + } +} diff --git a/acb_dirichlet/group_init.c b/acb_dirichlet/group_init.c index eba5b84a..ef144b60 100644 --- a/acb_dirichlet/group_init.c +++ b/acb_dirichlet/group_init.c @@ -41,6 +41,8 @@ acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q) n_factor_t fac; G->q = q; + nmod_init(&G->mod, q); + G->q_odd = q; G->q_even = 1; @@ -63,6 +65,7 @@ acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q) G->generators = flint_malloc(G->num * sizeof(ulong)); G->phi = flint_malloc(G->num * sizeof(ulong)); G->PHI = flint_malloc(G->num * sizeof(ulong)); + G->dlog = NULL; /* even part */ G->expo = G->phi_q = 1; diff --git a/acb_dirichlet/test/t-thetanull.c b/acb_dirichlet/test/t-thetanull.c new file mode 100644 index 00000000..dc7ca96a --- /dev/null +++ b/acb_dirichlet/test/t-thetanull.c @@ -0,0 +1,143 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2016 Pascal Molin + +******************************************************************************/ + +#include "acb_dirichlet.h" +#include + +#define PI 3.14159265358 +#define LOG2 0.69314718055 + +int main() +{ + slong prec = 80; + ulong q; + + flint_printf("theta null...."); + fflush(stdout); + + /* check the only theta functions + * theta(chi) = sum chi(k)* k^odd * exp(-Pi * k^2 / q) + * vanishing at 1 correspond to two specific + * characters of moduli 300 and 600 + conjugates + */ + + for (q = 3; q < 800; q ++) + { + acb_dirichlet_group_t G; + acb_dirichlet_conrey_t x; + acb_dirichlet_char_t chi; + ulong * v, nv, k; + + acb_t zeta, sum; + acb_ptr z; + + arb_t eq; + arb_ptr t, kt, tt; + + if (q % 4 == 2) + continue; + + acb_dirichlet_group_init(G, q); + acb_dirichlet_conrey_init(x, G); + acb_dirichlet_char_init(chi, G); + + acb_init(zeta); + acb_dirichlet_zeta(zeta, G->expo, prec); + + z = _acb_vec_init(G->expo); + _acb_vec_set_powers(z, zeta, G->expo, prec); + + nv = ceil(sqrt((double)q * prec * LOG2 / PI)) + 2; + v = flint_malloc(nv * sizeof(ulong)); + + arb_init(eq); + arb_const_pi(eq, prec); + arb_div_ui(eq, eq, q, prec); + arb_neg(eq, eq); + arb_exp(eq, eq, prec); + + t = _arb_vec_init(nv); + acb_dirichlet_arb_quadratic_powers(t, nv, eq, prec); + + kt = _arb_vec_init(nv); + for (k = 1; k < nv; k++) + arb_mul_ui(kt + k, t + k, k, prec); + + /* theta function on primitive characters */ + acb_init(sum); + acb_dirichlet_conrey_first_primitive(x, G); + + while (1) { + ulong m; + acb_zero(sum); + + acb_dirichlet_char_conrey(chi, G, x); + acb_dirichlet_chi_vec_loop(v, nv, G, chi); + + m = G->expo / chi->order; + /* + flint_printf("Theta(chi_%wu(%wu)) (m=%wu)\n", q, chi->n, m); + */ + tt = acb_dirichlet_char_is_odd(G, chi) ? kt : t; + for (k = 1; k < nv; k++) + if (v[k] != CHI_NULL) + acb_addmul_arb(sum, z + (v[k] * m), tt + k, prec); + + if ((q == 300 && (chi->n == 271 || chi->n == 131)) + || (q == 600 && (chi->n == 11 || chi->n == 91))) + { + if (!acb_contains_zero(sum)) + { + flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->n); + acb_printd(sum, 10); + flint_printf("\n"); + abort(); + } + } + else if (acb_contains_zero(sum)) + { + flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->n); + acb_printd(sum, 10); + flint_printf("\n"); + abort(); + } + + if (acb_dirichlet_conrey_next_primitive(x, G) == G->num) + break; + } + _acb_vec_clear(z, G->expo); + _arb_vec_clear(t, nv); + acb_clear(zeta); + acb_clear(sum); + arb_clear(eq); + acb_dirichlet_group_clear(G); + acb_dirichlet_char_clear(chi); + acb_dirichlet_conrey_clear(x); + } + + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/acb_dirichlet/zeta.c b/acb_dirichlet/zeta.c new file mode 100644 index 00000000..111f6c7b --- /dev/null +++ b/acb_dirichlet/zeta.c @@ -0,0 +1,36 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2016 Pascal Molin + +******************************************************************************/ + +#include "acb_dirichlet.h" + +void +acb_dirichlet_zeta(acb_t res, ulong order, slong prec) +{ + fmpq_t t; + fmpq_init(t); + fmpq_set_si(t, 2, order); + arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec); + fmpq_clear(t); +} diff --git a/dlog.h b/dlog.h index 16bdfd98..c35e5ff3 100644 --- a/dlog.h +++ b/dlog.h @@ -191,6 +191,7 @@ ulong dlog_bsgs(const dlog_bsgs_t t, ulong b); ulong dlog_rho(const dlog_rho_t t, ulong b); /*#define dlog_bsgs(t, b) n_discrete_log_bsgs_table(t, b)*/ +void dlog_precomp_modpe_init(dlog_precomp_t pre, ulong a, ulong p, ulong e, ulong pe, ulong num); void dlog_precomp_n_init(dlog_precomp_t pre, ulong a, ulong mod, ulong n, ulong num); void dlog_precomp_p_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong num); void dlog_precomp_pe_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong e, ulong pe, ulong num);