diff --git a/acb_dirichlet.h b/acb_dirichlet.h index da8d5498..fae55e15 100644 --- a/acb_dirichlet.h +++ b/acb_dirichlet.h @@ -24,18 +24,30 @@ typedef struct ulong q; /* modulus */ ulong q_even; /* even part of modulus */ ulong q_odd; /* odd part of modulus */ - ulong phi_q; /* phi(q) */ - ulong phi_q_odd; /* phi(q_odd) */ - slong num; /* number of odd prime factors */ - ulong * primes; /* odd primes p[k] */ + ulong phi_q; /* phi(q) = group size */ + ulong expo; /* group exponent = lcm(phi(q_even), phi(p[k]^e[k]) ) */ + slong neven; /* number of even components */ + slong num; /* number of prime components (even + odd) */ + ulong * primes; /* primes p[k] */ ulong * exponents; /* exponents e[k] */ - ulong * generators; /* generator for each odd prime p[k] */ - ulong * PHI; /* PHI(k) = phi(q_odd) / phi(p[k]^e[k]) */ + ulong * generators; /* generator for each prime p[k] */ + ulong * phi; /* phi(k) = phi(p[k]^e[k]) */ + ulong * PHI; /* PHI(k) = expo / phi(k) */ } acb_dirichlet_group_struct; typedef acb_dirichlet_group_struct acb_dirichlet_group_t[1]; +/* elements of the group, keep both number and log */ +typedef struct +{ + ulong n; /* number */ + ulong * log; /* s.t. prod generators[k]^log[k] = number */ +} +acb_conrey_struct; + +typedef acb_conrey_struct acb_conrey_t[1]; + void _acb_dirichlet_euler_product_real_ui(arb_t res, ulong s, const signed char * chi, int mod, int reciprocal, slong prec); @@ -45,8 +57,52 @@ void acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q); void acb_dirichlet_group_clear(acb_dirichlet_group_t G); +void acb_conrey_init(acb_conrey_t x, const acb_dirichlet_group_t G); + +void acb_conrey_clear(acb_conrey_t x); + +void acb_conrey_one(acb_conrey_t x, const acb_dirichlet_group_t G); + +void acb_conrey_log(acb_conrey_t x, const acb_dirichlet_group_t G, ulong m); + +int acb_conrey_next(acb_conrey_t x, const acb_dirichlet_group_t G); + +long n_dirichlet_chi_conrey(const acb_dirichlet_group_t G, const acb_conrey_t a, const acb_conrey_t b); + +long n_dirichlet_chi(const acb_dirichlet_group_t G, ulong m, ulong n); + +void acb_dirichlet_chi_conrey(acb_t res, const acb_dirichlet_group_t G, const acb_conrey_t a, const acb_conrey_t b, slong prec); + void acb_dirichlet_chi(acb_t res, const acb_dirichlet_group_t G, ulong m, ulong n, slong prec); +/* introducing character type */ + +/* character = reduced exponents, keep its order and number */ +typedef struct +{ + ulong q; /* modulus */ + ulong n; /* number */ + ulong order; /* order */ + ulong * expo; /* reduced exponents ( order * log[k] / gcd( ) ) */ +} +acb_dirichlet_char_struct; + +typedef acb_dirichlet_char_struct acb_dirichlet_char_t[1]; + +void acb_dirichlet_char_init(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G); + +void acb_dirichlet_char_clear(acb_dirichlet_char_t chi); + +void acb_dirichlet_char(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, ulong n); + +long n_dirichlet_char_eval(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n); + +void acb_dirichlet_char_eval(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n, slong prec); + +void n_dirichlet_char_vec(long *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong nv); + +void acb_dirichlet_char_vec(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n, slong prec); + #ifdef __cplusplus } #endif diff --git a/acb_dirichlet/char.c b/acb_dirichlet/char.c new file mode 100644 index 00000000..869409bc --- /dev/null +++ b/acb_dirichlet/char.c @@ -0,0 +1,66 @@ +/*============================================================================= + + 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) 2015 Jonathan Bober + Copyright (C) 2016 Fredrik Johansson + +******************************************************************************/ + +#include "acb_dirichlet.h" + +void +acb_dirichlet_char_init(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G) { + chi->expo = flint_malloc(G->num * sizeof(ulong)); +} + +void +acb_dirichlet_char_clear(acb_dirichlet_char_t chi) { + flint_free(chi->expo); +} + +void +acb_dirichlet_char_log(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, const acb_conrey_t x) +{ + ulong k, g; + g = G->expo; + /* modify exponents */ + for (k = 0; k < G->num; k++) + { + chi->expo[k] = (x->log[k] * G->PHI[k]) % G->expo; + g = n_gcd(g, chi->expo[k]); + } + for (k = 0; k < G->num; k++) + chi->expo[k] = chi->expo[k] / g; + + chi->q = G->q; + chi->order = G->expo / g; + chi->n = x->n; +} + +void +acb_dirichlet_char(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, ulong n) +{ + acb_conrey_t x; + x->log = chi->expo; + acb_conrey_log(x, G, n); + acb_dirichlet_char_log(chi, G, x); +} + diff --git a/acb_dirichlet/chareval.c b/acb_dirichlet/chareval.c new file mode 100644 index 00000000..bafdbbdc --- /dev/null +++ b/acb_dirichlet/chareval.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) 2015 Jonathan Bober + Copyright (C) 2016 Fredrik Johansson + +******************************************************************************/ + +#include "acb_dirichlet.h" + +long +n_dirichlet_char_eval(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n) +{ + ulong v = 0, k; + acb_conrey_t x; + acb_conrey_init(x, G); + acb_conrey_log(x, G, n); + for (k = 0; k < G->num; k++) + v = (v + chi->expo[k] * x->log[k]) % chi->order; + acb_conrey_clear(x); + return v; +} + +void +fmpq_dirichlet_char_eval(fmpq_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n) +{ + fmpq_set_si(res, n_dirichlet_char_eval(G, chi, n), chi->order); +} + +void +acb_dirichlet_char_eval(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n, slong prec) +{ + fmpq_t t; + fmpq_init(t); + fmpq_dirichlet_char_eval(t, G, chi, n); + arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec); + fmpq_clear(t); +} diff --git a/acb_dirichlet/charevalvec.c b/acb_dirichlet/charevalvec.c new file mode 100644 index 00000000..ddb16efe --- /dev/null +++ b/acb_dirichlet/charevalvec.c @@ -0,0 +1,67 @@ +/*============================================================================= + + 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) 2015 Jonathan Bober + Copyright (C) 2016 Fredrik Johansson + +******************************************************************************/ + +#include "acb_dirichlet.h" + +static void +set_non_invertible_values(long *v, const acb_dirichlet_group_t G, ulong nv) +{ + ulong k, l; + if (G->q_even > 1) + { + for (k = 2; k < nv; k += 2) + v[k] = -1; + } + for (l = 0; l < G->num; l++) + { + ulong p = G->primes[k]; + for (k = p; k < nv; k += p) + v[k] = -1; + } +} + +void +n_dirichlet_char_vec(long *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong nv) +{ + ulong t, k, j; + acb_conrey_t x; + acb_conrey_init(x, G); + acb_conrey_one(x, G); + t = v[1] = 0; + while( (j = acb_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 */ + set_non_invertible_values(v, G, nv); + /* copy outside modulus */ + for (k = G->q + 1; k < nv ; k++ ) + v[k] = v[k-G->q]; + acb_conrey_clear(x); +} diff --git a/acb_dirichlet/chi.c b/acb_dirichlet/chi.c index 182b0540..fb4223d6 100644 --- a/acb_dirichlet/chi.c +++ b/acb_dirichlet/chi.c @@ -14,114 +14,50 @@ /* todo: modular arithmetic */ -static ulong -chi_odd_exponent(const acb_dirichlet_group_t G, ulong m, ulong n) +long +n_dirichlet_chi_conrey(const acb_dirichlet_group_t G, const acb_conrey_t a, const acb_conrey_t b) { - ulong x, k, pk, gk, logm, logn; - - x = 0; - - for (k = 0; k < G->num; k++) - { - pk = n_pow(G->primes[k], G->exponents[k]); - gk = G->generators[k] % pk; - - logm = n_discrete_log_bsgs(m % pk, gk, pk); - logn = n_discrete_log_bsgs(n % pk, gk, pk); - - x = (x + G->PHI[k] * logm * logn) % G->phi_q_odd; - } - - return x; + ulong x, k; + x = 0; + for (k = 0; k < G->num; k++) + x = (x + G->PHI[k] * a->log[k] * b->log[k]) % G->expo; + return x; } -static ulong -chi_even_exponent(const acb_dirichlet_group_t G, ulong m, ulong n) +long +n_dirichlet_chi(const acb_dirichlet_group_t G, ulong m, ulong n) { - ulong x; - ulong q_even = G->q_even; + ulong x; + acb_conrey_t a, b; + acb_conrey_init(a, G); + acb_conrey_init(b, G); - if (q_even <= 2) - return 0; + acb_conrey_log(a, G, m); + acb_conrey_log(b, G, n); + x = n_dirichlet_chi_conrey(G, a, b); - x = 0; + acb_conrey_clear(a); + acb_conrey_clear(b); - if ((m % 4 == 3) && (n % 4 == 3)) - x = q_even / 8; + return x; +} - if (q_even > 4) - { - ulong g2, logm, logn; - - g2 = 5; - - if (m % 4 == 3) - { - m = n_negmod(m, q_even); - } - - if (n % 4 == 3) - { - n = n_negmod(n, q_even); - } - - logm = n_discrete_log_bsgs(m % q_even, g2, q_even); - logn = n_discrete_log_bsgs(n % q_even, g2, q_even); - x += logm * logn; - } - - return x % (q_even / 4); +void +acb_dirichlet_chi_conrey(acb_t res, const acb_dirichlet_group_t G, const acb_conrey_t a, const acb_conrey_t b, slong prec) +{ + fmpq_t t; + fmpq_init(t); + fmpq_set_si(t, n_dirichlet_chi_conrey(G, a, b), G->expo); + arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec); + fmpq_clear(t); } void acb_dirichlet_chi(acb_t res, const acb_dirichlet_group_t G, ulong m, ulong n, slong prec) { - fmpq_t t, u; - - ulong odd_part, even_part; - ulong q_even = G->q_even; - ulong q_odd = G->q_odd; - - if ((q_even > 1 && (n % 2 == 0)) || (q_odd > 1 && (n_gcd(q_odd, n) != 1))) - { - acb_zero(res); - return; - } - - if (q_even > 2) - { - if (q_even == 4) - { - if (m % 4 == 3 && n % 4 == 3) - even_part = q_even / 2; /* -1 */ - else - even_part = 0; /* 1 */ - } - else - { - even_part = 4 * chi_even_exponent(G, m % q_even, n % q_even); - } - } - else - { - even_part = 0; - } - - if (q_odd > 1) - odd_part = chi_odd_exponent(G, m % q_odd, n % q_odd); - else - odd_part = 0; - - fmpq_init(t); - fmpq_init(u); - - fmpq_set_si(t, 2 * even_part, q_even); - fmpq_set_si(u, 2 * odd_part, G->phi_q_odd); - fmpq_add(t, t, u); - - arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec); - - fmpq_clear(t); - fmpq_clear(u); + fmpq_t t; + fmpq_init(t); + fmpq_set_si(t, n_dirichlet_chi(G, m, n), G->expo); + arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec); + fmpq_clear(t); } - diff --git a/acb_dirichlet/conrey_log.c b/acb_dirichlet/conrey_log.c new file mode 100644 index 00000000..f7dac634 --- /dev/null +++ b/acb_dirichlet/conrey_log.c @@ -0,0 +1,91 @@ +/*============================================================================= + + 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) 2015 Jonathan Bober + Copyright (C) 2016 Fredrik Johansson + +******************************************************************************/ + +#include "acb_dirichlet.h" + +/* todo: modular arithmetic + discrete log can be computed along exponents or using p-adic log +*/ + +void +acb_conrey_init(acb_conrey_t x, const acb_dirichlet_group_t G) { + x->log = flint_malloc(G->num * sizeof(ulong)); +} + +void +acb_conrey_one(acb_conrey_t x, const acb_dirichlet_group_t G) { + ulong k; + for (k = 0; k < G->num ; k++) + x->log[k] = 0; + x->n = 1; +} + +void +acb_conrey_clear(acb_conrey_t x) { + flint_free(x->log); +} + +/* TODO: use precomputations in G if present */ +void +acb_conrey_log(acb_conrey_t x, const acb_dirichlet_group_t G, ulong m) +{ + 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) ? n_negmod(m, q_even) : m % q_even; + x->log[1] = n_discrete_log_bsgs(m2, g2, q_even); + } + /* odd part */ + for (k = G->neven; k < G->num; k++) + { + pk = n_pow(G->primes[k], G->exponents[k]); + gk = G->generators[k] % pk; + x->log[k] = n_discrete_log_bsgs(m % pk, gk, pk); + } + /* keep value m */ + x->n = m; +} + +int +acb_conrey_next(acb_conrey_t x, const acb_dirichlet_group_t G) +{ + /* update index */ + ulong k; + for (k=0; k < G->num ; k++) + { + x->n = x->n * G->generators[k]; + if(x->log[k]++ < G->phi[k]) + break; + x->log[k] = 0; + } + /* return last index modified */ + return k; +} diff --git a/acb_dirichlet/group_init.c b/acb_dirichlet/group_init.c index a57ff083..8820ed2f 100644 --- a/acb_dirichlet/group_init.c +++ b/acb_dirichlet/group_init.c @@ -36,6 +36,7 @@ void acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q) { slong k; + ulong e2 = 0; n_factor_t fac; G->q = q; @@ -46,39 +47,56 @@ acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q) { G->q_odd /= 2; G->q_even *= 2; + e2++; } n_factor_init(&fac); n_factor(&fac, G->q_odd, 1); - G->num = fac.num; + /* number of components at p=2 */ + G->neven = (e2 >= 3) ? 2 : (e2 == 2) ? 1 : 0; + G->num = G->neven + fac.num; G->primes = flint_malloc(G->num * sizeof(ulong)); G->exponents = flint_malloc(G->num * sizeof(ulong)); G->generators = flint_malloc(G->num * sizeof(ulong)); + G->phi = flint_malloc(G->num * sizeof(ulong)); G->PHI = flint_malloc(G->num * sizeof(ulong)); - for (k = 0; k < G->num; k++) + /* even part */ + if (G->q_even <= 2) { - G->primes[k] = fac.p[k]; - G->exponents[k] = fac.exp[k]; + G->expo = G->phi_q = 1; + } + else + { + G->phi_q = G->q_even / 2; + G->expo = G->phi_q / 2; } - G->phi_q_odd = 1; - for (k = 0; k < G->num; k++) - G->phi_q_odd *= (G->primes[k] - 1) * n_pow(G->primes[k], G->exponents[k]-1); + for (k = 0; k < G->neven; k++) + { + G->primes[k] = 2; + G->exponents[k] = (k==0) ? 1 : e2-2; + G->generators[k] = (k==0) ? -1 : 5; + G->phi[k] = (k==0) ? 1 : G->expo; + } - if (G->q_even == 1) - G->phi_q = G->phi_q_odd; - else - G->phi_q = G->phi_q_odd * (G->q_even / 2); + for (k = G->neven; k < G->num; k++) + { + ulong phik, p1; + G->primes[k] = fac.p[k - G->neven]; + G->exponents[k] = fac.exp[k - G->neven]; + G->generators[k] = primitive_root_p_and_p2(G->primes[k]); + p1 = G->primes[k] - 1; + phik = p1 * n_pow(G->primes[k], G->exponents[k]-1); + G->expo *= phik / n_gcd(G->expo, p1); + G->phi_q *= phik; + G->phi[k] = phik; + } for (k = 0; k < G->num; k++) { - ulong phi; - - G->generators[k] = primitive_root_p_and_p2(G->primes[k]); - phi = n_pow(G->primes[k], G->exponents[k] - 1) * (G->primes[k] - 1); - G->PHI[k] = G->phi_q_odd / phi; + G->PHI[k] = G->expo / G->phi[k]; + /* FIXME: generators[k] should be lifted mod q! */ } } - diff --git a/acb_dirichlet/l.c b/acb_dirichlet/l.c index 645c36e5..d86ab4b2 100644 --- a/acb_dirichlet/l.c +++ b/acb_dirichlet/l.c @@ -16,26 +16,30 @@ acb_dirichlet_l(acb_t res, const acb_t s, const acb_dirichlet_group_t G, ulong m, slong prec) { acb_t chi, t, u, a; - ulong k; + acb_conrey_t cm, cn; acb_init(chi); + acb_conrey_init(cm, G); + acb_conrey_init(cn, G); acb_init(t); acb_init(u); acb_init(a); + acb_conrey_log(cm, G, m); + acb_conrey_one(cn, G); acb_zero(t); - for (k = 1; k <= G->q; k++) - { - acb_dirichlet_chi(chi, G, m, k, prec); + while (1) { + /* todo: use n_dirichlet_chi and precomputed roots instead */ + acb_dirichlet_chi_conrey(chi, G, cm, cn, prec); - if (!acb_is_zero(chi)) - { - acb_set_ui(a, k); - acb_div_ui(a, a, G->q, prec); - acb_hurwitz_zeta(u, s, a, prec); - acb_addmul(t, chi, u, prec); - } + acb_set_ui(a, cn->n); + acb_div_ui(a, a, G->q, prec); + acb_hurwitz_zeta(u, s, a, prec); + acb_addmul(t, chi, u, prec); + + if (acb_conrey_next(cn, G) == G->num) + break; } acb_set_ui(u, G->q); @@ -43,6 +47,9 @@ acb_dirichlet_l(acb_t res, const acb_t s, acb_pow(u, u, a, prec); acb_mul(res, t, u, prec); + acb_conrey_clear(cm); + acb_conrey_clear(cn); + acb_clear(chi); acb_clear(t); acb_clear(u);