hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 17:01:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

remove all conrey

Browse source
This commit is contained in:
Pascal 2016-10-08 17:22:19 +02:00
parent e29ca5d50e
commit 00118eb55f
67 changed files with 757 additions and 1053 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
acb_dirichlet.h
View file

@ -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);

4
acb_dirichlet/chi.c
View file

@ -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);
}

7
acb_dirichlet/chi_vec.c
View file

@ -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++)

26
acb_dirichlet/gauss_sum.c
View file

@ -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);
}

12
acb_dirichlet/gauss_sum_factor.c
View file

@ -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);

8
acb_dirichlet/gauss_sum_order2.c
View file

@ -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);
}
}

7
acb_dirichlet/gauss_sum_theta.c
View file

@ -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);

20
acb_dirichlet/jacobi_sum.c
View file

@ -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);

17
acb_dirichlet/jacobi_sum_factor.c
View file

@ -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)

2
acb_dirichlet/jacobi_sum_gauss.c
View file

@ -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);

33
acb_dirichlet/jacobi_sum_naive.c
View file

@ -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);

25
acb_dirichlet/l_hurwitz.c
View file

@ -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);

2
acb_dirichlet/pairing.c
View file

@ -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

4
acb_dirichlet/pairing_conrey.c
View file

@ -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

4
acb_dirichlet/root_number.c
View file

@ -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);

12
acb_dirichlet/test/t-chi.c
View file

@ -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();
}

21
acb_dirichlet/test/t-gauss.c
View file

@ -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();

7
acb_dirichlet/test/t-jacobi.c
View file

@ -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();
}

2
acb_dirichlet/test/t-l.c
View file

@ -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++)
{

18
acb_dirichlet/test/t-thetanull.c
View file

@ -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);

27
acb_dirichlet/theta_arb.c
View file

@ -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);

2
acb_dirichlet/ui_theta_arb.c
View file

@ -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);

124
dirichlet.h
View file

@ -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
}

19
dirichlet/char.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

19
dirichlet/char_clear.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#include "dirichlet.h"
void
dirichlet_fullchar_clear(dirichlet_fullchar_t chi)
{
dirichlet_char_clear(chi->x);
flint_free(chi->expo);
}

32
dirichlet/char_conrey.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

38
dirichlet/char_eq_deep.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

20
dirichlet/char_first_primitive.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

26
dirichlet/char_init.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

19
dirichlet/char_mul.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

22
dirichlet/char_next.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

22
dirichlet/char_next_primitive.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

46
dirichlet/char_normalize.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

24
dirichlet/char_primitive.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

23
dirichlet/char_print.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

4
dirichlet/ui_chi.c → dirichlet/chi.c
View file

@ -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;

19
dirichlet/char_one.c → dirichlet/chi_char.c
View file

@ -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;
}

30
dirichlet/chi_vec.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

26
dirichlet/ui_chi_vec_loop.c → dirichlet/chi_vec_loop.c
View file

@ -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);
}

82
dirichlet/chi_vec_primeloop.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

2
dirichlet/conrey_conductor.c
View file

@ -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;

6
dirichlet/conrey_index.c
View file

@ -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);

36
dirichlet/conrey_lift.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

2
dirichlet/conrey_log.c
View file

@ -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
{

42
dirichlet/conrey_lower.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

2
dirichlet/conrey_mul.c
View file

@ -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);
}

2
dirichlet/conrey_next.c
View file

@ -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;
}

4
dirichlet/conrey_next_primitive.c
View file

@ -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

2
dirichlet/conrey_order.c
View file

@ -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;

2
dirichlet/conrey_parity.c
View file

@ -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++)

2
dirichlet/conrey_pow.c
View file

@ -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);
}

47
dirichlet/conrey_primitive.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}
}
}
}

22
dirichlet/group_init.c
View file

@ -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++;
}

2
dirichlet/index_conrey.c
View file

@ -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;
}

4
dirichlet/ui_pairing.c → dirichlet/pairing.c
View file

@ -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);

235
dirichlet/test/t-chars.c
View file

@ -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);

153
dirichlet/test/t-conrey.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

199
dirichlet/test/t-properties.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

36
acb_dirichlet/test/t-vec.c → dirichlet/test/t-vec.c
View file

@ -9,7 +9,7 @@
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#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();

24
dirichlet/ui_chi_conrey.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

21
dirichlet/ui_chi_vec.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

72
dirichlet/ui_chi_vec_primeloop.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

2
dirichlet/ui_conductor.c
View file

@ -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;

2
dirichlet/ui_order.c
View file

@ -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;

26
dirichlet/ui_pairing_conrey.c
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

2
dirichlet/ui_parity.c
View file

@ -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;

2
dirichlet/ui_vec_set_null.c
View file

@ -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)