From 00118eb55f2da102e43622ff8b8514078a2730fd Mon Sep 17 00:00:00 2001 From: Pascal Date: Sat, 8 Oct 2016 17:22:19 +0200 Subject: [PATCH] remove all conrey --- acb_dirichlet.h | 4 +- acb_dirichlet/chi.c | 4 +- acb_dirichlet/chi_vec.c | 7 +- acb_dirichlet/gauss_sum.c | 26 +- acb_dirichlet/gauss_sum_factor.c | 12 +- acb_dirichlet/gauss_sum_order2.c | 8 +- acb_dirichlet/gauss_sum_theta.c | 7 +- acb_dirichlet/jacobi_sum.c | 20 +- acb_dirichlet/jacobi_sum_factor.c | 17 +- acb_dirichlet/jacobi_sum_gauss.c | 2 +- acb_dirichlet/jacobi_sum_naive.c | 33 ++- acb_dirichlet/l_hurwitz.c | 25 +- acb_dirichlet/pairing.c | 2 +- acb_dirichlet/pairing_conrey.c | 4 +- acb_dirichlet/root_number.c | 4 +- acb_dirichlet/test/t-chi.c | 12 +- acb_dirichlet/test/t-gauss.c | 21 +- acb_dirichlet/test/t-jacobi.c | 7 +- acb_dirichlet/test/t-l.c | 2 +- acb_dirichlet/test/t-thetanull.c | 18 +- acb_dirichlet/theta_arb.c | 27 +- acb_dirichlet/ui_theta_arb.c | 2 +- dirichlet.h | 124 ++------- dirichlet/char.c | 19 -- dirichlet/char_clear.c | 19 -- dirichlet/char_conrey.c | 32 --- dirichlet/char_eq_deep.c | 38 --- dirichlet/char_first_primitive.c | 20 -- dirichlet/char_init.c | 26 -- dirichlet/char_mul.c | 19 -- dirichlet/char_next.c | 22 -- dirichlet/char_next_primitive.c | 22 -- dirichlet/char_normalize.c | 46 ---- dirichlet/char_primitive.c | 24 -- dirichlet/char_print.c | 23 -- dirichlet/{ui_chi.c => chi.c} | 4 +- dirichlet/{char_one.c => chi_char.c} | 19 +- dirichlet/chi_vec.c | 30 +++ .../{ui_chi_vec_loop.c => chi_vec_loop.c} | 26 +- dirichlet/chi_vec_primeloop.c | 82 ++++++ dirichlet/conrey_conductor.c | 2 +- dirichlet/conrey_index.c | 6 +- dirichlet/conrey_lift.c | 36 +++ dirichlet/conrey_log.c | 2 +- dirichlet/conrey_lower.c | 42 ++++ dirichlet/conrey_mul.c | 2 +- dirichlet/conrey_next.c | 2 +- dirichlet/conrey_next_primitive.c | 4 +- dirichlet/conrey_order.c | 2 +- dirichlet/conrey_parity.c | 2 +- dirichlet/conrey_pow.c | 2 +- dirichlet/conrey_primitive.c | 47 ---- dirichlet/group_init.c | 22 +- dirichlet/index_conrey.c | 2 +- dirichlet/{ui_pairing.c => pairing.c} | 4 +- dirichlet/test/t-chars.c | 235 +++++++++--------- dirichlet/test/t-conrey.c | 153 ------------ dirichlet/test/t-properties.c | 199 +++++++++++++++ {acb_dirichlet => dirichlet}/test/t-vec.c | 36 ++- dirichlet/ui_chi_conrey.c | 24 -- dirichlet/ui_chi_vec.c | 21 -- dirichlet/ui_chi_vec_primeloop.c | 72 ------ dirichlet/ui_conductor.c | 2 +- dirichlet/ui_order.c | 2 +- dirichlet/ui_pairing_conrey.c | 26 -- dirichlet/ui_parity.c | 2 +- dirichlet/ui_vec_set_null.c | 2 +- 67 files changed, 757 insertions(+), 1053 deletions(-) delete mode 100644 dirichlet/char.c delete mode 100644 dirichlet/char_clear.c delete mode 100644 dirichlet/char_conrey.c delete mode 100644 dirichlet/char_eq_deep.c delete mode 100644 dirichlet/char_first_primitive.c delete mode 100644 dirichlet/char_init.c delete mode 100644 dirichlet/char_mul.c delete mode 100644 dirichlet/char_next.c delete mode 100644 dirichlet/char_next_primitive.c delete mode 100644 dirichlet/char_normalize.c delete mode 100644 dirichlet/char_primitive.c delete mode 100644 dirichlet/char_print.c rename dirichlet/{ui_chi.c => chi.c} (83%) rename dirichlet/{char_one.c => chi_char.c} (60%) create mode 100644 dirichlet/chi_vec.c rename dirichlet/{ui_chi_vec_loop.c => chi_vec_loop.c} (57%) create mode 100644 dirichlet/chi_vec_primeloop.c create mode 100644 dirichlet/conrey_lift.c create mode 100644 dirichlet/conrey_lower.c delete mode 100644 dirichlet/conrey_primitive.c rename dirichlet/{ui_pairing.c => pairing.c} (87%) delete mode 100644 dirichlet/test/t-conrey.c create mode 100644 dirichlet/test/t-properties.c rename {acb_dirichlet => dirichlet}/test/t-vec.c (51%) delete mode 100644 dirichlet/ui_chi_conrey.c delete mode 100644 dirichlet/ui_chi_vec.c delete mode 100644 dirichlet/ui_chi_vec_primeloop.c delete mode 100644 dirichlet/ui_pairing_conrey.c diff --git a/acb_dirichlet.h b/acb_dirichlet.h index 607e5e44..7330f572 100644 --- a/acb_dirichlet.h +++ b/acb_dirichlet.h @@ -32,8 +32,8 @@ void _acb_dirichlet_euler_product_real_ui(arb_t res, ulong s, void acb_dirichlet_eta(acb_t res, const acb_t s, slong prec); -void acb_dirichlet_pairing_conrey(acb_t res, const dirichlet_group_t G, const dirichlet_conrey_t a, const dirichlet_conrey_t b, slong prec); void acb_dirichlet_pairing(acb_t res, const dirichlet_group_t G, ulong m, ulong n, slong prec); +void acb_dirichlet_chi_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 */ typedef struct @@ -71,7 +71,7 @@ void acb_dirichlet_ui_theta_arb(acb_t res, const dirichlet_group_t G, ulong a, c void acb_dirichlet_gauss_sum_naive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec); void acb_dirichlet_gauss_sum_factor(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec); -void acb_dirichlet_gauss_sum_order2(acb_t res, const dirichlet_char_t chi, slong prec); +void acb_dirichlet_gauss_sum_order2(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec); void acb_dirichlet_gauss_sum_theta(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec); void acb_dirichlet_gauss_sum(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec); diff --git a/acb_dirichlet/chi.c b/acb_dirichlet/chi.c index e8b1d6b2..3d67536c 100644 --- a/acb_dirichlet/chi.c +++ b/acb_dirichlet/chi.c @@ -17,14 +17,14 @@ void acb_dirichlet_chi(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, ulong n, slong prec) { ulong expo; - expo = dirichlet_ui_chi(G, chi, n); + expo = dirichlet_chi(G, chi, n); if (expo == DIRICHLET_CHI_NULL) acb_zero(res); else { fmpq_t t; fmpq_init(t); - fmpq_set_si(t, 2 * expo , chi->order.n); + fmpq_set_si(t, 2 * expo , G->expo); arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec); fmpq_clear(t); } diff --git a/acb_dirichlet/chi_vec.c b/acb_dirichlet/chi_vec.c index f5f45b18..c3462cdc 100644 --- a/acb_dirichlet/chi_vec.c +++ b/acb_dirichlet/chi_vec.c @@ -15,13 +15,14 @@ void acb_dirichlet_chi_vec(acb_ptr v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv, slong prec) { slong k; - ulong * a; + ulong * a, order; acb_dirichlet_powers_t t; a = flint_malloc(nv * sizeof(ulong)); - dirichlet_ui_chi_vec(a, G, chi, nv); + order = dirichlet_order_char(G, chi); + dirichlet_chi_vec_order(a, G, chi, order, nv); - acb_dirichlet_powers_init(t, chi->order.n, nv, prec); + acb_dirichlet_powers_init(t, order, nv, prec); acb_zero(v + 0); for (k = 0; k < nv; k++) diff --git a/acb_dirichlet/gauss_sum.c b/acb_dirichlet/gauss_sum.c index 3b8447ff..15816498 100644 --- a/acb_dirichlet/gauss_sum.c +++ b/acb_dirichlet/gauss_sum.c @@ -12,11 +12,11 @@ #include "acb_dirichlet.h" static void -gauss_sum_non_primitive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec) +gauss_sum_non_primitive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, ulong cond, slong prec) { slong k, mu = 1; - ulong NN0 = G->q / chi->conductor; + ulong NN0 = G->q / cond; /* G(chi) = mu(N/N0)chi0(N/N0)G(chi0) */ @@ -45,7 +45,7 @@ gauss_sum_non_primitive(acb_t res, const dirichlet_group_t G, const dirichlet_ch mu *= -1; } - if (chi->x->n == 1) + if (chi->n == 1) { acb_set_si(res, mu); } @@ -56,15 +56,15 @@ gauss_sum_non_primitive(acb_t res, const dirichlet_group_t G, const dirichlet_ch dirichlet_char_t chi0; acb_t z; - dirichlet_subgroup_init(G0, G, chi->conductor); - dirichlet_char_init(chi0, G); - dirichlet_char_primitive(chi0, G0, G, chi); + /* primitive char associated to chi */ + dirichlet_subgroup_init(G0, G, cond); + dirichlet_char_init(chi0, G0); + dirichlet_char_lower(chi0, G0, chi, G); acb_init(z); acb_dirichlet_gauss_sum(z, G0, chi0, prec); acb_dirichlet_chi(res, G0, chi0, NN0, prec); - acb_mul(res, res, z, prec); acb_mul_si(res, res, mu, prec); @@ -77,14 +77,15 @@ gauss_sum_non_primitive(acb_t res, const dirichlet_group_t G, const dirichlet_ch void acb_dirichlet_gauss_sum(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec) { + ulong cond = dirichlet_conductor_char(G, chi); /* TODO: no need, factor also does it... */ - if (chi->conductor != G->q) + if (cond != G->q) { - gauss_sum_non_primitive(res, G, chi, prec); + gauss_sum_non_primitive(res, G, chi, cond, prec); } - else if (chi->order.n <= 2) + else if (dirichlet_char_is_real(G, chi)) { - acb_dirichlet_gauss_sum_order2(res, chi, prec); + acb_dirichlet_gauss_sum_order2(res, G, chi, prec); } else if (G->num > 1 && G->num > G->neven) { @@ -92,6 +93,7 @@ acb_dirichlet_gauss_sum(acb_t res, const dirichlet_group_t G, const dirichlet_ch } else { + /* must be non primitive */ if (acb_dirichlet_theta_length_d(G->q, 1, prec) > G->q) acb_dirichlet_gauss_sum_naive(res, G, chi, prec); else @@ -104,7 +106,7 @@ acb_dirichlet_gauss_sum_ui(acb_t res, const dirichlet_group_t G, ulong a, slong { dirichlet_char_t chi; dirichlet_char_init(chi, G); - dirichlet_char(chi, G, a); + dirichlet_char_log(chi, G, a); acb_dirichlet_gauss_sum(res, G, chi, prec); dirichlet_char_clear(chi); } diff --git a/acb_dirichlet/gauss_sum_factor.c b/acb_dirichlet/gauss_sum_factor.c index 29fa0a32..fe96bd87 100644 --- a/acb_dirichlet/gauss_sum_factor.c +++ b/acb_dirichlet/gauss_sum_factor.c @@ -22,7 +22,7 @@ acb_dirichlet_gauss_sum_factor(acb_t res, const dirichlet_group_t G, const diric for (k = (G->neven == 2); k < G->num; k++) { /* if e > 1 and not primitive, 0 */ - if (chi->x->log[k] % G->P[k].p == 0 && G->P[k].e > 1) + if (chi->log[k] % G->P[k].p == 0 && G->P[k].e > 1) { acb_zero(res); return; @@ -42,17 +42,15 @@ acb_dirichlet_gauss_sum_factor(acb_t res, const dirichlet_group_t G, const diric dirichlet_subgroup_init(Gp, G, pe); dirichlet_char_init(chip, Gp); - chip->x->n = chi->x->n % pe; + chip->n = chi->n % pe; if (k == 1 && G->neven == 2) { - chip->x->log[0] = chi->x->log[0]; - chip->x->log[1] = chi->x->log[1]; + chip->log[0] = chi->log[0]; + chip->log[1] = chi->log[1]; } else - chip->x->log[0] = chi->x->log[k]; - - dirichlet_char_conrey(chip, Gp, NULL); + chip->log[0] = chi->log[k]; /* chi_pe(a, q/pe) * G_pe(a) */ acb_dirichlet_gauss_sum(tmp, Gp, chip, prec); diff --git a/acb_dirichlet/gauss_sum_order2.c b/acb_dirichlet/gauss_sum_order2.c index 39fed38b..863626d8 100644 --- a/acb_dirichlet/gauss_sum_order2.c +++ b/acb_dirichlet/gauss_sum_order2.c @@ -12,16 +12,16 @@ #include "acb_dirichlet.h" void -acb_dirichlet_gauss_sum_order2(acb_t res, const dirichlet_char_t chi, slong prec) +acb_dirichlet_gauss_sum_order2(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec) { - if (chi->parity) + if (dirichlet_parity_char(G, chi)) { arb_zero(acb_realref(res)); - arb_sqrt_ui(acb_imagref(res), chi->q, prec); + arb_sqrt_ui(acb_imagref(res), G->q, prec); } else { arb_zero(acb_imagref(res)); - arb_sqrt_ui(acb_realref(res), chi->q, prec); + arb_sqrt_ui(acb_realref(res), G->q, prec); } } diff --git a/acb_dirichlet/gauss_sum_theta.c b/acb_dirichlet/gauss_sum_theta.c index ab97f384..21a5a09a 100644 --- a/acb_dirichlet/gauss_sum_theta.c +++ b/acb_dirichlet/gauss_sum_theta.c @@ -14,8 +14,9 @@ void acb_dirichlet_gauss_sum_theta(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec) { - if (chi->conductor < G->q || (G->q == 300 && (chi->x->n == 71 || chi->x->n == 131)) - || (G->q == 600 && (chi->x->n == 11 || chi->x->n == 491))) + ulong cond = dirichlet_conductor_char(G, chi); + if (cond < G->q || (G->q == 300 && (chi->n == 71 || chi->n == 131)) + || (G->q == 600 && (chi->n == 11 || chi->n == 491))) { flint_printf("gauss_sum_theta: non available for non primitive character" "or exceptional characters chi_300(71,.), chi_300(131,.), " @@ -26,7 +27,7 @@ acb_dirichlet_gauss_sum_theta(acb_t res, const dirichlet_group_t G, const dirich { acb_t iq; acb_init(iq); - acb_dirichlet_gauss_sum_order2(iq, chi, prec); + acb_dirichlet_gauss_sum_order2(iq, G, chi, prec); acb_dirichlet_root_number_theta(res, G, chi, prec); acb_mul(res, res, iq, prec); acb_clear(iq); diff --git a/acb_dirichlet/jacobi_sum.c b/acb_dirichlet/jacobi_sum.c index ab67d3c4..8ce90fad 100644 --- a/acb_dirichlet/jacobi_sum.c +++ b/acb_dirichlet/jacobi_sum.c @@ -47,16 +47,16 @@ acb_dirichlet_jacobi_sum(acb_t res, const dirichlet_group_t G, const dirichlet_c { acb_zero(res); } - else if (chi1->x->n == 1 || chi2->x->n == 1) + else if (chi1->n == 1 || chi2->n == 1) { - ulong cond = (chi1->x->n == 1) ? chi2->conductor : chi1->conductor; + ulong cond = (chi1->n == 1) ? dirichlet_conductor_char(G, chi2) : dirichlet_conductor_char(G, chi1); acb_set_si(res, jacobi_one(G, cond)); } - else if (nmod_mul(chi1->x->n, chi2->x->n, G->mod) == 1) + else if (nmod_mul(chi1->n, chi2->n, G->mod) == 1) { ulong n; - n = jacobi_one(G, chi1->conductor); - if (chi1->parity) + n = jacobi_one(G, dirichlet_conductor_char(G, chi1)); + if (dirichlet_parity_char(G, chi1)) acb_set_si(res, -n); else acb_set_si(res, n); @@ -83,14 +83,14 @@ acb_dirichlet_jacobi_sum_ui(acb_t res, const dirichlet_group_t G, ulong a, ulong } else if (a == 1 || b == 1) { - ulong cond = (a == 1) ? dirichlet_ui_conductor(G, b) : dirichlet_ui_conductor(G, a); + ulong cond = (a == 1) ? dirichlet_conductor_ui(G, b) : dirichlet_conductor_ui(G, a); acb_set_si(res, jacobi_one(G, cond)); } else if (nmod_mul(a, b, G->mod) == 1) { ulong n; - n = jacobi_one(G, dirichlet_ui_conductor(G, a)); - if (dirichlet_ui_parity(G, a)) + n = jacobi_one(G, dirichlet_conductor_ui(G, a)); + if (dirichlet_parity_ui(G, a)) acb_set_si(res, -n); else acb_set_si(res, n); @@ -100,8 +100,8 @@ acb_dirichlet_jacobi_sum_ui(acb_t res, const dirichlet_group_t G, ulong a, ulong dirichlet_char_t chi1, chi2; dirichlet_char_init(chi1, G); dirichlet_char_init(chi2, G); - dirichlet_char(chi1, G, a); - dirichlet_char(chi2, G, b); + dirichlet_char_log(chi1, G, a); + dirichlet_char_log(chi2, G, b); acb_dirichlet_jacobi_sum(res, G, chi1, chi2, prec); dirichlet_char_clear(chi1); dirichlet_char_clear(chi2); diff --git a/acb_dirichlet/jacobi_sum_factor.c b/acb_dirichlet/jacobi_sum_factor.c index ddecf60a..90c9bf44 100644 --- a/acb_dirichlet/jacobi_sum_factor.c +++ b/acb_dirichlet/jacobi_sum_factor.c @@ -28,15 +28,15 @@ acb_dirichlet_jacobi_sum_factor(acb_t res, const dirichlet_group_t G, const dir p = G->P[k].p; e = G->P[k].e; pe = G->P[k].pe; - ap = chi1->x->n % pe.n; - bp = chi2->x->n % pe.n; + ap = chi1->n % pe.n; + bp = chi2->n % pe.n; if (ap == 1 || bp == 1 || nmod_mul(ap, bp, pe) == 1) { slong r; ulong cond; - cond = (ap == 1) ? chi2->conductor : chi1->conductor; + cond = (ap == 1) ? dirichlet_conductor_char(G, chi2) : dirichlet_conductor_char(G, chi1); r = jacobi_one_prime(p, e, pe.n, cond); /* chi(a,-1) if ap * bp = 1 */ @@ -54,13 +54,10 @@ acb_dirichlet_jacobi_sum_factor(acb_t res, const dirichlet_group_t G, const dir dirichlet_char_init(chi1p, Gp); dirichlet_char_init(chi2p, Gp); - chi1p->x->n = ap; - chi1p->x->log[0] = chi1->x->log[k]; - chi2p->x->n = ap; - chi2p->x->log[0] = chi2->x->log[k]; - - dirichlet_char_conrey(chi1p, Gp, NULL); - dirichlet_char_conrey(chi2p, Gp, NULL); + chi1p->n = ap; + chi1p->log[0] = chi1->log[k]; + chi2p->n = ap; + chi2p->log[0] = chi2->log[k]; /* TODO: work out gauss relations for e > 1 */ if (p <= 100 || e > 1) diff --git a/acb_dirichlet/jacobi_sum_gauss.c b/acb_dirichlet/jacobi_sum_gauss.c index 5d825ec6..7b2cdc2a 100644 --- a/acb_dirichlet/jacobi_sum_gauss.c +++ b/acb_dirichlet/jacobi_sum_gauss.c @@ -25,7 +25,7 @@ acb_dirichlet_jacobi_sum_gauss(acb_t res, const dirichlet_group_t G, const diric acb_init(tmp); acb_dirichlet_gauss_sum(res, G, chi1, prec); - if (chi2->x->n == chi1->x->n) + if (chi2->n == chi1->n) acb_set(tmp, res); else acb_dirichlet_gauss_sum(tmp, G, chi2, prec); diff --git a/acb_dirichlet/jacobi_sum_naive.c b/acb_dirichlet/jacobi_sum_naive.c index 5541e6f6..bbbbfc0e 100644 --- a/acb_dirichlet/jacobi_sum_naive.c +++ b/acb_dirichlet/jacobi_sum_naive.c @@ -15,31 +15,30 @@ void acb_dirichlet_jacobi_sum_naive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec) { - ulong k1, k2, m1, m2, g, e; + ulong k1, k2, m1, m2, g, e, m; ulong * v1, * v2; slong *v; - nmod_t order; + nmod_t expo; acb_t z; v1 = flint_malloc(G->q * sizeof(ulong)); v2 = flint_malloc(G->q * sizeof(ulong)); - dirichlet_ui_vec_set_null(v1, G, G->q); - dirichlet_ui_chi_vec_loop(v1, G, chi1, G->q); + dirichlet_vec_set_null(v1, G, G->q); + dirichlet_chi_vec_loop(v1, G, chi1, G->q); - dirichlet_ui_vec_set_null(v2, G, G->q); - dirichlet_ui_chi_vec_loop(v2, G, chi2, G->q); + dirichlet_vec_set_null(v2, G, G->q); + dirichlet_chi_vec_loop(v2, G, chi2, G->q); - m1 = chi1->order.n; - m2 = chi2->order.n; - g = n_gcd(m1, m2); - nmod_init(&order, m1 * (m2 / g)); - m1 = order.n / m1; - m2 = order.n / m2; + nmod_init(&expo, G->expo); + m1 = dirichlet_order_char(G, chi1); + m2 = dirichlet_order_char(G, chi2); + g = m1 * m2 / n_gcd(m1, m2); + m = G->expo / g; - v = flint_malloc(order.n * sizeof(slong)); + v = flint_malloc(g * sizeof(slong)); - for (e = 0; e < order.n; e++) + for (e = 0; e < g; e++) v[e] = 0; for (k1 = 2, k2 = G->q - 1; k2 > 1; k1++, k2--) @@ -47,13 +46,13 @@ acb_dirichlet_jacobi_sum_naive(acb_t res, const dirichlet_group_t G, const diric if (v1[k1] == DIRICHLET_CHI_NULL || v2[k2] == DIRICHLET_CHI_NULL) continue; - e = nmod_add(v1[k1] * m1, v2[k2] * m2, order); + e = nmod_add(v1[k1], v2[k2], expo) / m; v[e]++; } acb_init(z); - acb_nth_root(z, order.n, prec); - acb_dirichlet_si_poly_evaluate(res, v, order.n, z, prec); + acb_nth_root(z, g, prec); + acb_dirichlet_si_poly_evaluate(res, v, g, z, prec); acb_clear(z); flint_free(v); diff --git a/acb_dirichlet/l_hurwitz.c b/acb_dirichlet/l_hurwitz.c index f219e6e8..5d10e72a 100644 --- a/acb_dirichlet/l_hurwitz.c +++ b/acb_dirichlet/l_hurwitz.c @@ -18,18 +18,17 @@ void acb_dirichlet_l_hurwitz(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec) { - ulong chin; + ulong order, chin, mult; acb_t t, u, a; acb_ptr z; - dirichlet_conrey_t cn; + dirichlet_char_t cn; int deflate; /* remove pole in Hurwitz zeta at s = 1 */ deflate = 0; if (acb_is_one(s)) { - /* character is principal */ - if (chi->x->n == 1) + if (dirichlet_char_is_principal(chi)) { acb_indeterminate(res); return; @@ -37,21 +36,23 @@ acb_dirichlet_l_hurwitz(acb_t res, const acb_t s, deflate = 1; } - dirichlet_conrey_init(cn, G); + dirichlet_char_init(cn, G); acb_init(t); acb_init(u); acb_init(a); - dirichlet_conrey_one(cn, G); + dirichlet_char_one(cn, G); acb_zero(t); prec += n_clog(G->phi_q, 2); - z = _acb_vec_init(chi->order.n); - _acb_vec_nth_roots(z, chi->order.n, prec); + order = dirichlet_order_char(G, chi); + mult = G->expo / order; + z = _acb_vec_init(order); + _acb_vec_nth_roots(z, order, prec); do { - chin = dirichlet_ui_chi_conrey(G, chi, cn); + chin = dirichlet_chi_char(G, chi, cn) / mult; acb_set_ui(a, cn->n); acb_div_ui(a, a, G->q, prec); @@ -63,16 +64,16 @@ acb_dirichlet_l_hurwitz(acb_t res, const acb_t s, acb_addmul(t, z + chin, u, prec); - } while (dirichlet_conrey_next(cn, G) >= 0); + } while (dirichlet_char_next(cn, G) >= 0); acb_set_ui(u, G->q); acb_neg(a, s); acb_pow(u, u, a, prec); acb_mul(res, t, u, prec); - dirichlet_conrey_clear(cn); + dirichlet_char_clear(cn); - _acb_vec_clear(z, chi->order.n); + _acb_vec_clear(z, order); acb_clear(t); acb_clear(u); acb_clear(a); diff --git a/acb_dirichlet/pairing.c b/acb_dirichlet/pairing.c index c05cbeb1..871e405f 100644 --- a/acb_dirichlet/pairing.c +++ b/acb_dirichlet/pairing.c @@ -15,7 +15,7 @@ void acb_dirichlet_pairing(acb_t res, const dirichlet_group_t G, ulong m, ulong n, slong prec) { ulong expo; - expo = dirichlet_ui_pairing(G, m, n); + expo = dirichlet_pairing(G, m, n); if (expo == DIRICHLET_CHI_NULL) acb_zero(res); else diff --git a/acb_dirichlet/pairing_conrey.c b/acb_dirichlet/pairing_conrey.c index 18d8aa54..c9c1284a 100644 --- a/acb_dirichlet/pairing_conrey.c +++ b/acb_dirichlet/pairing_conrey.c @@ -12,10 +12,10 @@ #include "acb_dirichlet.h" void -acb_dirichlet_pairing_conrey(acb_t res, const dirichlet_group_t G, const dirichlet_conrey_t a, const dirichlet_conrey_t b, slong prec) +acb_dirichlet_chi_char(acb_t res, const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b, slong prec) { ulong expo; - expo = dirichlet_ui_pairing_conrey(G, a, b); + expo = dirichlet_chi_char(G, a, b); if (expo == DIRICHLET_CHI_NULL) acb_zero(res); else diff --git a/acb_dirichlet/root_number.c b/acb_dirichlet/root_number.c index 33324182..1e4869b9 100644 --- a/acb_dirichlet/root_number.c +++ b/acb_dirichlet/root_number.c @@ -31,7 +31,7 @@ acb_dirichlet_root_number_theta(acb_t res, const dirichlet_group_t G, const diri void acb_dirichlet_root_number(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec) { - if (chi->conductor < G->q) + if (dirichlet_conductor_char(G, chi) < G->q) { flint_printf("root number: need primitive character\n"); abort(); @@ -40,7 +40,7 @@ acb_dirichlet_root_number(acb_t res, const dirichlet_group_t G, const dirichlet_ { acb_t iq; acb_init(iq); - acb_dirichlet_gauss_sum_order2(iq, chi, prec); + acb_dirichlet_gauss_sum_order2(iq, G, chi, prec); acb_dirichlet_gauss_sum(res, G, chi, prec); acb_div(res, res, iq, prec); acb_clear(iq); diff --git a/acb_dirichlet/test/t-chi.c b/acb_dirichlet/test/t-chi.c index e12eb7c0..ba67142b 100644 --- a/acb_dirichlet/test/t-chi.c +++ b/acb_dirichlet/test/t-chi.c @@ -44,7 +44,7 @@ int main() m = 1 + n_randint(state, q); } while (n_gcd(q, m) != 1); - dirichlet_char(chi, G, m); + dirichlet_char_log(chi, G, m); n1 = n_randint(state, 1000); n2 = n_randint(state, 1000); @@ -53,7 +53,6 @@ int main() acb_dirichlet_pairing(zn2, G, m, n1, 53); if (!acb_overlaps(zn1, zn2)) { - dirichlet_conrey_t x; flint_printf("FAIL: overlap\n\n"); flint_printf("q = %wu\n\n", q); flint_printf("m = %wu\n\n", m); @@ -61,11 +60,10 @@ int main() flint_printf("char = "); acb_printd(zn1, 15); flint_printf("\n\n"); flint_printf("pairing = "); acb_printd(zn2, 15); flint_printf("\n\n"); dirichlet_char_print(G, chi); - dirichlet_conrey_init(x, G); - dirichlet_conrey_log(x, G, m); - flint_printf("log(m) = "); dirichlet_conrey_print(G, x); - dirichlet_conrey_log(x, G, n1); - flint_printf("log(n1) = "); dirichlet_conrey_print(G, x); + dirichlet_char_log(chi, G, m); + flint_printf("log(m) = "); dirichlet_char_print(G, chi); + dirichlet_char_log(chi, G, n1); + flint_printf("log(n1) = "); dirichlet_char_print(G, chi); abort(); } diff --git a/acb_dirichlet/test/t-gauss.c b/acb_dirichlet/test/t-gauss.c index 6b2afdfd..5cc2d6da 100644 --- a/acb_dirichlet/test/t-gauss.c +++ b/acb_dirichlet/test/t-gauss.c @@ -24,29 +24,25 @@ int main() for (q = 3; q < 250; q ++) { dirichlet_group_t G; - dirichlet_conrey_t x; dirichlet_char_t chi; acb_t s1, s2, s3, s4; dirichlet_group_init(G, q); - dirichlet_conrey_init(x, G); dirichlet_char_init(chi, G); acb_init(s1); acb_init(s2); acb_init(s3); acb_init(s4); - dirichlet_conrey_one(x, G); + dirichlet_char_one(chi, G); while (1) { - dirichlet_char_conrey(chi, G, x); - acb_dirichlet_gauss_sum_naive(s1, G, chi, prec); acb_dirichlet_gauss_sum(s2, G, chi, prec); acb_dirichlet_gauss_sum_factor(s3, G, chi, prec); - if (chi->conductor == G->q) + if (dirichlet_conductor_char(G, chi) == G->q) acb_dirichlet_gauss_sum_theta(s4, G, chi, prec); else acb_set(s4, s1); @@ -55,19 +51,19 @@ int main() || !acb_overlaps(s1, s3) || !acb_overlaps(s1, s4)) { - flint_printf("FAIL: G(chi_%wu(%wu))\n\n", q, chi->x->n); - flint_printf("\nnaive ", q, x->n); + flint_printf("FAIL: G(chi_%wu(%wu))\n\n", q, chi->n); + flint_printf("\nnaive "); acb_printd(s1, 25); - flint_printf("\ndefault ", q, x->n); + flint_printf("\ndefault "); acb_printd(s2, 25); - flint_printf("\nfactor ", q, x->n); + flint_printf("\nfactor "); acb_printd(s3, 25); - flint_printf("\ntheta ", q, x->n); + flint_printf("\ntheta "); acb_printd(s4, 25); abort(); } - if (dirichlet_conrey_next(x, G) < 0) + if (dirichlet_char_next(chi, G) < 0) break; } acb_clear(s1); @@ -77,7 +73,6 @@ int main() dirichlet_group_clear(G); dirichlet_char_clear(chi); - dirichlet_conrey_clear(x); } flint_cleanup(); diff --git a/acb_dirichlet/test/t-jacobi.c b/acb_dirichlet/test/t-jacobi.c index e7dfa520..b68b7a09 100644 --- a/acb_dirichlet/test/t-jacobi.c +++ b/acb_dirichlet/test/t-jacobi.c @@ -52,15 +52,16 @@ int main() if (!acb_overlaps(s1, s2)) { flint_printf("FAIL: J_%wu(%wu,%wu)", - q, chi1->x->n, chi2->x->n); + q, chi1->n, chi2->n); flint_printf("\nnaive "); acb_printd(s1, 25); flint_printf("\ndefault "); acb_printd(s2, 25); flint_printf("\n"); flint_printf("cond = %wu, %wu, %wu\n", - chi1->conductor, chi2->conductor, - dirichlet_ui_conductor(G, nmod_mul(chi1->x->n, chi2->x->n, G->mod)) + dirichlet_conductor_char(G, chi1), + dirichlet_conductor_char(G, chi2), + dirichlet_conductor_ui(G, nmod_mul(chi1->n, chi2->n, G->mod)) ); abort(); } diff --git a/acb_dirichlet/test/t-l.c b/acb_dirichlet/test/t-l.c index e37a5efb..5e472569 100644 --- a/acb_dirichlet/test/t-l.c +++ b/acb_dirichlet/test/t-l.c @@ -99,7 +99,7 @@ int main() dirichlet_group_init(G, q[i]); dirichlet_char_init(chi, G); - dirichlet_char(chi, G, m[i]); + dirichlet_char_log(chi, G, m[i]); for (j = 0; j < nx; j++) { diff --git a/acb_dirichlet/test/t-thetanull.c b/acb_dirichlet/test/t-thetanull.c index f40ac1e2..a5a37e09 100644 --- a/acb_dirichlet/test/t-thetanull.c +++ b/acb_dirichlet/test/t-thetanull.c @@ -68,24 +68,22 @@ int main() dirichlet_char_first_primitive(chi, G); do { - ulong m; + acb_zero(sum); + dirichlet_chi_vec(v, G, chi, nv); - dirichlet_ui_chi_vec(v, G, chi, nv); - - m = G->expo / chi->order.n; - tt = dirichlet_char_parity(chi) ? kt : t; + tt = dirichlet_parity_char(G, chi) ? kt : t; for (k = 1; k < nv; k++) if (v[k] != DIRICHLET_CHI_NULL) - acb_addmul_arb(sum, z + (v[k] * m), tt + k, prec); + acb_addmul_arb(sum, z + v[k], tt + k, prec); - if ((q == 300 && (chi->x->n == 71 || chi->x->n == 131)) - || (q == 600 && (chi->x->n == 11 || chi->x->n == 491))) + if ((q == 300 && (chi->n == 71 || chi->n == 131)) + || (q == 600 && (chi->n == 11 || chi->n == 491))) { if (!acb_contains_zero(sum)) { - flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->x->n); + flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->n); acb_printd(sum, 10); flint_printf("\n"); dirichlet_char_print(G, chi); @@ -95,7 +93,7 @@ int main() } else if (acb_contains_zero(sum)) { - flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->x->n); + flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->n); acb_printd(sum, 10); flint_printf("\n"); dirichlet_char_print(G, chi); diff --git a/acb_dirichlet/theta_arb.c b/acb_dirichlet/theta_arb.c index f6a643e4..4e7774ee 100644 --- a/acb_dirichlet/theta_arb.c +++ b/acb_dirichlet/theta_arb.c @@ -15,14 +15,18 @@ void _acb_dirichlet_theta_arb_smallorder(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, const arb_t xt, slong len, slong prec) { + ulong order; ulong * a; + int parity; acb_dirichlet_powers_t z; + parity = dirichlet_parity_char(G, chi); + order = dirichlet_order_char(G, chi); a = flint_malloc(len * sizeof(ulong)); - dirichlet_ui_chi_vec(a, G, chi, len); + dirichlet_chi_vec_order(a, G, chi, order, len); - _acb_dirichlet_powers_init(z, chi->order.n, 0, 0, prec); - acb_dirichlet_qseries_arb_powers_smallorder(res, xt, chi->parity, a, z, len, prec); + _acb_dirichlet_powers_init(z, order, 0, 0, prec); + acb_dirichlet_qseries_arb_powers_smallorder(res, xt, parity, a, z, len, prec); acb_dirichlet_powers_clear(z); flint_free(a); @@ -34,7 +38,7 @@ _acb_dirichlet_theta_arb_series(acb_t res, const dirichlet_group_t G, const diri acb_ptr a; a = _acb_vec_init(len); acb_dirichlet_chi_vec(a, G, chi, len, prec); - if (chi->parity) + if (dirichlet_parity_char(G, chi)) { slong k; for (k = 2; k < len; k++) @@ -47,15 +51,18 @@ _acb_dirichlet_theta_arb_series(acb_t res, const dirichlet_group_t G, const diri 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 * a; + ulong order, * a; acb_dirichlet_powers_t z; + int parity; + parity = dirichlet_parity_char(G, chi); + order = dirichlet_order_char(G, chi); a = flint_malloc(len * sizeof(ulong)); - dirichlet_ui_chi_vec(a, G, chi, len); + dirichlet_chi_vec_order(a, G, chi, order, len); - acb_dirichlet_powers_init(z, chi->order.n, len, prec); + acb_dirichlet_powers_init(z, order, len, prec); - acb_dirichlet_qseries_arb_powers_naive(res, xt, chi->parity, a, z, len, prec); + acb_dirichlet_qseries_arb_powers_naive(res, xt, parity, a, z, len, prec); acb_dirichlet_powers_clear(z); flint_free(a); @@ -65,6 +72,7 @@ void acb_dirichlet_theta_arb(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, const arb_t t, slong prec) { slong len; + ulong order; arb_t xt; mag_t e; @@ -80,7 +88,8 @@ acb_dirichlet_theta_arb(acb_t res, const dirichlet_group_t G, const dirichlet_ch arb_exp(xt, xt, prec); /* TODO: tune this limit */ - if (chi->order.n < 30) + order = dirichlet_order_char(G, chi); + if (order < 30) _acb_dirichlet_theta_arb_smallorder(res, G, chi, xt, len, prec); else _acb_dirichlet_theta_arb_naive(res, G, chi, xt, len, prec); diff --git a/acb_dirichlet/ui_theta_arb.c b/acb_dirichlet/ui_theta_arb.c index f98d2e42..51f7043c 100644 --- a/acb_dirichlet/ui_theta_arb.c +++ b/acb_dirichlet/ui_theta_arb.c @@ -17,7 +17,7 @@ acb_dirichlet_ui_theta_arb(acb_t res, const dirichlet_group_t G, ulong a, const dirichlet_char_t chi; dirichlet_char_init(chi, G); - dirichlet_char(chi, G, a); + dirichlet_char_log(chi, G, a); acb_dirichlet_theta_arb(res, G, chi, t, prec); diff --git a/dirichlet.h b/dirichlet.h index 01d0cfca..673bde66 100644 --- a/dirichlet.h +++ b/dirichlet.h @@ -24,13 +24,15 @@ extern "C" { #endif +#define MAX_FACTORS 15 + /* should this dlog pointer be in the prime or the global group? */ typedef struct { ulong p; /* underlying prime */ int e; /* exponent */ nmod_t pe; /* modulus */ - ulong phi; /* phi(p^e) */ + nmod_t phi; /* phi(p^e) */ ulong g; /* conrey generator */ dlog_precomp_struct * dlog; /* precomputed data for discrete log mod p^e */ } @@ -69,9 +71,9 @@ void dirichlet_group_dlog_clear(dirichlet_group_t G); /* properties of elements without log */ ulong dirichlet_number_primitive(const dirichlet_group_t G); -ulong dirichlet_ui_conductor(const dirichlet_group_t G, ulong a); -int dirichlet_ui_parity(const dirichlet_group_t G, ulong a); -ulong dirichlet_ui_order(const dirichlet_group_t G, ulong a); +ulong dirichlet_conductor_ui(const dirichlet_group_t G, ulong a); +int dirichlet_parity_ui(const dirichlet_group_t G, ulong a); +ulong dirichlet_order_ui(const dirichlet_group_t G, ulong a); /* characters, keep both number and log */ typedef struct @@ -103,9 +105,9 @@ dirichlet_char_eq(const dirichlet_char_t x, const dirichlet_char_t y) } int dirichlet_char_eq_deep(const dirichlet_group_t G, const dirichlet_char_t x, const dirichlet_char_t y); -int dirichlet_char_parity(const dirichlet_group_t G, const dirichlet_char_t x); -ulong dirichlet_char_conductor(const dirichlet_group_t G, const dirichlet_char_t x); -ulong dirichlet_char_order(const dirichlet_group_t G, const dirichlet_char_t x); +int dirichlet_parity_char(const dirichlet_group_t G, const dirichlet_char_t x); +ulong dirichlet_conductor_char(const dirichlet_group_t G, const dirichlet_char_t x); +ulong dirichlet_order_char(const dirichlet_group_t G, const dirichlet_char_t x); void dirichlet_char_log(dirichlet_char_t x, const dirichlet_group_t G, ulong m); ulong dirichlet_char_exp(dirichlet_char_t x, const dirichlet_group_t G); @@ -121,108 +123,36 @@ int dirichlet_char_next_primitive(dirichlet_char_t x, const dirichlet_group_t G) void dirichlet_char_mul(dirichlet_char_t c, const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b); void dirichlet_char_pow(dirichlet_char_t c, const dirichlet_group_t G, const dirichlet_char_t a, ulong n); -void dirichlet_char_primitive(dirichlet_char_t y, const dirichlet_group_t G, const dirichlet_char_t x, ulong cond); + +void dirichlet_char_lower(dirichlet_char_t y, const dirichlet_group_t H, const dirichlet_char_t x, const dirichlet_group_t G); +void dirichlet_char_lift(dirichlet_char_t y, const dirichlet_group_t G, const dirichlet_char_t x, const dirichlet_group_t H); #define DIRICHLET_CHI_NULL UWORD_MAX -ulong dirichlet_ui_pairing_char(const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b); -ulong dirichlet_ui_pairing(const dirichlet_group_t G, ulong m, ulong n); +ulong dirichlet_pairing(const dirichlet_group_t G, ulong m, ulong n); -void dirichlet_pairing_char(acb_t res, const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b, slong prec); -void dirichlet_pairing(acb_t res, const dirichlet_group_t G, ulong m, ulong n, slong prec); - -/* introducing character type */ - -/* character = reduced exponents, keep order, number and conductor */ -typedef struct +DIRICHLET_INLINE int +dirichlet_char_is_principal(const dirichlet_char_t chi) { - ulong q; /* modulus */ - nmod_t order; /* order */ - dirichlet_char_t x; - ulong * expo; /* reduced exponents ( log[k] * PHI[k] / gcd( ) ) */ - int parity; /* 0 for even char, 1 for odd */ - ulong conductor; -} -dirichlet_fullchar_struct; - -typedef dirichlet_fullchar_struct dirichlet_fullchar_t[1]; - -DIRICHLET_INLINE ulong -dirichlet_fullchar_order(const dirichlet_fullchar_t chi) -{ - return chi->order.n; -} - -DIRICHLET_INLINE ulong -dirichlet_fullchar_conductor(const dirichlet_fullchar_t chi) -{ - return chi->conductor; + return (chi->n == 1); } DIRICHLET_INLINE int -dirichlet_fullchar_parity(const dirichlet_fullchar_t chi) +dirichlet_char_is_real(const dirichlet_group_t G, const dirichlet_char_t chi) { - return chi->parity; + return (nmod_mul(chi->n, chi->n, G->mod) == 1); } -void dirichlet_fullchar_init(dirichlet_fullchar_t chi, const dirichlet_group_t G); -void dirichlet_fullchar_clear(dirichlet_fullchar_t chi); -void dirichlet_fullchar_print(const dirichlet_group_t G, const dirichlet_fullchar_t chi); +ulong dirichlet_chi_char(const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b); +ulong dirichlet_chi(const dirichlet_group_t G, const dirichlet_char_t chi, ulong n); -DIRICHLET_INLINE void -dirichlet_fullchar_set(dirichlet_fullchar_t chi1, const dirichlet_group_t G, const dirichlet_fullchar_t chi2) -{ - slong k; - - chi1->q = chi2->q; - chi1->conductor = chi2->conductor; - chi1->order = chi2->order; - chi1->parity = chi2->parity; - dirichlet_char_set(chi1->x, G, chi2->x); - for (k = 0; k < G->num; k++) - chi1->expo[k] = chi2->expo[k]; -} - -DIRICHLET_INLINE int -dirichlet_fullchar_eq(const dirichlet_fullchar_t chi1, const dirichlet_fullchar_t chi2) -{ - return (chi1->q == chi2->q && chi1->x->n == chi2->x->n); -} - -int dirichlet_fullchar_eq_deep(const dirichlet_group_t G, const dirichlet_fullchar_t chi1, const dirichlet_fullchar_t chi2); -DIRICHLET_INLINE int -dirichlet_fullchar_is_principal(const dirichlet_fullchar_t chi) -{ - return (chi->x->n == 1); -} -DIRICHLET_INLINE int -dirichlet_fullchar_is_real(const dirichlet_fullchar_t chi) -{ - return (chi->order.n <= 2); -} - -void dirichlet_fullchar(dirichlet_fullchar_t chi, const dirichlet_group_t G, ulong n); -void dirichlet_fullchar_char(dirichlet_fullchar_t chi, const dirichlet_group_t G, const dirichlet_char_t x); -void dirichlet_fullchar_set_expo(dirichlet_fullchar_t chi, const dirichlet_group_t G); -void dirichlet_fullchar_normalize(dirichlet_fullchar_t chi, const dirichlet_group_t G); -void dirichlet_fullchar_denormalize(dirichlet_fullchar_t chi, const dirichlet_group_t G); - -void dirichlet_fullchar_mul(dirichlet_fullchar_t chi12, const dirichlet_group_t G, const dirichlet_fullchar_t chi1, const dirichlet_fullchar_t chi2); -void dirichlet_fullchar_primitive(dirichlet_fullchar_t chi0, const dirichlet_group_t G0, const dirichlet_group_t G, const dirichlet_fullchar_t chi); - -void dirichlet_fullchar_one(dirichlet_fullchar_t chi, const dirichlet_group_t G); -void dirichlet_fullchar_first_primitive(dirichlet_fullchar_t chi, const dirichlet_group_t G); - -int dirichlet_fullchar_next(dirichlet_fullchar_t chi, const dirichlet_group_t G); -int dirichlet_fullchar_next_primitive(dirichlet_fullchar_t chi, const dirichlet_group_t G); - -ulong dirichlet_ui_chi_char(const dirichlet_group_t G, const dirichlet_fullchar_t chi, const dirichlet_char_t x); -ulong dirichlet_ui_chi(const dirichlet_group_t G, const dirichlet_fullchar_t chi, ulong n); - -void dirichlet_ui_vec_set_null(ulong *v, const dirichlet_group_t G, slong nv); -void dirichlet_ui_chi_vec_loop(ulong *v, const dirichlet_group_t G, const dirichlet_fullchar_t chi, slong nv); -void dirichlet_ui_chi_vec_primeloop(ulong *v, const dirichlet_group_t G, const dirichlet_fullchar_t chi, slong nv); -void dirichlet_ui_chi_vec(ulong *v, const dirichlet_group_t G, const dirichlet_fullchar_t chi, slong nv); +void dirichlet_vec_set_null(ulong *v, const dirichlet_group_t G, slong nv); +void dirichlet_chi_vec_loop(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv); +void dirichlet_chi_vec_primeloop(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv); +void dirichlet_chi_vec(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv); +void dirichlet_chi_vec_loop_order(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, ulong order, slong nv); +void dirichlet_chi_vec_primeloop_order(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, ulong order, slong nv); +void dirichlet_chi_vec_order(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, ulong order, slong nv); #ifdef __cplusplus } diff --git a/dirichlet/char.c b/dirichlet/char.c deleted file mode 100644 index c8f57ad9..00000000 --- a/dirichlet/char.c +++ /dev/null @@ -1,19 +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 . -*/ - -#include "dirichlet.h" - -void -dirichlet_fullchar(dirichlet_fullchar_t chi, const dirichlet_group_t G, ulong n) -{ - dirichlet_char_log(chi->x, G, n); - dirichlet_fullchar_char(chi, G, NULL); -} diff --git a/dirichlet/char_clear.c b/dirichlet/char_clear.c deleted file mode 100644 index f4432c79..00000000 --- a/dirichlet/char_clear.c +++ /dev/null @@ -1,19 +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 . -*/ - -#include "dirichlet.h" - -void -dirichlet_fullchar_clear(dirichlet_fullchar_t chi) -{ - dirichlet_char_clear(chi->x); - flint_free(chi->expo); -} diff --git a/dirichlet/char_conrey.c b/dirichlet/char_conrey.c deleted file mode 100644 index 4f82fbfa..00000000 --- a/dirichlet/char_conrey.c +++ /dev/null @@ -1,32 +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 . -*/ - -#include "dirichlet.h" - -/* char n has exponents = log[k]*PHI[k] / gcd and order expo / gcd - * so that log = expo[k] */ -void -dirichlet_fullchar_char(dirichlet_fullchar_t chi, const dirichlet_group_t G, const dirichlet_char_t x) -{ - /* assume chi->x already set if x == NULL */ - if (x == NULL) - x = chi->x; - else - dirichlet_char_set(chi->x, G, x); - - chi->q = G->q; - chi->parity = dirichlet_char_parity(G, x); - chi->conductor = dirichlet_char_conductor(G, x); - - dirichlet_fullchar_set_expo(chi, G); - /* optional: divide by gcd to obtain true order */ - dirichlet_fullchar_normalize(chi, G); -} diff --git a/dirichlet/char_eq_deep.c b/dirichlet/char_eq_deep.c deleted file mode 100644 index d4e845f6..00000000 --- a/dirichlet/char_eq_deep.c +++ /dev/null @@ -1,38 +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 . -*/ - -#include "dirichlet.h" - -int -dirichlet_fullchar_eq_deep(const dirichlet_group_t G, const dirichlet_fullchar_t chi1, const dirichlet_fullchar_t chi2) -{ - dirichlet_char_t x, y; - - if (chi1->q != chi2->q) - return 0; - - if (chi1->order.n != chi2->order.n) - return 0; - - if (chi1->conductor != chi2->conductor) - return 0; - - if (!dirichlet_char_eq_deep(G, chi1->x, chi2->x)) - return 0; - - x->n = y->n = 1; - x->log = chi1->expo; - y->log = chi2->expo; - if (!dirichlet_char_eq_deep(G, x, y)) - return 0; - - return 1; -} diff --git a/dirichlet/char_first_primitive.c b/dirichlet/char_first_primitive.c deleted file mode 100644 index acf378a3..00000000 --- a/dirichlet/char_first_primitive.c +++ /dev/null @@ -1,20 +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 . -*/ - -#include "dirichlet.h" - -void -dirichlet_fullchar_first_primitive(dirichlet_fullchar_t chi, const dirichlet_group_t G) -{ - dirichlet_char_first_primitive(chi->x, G); - dirichlet_fullchar_char(chi, G, NULL); - dirichlet_fullchar_normalize(chi, G); -} diff --git a/dirichlet/char_init.c b/dirichlet/char_init.c deleted file mode 100644 index 67d5c9a6..00000000 --- a/dirichlet/char_init.c +++ /dev/null @@ -1,26 +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 . -*/ - -#include "dirichlet.h" - -void -dirichlet_fullchar_init(dirichlet_fullchar_t chi, const dirichlet_group_t G) -{ - slong k; - dirichlet_char_init(chi->x, G); - chi->expo = flint_malloc(G->num * sizeof(ulong)); - chi->q = G->q; - chi->conductor = 1; - chi->parity = 0; - for (k = 0; k < G->num; k++) - chi->expo[k] = 0; - nmod_init(&(chi->order), 1); -} diff --git a/dirichlet/char_mul.c b/dirichlet/char_mul.c deleted file mode 100644 index 28ceee77..00000000 --- a/dirichlet/char_mul.c +++ /dev/null @@ -1,19 +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 . -*/ - -#include "dirichlet.h" - -void -dirichlet_fullchar_mul(dirichlet_fullchar_t chi12, const dirichlet_group_t G, const dirichlet_fullchar_t chi1, const dirichlet_fullchar_t chi2) -{ - dirichlet_char_mul(chi12->x, G, chi1->x, chi2->x); - dirichlet_fullchar_char(chi12, G, NULL); -} diff --git a/dirichlet/char_next.c b/dirichlet/char_next.c deleted file mode 100644 index f00206a0..00000000 --- a/dirichlet/char_next.c +++ /dev/null @@ -1,22 +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 . -*/ - -#include "dirichlet.h" - -int -dirichlet_fullchar_next(dirichlet_fullchar_t chi, const dirichlet_group_t G) -{ - int k; - k = dirichlet_char_next(chi->x, G); - dirichlet_fullchar_char(chi, G, NULL); - /* return last index modified */ - return k; -} diff --git a/dirichlet/char_next_primitive.c b/dirichlet/char_next_primitive.c deleted file mode 100644 index eba33457..00000000 --- a/dirichlet/char_next_primitive.c +++ /dev/null @@ -1,22 +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 . -*/ - -#include "dirichlet.h" - -int -dirichlet_fullchar_next_primitive(dirichlet_fullchar_t chi, const dirichlet_group_t G) -{ - int k; - k = dirichlet_char_next_primitive(chi->x, G); - dirichlet_fullchar_char(chi, G, NULL); - /* return last index modified */ - return k; -} diff --git a/dirichlet/char_normalize.c b/dirichlet/char_normalize.c deleted file mode 100644 index a7287c74..00000000 --- a/dirichlet/char_normalize.c +++ /dev/null @@ -1,46 +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 . -*/ - -#include "dirichlet.h" - -void -dirichlet_fullchar_set_expo(dirichlet_fullchar_t chi, const dirichlet_group_t G) -{ - slong k; - for (k = 0; k < G->num; k++) - /* no overflow: log[k] < phi[k] and G->expo = phi[k] * PHI[k] */ - chi->expo[k] = chi->x->log[k] * G->PHI[k]; -} - -void -dirichlet_fullchar_normalize(dirichlet_fullchar_t chi, const dirichlet_group_t G) -{ - ulong k, g; - g = G->expo; - - for (k = 0; k < G->num; k++) - g = n_gcd(g, chi->expo[k]); - - for (k = 0; k < G->num; k++) - chi->expo[k] = chi->expo[k] / g; - - nmod_init(&chi->order, G->expo / g); -} - -void -dirichlet_fullchar_denormalize(dirichlet_fullchar_t chi, const dirichlet_group_t G) -{ - ulong k, g; - g = G->expo / chi->order.n; - - for (k = 0; k < G->num; k++) - chi->expo[k] *= g; -} diff --git a/dirichlet/char_primitive.c b/dirichlet/char_primitive.c deleted file mode 100644 index c8fc771f..00000000 --- a/dirichlet/char_primitive.c +++ /dev/null @@ -1,24 +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 . -*/ - -#include "dirichlet.h" - -void -dirichlet_fullchar_primitive(dirichlet_fullchar_t chi0, const dirichlet_group_t G0, const dirichlet_group_t G, const dirichlet_fullchar_t chi) -{ - chi0->q = chi->conductor; - chi0->parity = chi->parity; - chi0->conductor = chi->conductor; - dirichlet_char_primitive(chi0->x, G, chi->x, chi->conductor); - dirichlet_fullchar_set_expo(chi0, G0); - /* optional: divide by gcd to obtain true order */ - dirichlet_fullchar_normalize(chi0, G0); -} diff --git a/dirichlet/char_print.c b/dirichlet/char_print.c deleted file mode 100644 index 6898bfe2..00000000 --- a/dirichlet/char_print.c +++ /dev/null @@ -1,23 +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 . -*/ - -#include "dirichlet.h" - -void -dirichlet_fullchar_print(const dirichlet_group_t G, const dirichlet_fullchar_t chi) -{ - dirichlet_char_t x; - flint_printf("chi_%wu(%wu,.) of order %wu, parity %wd, index ", G->q, chi->x->n, chi->order, chi->parity); - dirichlet_char_print(G, chi->x); - flint_printf(" and exponents "); - x->log = chi->expo; - dirichlet_char_print(G, x); -} diff --git a/dirichlet/ui_chi.c b/dirichlet/chi.c similarity index 83% rename from dirichlet/ui_chi.c rename to dirichlet/chi.c index 7133e614..43901218 100644 --- a/dirichlet/ui_chi.c +++ b/dirichlet/chi.c @@ -12,7 +12,7 @@ #include "dirichlet.h" ulong -dirichlet_ui_chi(const dirichlet_group_t G, const dirichlet_fullchar_t chi, ulong n) +dirichlet_chi(const dirichlet_group_t G, const dirichlet_char_t chi, ulong n) { if (n_gcd(G->q, n) > 1) { @@ -25,7 +25,7 @@ dirichlet_ui_chi(const dirichlet_group_t G, const dirichlet_fullchar_t chi, ulon dirichlet_char_init(x, G); dirichlet_char_log(x, G, n); - v = dirichlet_ui_chi_char(G, chi, x); + v = dirichlet_chi_char(G, chi, x); dirichlet_char_clear(x); return v; diff --git a/dirichlet/char_one.c b/dirichlet/chi_char.c similarity index 60% rename from dirichlet/char_one.c rename to dirichlet/chi_char.c index 7f872555..951b4cb2 100644 --- a/dirichlet/char_one.c +++ b/dirichlet/chi_char.c @@ -11,15 +11,16 @@ #include "dirichlet.h" -void -dirichlet_fullchar_one(dirichlet_fullchar_t chi, const dirichlet_group_t G) +ulong +dirichlet_chi_char(const dirichlet_group_t G, const dirichlet_char_t chi, const dirichlet_char_t x) { - slong k; - dirichlet_char_one(chi->x, G); - chi->q = G->q; - chi->conductor = 1; - chi->parity = 0; + ulong v = 0, k; + nmod_t order; + + nmod_init(&order, G->expo); + for (k = 0; k < G->num; k++) - chi->expo[k] = 0; - nmod_init(&(chi->order), 1); + v = nmod_add(v, G->PHI[k] * nmod_mul(chi->log[k], x->log[k], G->P[k].phi), order); + + return v; } diff --git a/dirichlet/chi_vec.c b/dirichlet/chi_vec.c new file mode 100644 index 00000000..f56fc293 --- /dev/null +++ b/dirichlet/chi_vec.c @@ -0,0 +1,30 @@ +/* + 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 . +*/ + +#include "dirichlet.h" + +void +dirichlet_chi_vec(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv) +{ + if (2 * nv > G->phi_q) + dirichlet_chi_vec_loop(v, G, chi, nv); + else + dirichlet_chi_vec_primeloop(v, G, chi, nv); +} + +void +dirichlet_chi_vec_order(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, ulong order, slong nv) +{ + if (2 * nv > G->phi_q) + dirichlet_chi_vec_loop_order(v, G, chi, order, nv); + else + dirichlet_chi_vec_primeloop_order(v, G, chi, order, nv); +} diff --git a/dirichlet/ui_chi_vec_loop.c b/dirichlet/chi_vec_loop.c similarity index 57% rename from dirichlet/ui_chi_vec_loop.c rename to dirichlet/chi_vec_loop.c index 8ff81d76..d58341cf 100644 --- a/dirichlet/ui_chi_vec_loop.c +++ b/dirichlet/chi_vec_loop.c @@ -11,17 +11,33 @@ #include "dirichlet.h" +static void +dirichlet_exponents_char(ulong * expo, const dirichlet_group_t G, const dirichlet_char_t chi, ulong order) +{ + slong k; + ulong factor = G->expo / order; + for (k = 0; k < G->num; k++) + /* no overflow: log[k] < phi[k] and G->expo = phi[k] * PHI[k] */ + expo[k] = (chi->log[k] * G->PHI[k]) / factor; +} + /* loop over whole group */ void -dirichlet_ui_chi_vec_loop(ulong *v, const dirichlet_group_t G, const dirichlet_fullchar_t chi, slong nv) +dirichlet_chi_vec_loop_order(ulong * v, const dirichlet_group_t G, const dirichlet_char_t chi, ulong order, slong nv) { int j; ulong t; slong k; + ulong expo[MAX_FACTORS]; dirichlet_char_t x; + nmod_t o; + dirichlet_char_init(x, G); dirichlet_char_one(x, G); + dirichlet_exponents_char(expo, G, chi, order); + nmod_init(&o, order); + for (k = 0; k < nv; k++) v[k] = DIRICHLET_CHI_NULL; @@ -31,7 +47,7 @@ dirichlet_ui_chi_vec_loop(ulong *v, const dirichlet_group_t G, const dirichlet_f { /* exponents were modified up to j */ for (k = G->num - 1; k >= j; k--) - t = nmod_add(t, chi->expo[k], chi->order); + t = nmod_add(t, expo[k], o); if (x->n < nv) v[x->n] = t; @@ -46,3 +62,9 @@ dirichlet_ui_chi_vec_loop(ulong *v, const dirichlet_group_t G, const dirichlet_f dirichlet_char_clear(x); } + +void +dirichlet_chi_vec_loop(ulong * v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv) +{ + dirichlet_chi_vec_loop_order(v, G, chi, G->expo, nv); +} diff --git a/dirichlet/chi_vec_primeloop.c b/dirichlet/chi_vec_primeloop.c new file mode 100644 index 00000000..90536a07 --- /dev/null +++ b/dirichlet/chi_vec_primeloop.c @@ -0,0 +1,82 @@ +/* + 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 . +*/ + +#include "dirichlet.h" + +static void +dirichlet_chi_vec_evenpart(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, ulong order, slong nv) +{ + ulong mult = G->expo / order; + if (G->neven >= 1 && chi->log[0]) + { + ulong x, c3 = G->PHI[0] / mult; + for (x = 3; x < nv; x += 4) + v[x] = c3; + } + if (G->neven == 2 && chi->log[1]) + { + ulong x, g, vx, xp, c4; + nmod_t pe, o; + nmod_init(&o, order); + + pe = G->P[1].pe; + g = G->P[1].g; + + vx = c4 = (chi->log[1] * G->PHI[1]) / mult; + for (x = g; x > 1;) + { + + for (xp = x; xp < nv; xp += pe.n) + v[xp] = nmod_add(v[xp], vx, o); + + for (xp = pe.n - x; xp < nv; xp += pe.n) + v[xp] = nmod_add(v[xp], vx, o); + + x = nmod_mul(x, g, pe); + vx = nmod_add(vx, c4, o); + } + } +} + +/* loop over primary components */ +void +dirichlet_chi_vec_primeloop_order(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, ulong order, slong nv) +{ + slong k, l; + ulong mult = G->expo / order; + nmod_t o; + nmod_init(&o, order); + + for (k = 0; k < nv; k++) + v[k] = 0; + + if (G->neven) + dirichlet_chi_vec_evenpart(v, G, chi, order, nv); + + for (l = G->neven; l < G->num; l++) + { + dirichlet_prime_group_struct P = G->P[l]; + + /* FIXME: there may be some precomputed dlog in P if needed */ + if (P.dlog == NULL) + dlog_vec_add(v, nv, P.g, (chi->log[l] * G->PHI[l]) / mult, P.pe, P.phi.n, o); + else + dlog_vec_add_precomp(v, nv, P.dlog, P.g, (chi->log[l] * G->PHI[l]) / mult, P.pe, P.phi.n, o); + + } + dirichlet_vec_set_null(v, G, nv); +} + +void +dirichlet_chi_vec_primeloop(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv) +{ + dirichlet_chi_vec_primeloop_order(v, G, chi, G->expo, nv); +} diff --git a/dirichlet/conrey_conductor.c b/dirichlet/conrey_conductor.c index 24afc2b9..281c4ccb 100644 --- a/dirichlet/conrey_conductor.c +++ b/dirichlet/conrey_conductor.c @@ -12,7 +12,7 @@ #include "dirichlet.h" ulong -dirichlet_char_conductor(const dirichlet_group_t G, const dirichlet_char_t x) +dirichlet_conductor_char(const dirichlet_group_t G, const dirichlet_char_t x) { int k, f; ulong cond = 1; diff --git a/dirichlet/conrey_index.c b/dirichlet/conrey_index.c index ddea9463..50d14cfb 100644 --- a/dirichlet/conrey_index.c +++ b/dirichlet/conrey_index.c @@ -30,10 +30,10 @@ dirichlet_char_index(dirichlet_char_t x, const dirichlet_group_t G, ulong j) { slong k; - for (k = 0; k < G->num; k++) + for (k = G->num - 1; k >= 0; k--) { - x->log[k] = j % G->P[k].phi; - j = j / G->P[k].phi; + x->log[k] = j % G->P[k].phi.n; + j = j / G->P[k].phi.n; } dirichlet_char_exp(x, G); diff --git a/dirichlet/conrey_lift.c b/dirichlet/conrey_lift.c new file mode 100644 index 00000000..36854a70 --- /dev/null +++ b/dirichlet/conrey_lift.c @@ -0,0 +1,36 @@ +/* + 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 . +*/ + +#include "dirichlet.h" + +void +dirichlet_char_lift(dirichlet_char_t y, const dirichlet_group_t G, const dirichlet_char_t x, const dirichlet_group_t H) +{ + slong k, l; + + if (G->q % H->q != 0) + { + flint_printf("conrey_lift: lower modulus %wu does not divide %wu\n", H->q, G->q); + abort(); + } + + for (k = 0, l = 0; k < G->num && l < H->num; k++) + { + if (G->P[k].p == H->P[l].p) + { + slong e = n_pow(G->P[k].p, G->P[k].e - H->P[l].e); + y->log[k] = x->log[l] * e; + l++; + } + } + + dirichlet_char_exp(y, G); +} diff --git a/dirichlet/conrey_log.c b/dirichlet/conrey_log.c index e53245ad..a989627f 100644 --- a/dirichlet/conrey_log.c +++ b/dirichlet/conrey_log.c @@ -39,7 +39,7 @@ dirichlet_char_log(dirichlet_char_t x, const dirichlet_group_t G, ulong m) dirichlet_prime_group_struct P = G->P[k]; if (P.dlog == NULL) { - x->log[k] = dlog_once(m % P.pe.n, P.g, P.pe, P.phi); + x->log[k] = dlog_once(m % P.pe.n, P.g, P.pe, P.phi.n); } else { diff --git a/dirichlet/conrey_lower.c b/dirichlet/conrey_lower.c new file mode 100644 index 00000000..16bdcc8d --- /dev/null +++ b/dirichlet/conrey_lower.c @@ -0,0 +1,42 @@ +/* + 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 . +*/ + +#include "dirichlet.h" + +void +dirichlet_char_lower(dirichlet_char_t y, const dirichlet_group_t H, const dirichlet_char_t x, const dirichlet_group_t G) +{ + slong k, l; + + if (G->q % H->q != 0) + { + flint_printf("conrey_lower: lower modulus %wu does not divide %wu\n", H->q, G->q); + abort(); + } + + for (k = 0, l = 0; k < G->num && l < H->num; k++) + { + ulong p = G->P[k].p; + if (p == H->P[l].p) + { + ulong pef = n_pow(G->P[k].p, G->P[k].e - H->P[l].e); + ulong a = x->log[k]; + if (a % pef) + { + flint_printf("conrey_lower: conductor does not divide lower modulus %wu", H->q); + abort(); + } + y->log[l] = a / pef; + l++; + } + } + dirichlet_char_exp(y, H); +} diff --git a/dirichlet/conrey_mul.c b/dirichlet/conrey_mul.c index a8071405..6194455d 100644 --- a/dirichlet/conrey_mul.c +++ b/dirichlet/conrey_mul.c @@ -16,6 +16,6 @@ dirichlet_char_mul(dirichlet_char_t c, const dirichlet_group_t G, const dirichle { ulong k; for (k = 0; k < G->num ; k++) - c->log[k] = (a->log[k] + b->log[k]) % G->P[k].phi; + c->log[k] = nmod_add(a->log[k], b->log[k], G->P[k].phi); c->n = nmod_mul(a->n, b->n, G->mod); } diff --git a/dirichlet/conrey_next.c b/dirichlet/conrey_next.c index 4a1db3c1..140ecbf0 100644 --- a/dirichlet/conrey_next.c +++ b/dirichlet/conrey_next.c @@ -20,7 +20,7 @@ dirichlet_char_next(dirichlet_char_t x, const dirichlet_group_t G) { x->n = nmod_mul(x->n, G->generators[k], G->mod); x->log[k] += 1; - if (x->log[k] < G->P[k].phi) + if (x->log[k] < G->P[k].phi.n) break; x->log[k] = 0; } diff --git a/dirichlet/conrey_next_primitive.c b/dirichlet/conrey_next_primitive.c index a39b6fc1..9831d785 100644 --- a/dirichlet/conrey_next_primitive.c +++ b/dirichlet/conrey_next_primitive.c @@ -38,9 +38,9 @@ dirichlet_char_next_primitive(dirichlet_char_t x, const dirichlet_group_t G) x->n = nmod_mul(x->n, G->generators[k], G->mod); x->log[k]++; } - if (x->log[k] < G->P[k].phi) + if (x->log[k] < G->P[k].phi.n) break; - if (x->log[k] == G->P[k].phi) + if (x->log[k] == G->P[k].phi.n) x->n = nmod_mul(x->n, G->generators[k], G->mod); x->log[k] = 1; #else diff --git a/dirichlet/conrey_order.c b/dirichlet/conrey_order.c index d8da6fc7..51875f2b 100644 --- a/dirichlet/conrey_order.c +++ b/dirichlet/conrey_order.c @@ -12,7 +12,7 @@ #include "dirichlet.h" ulong -dirichlet_char_order(const dirichlet_group_t G, const dirichlet_char_t x) +dirichlet_order_char(const dirichlet_group_t G, const dirichlet_char_t x) { ulong k, g; g = G->expo; diff --git a/dirichlet/conrey_parity.c b/dirichlet/conrey_parity.c index 3ba6d181..dcaba5a1 100644 --- a/dirichlet/conrey_parity.c +++ b/dirichlet/conrey_parity.c @@ -12,7 +12,7 @@ #include "dirichlet.h" int -dirichlet_char_parity(const dirichlet_group_t G, const dirichlet_char_t x) +dirichlet_parity_char(const dirichlet_group_t G, const dirichlet_char_t x) { int k, odd = 0; for (k = 0; k < G->num; k++) diff --git a/dirichlet/conrey_pow.c b/dirichlet/conrey_pow.c index 3a682bd7..804b9767 100644 --- a/dirichlet/conrey_pow.c +++ b/dirichlet/conrey_pow.c @@ -16,6 +16,6 @@ dirichlet_char_pow(dirichlet_char_t c, const dirichlet_group_t G, const dirichle { ulong k; for (k = 0; k < G->num ; k++) - c->log[k] = n_mulmod2(a->log[k], n, G->P[k].phi); + c->log[k] = nmod_mul(a->log[k], n, G->P[k].phi); c->n = nmod_pow_ui(a->n, n, G->mod); } diff --git a/dirichlet/conrey_primitive.c b/dirichlet/conrey_primitive.c deleted file mode 100644 index 2a73c64b..00000000 --- a/dirichlet/conrey_primitive.c +++ /dev/null @@ -1,47 +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 . -*/ - -#include "dirichlet.h" - -void -dirichlet_char_primitive(dirichlet_char_t y, const dirichlet_group_t G, const dirichlet_char_t x, ulong cond) -{ - slong k, l; - - l = 0; - if (cond % 4 == 0) - { - y->log[l++] = x->log[0]; - - if (cond % 8 == 0) - { - ulong l2 = x->log[1]; - n_remove(&l2, 2); - y->log[l++] = l2; - } - } - - for (k = G->neven; k < G->num; k++) - { - if (x->log[k]) - { - ulong p, lp; - p = G->P[k].p; - if (cond % p == 0) - { - lp = x->log[k]; - n_remove(&lp, p); - y->log[l++] = lp; - } - } - } - -} diff --git a/dirichlet/group_init.c b/dirichlet/group_init.c index 4a98121f..7c6c9523 100644 --- a/dirichlet/group_init.c +++ b/dirichlet/group_init.c @@ -45,12 +45,12 @@ dirichlet_prime_group_init(dirichlet_prime_group_struct * P, ulong p, int e) if (p == 2) { P->e = 2; - P->phi = 2; + nmod_init(&P->phi, 2); P->g = (1 << e) - 1; } else { - P->phi = 1 << (e - 2); + nmod_init(&P->phi, 1 << (e - 2)); P->g = 5; } } @@ -58,7 +58,7 @@ dirichlet_prime_group_init(dirichlet_prime_group_struct * P, ulong p, int e) { ulong pe1; pe1 = n_pow(p, e - 1); - P->phi = (p-1) * pe1; + nmod_init(&P->phi, (p-1) * pe1); nmod_init(&P->pe, p * pe1); P->g = primitive_root_p_and_p2(p); } @@ -75,19 +75,19 @@ dirichlet_group_lift_generators(dirichlet_group_t G) if (G->neven) { G->phi_q = G->q_even / 2; - G->expo = P[G->neven - 1].phi; + G->expo = P[G->neven - 1].phi.n; } for (k = G->neven; k < G->num; k++) { - G->phi_q *= P[k].phi; - G->expo *= P[k].phi / n_gcd(G->expo, P[k].p - 1); + G->phi_q *= P[k].phi.n; + G->expo *= P[k].phi.n / n_gcd(G->expo, P[k].p - 1); } for (k = 0; k < G->num; k++) { nmod_t pe; ulong qpe, v; - G->PHI[k] = G->expo / G->P[k].phi; + G->PHI[k] = G->expo / G->P[k].phi.n; /* lift generators mod q */ /* u * p^e + v * q/p^e = 1 -> g mod q = 1 + (g-1) * v*(q/p^e) */ pe = G->P[k].pe; @@ -199,6 +199,8 @@ dirichlet_subgroup_init(dirichlet_group_t H, const dirichlet_group_t G, ulong h) H->P[j].e = s[k]; if (k == 0) H->P[j].g = H->q_even - 1; + else + nmod_init(&H->P[j].phi, H->q_even / 4); } j++; } @@ -210,9 +212,11 @@ dirichlet_subgroup_init(dirichlet_group_t H, const dirichlet_group_t G, ulong h) H->P[j] = G->P[k]; if (s[k] < G->P[k].e) { - ulong p = H->P[j].p; + ulong pf, p = H->P[j].p; H->P[j].e = s[k]; - nmod_init(&H->P[j].pe, n_pow(p, s[k])); + pf = n_pow(p, s[k]); + nmod_init(&H->P[j].pe, pf); + nmod_init(&H->P[j].phi, (p-1) * pf / p); } j++; } diff --git a/dirichlet/index_conrey.c b/dirichlet/index_conrey.c index 4747d4e4..18993f32 100644 --- a/dirichlet/index_conrey.c +++ b/dirichlet/index_conrey.c @@ -32,7 +32,7 @@ dirichlet_index_char(const dirichlet_group_t G, const dirichlet_char_t x) ulong j = 0; for (k = 0; k < G->num; k++) - j = j * G->P[k].phi + x->log[k]; + j = j * G->P[k].phi.n + x->log[k]; return j; } diff --git a/dirichlet/ui_pairing.c b/dirichlet/pairing.c similarity index 87% rename from dirichlet/ui_pairing.c rename to dirichlet/pairing.c index 4d572ff4..cf4fd648 100644 --- a/dirichlet/ui_pairing.c +++ b/dirichlet/pairing.c @@ -12,7 +12,7 @@ #include "dirichlet.h" ulong -dirichlet_ui_pairing(const dirichlet_group_t G, ulong m, ulong n) +dirichlet_pairing(const dirichlet_group_t G, ulong m, ulong n) { ulong x; dirichlet_char_t a, b; @@ -25,7 +25,7 @@ dirichlet_ui_pairing(const dirichlet_group_t G, ulong m, ulong n) dirichlet_char_log(a, G, m); dirichlet_char_log(b, G, n); - x = dirichlet_ui_pairing_char(G, a, b); + x = dirichlet_chi_char(G, a, b); dirichlet_char_clear(a); dirichlet_char_clear(b); diff --git a/dirichlet/test/t-chars.c b/dirichlet/test/t-chars.c index d2a79a30..06b849a8 100644 --- a/dirichlet/test/t-chars.c +++ b/dirichlet/test/t-chars.c @@ -13,140 +13,137 @@ int main() { - slong iter, bits; + slong iter; flint_rand_t state; - flint_printf("chars...."); + flint_printf("char...."); fflush(stdout); flint_randinit(state); - for (bits = 5; bits <= 30; bits += 5) + + for (iter = 0; iter < 3000 * arb_test_multiplier(); iter++) { + dirichlet_group_t G; + dirichlet_char_t x, y; + ulong q, n, k, sum; + slong ref; - for (iter = 0; iter < 50; iter++) + q = 1 + n_randint(state, 1000 * (1 + iter / 100)); + + dirichlet_group_init(G, q); + + dirichlet_char_init(x, G); + dirichlet_char_init(y, G); + + /* check group size and elements */ + dirichlet_char_one(x, G); + sum = 1; + + for (n = 1; dirichlet_char_next(x, G) >= 0; n++) + sum += x->n * x->n; + + if (FLINT_BITS == 64 || q < 1024) { - dirichlet_group_t G; - dirichlet_conrey_t x; - dirichlet_char_t chi, chi2; - ulong q, iter2; + /* use http://oeis.org/A053818 to check all elements + * are gone through */ + ref = (q % 4 == 2) ? -2 : 1; + for (k = (G->neven == 2); k < G->num; k++) + ref = - ref * G->P[k].p; + ref = ( G->phi_q * (2 * q * q + ref) ) / 6; - q = 2 + n_randint(state, 1 << bits); - - dirichlet_group_init(G, q); - dirichlet_conrey_init(x, G); - dirichlet_char_init(chi, G); - dirichlet_char_init(chi2, G); - - dirichlet_group_dlog_precompute(G, 50); - - /* check number char properties */ - for (iter2 = 0; iter2 < 100; iter2++) + if (n != G->phi_q) { - int par; - ulong m, n; - ulong order, chim1, pairing, cn, cm, cond; + flint_printf("FAIL: group size\n\n"); + flint_printf("q = %wu\n\n", q); + flint_printf("phi(q) = %wu\n\n", G->phi_q); + flint_printf("loop index = %wu\n\n", n); + abort(); + } + if (sum != ref && q > 1) + { + flint_printf("FAIL: sum test\n\n"); + flint_printf("q = %wu\n\n", q); + flint_printf("sum k^2 = %wu\n\n", ref); + flint_printf("sum obtained = %wu\n\n", sum); + abort(); + } + } - do - m = n_randint(state, q); - while (n_gcd(q, m) > 1); - - dirichlet_char(chi, G, m); - dirichlet_conrey_log(x, G, m); - dirichlet_char_conrey(chi2, G, x); - - if (!dirichlet_char_eq_deep(G, chi, chi2)) - { - flint_printf("FAIL: init char\n\n"); - flint_printf("q = %wu\n\n", q); - flint_printf("m = %wu\n\n", m); - dirichlet_char_print(G, chi); - flint_printf("\n"); - dirichlet_char_print(G, chi2); - flint_printf("\n"); - abort(); - } - - order = dirichlet_ui_order(G, m); - if (order != chi->order.n) - { - flint_printf("FAIL: order\n\n"); - flint_printf("q = %wu\n\n", q); - flint_printf("m = %wu\n\n", m); - flint_printf("order(m) = %wu\n\n", order); - flint_printf("chi->order = %wu\n\n", chi->order); - abort(); - } - - cond = dirichlet_ui_conductor(G, m); - if (cond != chi->conductor) - { - flint_printf("FAIL: conductor\n\n"); - flint_printf("q = %wu\n\n", q); - flint_printf("m = %wu\n\n", m); - flint_printf("conductor(m) = %wu\n\n", cond); - flint_printf("chi->conductor = %wu\n\n", chi->conductor); - abort(); - } - - par = dirichlet_ui_parity(G, m); - chim1 = dirichlet_ui_chi(G, chi, q - 1); - if (dirichlet_char_parity(chi) != par || par != (chim1 != 0)) - { - flint_printf("FAIL: parity\n\n"); - flint_printf("q = %wu\n\n", q); - flint_printf("m = %wu\n\n", m); - flint_printf("chi(-1) = %wu\n\n", chim1); - flint_printf("char_parity = %d", dirichlet_char_parity(chi)); - flint_printf("parity_ui = %d", par); - dirichlet_char_print(G, chi); - abort(); - } - - do - n = n_randint(state, q); - while (n_gcd(q, n) > 1); - - dirichlet_char(chi2, G, n); - pairing = dirichlet_ui_pairing(G, m, n); - cn = dirichlet_ui_chi(G, chi, n) * (G->expo / chi->order.n); - cm = dirichlet_ui_chi(G, chi2, m) * (G->expo / chi2->order.n); - - if (pairing != cn || pairing != cm) - { - flint_printf("FAIL: pairing\n\n"); - flint_printf("q = %wu\n\n", q); - flint_printf("m = %wu\n\n", m); - flint_printf("n = %wu\n\n", n); - flint_printf("chi(m,n) = %wu\n\n", pairing); - flint_printf("chi(m)(n) = %wu\n\n", cn); - flint_printf("chi(n)(m) = %wu\n\n", cm); - abort(); - } - - dirichlet_conrey_next(x, G); - dirichlet_char_next(chi, G); - dirichlet_char_conrey(chi2, G, x); - - if (!dirichlet_char_eq_deep(G, chi, chi2)) - { - flint_printf("FAIL: next char\n\n"); - flint_printf("q = %wu\n\n", q); - flint_printf("m = %wu\n\n", m); - dirichlet_char_print(G, chi); - flint_printf("\n"); - dirichlet_char_print(G, chi2); - flint_printf("\n"); - abort(); - } + if (q % 4 != 2) + { + dirichlet_char_first_primitive(x, G); + for (n = 1; dirichlet_char_next_primitive(x, G) >= 0; n++); + ref = dirichlet_number_primitive(G); + if (n != ref) + { + flint_printf("FAIL: number of primitive elements\n\n"); + flint_printf("q = %wu\n\n", q); + flint_printf("# primitive = %wu\n\n", ref); + flint_printf("loop index = %wu\n\n", n); + abort(); } - dirichlet_group_dlog_clear(G); + /* some random elements, check log and exp */ + for (n = 0; n < 30; n++) + { + slong k; + ulong m; - dirichlet_char_clear(chi); - dirichlet_char_clear(chi2); - dirichlet_conrey_clear(x); - dirichlet_group_clear(G); + for (m = 1; n_gcd(m, q) > 1; m = n_randint(state, q)); + dirichlet_char_log(x, G, m); + + if (m != dirichlet_char_exp(x, G)) + { + flint_printf("FAIL: char log and exp\n\n"); + flint_printf("q = %wu\n\n", q); + flint_printf("m = %wu\n\n", m); + flint_printf("char = "); + dirichlet_char_print(G, x); + flint_printf("\n\nnumber = %wu\n\n", x->n); + abort(); + } + + for (k = 0; k < G->num; k++) + x->log[k] = n_randint(state, G->P[k].phi.n); + + m = dirichlet_char_exp(x, G); + dirichlet_char_log(y, G, m); + + if (!dirichlet_char_eq_deep(G, x, y)) + { + flint_printf("FAIL: char exp and log\n\n"); + flint_printf("q = %wu\n\n", q); + flint_printf("char = "); + dirichlet_char_print(G, x); + flint_printf("\n\nm = %wu\n\n", m); + flint_printf("log = "); + dirichlet_char_print(G, y); + flint_printf("\n\nnumber = %wu\n\n", y->n); + abort(); + } + + dirichlet_char_next_primitive(x, G); + m = x->n; + + if (m != dirichlet_char_exp(x, G)) + { + flint_printf("FAIL: char number next primitive\n\n"); + flint_printf("q = %wu\n\n", q); + flint_printf("char = "); + dirichlet_char_print(G, y); + flint_printf(", m = %wu\n\n", y->n); + flint_printf("next primitive = "); + dirichlet_char_print(G, x); + flint_printf(", m = %wu\n\n", m); + flint_printf("exp = %wu\n\n", x->n); + abort(); + } + } } + + dirichlet_char_clear(x); + dirichlet_char_clear(y); + dirichlet_group_clear(G); } flint_randclear(state); diff --git a/dirichlet/test/t-conrey.c b/dirichlet/test/t-conrey.c deleted file mode 100644 index 525d7c00..00000000 --- a/dirichlet/test/t-conrey.c +++ /dev/null @@ -1,153 +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 . -*/ - -#include "dirichlet.h" - -int main() -{ - slong iter; - flint_rand_t state; - - flint_printf("conrey...."); - fflush(stdout); - flint_randinit(state); - - for (iter = 0; iter < 3000 * arb_test_multiplier(); iter++) - { - dirichlet_group_t G; - dirichlet_conrey_t x, y; - ulong q, n, k, sum; - slong ref; - - q = 1 + n_randint(state, 1000 * (1 + iter / 100)); - - dirichlet_group_init(G, q); - - dirichlet_conrey_init(x, G); - dirichlet_conrey_init(y, G); - - /* check group size and elements */ - dirichlet_conrey_one(x, G); - sum = 1; - - for (n = 1; dirichlet_conrey_next(x, G) >= 0; n++) - sum += x->n * x->n; - - if (FLINT_BITS == 64 || q < 1024) - { - /* use http://oeis.org/A053818 to check all elements - * are gone through */ - ref = (q % 4 == 2) ? -2 : 1; - for (k = (G->neven == 2); k < G->num; k++) - ref = - ref * G->P[k].p; - ref = ( G->phi_q * (2 * q * q + ref) ) / 6; - - if (n != G->phi_q) - { - flint_printf("FAIL: group size\n\n"); - flint_printf("q = %wu\n\n", q); - flint_printf("phi(q) = %wu\n\n", G->phi_q); - flint_printf("loop index = %wu\n\n", n); - abort(); - } - if (sum != ref && q > 1) - { - flint_printf("FAIL: sum test\n\n"); - flint_printf("q = %wu\n\n", q); - flint_printf("sum k^2 = %wu\n\n", ref); - flint_printf("sum obtained = %wu\n\n", sum); - abort(); - } - } - - if (q % 4 != 2) - { - dirichlet_conrey_first_primitive(x, G); - for (n = 1; dirichlet_conrey_next_primitive(x, G) >= 0; n++); - - ref = dirichlet_number_primitive(G); - if (n != ref) - { - flint_printf("FAIL: number of primitive elements\n\n"); - flint_printf("q = %wu\n\n", q); - flint_printf("# primitive = %wu\n\n", ref); - flint_printf("loop index = %wu\n\n", n); - abort(); - } - - /* some random elements, check log and exp */ - for (n = 0; n < 30; n++) - { - slong k; - ulong m; - - for (m = 1; n_gcd(m, q) > 1; m = n_randint(state, q)); - dirichlet_conrey_log(x, G, m); - - if (m != dirichlet_conrey_exp(x, G)) - { - flint_printf("FAIL: conrey log and exp\n\n"); - flint_printf("q = %wu\n\n", q); - flint_printf("m = %wu\n\n", m); - flint_printf("conrey = "); - dirichlet_conrey_print(G, x); - flint_printf("\n\nnumber = %wu\n\n", x->n); - abort(); - } - - for (k = 0; k < G->num; k++) - x->log[k] = n_randint(state, G->P[k].phi); - - m = dirichlet_conrey_exp(x, G); - dirichlet_conrey_log(y, G, m); - - if (!dirichlet_conrey_eq_deep(G, x, y)) - { - flint_printf("FAIL: conrey exp and log\n\n"); - flint_printf("q = %wu\n\n", q); - flint_printf("conrey = "); - dirichlet_conrey_print(G, x); - flint_printf("\n\nm = %wu\n\n", m); - flint_printf("log = "); - dirichlet_conrey_print(G, y); - flint_printf("\n\nnumber = %wu\n\n", y->n); - abort(); - } - - dirichlet_conrey_next_primitive(x, G); - m = x->n; - - if (m != dirichlet_conrey_exp(x, G)) - { - flint_printf("FAIL: conrey number next primitive\n\n"); - flint_printf("q = %wu\n\n", q); - flint_printf("conrey = "); - dirichlet_conrey_print(G, y); - flint_printf(", m = %wu\n\n", y->n); - flint_printf("next primitive = "); - dirichlet_conrey_print(G, x); - flint_printf(", m = %wu\n\n", m); - flint_printf("exp = %wu\n\n", x->n); - abort(); - } - } - } - - dirichlet_conrey_clear(x); - dirichlet_conrey_clear(y); - dirichlet_group_clear(G); - } - - flint_randclear(state); - flint_cleanup(); - flint_printf("PASS\n"); - return EXIT_SUCCESS; -} diff --git a/dirichlet/test/t-properties.c b/dirichlet/test/t-properties.c new file mode 100644 index 00000000..11cd1b79 --- /dev/null +++ b/dirichlet/test/t-properties.c @@ -0,0 +1,199 @@ +/* + 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 . +*/ + +#include "dirichlet.h" + +static void +check_eq(ulong p1, ulong p2, ulong q, ulong m, char * fct1, char * fct2) +{ + if (p1 != p2) + { + flint_printf("FAIL:\n\n"); + flint_printf("chi_%wu(%wu,.)\n\n", q, m); + flint_printf("%s = %wu\n\n", fct1, p1); + flint_printf("%s = %wu\n\n", fct2, p2); + abort(); + } +} + +static ulong +random_divisor(flint_rand_t state, const dirichlet_group_t G) +{ + ulong d; + slong k; + d = 1; + for (k = (G->neven == 2); k < G->num; k++) + d *= n_pow(G->P[k].p, n_randint(state, G->P[k].e)); + return d; +} + +int main() +{ + slong iter, bits; + flint_rand_t state; + + flint_printf("properties...."); + fflush(stdout); + flint_randinit(state); + for (bits = 5; bits <= 30; bits += 5) + { + + for (iter = 0; iter < 50; iter++) + { + dirichlet_group_t G; + dirichlet_char_t chi, psi; + ulong q, iter2; + + q = 2 + n_randint(state, 1 << bits); + + dirichlet_group_init(G, q); + dirichlet_char_init(chi, G); + dirichlet_char_init(psi, G); + + /* check number char properties */ + for (iter2 = 0; iter2 < 100; iter2++) + { + ulong m, n; + ulong p1, p2, pairing, cm, cn, q2, q3; + dirichlet_group_t G2, G3; + dirichlet_char_t chi2, chi3; + + if (iter2 == 50) + dirichlet_group_dlog_precompute(G, 5); + + /* one random character */ + do + m = n_randint(state, q); + while (n_gcd(q, m) > 1); + + dirichlet_char_log(chi, G, m); + + p1 = dirichlet_order_ui(G, m); + p2 = dirichlet_order_char(G, chi); + check_eq(p1, p2, q, m, "order m", "order chi"); + + p1 = dirichlet_conductor_ui(G, m); + p2 = dirichlet_conductor_char(G, chi); + check_eq(p1, p2, q, m, "conductor m", "conductor chi"); + + p1 = dirichlet_parity_ui(G, m); + p2 = dirichlet_parity_char(G, chi); + check_eq(p1, p2, q, m, "parity m", "parity chi"); + + p1 = dirichlet_char_is_real(G, chi); + p2 = (dirichlet_order_char(G, chi) <= 2); + check_eq(p1, p2, q, m, "is_real", "(order <= 2)"); + + /* check index */ + p1 = dirichlet_index_char(G, chi); + dirichlet_char_index(psi, G, p1); + + if (!dirichlet_char_eq_deep(G, chi, psi)) + { + flint_printf("FAIL: index\n\n"); + flint_printf("q = %wu\n\n", q); + flint_printf("m = %wu\n\n", m); + flint_printf("chi = "); dirichlet_char_print(G, chi); + flint_printf("\n\nindex(chi) = %wu\n\n", p1); + flint_printf("psi(index) = %wu\n\n", psi->n); + flint_printf("psi = "); dirichlet_char_print(G, psi); + flint_printf("\n\n"); + abort(); + } + + /* lift to higher modulus */ + q2 = q * (1 + n_randint(state, 100)); + + dirichlet_group_init(G2, q2); + dirichlet_char_init(chi2, G2); + dirichlet_char_lift(chi2, G2, chi, G); + + p1 = dirichlet_conductor_char(G, chi); + p2 = dirichlet_conductor_char(G2, chi2); + check_eq(p1, p2, q, m, "conductor chi", "conductor lift"); + + p1 = dirichlet_order_char(G, chi); + p2 = dirichlet_order_char(G2, chi2); + check_eq(p1, p2, q, m, "order chi", "order lift"); + + /* and lower */ + + dirichlet_char_lower(psi, G, chi2, G2); + if (!dirichlet_char_eq_deep(G, chi, psi)) + { + flint_printf("FAIL: lift and lower back\n\n"); + flint_printf("q = %wu\n\nchi = ", q); + dirichlet_char_print(G, chi); + flint_printf("\n\nq2 = %wu\n\nchi2 = ", q2); + dirichlet_char_print(G2, chi2); + flint_printf("\n\nq = %wu\n\npsi = ", q); + dirichlet_char_print(G, psi); + flint_printf("\n\n"); + abort(); + } + + q3 = dirichlet_conductor_char(G, chi) * random_divisor(state, G); + q3 = n_gcd(q, q3); + + dirichlet_group_init(G3, q3); + dirichlet_char_init(chi3, G3); + dirichlet_char_lower(chi3, G3, chi2, G2); + + p1 = dirichlet_conductor_char(G, chi); + p2 = dirichlet_conductor_char(G3, chi3); + check_eq(p1, p2, q, m, "conductor chi", "conductor lower"); + + p1 = dirichlet_order_char(G, chi); + p2 = dirichlet_order_char(G3, chi3); + check_eq(p1, p2, q, m, "order chi", "order lower"); + + dirichlet_char_clear(chi3); + dirichlet_group_clear(G3); + dirichlet_char_clear(chi2); + dirichlet_group_clear(G2); + + /* another random character */ + do + n = n_randint(state, q); + while (n_gcd(q, n) > 1); + + dirichlet_char_log(psi, G, n); + pairing = dirichlet_pairing(G, m, n); + cn = dirichlet_chi(G, chi, n); + cm = dirichlet_chi(G, psi, m); + + if (pairing != cn || pairing != cm) + { + flint_printf("FAIL: pairing\n\n"); + flint_printf("q = %wu\n\n", q); + flint_printf("m = %wu\n\n", m); + flint_printf("n = %wu\n\n", n); + flint_printf("chi(m,n) = %wu\n\n", pairing); + flint_printf("chi(m)(n) = %wu\n\n", cn); + flint_printf("chi(n)(m) = %wu\n\n", cm); + abort(); + } + + } + + dirichlet_group_dlog_clear(G); + + dirichlet_char_clear(chi); + dirichlet_char_clear(psi); + dirichlet_group_clear(G); + } + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/acb_dirichlet/test/t-vec.c b/dirichlet/test/t-vec.c similarity index 51% rename from acb_dirichlet/test/t-vec.c rename to dirichlet/test/t-vec.c index d3cea907..edb2ad08 100644 --- a/acb_dirichlet/test/t-vec.c +++ b/dirichlet/test/t-vec.c @@ -9,7 +9,7 @@ (at your option) any later version. See . */ -#include "acb_dirichlet.h" +#include "dirichlet.h" static ulong vec_diff(ulong * v, ulong * ref, ulong nv) @@ -31,41 +31,55 @@ int main() for (q = 2; q < 600; q ++) { dirichlet_group_t G; - dirichlet_conrey_t x; dirichlet_char_t chi; ulong * v1, * v2, nv, k; dirichlet_group_init(G, q); - dirichlet_conrey_init(x, G); dirichlet_char_init(chi, G); nv = 100; v1 = flint_malloc(nv * sizeof(ulong)); v2 = flint_malloc(nv * sizeof(ulong)); - dirichlet_conrey_one(x, G); + dirichlet_char_one(chi, G); do { - dirichlet_char_conrey(chi, G, x); + ulong order; - dirichlet_ui_chi_vec_loop(v1, G, chi, nv); - dirichlet_ui_chi_vec_primeloop(v2, G, chi, nv); + dirichlet_chi_vec_loop(v1, G, chi, nv); + dirichlet_chi_vec_primeloop(v2, G, chi, nv); if ((k = vec_diff(v1, v2, nv))) { - flint_printf("FAIL: chi(%wu,%wu) = %wu != chi(%wu,%wu) = %wu mod %wu (modulus = %wu)\n", - chi->x->n, k, v1[k], chi->x->n, k, v2[k], chi->order, q); + flint_printf("FAIL: chi_%wu(%wu,%wu) [mod %wu]\n", q, chi->n, k, G->expo); + flint_printf("vec_loop -> %wu\n", v1[k]); + flint_printf("vec_primeloop -> %wu\n", v2[k]); + flint_printf("pairing = %wu\n", dirichlet_pairing(G, chi->n, k)); abort(); } - } while (dirichlet_conrey_next(x, G) >= 0); + order = dirichlet_order_char(G, chi); + dirichlet_chi_vec_loop_order(v1, G, chi, order, nv); + dirichlet_chi_vec_primeloop_order(v2, G, chi, order, nv); + + if ((k = vec_diff(v1, v2, nv))) + { + flint_printf("FAIL: chi_%wu(%wu,%wu) [mod %wu]\n", q, chi->n, k, order); + flint_printf("vec_loop -> %wu\n", v1[k]); + flint_printf("vec_primeloop -> %wu\n", v2[k]); + flint_printf("pairing = %wu mod %wu\n", dirichlet_pairing(G, chi->n, k), G->expo); + abort(); + } + + + + } while (dirichlet_char_next(chi, G) >= 0); flint_free(v1); flint_free(v2); dirichlet_group_clear(G); dirichlet_char_clear(chi); - dirichlet_conrey_clear(x); } flint_cleanup(); diff --git a/dirichlet/ui_chi_conrey.c b/dirichlet/ui_chi_conrey.c deleted file mode 100644 index 44f93fcc..00000000 --- a/dirichlet/ui_chi_conrey.c +++ /dev/null @@ -1,24 +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 . -*/ - -#include "dirichlet.h" - -ulong -dirichlet_ui_chi_char(const dirichlet_group_t G, const dirichlet_fullchar_t chi, const dirichlet_char_t x) -{ - ulong v = 0, k; - - /* TODO: nmod_addmul? */ - for (k = 0; k < G->num; k++) - v = nmod_add(v, nmod_mul(chi->expo[k], x->log[k], chi->order), chi->order); - - return v; -} diff --git a/dirichlet/ui_chi_vec.c b/dirichlet/ui_chi_vec.c deleted file mode 100644 index 92ca6817..00000000 --- a/dirichlet/ui_chi_vec.c +++ /dev/null @@ -1,21 +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 . -*/ - -#include "dirichlet.h" - -void -dirichlet_ui_chi_vec(ulong *v, const dirichlet_group_t G, const dirichlet_fullchar_t chi, slong nv) -{ - if (2 * nv > G->q) - dirichlet_ui_chi_vec_loop(v, G, chi, nv); - else - dirichlet_ui_chi_vec_primeloop(v, G, chi, nv); -} diff --git a/dirichlet/ui_chi_vec_primeloop.c b/dirichlet/ui_chi_vec_primeloop.c deleted file mode 100644 index 31dd4ca7..00000000 --- a/dirichlet/ui_chi_vec_primeloop.c +++ /dev/null @@ -1,72 +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 . -*/ - -#include "dirichlet.h" - -static void -chi_vec_evenpart(ulong *v, const dirichlet_group_t G, const dirichlet_fullchar_t chi, slong nv) -{ - - ulong c3, c4, x; - - if (G->neven >= 1 && (c3 = chi->expo[0])) - { - for (x = 3; x < nv; x += 4) - v[x] = c3; - } - if (G->neven == 2 && (c4 = chi->expo[1])) - { - ulong g, vx, xp; - nmod_t pe; - vx = c4; - pe = G->P[1].pe; - g = G->P[1].g; - - for (x = g; x > 1;) - { - - for (xp = x; xp < nv; xp += pe.n) - v[xp] = nmod_add(v[xp], vx, chi->order); - - for (xp = pe.n - x; xp < nv; xp += pe.n) - v[xp] = nmod_add(v[xp], vx, chi->order); - - x = nmod_mul(x, g, pe); - vx = nmod_add(vx, c4, chi->order); - } - } -} - -/* loop over primary components */ -void -dirichlet_ui_chi_vec_primeloop(ulong *v, const dirichlet_group_t G, const dirichlet_fullchar_t chi, slong nv) -{ - slong k, l; - - for (k = 0; k < nv; k++) - v[k] = 0; - - if (G->neven) - chi_vec_evenpart(v, G, chi, nv); - - for (l = G->neven; l < G->num; l++) - { - dirichlet_prime_group_struct P = G->P[l]; - - /* FIXME: there may be some precomputed dlog in P if needed */ - if (P.dlog == NULL) - dlog_vec_add(v, nv, P.g, chi->expo[l], P.pe, P.phi, chi->order); - else - dlog_vec_add_precomp(v, nv, P.dlog, P.g, chi->expo[l], P.pe, P.phi, chi->order); - - } - dirichlet_ui_vec_set_null(v, G, nv); -} diff --git a/dirichlet/ui_conductor.c b/dirichlet/ui_conductor.c index 4c2d1295..63d68273 100644 --- a/dirichlet/ui_conductor.c +++ b/dirichlet/ui_conductor.c @@ -12,7 +12,7 @@ #include "dirichlet.h" ulong -dirichlet_ui_conductor(const dirichlet_group_t G, ulong a) +dirichlet_conductor_ui(const dirichlet_group_t G, ulong a) { slong k; ulong ap, cond; diff --git a/dirichlet/ui_order.c b/dirichlet/ui_order.c index 6c58db89..c812d015 100644 --- a/dirichlet/ui_order.c +++ b/dirichlet/ui_order.c @@ -31,7 +31,7 @@ nmod_order_precomp(ulong a, nmod_t mod, ulong expo, n_factor_t fac) } ulong -dirichlet_ui_order(const dirichlet_group_t G, ulong a) +dirichlet_order_ui(const dirichlet_group_t G, ulong a) { n_factor_t fac; diff --git a/dirichlet/ui_pairing_conrey.c b/dirichlet/ui_pairing_conrey.c deleted file mode 100644 index 3488f873..00000000 --- a/dirichlet/ui_pairing_conrey.c +++ /dev/null @@ -1,26 +0,0 @@ -/* - Copyright (C) 2015 Jonathan Bober - Copyright (C) 2016 Fredrik Johansson - 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 . -*/ - -#include "dirichlet.h" - -ulong -dirichlet_ui_pairing_char(const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b) -{ - ulong x, k; - x = 0; - - for (k = 0; k < G->num; k++) - x = n_addmod(x, G->PHI[k] * n_mulmod2(a->log[k], b->log[k], G->P[k].phi), G->expo); - - return x; -} diff --git a/dirichlet/ui_parity.c b/dirichlet/ui_parity.c index 6312dc81..44191e25 100644 --- a/dirichlet/ui_parity.c +++ b/dirichlet/ui_parity.c @@ -12,7 +12,7 @@ #include "dirichlet.h" int -dirichlet_ui_parity(const dirichlet_group_t G, ulong a) +dirichlet_parity_ui(const dirichlet_group_t G, ulong a) { int par; diff --git a/dirichlet/ui_vec_set_null.c b/dirichlet/ui_vec_set_null.c index 834a5c0b..65a05aef 100644 --- a/dirichlet/ui_vec_set_null.c +++ b/dirichlet/ui_vec_set_null.c @@ -12,7 +12,7 @@ #include "dirichlet.h" void -dirichlet_ui_vec_set_null(ulong *v, const dirichlet_group_t G, slong nv) +dirichlet_vec_set_null(ulong *v, const dirichlet_group_t G, slong nv) { slong k, l; if (G->q_even > 1)