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

Merge remote-tracking branch 'pascal/dirichlet' into dirichlet

Browse source
This commit is contained in:
Fredrik Johansson 2016-09-06 20:40:29 +02:00
commit 2e8ab6f3d8
136 changed files with 8197 additions and 242 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
Makefile.in
View file

@ -14,7 +14,7 @@ AT=@
BUILD_DIRS = fmpr arf mag arb arb_mat arb_poly arb_calc acb acb_mat acb_poly \
acb_calc acb_hypgeom acb_modular acb_dirichlet arb_hypgeom bernoulli hypgeom \
fmpz_extras bool_mat partitions \
fmpz_extras bool_mat partitions dlog \
$(EXTRA_BUILD_DIRS)
TEMPLATE_DIRS =

283
acb_dirichlet.h
View file

@ -1,55 +1,282 @@
/*
/*=============================================================================
*
*
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2015 Jonathan Bober
Copyright (C) 2016 Fredrik Johansson
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/>.
*/
******************************************************************************/
#ifndef ACB_DIRICHLET_H
#define ACB_DIRICHLET_H
#ifdef ACB_DIRICHLET_INLINES_C
#define ACB_DIRICHLET_INLINE
#else
#define ACB_DIRICHLET_INLINE static __inline__
#endif
#include "acb.h"
#include "dlog.h"
#ifdef __cplusplus
extern "C" {
#endif
typedef struct
{
ulong q; /* modulus */
ulong q_even; /* even part of modulus */
ulong q_odd; /* odd part of modulus */
ulong phi_q; /* phi(q) */
ulong phi_q_odd; /* phi(q_odd) */
slong num; /* number of odd prime factors */
ulong * primes; /* odd primes p[k] */
ulong * exponents; /* exponents e[k] */
ulong * generators; /* generator for each odd prime p[k] */
ulong * PHI; /* PHI(k) = phi(q_odd) / phi(p[k]^e[k]) */
}
acb_dirichlet_group_struct;
typedef acb_dirichlet_group_struct acb_dirichlet_group_t[1];
void _acb_dirichlet_euler_product_real_ui(arb_t res, ulong s,
const signed char * chi, int mod, int reciprocal, slong prec);
void acb_dirichlet_eta(acb_t res, const acb_t s, slong prec);
/* 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) */
ulong g; /* conrey generator */
dlog_precomp_struct * dlog; /* precomputed data for discrete log mod p^e */
}
acb_dirichlet_prime_group_struct;
typedef struct
{
ulong q; /* modulus */
ulong q_even; /* even part of modulus */
nmod_t mod; /* modulus with precomputed inverse */
ulong rad_q; /* radical = product of odd primes */
ulong phi_q; /* phi(q) = group size */
slong neven; /* number of even components (in 0,1,2)*/
slong num; /* number of prime components (even + odd) */
ulong expo; /* exponent = largest order in G */
acb_dirichlet_prime_group_struct * P;
ulong * generators; /* generators lifted mod q */
ulong * PHI; /* PHI(k) = expo / phi(k) */
}
acb_dirichlet_group_struct;
typedef acb_dirichlet_group_struct acb_dirichlet_group_t[1];
ACB_DIRICHLET_INLINE ulong
acb_dirichlet_group_size(const acb_dirichlet_group_t G)
{
return G->phi_q;
}
void acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q);
void acb_dirichlet_subgroup_init(acb_dirichlet_group_t H, const acb_dirichlet_group_t G, ulong h);
void acb_dirichlet_group_clear(acb_dirichlet_group_t G);
void acb_dirichlet_group_dlog_precompute(acb_dirichlet_group_t G, ulong num);
void acb_dirichlet_chi(acb_t res, const acb_dirichlet_group_t G, ulong m, ulong n, slong prec);
/* properties of elements without log */
ulong acb_dirichlet_number_primitive(const acb_dirichlet_group_t G);
ulong acb_dirichlet_ui_conductor(const acb_dirichlet_group_t G, ulong a);
int acb_dirichlet_ui_parity(const acb_dirichlet_group_t G, ulong a);
ulong acb_dirichlet_ui_order(const acb_dirichlet_group_t G, ulong a);
/* elements of the group, keep both number and log */
typedef struct
{
ulong n; /* number */
ulong * log; /* s.t. prod generators[k]^log[k] = number */
}
acb_dirichlet_conrey_struct;
typedef acb_dirichlet_conrey_struct acb_dirichlet_conrey_t[1];
void acb_dirichlet_conrey_init(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G);
void acb_dirichlet_conrey_clear(acb_dirichlet_conrey_t x);
void acb_dirichlet_conrey_print(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x);
ACB_DIRICHLET_INLINE void
acb_dirichlet_conrey_copy(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t y)
{
slong k;
x->n = y->n;
for (k = 0; k < G->num; k++)
x->log[k] = y->log[k];
}
int acb_dirichlet_conrey_eq(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x, const acb_dirichlet_conrey_t y);
int acb_dirichlet_conrey_parity(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x);
ulong acb_dirichlet_conrey_conductor(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x);
void acb_dirichlet_conrey_log(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G, ulong m);
ulong acb_dirichlet_conrey_exp(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G);
void acb_dirichlet_conrey_one(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G);
void acb_dirichlet_conrey_first_primitive(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G);
int acb_dirichlet_conrey_next(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G);
int acb_dirichlet_conrey_next_primitive(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G);
void acb_dirichlet_conrey_mul(acb_dirichlet_conrey_t c, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t a, const acb_dirichlet_conrey_t b);
void acb_dirichlet_conrey_pow(acb_dirichlet_conrey_t c, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t a, ulong n);
void acb_dirichlet_conrey_primitive(acb_dirichlet_conrey_t y, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x, ulong cond);
#define ACB_DIRICHLET_CHI_NULL UWORD_MAX
ulong acb_dirichlet_ui_pairing_conrey(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t a, const acb_dirichlet_conrey_t b);
ulong acb_dirichlet_ui_pairing(const acb_dirichlet_group_t G, ulong m, ulong n);
void acb_dirichlet_pairing_conrey(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t a, const acb_dirichlet_conrey_t b, slong prec);
void acb_dirichlet_pairing(acb_t res, const acb_dirichlet_group_t G, ulong m, ulong n, slong prec);
/* introducing character type */
/* character = reduced exponents, keep order, number and conductor */
typedef struct
{
ulong q; /* modulus */
nmod_t order; /* order */
acb_dirichlet_conrey_t x;
ulong * expo; /* reduced exponents ( log[k] * PHI[k] / gcd( ) ) */
int parity; /* 0 for even char, 1 for odd */
ulong conductor;
}
acb_dirichlet_char_struct;
typedef acb_dirichlet_char_struct acb_dirichlet_char_t[1];
ACB_DIRICHLET_INLINE ulong
acb_dirichlet_char_order(const acb_dirichlet_char_t chi)
{
return chi->order.n;
}
ACB_DIRICHLET_INLINE ulong
acb_dirichlet_char_conductor(const acb_dirichlet_char_t chi)
{
return chi->conductor;
}
ACB_DIRICHLET_INLINE int
acb_dirichlet_char_parity(const acb_dirichlet_char_t chi)
{
return chi->parity;
}
void acb_dirichlet_char_init(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
void acb_dirichlet_char_clear(acb_dirichlet_char_t chi);
void acb_dirichlet_char_print(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi);
int acb_dirichlet_char_eq(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2);
void acb_dirichlet_char(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, ulong n);
void acb_dirichlet_char_conrey(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x);
void acb_dirichlet_char_set_expo(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
void acb_dirichlet_char_normalize(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
void acb_dirichlet_char_denormalize(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
void acb_dirichlet_char_mul(acb_dirichlet_char_t chi12, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2);
void acb_dirichlet_char_primitive(acb_dirichlet_char_t chi0, const acb_dirichlet_group_t G0, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi);
/*ulong acb_dirichlet_char_conductor(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi);*/
void acb_dirichlet_char_one(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
void acb_dirichlet_char_first_primitive(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
int acb_dirichlet_char_next(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
int acb_dirichlet_char_next_primitive(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
ulong acb_dirichlet_ui_chi_conrey(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const acb_dirichlet_conrey_t x);
ulong acb_dirichlet_ui_chi(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n);
void acb_dirichlet_ui_vec_set_null(ulong *v, const acb_dirichlet_group_t G, slong nv);
void acb_dirichlet_ui_chi_vec_loop(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv);
void acb_dirichlet_ui_chi_vec_primeloop(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv);
void acb_dirichlet_ui_chi_vec(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv);
/* precompute powers of a root of unity */
typedef struct
{
ulong order;
acb_t z;
slong size;
slong depth;
acb_ptr * Z;
}
acb_dirichlet_powers_struct;
typedef acb_dirichlet_powers_struct acb_dirichlet_powers_t[1];
void _acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong size, slong depth, slong prec);
void acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong num, slong prec);
void acb_dirichlet_powers_clear(acb_dirichlet_powers_t t);
void acb_dirichlet_power(acb_t z, const acb_dirichlet_powers_t t, ulong n, slong prec);
void acb_dirichlet_nth_root(acb_t res, ulong order, slong prec);
void acb_dirichlet_chi(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n, slong prec);
void acb_dirichlet_chi_vec(acb_ptr v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv, slong prec);
void acb_dirichlet_arb_quadratic_powers(arb_ptr v, slong nv, const arb_t x, slong prec);
void acb_dirichlet_qseries_eval_arb(acb_t res, acb_srcptr a, const arb_t x, slong len, slong prec);
ulong acb_dirichlet_theta_length_d(ulong q, double x, slong prec);
ulong acb_dirichlet_theta_length(ulong q, const arb_t x, slong prec);
void acb_dirichlet_arb_theta_naive(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec);
void acb_dirichlet_arb_theta_smallorder(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec);
void acb_dirichlet_chi_theta_arb(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t t, slong prec);
void acb_dirichlet_ui_theta_arb(acb_t res, const acb_dirichlet_group_t G, ulong a, const arb_t t, slong prec);
void acb_dirichlet_gauss_sum_naive(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec);
void acb_dirichlet_gauss_sum_theta(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec);
void acb_dirichlet_gauss_sum_factor(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec);
void acb_dirichlet_gauss_sum(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec);
void acb_dirichlet_si_poly_evaluate(acb_t res, slong * v, slong len, const acb_t z, slong prec);
void acb_dirichlet_jacobi_sum_naive(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2, slong prec);
void acb_dirichlet_jacobi_sum(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2, slong prec);
void acb_dirichlet_l_hurwitz(acb_t res, const acb_t s, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec);
void acb_dirichlet_l_vec_hurwitz(acb_ptr res, const acb_t s, const acb_dirichlet_group_t G, slong prec);
/* Discrete Fourier Transform */
void acb_dirichlet_vec_nth_roots(acb_ptr z, slong len, slong prec);
void _acb_dirichlet_dft_pol(acb_ptr w, acb_srcptr v, acb_srcptr z, slong len, slong prec);
void acb_dirichlet_dft_pol(acb_ptr w, acb_srcptr v, slong len, slong prec);
void acb_dirichlet_dft_fast(acb_ptr w, acb_srcptr v, slong len, slong prec);
void acb_dirichlet_dft_prod(acb_ptr w, acb_srcptr v, slong * cyc, slong num, slong prec);
void acb_dirichlet_dft_conrey(acb_ptr w, acb_srcptr v, const acb_dirichlet_group_t G, slong prec);
void acb_dirichlet_dft(acb_ptr w, acb_srcptr v, const acb_dirichlet_group_t G, slong prec);
ACB_DIRICHLET_INLINE void
acb_vec_printd(acb_srcptr vec, slong len, slong digits)
{
slong i;
for (i = 0; i < len; i++)
acb_printd(vec + i, digits), flint_printf("\n");
}
#ifdef __cplusplus
}
#endif
#endif

38
acb_dirichlet/arb_quadratic_powers.c Normal file
View file

@ -0,0 +1,38 @@
/*
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 "acb_dirichlet.h"
/* TODO: BSGS can reduce to nv mul */
void
acb_dirichlet_arb_quadratic_powers(arb_ptr v, slong nv, const arb_t x, slong prec)
{
slong i;
arb_t dx, x2;
arb_init(dx);
arb_init(x2);
arb_set(dx, x);
arb_mul(x2, x, x, prec);
for (i = 0; i < nv; i++)
{
if (i == 0)
arb_one(v + i);
else if (i == 1)
arb_set_round(v + i, x, prec);
else
{
arb_mul(dx, dx, x2, prec);
arb_mul(v + i, v + i - 1, dx, prec);
}
}
arb_clear(dx);
arb_clear(x2);
}

50
acb_dirichlet/arb_theta_naive.c Normal file
View file

@ -0,0 +1,50 @@
/*
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 "acb_dirichlet.h"
#include "acb_poly.h"
void
acb_dirichlet_arb_theta_naive(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec)
{
slong k;
arb_t xk2, dx, x2;
acb_t zk;
arb_init(xk2);
arb_init(dx);
arb_init(x2);
acb_init(zk);
arb_set(dx, x);
arb_set(xk2, dx);
arb_mul(x2, dx, dx, prec);
acb_set_arb(res, xk2);
/* TODO: reduce prec */
for (k = 2; k < len; k++)
{
arb_mul(dx, dx, x2, prec);
arb_mul(xk2, xk2, dx, prec);
if (a[k] != ACB_DIRICHLET_CHI_NULL)
{
acb_dirichlet_power(zk, z, a[k], prec);
if (parity)
acb_mul_si(zk, zk, k, prec);
acb_addmul_arb(res, zk, xk2, prec);
}
}
arb_clear(xk2);
arb_clear(x2);
arb_clear(dx);
acb_clear(zk);
}

63
acb_dirichlet/arb_theta_smallorder.c Normal file
View file

@ -0,0 +1,63 @@
/*
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 "acb_dirichlet.h"
#include "acb_poly.h"
/* small order, multiply by chi at the end */
void
acb_dirichlet_arb_theta_smallorder(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec)
{
slong k;
ulong order = z->order;
arb_t xk2, kxk2, dx, x2;
acb_ptr t;
arb_init(xk2);
arb_init(dx);
arb_init(x2);
arb_init(kxk2);
arb_set(dx, x);
arb_set(xk2, x);
arb_mul(x2, x, x, prec);
t = _acb_vec_init(order);
_acb_vec_zero(t, order);
arb_set(acb_realref(t + 0), xk2);
/* TODO: reduce precision at each step */
for (k = 2; k < len; k++)
{
arb_mul(dx, dx, x2, prec);
arb_mul(xk2, xk2, dx, prec);
if (a[k] != ACB_DIRICHLET_CHI_NULL)
{
if (parity)
{
arb_mul_si(kxk2, xk2, k, prec);
arb_add(acb_realref(t + a[k]), acb_realref(t + a[k]), kxk2, prec);
}
else
{
arb_add(acb_realref(t + a[k]), acb_realref(t + a[k]), xk2, prec);
}
}
}
/* now Hörner */
_acb_poly_evaluate(res, t, order, z->z, prec);
_acb_vec_clear(t, order);
arb_clear(xk2);
arb_clear(x2);
arb_clear(dx);
}

19
acb_dirichlet/char.c Normal file
View file

@ -0,0 +1,19 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_char(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, ulong n)
{
acb_dirichlet_conrey_log(chi->x, G, n);
acb_dirichlet_char_conrey(chi, G, NULL);
}

19
acb_dirichlet/char_clear.c Normal file
View file

@ -0,0 +1,19 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_char_clear(acb_dirichlet_char_t chi)
{
acb_dirichlet_conrey_clear(chi->x);
flint_free(chi->expo);
}

32
acb_dirichlet/char_conrey.c Normal file
View file

@ -0,0 +1,32 @@
/*
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 "acb_dirichlet.h"
/* char n has exponents = log[k]*PHI[k] / gcd and order expo / gcd
* so that log = expo[k] */
void
acb_dirichlet_char_conrey(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x)
{
/* assume chi->x already set if x == NULL */
if (x == NULL)
x = chi->x;
else
acb_dirichlet_conrey_copy(chi->x, G, x);
chi->q = G->q;
chi->parity = acb_dirichlet_conrey_parity(G, x);
chi->conductor = acb_dirichlet_conrey_conductor(G, x);
acb_dirichlet_char_set_expo(chi, G);
/* optional: divide by gcd to obtain true order */
acb_dirichlet_char_normalize(chi, G);
}

38
acb_dirichlet/char_eq.c Normal file
View file

@ -0,0 +1,38 @@
/*
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 "acb_dirichlet.h"
int
acb_dirichlet_char_eq(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2)
{
acb_dirichlet_conrey_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 (!acb_dirichlet_conrey_eq(G, chi1->x, chi2->x))
return 0;
x->n = y->n = 1;
x->log = chi1->expo;
y->log = chi2->expo;
if (!acb_dirichlet_conrey_eq(G, x, y))
return 0;
return 1;
}

20
acb_dirichlet/char_first_primitive.c Normal file
View file

@ -0,0 +1,20 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_char_first_primitive(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
{
acb_dirichlet_conrey_first_primitive(chi->x, G);
acb_dirichlet_char_conrey(chi, G, NULL);
acb_dirichlet_char_normalize(chi, G);
}

19
acb_dirichlet/char_init.c Normal file
View file

@ -0,0 +1,19 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_char_init(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
{
acb_dirichlet_conrey_init(chi->x, G);
chi->expo = flint_malloc(G->num * sizeof(ulong));
}

19
acb_dirichlet/char_mul.c Normal file
View file

@ -0,0 +1,19 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_char_mul(acb_dirichlet_char_t chi12, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2)
{
acb_dirichlet_conrey_mul(chi12->x, G, chi1->x, chi2->x);
acb_dirichlet_char_conrey(chi12, G, NULL);
}

22
acb_dirichlet/char_next.c Normal file
View file

@ -0,0 +1,22 @@
/*
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 "acb_dirichlet.h"
int
acb_dirichlet_char_next(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
{
int k;
k = acb_dirichlet_conrey_next(chi->x, G);
acb_dirichlet_char_conrey(chi, G, NULL);
/* return last index modified */
return k;
}

22
acb_dirichlet/char_next_primitive.c Normal file
View file

@ -0,0 +1,22 @@
/*
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 "acb_dirichlet.h"
int
acb_dirichlet_char_next_primitive(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
{
int k;
k = acb_dirichlet_conrey_next_primitive(chi->x, G);
acb_dirichlet_char_conrey(chi, G, NULL);
/* return last index modified */
return k;
}

45
acb_dirichlet/char_normalize.c Normal file
View file

@ -0,0 +1,45 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_char_set_expo(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
{
slong k;
for (k = 0; k < G->num; k++)
chi->expo[k] = (chi->x->log[k] * G->PHI[k]) % G->expo;
}
void
acb_dirichlet_char_normalize(acb_dirichlet_char_t chi, const acb_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
acb_dirichlet_char_denormalize(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
{
ulong k, g;
g = G->expo / chi->order.n;
for (k = 0; k < G->num; k++)
chi->expo[k] *= g;
}

23
acb_dirichlet/char_one.c Normal file
View file

@ -0,0 +1,23 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_char_one(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
{
acb_dirichlet_conrey_one(chi->x, G);
chi->q = G->q;
chi->conductor = 1;
chi->parity = 0;
acb_dirichlet_char_set_expo(chi, G);
acb_dirichlet_char_normalize(chi, G);
}

24
acb_dirichlet/char_primitive.c Normal file
View file

@ -0,0 +1,24 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_char_primitive(acb_dirichlet_char_t chi0, const acb_dirichlet_group_t G0, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi)
{
chi0->q = chi->conductor;
chi0->parity = chi->parity;
chi0->conductor = chi->conductor;
acb_dirichlet_conrey_primitive(chi0->x, G, chi->x, chi->conductor);
acb_dirichlet_char_set_expo(chi0, G0);
/* optional: divide by gcd to obtain true order */
acb_dirichlet_char_normalize(chi0, G0);
}

23
acb_dirichlet/char_print.c Normal file
View file

@ -0,0 +1,23 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_char_print(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi)
{
acb_dirichlet_conrey_t x;
flint_printf("chi_%wu(%wu,.) of order %wu, parity %wd, index ", G->q, chi->x->n, chi->order, chi->parity);
acb_dirichlet_conrey_print(G, chi->x);
flint_printf(" and exponents ");
x->log = chi->expo;
acb_dirichlet_conrey_print(G, x);
}

116
acb_dirichlet/chi.c
View file

@ -1,6 +1,7 @@
/*
Copyright (C) 2015 Jonathan Bober
Copyright (C) 2016 Fredrik Johansson
Copyright (C) 2016 Pascal Molin
This file is part of Arb.
@ -12,116 +13,19 @@
#include "acb_dirichlet.h"
/* todo: modular arithmetic */
static ulong
chi_odd_exponent(const acb_dirichlet_group_t G, ulong m, ulong n)
{
ulong x, k, pk, gk, logm, logn;
x = 0;
for (k = 0; k < G->num; k++)
{
pk = n_pow(G->primes[k], G->exponents[k]);
gk = G->generators[k] % pk;
logm = n_discrete_log_bsgs(m % pk, gk, pk);
logn = n_discrete_log_bsgs(n % pk, gk, pk);
x = (x + G->PHI[k] * logm * logn) % G->phi_q_odd;
}
return x;
}
static ulong
chi_even_exponent(const acb_dirichlet_group_t G, ulong m, ulong n)
{
ulong x;
ulong q_even = G->q_even;
if (q_even <= 2)
return 0;
x = 0;
if ((m % 4 == 3) && (n % 4 == 3))
x = q_even / 8;
if (q_even > 4)
{
ulong g2, logm, logn;
g2 = 5;
if (m % 4 == 3)
{
m = n_negmod(m, q_even);
}
if (n % 4 == 3)
{
n = n_negmod(n, q_even);
}
logm = n_discrete_log_bsgs(m % q_even, g2, q_even);
logn = n_discrete_log_bsgs(n % q_even, g2, q_even);
x += logm * logn;
}
return x % (q_even / 4);
}
void
acb_dirichlet_chi(acb_t res, const acb_dirichlet_group_t G, ulong m, ulong n, slong prec)
acb_dirichlet_chi(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n, slong prec)
{
fmpq_t t, u;
ulong odd_part, even_part;
ulong q_even = G->q_even;
ulong q_odd = G->q_odd;
if ((q_even > 1 && (n % 2 == 0)) || (q_odd > 1 && (n_gcd(q_odd, n) != 1)))
{
ulong expo;
expo = acb_dirichlet_ui_chi(G, chi, n);
if (expo == ACB_DIRICHLET_CHI_NULL)
acb_zero(res);
return;
}
if (q_even > 2)
{
if (q_even == 4)
{
if (m % 4 == 3 && n % 4 == 3)
even_part = q_even / 2; /* -1 */
else
even_part = 0; /* 1 */
}
else
{
even_part = 4 * chi_even_exponent(G, m % q_even, n % q_even);
}
}
else
{
even_part = 0;
fmpq_t t;
fmpq_init(t);
fmpq_set_si(t, 2 * expo , chi->order.n);
arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
fmpq_clear(t);
}
if (q_odd > 1)
odd_part = chi_odd_exponent(G, m % q_odd, n % q_odd);
else
odd_part = 0;
fmpq_init(t);
fmpq_init(u);
fmpq_set_si(t, 2 * even_part, q_even);
fmpq_set_si(u, 2 * odd_part, G->phi_q_odd);
fmpq_add(t, t, u);
arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
fmpq_clear(t);
fmpq_clear(u);
}

131
acb_dirichlet/chi_theta_arb.c Normal file
View file

@ -0,0 +1,131 @@
/*
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 "acb_dirichlet.h"
#include "acb_poly.h"
/* x(t) = exp(-Pi / q * t^2) */
static void
acb_dirichlet_arb_theta_xt(arb_t xt, ulong q, const arb_t t, slong prec)
{
arb_const_pi(xt, prec);
arb_div_ui(xt, xt, q, prec);
arb_mul(xt, xt, t, prec);
arb_mul(xt, xt, t, prec);
arb_neg(xt, xt);
arb_exp(xt, xt, prec);
}
void
acb_dirichlet_chi_theta_arb_smallorder(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t xt, slong len, slong prec)
{
ulong * a;
acb_dirichlet_powers_t z;
a = flint_malloc(len * sizeof(ulong));
acb_dirichlet_ui_chi_vec(a, G, chi, len);
_acb_dirichlet_powers_init(z, chi->order.n, 0, 0, prec);
acb_dirichlet_arb_theta_smallorder(res, xt, chi->parity, a, z, len, prec);
acb_dirichlet_powers_clear(z);
flint_free(a);
}
void
acb_dirichlet_chi_theta_arb_series(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t xt, slong len, slong prec)
{
acb_ptr a;
a = _acb_vec_init(len);
acb_dirichlet_chi_vec(a, G, chi, len, prec);
if (chi->parity)
{
slong k;
for (k = 2; k < len; k++)
acb_mul_si(a + k, a + k, k, prec);
}
acb_dirichlet_qseries_eval_arb(res, a, xt, len, prec);
_acb_vec_clear(a, len);
}
void
acb_dirichlet_vec_mellin_series(acb_ptr res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong len, const arb_t t, slong n, slong prec)
{
slong k;
arb_t tk, xt, stk, st;
acb_ptr a;
a = _acb_vec_init(len);
acb_dirichlet_chi_vec(a, G, chi, len, prec);
if (chi->parity)
{
for (k = 2; k < len; k++)
acb_mul_si(a + k, a + k, k, prec);
}
arb_init(tk);
arb_init(xt);
arb_init(st);
arb_init(stk);
arb_sqrt(st, t, prec);
arb_one(tk);
arb_one(stk);
for (k = 0; k < n; k++)
{
acb_dirichlet_arb_theta_xt(xt, G->q, tk, prec);
/* TODO: reduce len */
acb_dirichlet_qseries_eval_arb(res + k, a, xt, len, prec);
acb_mul_arb(res + k, res + k, stk, prec);
arb_mul(tk, tk, t, prec);
arb_mul(stk, stk, st, prec);
}
arb_clear(xt);
arb_clear(tk);
arb_clear(stk);
arb_clear(st);
_acb_vec_clear(a, len);
}
void
acb_dirichlet_chi_theta_arb_naive(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t xt, slong len, slong prec)
{
ulong * a;
acb_dirichlet_powers_t z;
a = flint_malloc(len * sizeof(ulong));
acb_dirichlet_ui_chi_vec(a, G, chi, len);
acb_dirichlet_powers_init(z, chi->order.n, len, prec);
acb_dirichlet_arb_theta_naive(res, xt, chi->parity, a, z, len, prec);
acb_dirichlet_powers_clear(z);
flint_free(a);
}
void
acb_dirichlet_chi_theta_arb(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t t, slong prec)
{
slong len;
arb_t xt;
len = acb_dirichlet_theta_length(G->q, t, prec);
arb_init(xt);
acb_dirichlet_arb_theta_xt(xt, G->q, t, prec);
/* TODO: tune this limit */
if (chi->order.n < 30)
acb_dirichlet_chi_theta_arb_smallorder(res, G, chi, xt, len, prec);
else
acb_dirichlet_chi_theta_arb_naive(res, G, chi, xt, len, prec);
arb_clear(xt);
}

37
acb_dirichlet/chi_vec.c Normal file
View file

@ -0,0 +1,37 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_chi_vec(acb_ptr v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv, slong prec)
{
slong k;
ulong * a;
acb_dirichlet_powers_t t;
a = flint_malloc(nv * sizeof(ulong));
acb_dirichlet_ui_chi_vec(a, G, chi, nv);
acb_dirichlet_powers_init(t, chi->order.n, nv, prec);
acb_zero(v + 0);
for (k = 0; k < nv; k++)
{
if (a[k] != ACB_DIRICHLET_CHI_NULL)
acb_dirichlet_power(v + k, t, a[k], prec);
else
*(v + k) = *(v + 0);
}
acb_dirichlet_powers_clear(t);
flint_free(a);
}

22
acb_dirichlet/conrey.c Normal file
View file

@ -0,0 +1,22 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_conrey_init(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G) {
x->log = flint_malloc(G->num * sizeof(ulong));
}
void
acb_dirichlet_conrey_clear(acb_dirichlet_conrey_t x) {
flint_free(x->log);
}

45
acb_dirichlet/conrey_conductor.c Normal file
View file

@ -0,0 +1,45 @@
/*
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 "acb_dirichlet.h"
ulong
acb_dirichlet_conrey_conductor(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x)
{
int k, f;
ulong cond = 1;
if (G->neven >= 1 && x->log[0] == 1)
cond = 4;
if (G->neven == 2 && x->log[1])
{
ulong l2 = x->log[1];
f = n_remove(&l2, 2);
cond = 1 << (G->P[1].e - f);
}
for (k = G->neven; k < G->num; k++)
{
if (x->log[k])
{
ulong p, lp;
p = G->P[k].p;
lp = x->log[k];
f = n_remove(&lp, p);
if (f)
cond *= n_pow(p, G->P[k].e - f);
else
cond *= G->P[k].pe.n;
}
}
return cond;
}

27
acb_dirichlet/conrey_eq.c Normal file
View file

@ -0,0 +1,27 @@
/*
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 "acb_dirichlet.h"
int
acb_dirichlet_conrey_eq(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x, const acb_dirichlet_conrey_t y)
{
slong k;
if (x->n != y->n)
return 0;
for (k = 0; k < G->num; k++)
if (x->log[k] != y->log[k])
return 0;
return 1;
}

22
acb_dirichlet/conrey_exp.c Normal file
View file

@ -0,0 +1,22 @@
/*
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 "acb_dirichlet.h"
ulong
acb_dirichlet_conrey_exp(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G)
{
ulong k, n = 1;
for (k = 0; k < G->num; k++)
n = nmod_mul(n, nmod_pow_ui(G->generators[k], x->log[k], G->mod), G->mod);
x->n = n;
return n;
}

34
acb_dirichlet/conrey_first_primitive.c Normal file
View file

@ -0,0 +1,34 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_conrey_first_primitive(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G)
{
ulong k;
if (G->q % 4 == 2)
{
flint_printf("Exception (acb_dirichlet_conrey_first_primitive). No primitive element mod %wu.\n",G->q);
abort();
}
x->n = 1;
for (k = 0; k < G->num ; k++)
{
if (k == 0 && G->neven == 2)
x->log[k] = 0;
else
{
x->n = nmod_mul(x->n, G->generators[k], G->mod);
x->log[k] = 1;
}
}
}

51
acb_dirichlet/conrey_log.c Normal file
View file

@ -0,0 +1,51 @@
/*
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 "acb_dirichlet.h"
/* TODO: use dlog module instead of n_discrete_log_bsgs */
/* assume m is invertible */
void
acb_dirichlet_conrey_log(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G, ulong m)
{
ulong k, pk, gk;
/* even part */
if (G->neven >= 1)
{
x->log[0] = (m % 4 == 3);
if (G->neven == 2)
{
ulong m2 = (x->log[0]) ? -m % G->q_even : m % G->q_even;
if (G->P[1].dlog == NULL)
x->log[1] = n_discrete_log_bsgs(m2, 5, G->q_even);
else
x->log[1] = dlog_precomp(G->P[1].dlog, m2);
}
}
/* odd part */
for (k = G->neven; k < G->num; k++)
{
if (G->P[k].dlog == NULL)
{
pk = G->P[k].pe.n;
gk = G->P[k].g;
x->log[k] = n_discrete_log_bsgs(m % pk, gk, pk);
}
else
{
x->log[k] = dlog_precomp(G->P[k].dlog, m % G->P[k].pe.n);
}
}
/* keep value m */
x->n = m;
}

21
acb_dirichlet/conrey_mul.c Normal file
View file

@ -0,0 +1,21 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_conrey_mul(acb_dirichlet_conrey_t c, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t a, const acb_dirichlet_conrey_t b)
{
ulong k;
for (k = 0; k < G->num ; k++)
c->log[k] = (a->log[k] + b->log[k]) % G->P[k].phi;
c->n = nmod_mul(a->n, b->n, G->mod);
}

29
acb_dirichlet/conrey_next.c Normal file
View file

@ -0,0 +1,29 @@
/*
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 "acb_dirichlet.h"
int
acb_dirichlet_conrey_next(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G)
{
/* update index */
int k;
for (k = G->num - 1; k >= 0; k--)
{
x->n = nmod_mul(x->n, G->generators[k], G->mod);
x->log[k] += 1;
if (x->log[k] < G->P[k].phi)
break;
x->log[k] = 0;
}
/* return last index modified */
return k;
}

60
acb_dirichlet/conrey_next_primitive.c Normal file
View file

@ -0,0 +1,60 @@
/*
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 "acb_dirichlet.h"
int
acb_dirichlet_conrey_next_primitive(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G)
{
/* update index avoiding multiples of p except for first component
if 8|q */
slong k;
/*
if (G->neven == 2)
{
x->n = nmod_mul(x->n, G->generators[0], G->mod);
x->log[0]++;
if (x->log[0] == 1)
return 0;
x->log[0] = 0;
k = 1;
}
*/
for (k = G->num - 1; k >= 0; k--)
{
#if 1
x->n = nmod_mul(x->n, G->generators[k], G->mod);
x->log[k]++;
if (x->log[k] % G->P[k].p == 0)
{
x->n = nmod_mul(x->n, G->generators[k], G->mod);
x->log[k]++;
}
if (x->log[k] < G->P[k].phi)
break;
if (x->log[k] == G->P[k].phi)
x->n = nmod_mul(x->n, G->generators[k], G->mod);
x->log[k] = 1;
#else
do {
x->n = nmod_mul(x->n, G->generators[k], G->mod);
x->log[k]++;
} while (x->log[k] % G->P[k].p == 0);
if (x->log[k] < G->P[k].phi)
break;
if (x->log[k] == G->P[k].phi)
x->n = nmod_mul(x->n, G->generators[k], G->mod);
x->log[k] = 1;
#endif
}
/* return last index modified */
return k;
}

21
acb_dirichlet/conrey_one.c Normal file
View file

@ -0,0 +1,21 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_conrey_one(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G)
{
ulong k;
for (k = 0; k < G->num ; k++)
x->log[k] = 0;
x->n = 1;
}

24
acb_dirichlet/conrey_order.c Normal file
View file

@ -0,0 +1,24 @@
/*
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 "acb_dirichlet.h"
ulong
acb_dirichlet_conrey_order(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x)
{
ulong k, g;
g = G->expo;
for (k = 0; k < G->num; k++)
g = n_gcd(g, G->PHI[k] * x->log[k]);
return G->expo / g;
}

26
acb_dirichlet/conrey_parity.c Normal file
View file

@ -0,0 +1,26 @@
/*
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 "acb_dirichlet.h"
int
acb_dirichlet_conrey_parity(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x)
{
int k, odd = 0;
for (k = 0; k < G->num; k++)
{
if (k == 1 && G->neven == 2)
continue;
if (x->log[k] % 2)
odd = 1 - odd;
}
return odd;
}

21
acb_dirichlet/conrey_pow.c Normal file
View file

@ -0,0 +1,21 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_conrey_pow(acb_dirichlet_conrey_t c, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t a, ulong n)
{
ulong k;
for (k = 0; k < G->num ; k++)
c->log[k] = (a->log[k] * n) % G->P[k].phi;
c->n = nmod_pow_ui(a->n, n, G->mod);
}

47
acb_dirichlet/conrey_primitive.c Normal file
View file

@ -0,0 +1,47 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_conrey_primitive(acb_dirichlet_conrey_t y, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x, ulong cond)
{
int k, l, f;
l = 0;
if (cond % 4 == 0)
{
y->log[l++] = x->log[0];
if (cond % 8 == 0)
{
ulong l2 = x->log[1];
f = 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];
f = n_remove(&lp, p);
y->log[l++] = lp;
}
}
}
}

25
acb_dirichlet/conrey_print.c Normal file
View file

@ -0,0 +1,25 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_conrey_print(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x)
{
slong k;
if (G->num)
flint_printf("[%wu", x->log[0]);
else
flint_printf("[");
for (k = 1; k < G->num; k++)
flint_printf(", %wu", x->log[k]);
flint_printf("]");
}

59
acb_dirichlet/dft.c Normal file
View file

@ -0,0 +1,59 @@
/*
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 "acb_dirichlet.h"
/* dft, lexicographic conrey indexing, array size G->phi_q */
void
acb_dirichlet_dft_conrey(acb_ptr w, acb_srcptr v, const acb_dirichlet_group_t G, slong prec)
{
slong k, l, * cyc;
cyc = flint_malloc(G->num * sizeof(slong));
for (k = 0, l = G->num - 1; l >= 0; k++, l--)
cyc[k] = G->P[k].phi;
acb_dirichlet_dft_prod(w, v, cyc, G->num, prec);
flint_free(cyc);
}
/* dft, number indexing, array size G->q */
void
acb_dirichlet_dft(acb_ptr w, acb_srcptr v, const acb_dirichlet_group_t G, slong prec)
{
ulong i, len;
acb_ptr t1, t2;
acb_dirichlet_conrey_t x;
len = G->phi_q;
t1 = flint_malloc(len * sizeof(acb_struct));
acb_dirichlet_conrey_init(x, G);
acb_dirichlet_conrey_one(x, G);
for (i = 0; i < len; i++)
{
t1[i] = v[x->n];
acb_dirichlet_conrey_next(x, G);
};
t2 = _acb_vec_init(len);
acb_dirichlet_dft_conrey(t2, t1, G, prec);
acb_dirichlet_conrey_one(x, G);
for (i = 0; i < len; i++)
{
acb_set(w + x->n, t2 + i);
acb_dirichlet_conrey_next(x, G);
};
_acb_vec_clear(t2, len);
acb_dirichlet_conrey_clear(x);
flint_free(t1);
}

277
acb_dirichlet/dft_fast.c Normal file
View file

@ -0,0 +1,277 @@
/*
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 "acb_dirichlet.h"
/* shallow extract */
static acb_ptr
vec_extract(acb_srcptr v, slong step, slong len)
{
slong k, l;
acb_ptr res;
res = flint_malloc(len * sizeof(acb_struct));
for (k = 0, l = 0; k < len; k++, l+=step)
res[k] = v[l];
return res;
}
void
_acb_dft_base(acb_ptr w, acb_srcptr v, slong dv, acb_srcptr z, slong dz, slong n, slong prec)
{
acb_ptr v1, z1;
v1 = vec_extract(v, dv, n);
z1 = vec_extract(z, dz, n);
_acb_dirichlet_dft_pol(w, v1, z1, n, prec);
flint_free(v1);
flint_free(z1);
}
typedef struct
{
slong m;
slong M;
slong dv;
acb_srcptr z;
slong dz;
}
dft_cyc_step;
typedef struct
{
slong n;
acb_ptr z;
slong num;
dft_cyc_step * cyc;
}
acb_dft_cyc_struct;
typedef acb_dft_cyc_struct acb_dft_cyc_t[1];
typedef struct
{
slong m;
slong M;
slong dv;
acb_dft_cyc_t pre;
}
dft_prod_step;
typedef struct
{
slong n;
slong num;
dft_prod_step * cyc;
}
acb_dft_prod_struct;
typedef acb_dft_prod_struct acb_dft_prod_t[1];
void
_acb_dirichlet_dft_cyc_init(acb_dft_cyc_t t, slong len, slong dv, slong prec)
{
slong i, j, num, dz;
n_factor_t fac;
acb_ptr z;
t->n = len;
n_factor_init(&fac);
n_factor(&fac, len, 0);
num = 0;
for (i = 0; i < fac.num; i++)
num += fac.exp[i];
t->num = num;
t->cyc = flint_malloc(num * sizeof(dft_cyc_step));
t->z = z = _acb_vec_init(t->n);
acb_dirichlet_vec_nth_roots(z, t->n, prec);
num = 0;
dz = 1;
for (i = 0; i < fac.num; i++)
{
for (j = 0; j < fac.exp[i]; j++)
{
slong m = fac.p[i];
t->cyc[num].m = m;
t->cyc[num].M = (len /= m);
t->cyc[num].dv = dv;
t->cyc[num].z = z;
t->cyc[num].dz = dz;
dv *= m;
dz *= m;
num++;
}
}
}
static void
acb_dirichlet_dft_cyc_init(acb_dft_cyc_t t, slong len, slong prec)
{
_acb_dirichlet_dft_cyc_init(t, len, 1, prec);
}
void
acb_dirichlet_dft_cyc_clear(acb_dft_cyc_t t)
{
_acb_vec_clear(t->z, t->n);
flint_free(t->cyc);
}
void
_acb_dft_cyc(acb_ptr w, acb_srcptr v, dft_cyc_step * cyc, slong num, slong prec)
{
dft_cyc_step c;
if (num == 0)
abort(); /* or just copy v to w */
c = cyc[0];
if (num == 1)
{
_acb_dft_base(w, v, c.dv, c.z, c.dz, c.m, prec);
}
else
{
slong i, j;
slong m = c.m, M = c.M, dv = c.dv, dz = c.dz;
acb_srcptr z = c.z, vi;
acb_ptr t, wi;
t = _acb_vec_init(m * M);
wi = w; vi = v;
for (i = 0; i < m; i++)
{
_acb_dft_cyc(wi, vi, cyc + 1, num - 1, prec);
if (i)
{
for (j = 1; j < M; j++)
acb_mul(wi + j, wi + j, z + dz * i * j, prec);
}
wi += M;
vi += dv;
}
/* after first pass */
for (j = 0; j < M; j++)
_acb_dft_base(t + m * j, w + j, M, z, dz * M, m, prec);
/* reorder */
for (i = 0; i < m; i++)
for (j = 0; j < M; j++)
w[j + M * i] = t[i + m * j];
_acb_vec_clear(t, m * M);
}
}
void
acb_dirichlet_dft_cyc_precomp(acb_ptr w, acb_srcptr v, acb_dft_cyc_t t, slong prec)
{
_acb_dft_cyc(w, v, t->cyc, t->num, prec);
}
void
acb_dirichlet_dft_fast(acb_ptr w, acb_srcptr v, slong len, slong prec)
{
acb_dft_cyc_t t;
acb_dirichlet_dft_cyc_init(t, len, prec);
acb_dirichlet_dft_cyc_precomp(w, v, t, prec);
acb_dirichlet_dft_cyc_clear(t);
}
void
_acb_dft_prod(acb_ptr w, acb_srcptr v, dft_prod_step * cyc, slong num, slong prec)
{
dft_prod_step c;
if (num == 0)
{
flint_printf("error: reached num = 0 in dft_prod\n");
abort(); /* or just copy v to w */
}
c = cyc[0];
if (num == 1)
{
acb_dirichlet_dft_cyc_precomp(w, v, c.pre, prec);
}
else
{
slong i, j;
slong m = c.m, M = c.M, dv = c.dv;
acb_srcptr vi;
acb_ptr t, wi;
t = _acb_vec_init(m * M);
wi = w; vi = v;
for (i = 0; i < m; i++)
{
_acb_dft_prod(wi, vi, cyc + 1, num - 1, prec);
wi += M;
vi += dv; /* here = M */
}
/* after first pass */
for (j = 0; j < M; j++)
acb_dirichlet_dft_cyc_precomp(t + m * j, w + j, c.pre, prec);
/* reorder */
for (i = 0; i < m; i++)
for (j = 0; j < M; j++)
w[j + M * i] = t[i + m * j];
_acb_vec_clear(t, m * M);
}
}
void
acb_dirichlet_dft_prod_precomp(acb_ptr w, acb_srcptr v, acb_dft_prod_t t, slong prec)
{
if (t->num == 0)
{
acb_set(w + 0, v + 0);
}
else
{
_acb_dft_prod(w, v, t->cyc, t->num, prec);
}
}
void
acb_dirichlet_dft_prod_init(acb_dft_prod_t t, slong * cyc, slong num, slong prec)
{
slong i, len, dv;
t->num = num;
t->cyc = flint_malloc(num * sizeof(dft_prod_step));
len = 1; dv = 1;
for (i = 0; i < num; i++)
len *= cyc[i];
for (i = 0; i < num; i++)
{
slong m = cyc[i];
len /= m;
t->cyc[i].m = m;
t->cyc[i].M = len;
t->cyc[i].dv = len;
_acb_dirichlet_dft_cyc_init(t->cyc[i].pre, m, len, prec);
dv *= m;
}
}
void
acb_dirichlet_dft_prod_clear(acb_dft_prod_t t)
{
slong i;
for (i = 0; i < t->num; i++)
acb_dirichlet_dft_cyc_clear(t->cyc[i].pre);
flint_free(t->cyc);
}
void
acb_dirichlet_dft_prod(acb_ptr w, acb_srcptr v, slong * cyc, slong num, slong prec)
{
acb_dft_prod_t t;
acb_dirichlet_dft_prod_init(t, cyc, num, prec);
acb_dirichlet_dft_prod_precomp(w, v, t, prec);
acb_dirichlet_dft_prod_clear(t);
}

44
acb_dirichlet/dft_pol.c Normal file
View file

@ -0,0 +1,44 @@
/*
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 "acb_poly.h"
#include "acb_dirichlet.h"
/* all roots are already computed */
void
_acb_dirichlet_dft_pol(acb_ptr w, acb_srcptr v, acb_srcptr z, slong len, slong prec)
{
/* FIXME: huge accuracy loss */
#if 0
_acb_poly_evaluate_vec_fast(w, v, len, z, len, prec);
#elif 0
_acb_poly_evaluate_vec_iter(w, v, len, z, len, prec);
#else
slong i, j;
for (i = 0; i < len; i++)
{
acb_zero(w + i);
for (j = 0; j < len; j++)
acb_addmul(w + i, v + j, z + (i * j % len), prec);
}
#endif
}
void
acb_dirichlet_dft_pol(acb_ptr w, acb_srcptr v, slong len, slong prec)
{
acb_ptr z;
z = _acb_vec_init(len);
acb_dirichlet_vec_nth_roots(z, len, prec);
_acb_dirichlet_dft_pol(w, v, z, len, prec);
_acb_vec_clear(z, len);
}

1
acb_dirichlet/eta.c
View file

@ -34,4 +34,3 @@ acb_dirichlet_eta(acb_t res, const acb_t s, slong prec)
acb_clear(t);
}
}

109
acb_dirichlet/gauss_sum.c Normal file
View file

@ -0,0 +1,109 @@
/*
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 "acb_dirichlet.h"
static void
gauss_sum_non_primitive(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
{
slong k, mu = 1;
ulong NN0 = G->q / chi->conductor;
/* G(chi) = mu(N/N0)chi0(N/N0)G(chi0) */
if (NN0 % 2 == 0)
{
if (G->q % 4)
mu = -1;
else
{
acb_zero(res);
return;
}
}
for (k = G->neven; k < G->num; k++)
{
ulong p = G->P[k].p;
if (G->P[k].e > 1 && NN0 % (p*p) == 0)
{
acb_zero(res);
return;
}
if (NN0 % p == 0)
mu *= -1;
}
if (chi->x->n == 1)
{
acb_set_si(res, mu);
}
else
{
acb_dirichlet_group_t G0;
acb_dirichlet_char_t chi0;
acb_t z;
acb_dirichlet_subgroup_init(G0, G, chi->conductor);
acb_dirichlet_char_init(chi0, G);
acb_dirichlet_char_primitive(chi0, G0, G, chi);
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);
acb_dirichlet_group_clear(G0);
acb_dirichlet_char_clear(chi0);
acb_clear(z);
}
}
void
acb_dirichlet_gauss_sum(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
{
/* TODO: no need, factor also does it... */
if (chi->conductor != G->q)
{
gauss_sum_non_primitive(res, G, chi, prec);
}
else if (chi->order.n <= 2)
{
if (chi->parity)
{
arb_sqrt_ui(acb_imagref(res), G->q, prec);
arb_zero(acb_realref(res));
}
else
{
arb_sqrt_ui(acb_realref(res), G->q, prec);
arb_zero(acb_imagref(res));
}
}
else if (G->num > 1 && G->num > G->neven)
{
acb_dirichlet_gauss_sum_factor(res, G, chi, prec);
}
else
{
if (acb_dirichlet_theta_length_d(G->q, 1, prec) > G->q)
acb_dirichlet_gauss_sum_naive(res, G, chi, prec);
else
acb_dirichlet_gauss_sum_theta(res, G, chi, prec);
}
}

71
acb_dirichlet/gauss_sum_factor.c Normal file
View file

@ -0,0 +1,71 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_gauss_sum_factor(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
{
slong k;
acb_t tmp;
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)
{
acb_zero(res);
return;
}
}
/* factor */
acb_one(res);
acb_init(tmp);
for (k = (G->neven == 2); k < G->num; k++)
{
ulong pe = G->P[k].pe.n;
acb_dirichlet_group_t Gp;
acb_dirichlet_char_t chip;
acb_dirichlet_subgroup_init(Gp, G, pe);
acb_dirichlet_char_init(chip, Gp);
chip->x->n = chi->x->n % pe;
if (k == 1 && G->neven == 2)
{
chip->x->log[0] = chi->x->log[0];
chip->x->log[1] = chi->x->log[1];
}
else
chip->x->log[0] = chi->x->log[k];
acb_dirichlet_char_conrey(chip, Gp, NULL);
/* chi_pe(a, q/pe) * G_pe(a) */
acb_dirichlet_gauss_sum(tmp, Gp, chip, prec);
acb_mul(res, res, tmp, prec);
acb_dirichlet_chi(tmp, Gp, chip, (G->q / pe) % pe, prec);
acb_mul(res, res, tmp, prec);
acb_dirichlet_char_clear(chip);
acb_dirichlet_group_clear(Gp);
}
if (G->q_even == 2)
acb_neg(res, res);
acb_clear(tmp);
}

31
acb_dirichlet/gauss_sum_naive.c Normal file
View file

@ -0,0 +1,31 @@
/*
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 "acb_dirichlet.h"
#include "acb_poly.h"
void
acb_dirichlet_gauss_sum_naive(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
{
acb_t z;
acb_ptr v;
v = _acb_vec_init(G->q);
acb_dirichlet_chi_vec(v, G, chi, G->q, prec);
acb_init(z);
acb_dirichlet_nth_root(z, G->q, prec);
_acb_poly_evaluate(res, v, G->q, z, prec);
acb_clear(z);
_acb_vec_clear(v, G->q);
}

51
acb_dirichlet/gauss_sum_theta.c Normal file
View file

@ -0,0 +1,51 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_gauss_sum_theta(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
{
arb_t x;
acb_t eps;
arb_init(x);
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)))
{
/* or could use l'Hopital rule */
acb_dirichlet_gauss_sum_naive(res, G, chi, prec);
return;
}
arb_one(x);
acb_dirichlet_chi_theta_arb(res, G, chi, x, prec);
acb_init(eps);
acb_conj(eps, res);
acb_div(eps, res, eps, prec);
if (chi->parity)
{
arb_zero(acb_realref(res));
arb_sqrt_ui(acb_imagref(res), G->q, prec);
}
else
{
arb_zero(acb_imagref(res));
arb_sqrt_ui(acb_realref(res), G->q, prec);
}
acb_mul(res, res, eps, prec);
arb_clear(x);
acb_clear(eps);
}

4
acb_dirichlet/group_clear.c
View file

@ -15,9 +15,7 @@
void
acb_dirichlet_group_clear(acb_dirichlet_group_t G)
{
flint_free(G->primes);
flint_free(G->exponents);
flint_free(G->P);
flint_free(G->generators);
flint_free(G->PHI);
}

28
acb_dirichlet/group_dlog_precompute.c Normal file
View file

@ -0,0 +1,28 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_prime_group_dlog_precompute(acb_dirichlet_prime_group_struct * P, ulong num)
{
P->dlog = flint_malloc(sizeof(dlog_precomp_t));
dlog_precomp_modpe_init(P->dlog, P->g, P->p, P->e, P->pe.n, num);
}
void
acb_dirichlet_group_dlog_precompute(acb_dirichlet_group_t G, ulong num)
{
slong k;
for (k = 0; k < G->num; k++)
acb_dirichlet_prime_group_dlog_precompute(&G->P[k], num);
}

214
acb_dirichlet/group_init.c
View file

@ -1,6 +1,7 @@
/*
Copyright (C) 2015 Jonathan Bober
Copyright (C) 2016 Fredrik Johansson
Copyright (C) 2016 Pascal Molin
This file is part of Arb.
@ -33,52 +34,189 @@ primitive_root_p_and_p2(ulong p)
}
void
acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q)
acb_dirichlet_prime_group_init(acb_dirichlet_prime_group_struct * P, ulong p, int e)
{
P->p = p;
P->e = e;
if (p == 2 || p == 4)
{
P->p = 2;
nmod_init(&P->pe, 1 << e);
if (p == 2)
{
P->e = 2;
P->phi = 2;
P->g = (1 << e) - 1;
}
else
{
P->phi = 1 << (e - 2);
P->g = 5;
}
}
else
{
ulong pe1;
pe1 = n_pow(p, e - 1);
P->phi = (p-1) * pe1;
nmod_init(&P->pe, p * pe1);
P->g = primitive_root_p_and_p2(p);
}
P->dlog = NULL;
}
static void
acb_dirichlet_group_lift_generators(acb_dirichlet_group_t G)
{
slong k;
n_factor_t fac;
acb_dirichlet_prime_group_struct * P = G->P;
G->q = q;
G->q_odd = q;
G->q_even = 1;
while (G->q_odd % 2 == 0)
G->expo = G->phi_q = 1;
if (G->neven)
{
G->q_odd /= 2;
G->q_even *= 2;
G->phi_q = G->q_even / 2;
G->expo = P[G->neven - 1].phi;
}
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);
}
n_factor_init(&fac);
n_factor(&fac, G->q_odd, 1);
G->num = fac.num;
G->primes = flint_malloc(G->num * sizeof(ulong));
G->exponents = flint_malloc(G->num * sizeof(ulong));
G->generators = flint_malloc(G->num * sizeof(ulong));
G->PHI = flint_malloc(G->num * sizeof(ulong));
for (k = 0; k < G->num; k++)
{
G->primes[k] = fac.p[k];
G->exponents[k] = fac.exp[k];
}
G->phi_q_odd = 1;
for (k = 0; k < G->num; k++)
G->phi_q_odd *= (G->primes[k] - 1) * n_pow(G->primes[k], G->exponents[k]-1);
if (G->q_even == 1)
G->phi_q = G->phi_q_odd;
else
G->phi_q = G->phi_q_odd * (G->q_even / 2);
for (k = 0; k < G->num; k++)
{
ulong phi;
G->generators[k] = primitive_root_p_and_p2(G->primes[k]);
phi = n_pow(G->primes[k], G->exponents[k] - 1) * (G->primes[k] - 1);
G->PHI[k] = G->phi_q_odd / phi;
nmod_t pe;
ulong qpe, v;
G->PHI[k] = G->expo / G->P[k].phi;
/* 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;
qpe = G->q / pe.n;
if (G->q < G->P[k].pe.n)
{
flint_printf("lift generator %wu from %wu to %wu e=%wu\n",
G->P[k].g, G->P[k].pe.n, G->q, G->P[k].e);
}
v = nmod_inv(qpe % pe.n, pe);
/* no overflow since v * qpe < q */
G->generators[k] = (1 + (G->P[k].g-1) * v * qpe) % G->q;
}
}
void
acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q)
{
slong k;
ulong e2;
n_factor_t fac;
G->q = q;
nmod_init(&G->mod, q);
e2 = n_remove(&q, 2);
G->q_even = 1 << e2;
/* number of components at p=2 */
G->neven = (e2 >= 3) ? 2 : (e2 == 2) ? 1 : 0;
/* warning: only factor odd part */
n_factor_init(&fac);
n_factor(&fac, q, 1);
G->num = fac.num + G->neven;
G->P = flint_malloc(G->num * sizeof(acb_dirichlet_prime_group_struct));
G->generators = flint_malloc(G->num * sizeof(ulong));
G->PHI = flint_malloc(G->num * sizeof(ulong));
/* even part */
if (G->neven >= 1)
acb_dirichlet_prime_group_init(&G->P[0], 2, e2);
if (G->neven == 2)
acb_dirichlet_prime_group_init(&G->P[1], 4, e2);
/* odd part */
G->rad_q = 1;
for (k = G->neven; k < G->num; k++)
{
ulong p, e;
p = fac.p[k - G->neven];
e = fac.exp[k - G->neven];
G->rad_q *= p;
acb_dirichlet_prime_group_init(&G->P[k], p, e);
}
acb_dirichlet_group_lift_generators(G);
}
void
acb_dirichlet_subgroup_init(acb_dirichlet_group_t H, const acb_dirichlet_group_t G, ulong h)
{
int s[15]; /* selected components */
slong k, j, e2;
H->q = h;
nmod_init(&H->mod, h);
/* even components */
e2 = n_remove(&h, 2);
H->q_even = 1 << e2;
s[0] = s[1] = 0;
j = 0;
if (e2 >= 2)
s[j++] = 2;
if (e2 >= 3)
s[j++] = e2;
H->neven = j;
/* odd components */
for (k = G->neven; k < G->num; k++)
{
ulong p = G->P[k].p;
s[k] = n_remove(&h, p);
if (s[k] > 0)
{
j++;
H->rad_q *= p;
}
}
H->num = j;
H->P = flint_malloc(j * sizeof(acb_dirichlet_prime_group_struct));
H->generators = flint_malloc(j * sizeof(ulong));
H->PHI = flint_malloc(j * sizeof(ulong));
j = 0;
for (k = 0; k < H->neven; k++)
{
H->P[j] = G->P[k];
if (H->q_even < G->q_even)
{
nmod_init(&H->P[j].pe, H->q_even);
H->P[j].e = s[k];
if (k == 0)
H->P[j].g = H->q_even - 1;
}
j++;
}
for (k = G->neven; k < G->num; k++)
{
if (s[k])
{
H->P[j] = G->P[k];
if (s[k] < G->P[k].e)
{
ulong p = H->P[j].p;
H->P[j].e = s[k];
nmod_init(&H->P[j].pe, n_pow(p, s[k]));
}
j++;
}
}
acb_dirichlet_group_lift_generators(H);
}

169
acb_dirichlet/jacobi_sum.c Normal file
View file

@ -0,0 +1,169 @@
/*
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 "acb_dirichlet.h"
/* J_N(1,a) = sum on x = 1 mod some p | q */
static ulong
jacobi_one_prime(ulong p, ulong e, ulong pe, ulong cond)
{
if (e > 1 && cond % (p*p) == 0)
{
return 0;
}
else
{
slong r = (e > 1) ? pe / p : 1;
if (cond % p)
return r * (p - 2);
else
return -r;
}
}
static ulong
jacobi_one(const acb_dirichlet_group_t G, ulong cond)
{
slong k, r = 1;
for (k = 0; k < G->num; k++)
r *= jacobi_one_prime(G->P[k].p, G->P[k].e,
G->P[k].pe.n, cond);
return r;
}
/* should use only for prime power modulus */
static void
acb_dirichlet_jacobi_sum_gauss(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2, slong prec)
{
/* J_q(a,b)G_q(ab) = G_q(a)G_q(b) */
acb_t tmp;
acb_dirichlet_char_t chi12;
acb_dirichlet_char_init(chi12, G);
acb_dirichlet_char_mul(chi12, G, chi1, chi2);
acb_init(tmp);
acb_dirichlet_gauss_sum(res, G, chi1, prec);
if (chi2->x->n == chi1->x->n)
acb_set(tmp, res);
else
acb_dirichlet_gauss_sum(tmp, G, chi2, prec);
acb_mul(res, res, tmp, prec);
acb_dirichlet_gauss_sum(tmp, G, chi12, prec);
acb_div(res, res, tmp, prec);
acb_dirichlet_char_clear(chi12);
acb_clear(tmp);
}
static void
acb_dirichlet_jacobi_sum_primes(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2, slong prec)
{
slong k;
acb_t tmp;
acb_init(tmp);
acb_one(res);
/* TODO: efficient subgroup */
for (k = 0; k < G->num; k++)
{
nmod_t pe;
ulong p, e, ap, bp;
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;
if (ap == 1 || bp == 1 || nmod_mul(ap, bp, pe) == 1)
{
slong r;
ulong cond;
cond = (ap == 1) ? chi2->conductor : chi1->conductor;
r = jacobi_one_prime(p, e, pe.n, cond);
/* chi(a,-1) if ap * bp = 1 */
if (ap != 1 && bp != 1)
r *= n_jacobi_unsigned(ap, p);
acb_mul_si(res, res, r, prec);
}
else
{
acb_dirichlet_group_t Gp;
acb_dirichlet_char_t chi1p, chi2p;
acb_dirichlet_group_init(Gp, pe.n);
acb_dirichlet_char_init(chi1p, Gp);
acb_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];
acb_dirichlet_char_conrey(chi1p, Gp, NULL);
acb_dirichlet_char_conrey(chi2p, Gp, NULL);
/* TODO: work out gauss relations for e > 1 */
if (p <= 100 || e > 1)
acb_dirichlet_jacobi_sum_naive(tmp, Gp, chi1p, chi2p, prec);
else
acb_dirichlet_jacobi_sum_gauss(tmp, Gp, chi1p, chi2p, prec);
acb_mul(res, res, tmp, prec);
acb_dirichlet_char_clear(chi1p);
acb_dirichlet_char_clear(chi2p);
acb_dirichlet_group_clear(Gp);
}
}
acb_clear(tmp);
}
void
acb_dirichlet_jacobi_sum(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2, slong prec)
{
if (G->q_even > 1)
{
acb_zero(res);
}
else if (chi1->x->n == 1 || chi2->x->n == 1)
{
ulong cond = (chi1->x->n == 1) ? chi2->conductor : chi1->conductor;
acb_set_si(res, jacobi_one(G, cond));
}
else if (nmod_mul(chi1->x->n, chi2->x->n, G->mod) == 1)
{
ulong n;
n = jacobi_one(G, chi1->conductor);
if (chi1->parity)
acb_set_si(res, -n);
else
acb_set_si(res, n);
}
else
{
if (G->q <= 150)
acb_dirichlet_jacobi_sum_naive(res, G, chi1, chi2, prec);
else if (G->num > 1)
acb_dirichlet_jacobi_sum_primes(res, G, chi1, chi2, prec);
else if (G->P[0].e > 1)
acb_dirichlet_jacobi_sum_naive(res, G, chi1, chi2, prec);
else
acb_dirichlet_jacobi_sum_gauss(res, G, chi1, chi2, prec);
}
}

62
acb_dirichlet/jacobi_sum_naive.c Normal file
View file

@ -0,0 +1,62 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_jacobi_sum_naive(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2, slong prec)
{
ulong k1, k2, m1, m2, g, e;
ulong * v1, * v2;
slong *v;
nmod_t order;
acb_t z;
v1 = flint_malloc(G->q * sizeof(ulong));
v2 = flint_malloc(G->q * sizeof(ulong));
acb_dirichlet_ui_vec_set_null(v1, G, G->q);
acb_dirichlet_ui_chi_vec_loop(v1, G, chi1, G->q);
acb_dirichlet_ui_vec_set_null(v2, G, G->q);
acb_dirichlet_ui_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;
v = flint_malloc(order.n * sizeof(slong));
for (e = 0; e < order.n; e++)
v[e] = 0;
for (k1 = 2, k2 = G->q - 1; k2 > 1; k1++, k2--)
{
if (v1[k1] == ACB_DIRICHLET_CHI_NULL ||
v2[k2] == ACB_DIRICHLET_CHI_NULL)
continue;
e = nmod_add(v1[k1] * m1, v2[k2] * m2, order);
v[e]++;
}
acb_init(z);
acb_dirichlet_nth_root(z, order.n, prec);
acb_dirichlet_si_poly_evaluate(res, v, order.n, z, prec);
acb_clear(z);
flint_free(v);
flint_free(v2);
flint_free(v1);
}

51
acb_dirichlet/l.c
View file

@ -1,51 +0,0 @@
/*
Copyright (C) 2016 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_dirichlet.h"
void
acb_dirichlet_l(acb_t res, const acb_t s,
const acb_dirichlet_group_t G, ulong m, slong prec)
{
acb_t chi, t, u, a;
ulong k;
acb_init(chi);
acb_init(t);
acb_init(u);
acb_init(a);
acb_zero(t);
for (k = 1; k <= G->q; k++)
{
acb_dirichlet_chi(chi, G, m, k, prec);
if (!acb_is_zero(chi))
{
acb_set_ui(a, k);
acb_div_ui(a, a, G->q, prec);
acb_hurwitz_zeta(u, s, a, prec);
acb_addmul(t, chi, u, prec);
}
}
acb_set_ui(u, G->q);
acb_neg(a, s);
acb_pow(u, u, a, prec);
acb_mul(res, t, u, prec);
acb_clear(chi);
acb_clear(t);
acb_clear(u);
acb_clear(a);
}

59
acb_dirichlet/l_hurwitz.c Normal file
View file

@ -0,0 +1,59 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_l_hurwitz(acb_t res, const acb_t s,
const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
{
ulong chin;
acb_t t, u, a;
acb_ptr z;
acb_dirichlet_conrey_t cn;
acb_dirichlet_conrey_init(cn, G);
acb_init(t);
acb_init(u);
acb_init(a);
acb_dirichlet_conrey_one(cn, G);
acb_zero(t);
prec += n_clog(G->phi_q, 2);
z = _acb_vec_init(chi->order.n);
acb_dirichlet_vec_nth_roots(z, chi->order.n, prec);
do {
chin = acb_dirichlet_ui_chi_conrey(G, chi, cn);
acb_set_ui(a, cn->n);
acb_div_ui(a, a, G->q, prec);
acb_hurwitz_zeta(u, s, a, prec);
acb_addmul(t, z + chin, u, prec);
} while (acb_dirichlet_conrey_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);
acb_dirichlet_conrey_clear(cn);
_acb_vec_clear(z, chi->order.n);
acb_clear(t);
acb_clear(u);
acb_clear(a);
}

50
acb_dirichlet/l_vec_hurwitz.c Normal file
View file

@ -0,0 +1,50 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_l_vec_hurwitz(acb_ptr res, const acb_t s,
const acb_dirichlet_group_t G, slong prec)
{
acb_t a, qs;
acb_ptr zeta, z;
acb_dirichlet_conrey_t cn;
acb_dirichlet_conrey_init(cn, G);
acb_init(qs);
acb_init(a);
prec += n_clog(G->phi_q, 2);
acb_set_ui(qs, G->q);
acb_neg(a, s);
acb_pow(qs, qs, a, prec);
zeta = z = _acb_vec_init(G->phi_q);
acb_dirichlet_conrey_one(cn, G);
do {
acb_set_ui(a, cn->n);
acb_div_ui(a, a, G->q, prec);
acb_hurwitz_zeta(z, s, a, prec);
acb_mul(z, z, qs, prec);
z++;
} while (acb_dirichlet_conrey_next(cn, G) >= 0);
acb_dirichlet_dft_conrey(res, zeta, G, prec);
acb_dirichlet_conrey_clear(cn);
_acb_vec_clear(zeta, G->phi_q);
acb_clear(qs);
acb_clear(a);
}

42
acb_dirichlet/nth_root.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 "acb_dirichlet.h"
static void
_acb_dirichlet_nth_root(acb_t res, ulong order, slong prec)
{
fmpq_t t;
fmpq_init(t);
fmpq_set_si(t, 2, order);
arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
fmpq_clear(t);
}
void
acb_dirichlet_nth_root(acb_t res, ulong order, slong prec)
{
switch (order)
{
case 1:
acb_one(res);
break;
case 2:
acb_set_si(res, -1);
break;
case 4:
acb_onei(res);
break;
default:
_acb_dirichlet_nth_root(res, order, prec);
break;
}
}

37
acb_dirichlet/number_primitive.c Normal file
View file

@ -0,0 +1,37 @@
/*
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 "acb_dirichlet.h"
ulong
acb_dirichlet_number_primitive(const acb_dirichlet_group_t G)
{
if (G->q % 4 == 2)
return 0;
else
{
slong k;
ulong n = 1;
/* no overflow since result < G->q */
for (k = (G->neven == 2); k < G->num; k++)
{
ulong p = G->P[k].p, e = G->P[k].e;
if (e == 1)
n *= p - 2;
else
n *= (p * (p - 2) + 1) * n_pow(p, e - 2);
}
return n;
}
}

29
acb_dirichlet/pairing.c Normal file
View file

@ -0,0 +1,29 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_pairing(acb_t res, const acb_dirichlet_group_t G, ulong m, ulong n, slong prec)
{
ulong expo;
expo = acb_dirichlet_ui_pairing(G, m, n);
if (expo == ACB_DIRICHLET_CHI_NULL)
acb_zero(res);
else
{
fmpq_t t;
fmpq_init(t);
fmpq_set_si(t, 2 * expo, G->expo);
arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
fmpq_clear(t);
}
}

29
acb_dirichlet/pairing_conrey.c Normal file
View file

@ -0,0 +1,29 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_pairing_conrey(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t a, const acb_dirichlet_conrey_t b, slong prec)
{
ulong expo;
expo = acb_dirichlet_ui_pairing_conrey(G, a, b);
if (expo == ACB_DIRICHLET_CHI_NULL)
acb_zero(res);
else
{
fmpq_t t;
fmpq_init(t);
fmpq_set_si(t, 2 * expo, G->expo);
arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
fmpq_clear(t);
}
}

36
acb_dirichlet/power.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 "acb_dirichlet.h"
#include "acb_poly.h"
void
acb_dirichlet_power(acb_t z, const acb_dirichlet_powers_t t, ulong n, slong prec)
{
if (!t->depth)
{
acb_pow_ui(z, t->z, n, prec);
}
else
{
slong k;
ulong r;
r = n % t->size;
n = n / t->size;
acb_set(z, t->Z[0] + r);
for (k = 1; k < t->depth && n; k++)
{
r = n % t->size;
n = n / t->size;
acb_mul(z, z, t->Z[k] + r, prec);
}
}
}

24
acb_dirichlet/powers_clear.c Normal file
View file

@ -0,0 +1,24 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_powers_clear(acb_dirichlet_powers_t t)
{
slong k;
for (k = 0; k < t->depth; k++)
_acb_vec_clear(t->Z[k], t->size);
flint_free(t->Z);
acb_clear(t->z);
}

51
acb_dirichlet/powers_init.c Normal file
View file

@ -0,0 +1,51 @@
/*
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 "acb_dirichlet.h"
#include "acb_poly.h"
void
_acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong size, slong depth, slong prec)
{
slong k;
t->order = order;
acb_init(t->z);
acb_dirichlet_nth_root(t->z, order, prec);
t->size = size;
t->depth = depth;
if (depth)
{
acb_struct * z;
z = t->z + 0;
t->Z = flint_malloc(depth * sizeof(acb_ptr));
for (k = 0; k < depth; k++)
{
t->Z[k] = _acb_vec_init(size);
_acb_vec_set_powers(t->Z[k], z, size, prec);
z = t->Z[k] + 1;
}
}
else
{
t->Z = NULL;
}
}
void
acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong num, slong prec)
{
slong size, depth;
depth = (num > order) ? 1 : n_flog(order, num);
size = n_root(order, depth) + 1;
_acb_dirichlet_powers_init(t, order, size, depth, prec);
}

185
acb_dirichlet/profile/p-conrey.c Normal file
View file

@ -0,0 +1,185 @@
/*
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 <string.h>
#include "acb_dirichlet.h"
#include "profiler.h"
#define LOG 0
#define CSV 1
#define JSON 2
typedef ulong (*do_f) (ulong q1, ulong q2);
static ulong
do_gcd(ulong q1, ulong q2)
{
ulong n, q, k;
for (n = 0, q = q1; q <= q2; q++)
for (k = 1; k < q; k++)
n += (n_gcd(k, q) == 1);
return n;
}
static ulong
do_conrey(ulong q1, ulong q2)
{
ulong n, q;
for (n = 0, q = q1; q <= q2; q++)
{
acb_dirichlet_group_t G;
acb_dirichlet_conrey_t x;
acb_dirichlet_group_init(G, q);
acb_dirichlet_conrey_init(x, G);
acb_dirichlet_conrey_one(x, G);
n++;
for (; acb_dirichlet_conrey_next(x, G) >= 0; n++);
acb_dirichlet_conrey_clear(x);
acb_dirichlet_group_clear(G);
}
return n;
}
static ulong
do_chars(ulong q1, ulong q2)
{
ulong n, q;
for (n = 0, q = q1; q <= q2; q++)
{
acb_dirichlet_group_t G;
acb_dirichlet_char_t chi;
acb_dirichlet_group_init(G, q);
acb_dirichlet_char_init(chi, G);
acb_dirichlet_char_one(chi, G);
n++;
for (; acb_dirichlet_char_next(chi, G) >= 0; n++);
acb_dirichlet_char_clear(chi);
acb_dirichlet_group_clear(G);
}
return n;
}
static ulong
do_gcdpluscond(ulong q1, ulong q2)
{
ulong n, q, k;
for (n = 0, q = q1; q <= q2; q++)
{
acb_dirichlet_group_t G;
acb_dirichlet_group_init(G, q);
for (k = 1; k < q; k++)
{
/* known factors -> faster gcd */
slong i;
if (G->q_even > 1 && k % 2 == 0)
continue;
for (i = G->neven; i < G->num; i++)
if (k % G->P[i].p == 0)
break;
if (i == G->num)
acb_dirichlet_ui_conductor(G, k);
}
n += G->phi_q;
acb_dirichlet_group_clear(G);
}
return n;
}
int main(int argc, char *argv[])
{
int out;
ulong n, nref, maxq = 5000;
int l, nf = 4;
do_f func[4] = { do_gcd, do_conrey, do_chars, do_gcdpluscond };
char * name[4] = { "gcd", "conrey", "chars", "gcd+cond" };
int i, ni = 5;
ulong qmin[5] = { 2, 1000, 10000, 100000, 1000000 };
ulong qmax[5] = { 500, 3000, 11000, 100100, 1000010 };
if (argc < 2)
out = LOG;
else if (!strcmp(argv[1], "json"))
out = JSON;
else if (!strcmp(argv[1], "csv"))
out = CSV;
else if (!strcmp(argv[1], "log"))
out = LOG;
else
{
printf("usage: %s [log|csv|json]\n", argv[0]);
abort();
}
if (out == CSV)
flint_printf("# %-6s, %7s, %7s, %7s\n","name", "qmin", "qmax", "time");
for (i = 0; i < ni; i++)
{
if (out == LOG)
{
flint_printf("loop over all (Z/q)* for %wu<=q<=%wu\n", qmin[i], qmax[i]);
}
for (l = 0; l < nf; l++)
{
if (out == LOG)
flint_printf("%-8s ... ",name[l]);
else if (out == CSV)
flint_printf("%-8s, %7d, %7d, ",name[l],qmin[i],qmax[i]);
else if (out == JSON)
flint_printf("{ \"name\": \"%s\", \"qmin\": %d, \"qmax\": %d, \"time\": ",
name[l],qmin[i],qmax[i]);
TIMEIT_ONCE_START
(func[l])(qmin[i], qmax[i]);
TIMEIT_ONCE_STOP
if (l == 0)
nref = n;
else if (n != nref)
{
flint_printf("FAIL: wrong number of elements %wu != %wu\n\n",n, nref);
abort();
}
if (out == JSON)
flint_printf("}\n");
else
flint_printf("\n");
}
}
flint_cleanup();
return EXIT_SUCCESS;
}

98
acb_dirichlet/profile/p-thetanull.c Normal file
View file

@ -0,0 +1,98 @@
/*
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 "acb_dirichlet.h"
int main()
{
slong prec = 64;
n_primes_t iter;
ulong q;
flint_printf("thetanull....");
fflush(stdout);
/* look for vanishing theta values for prime power moduli */
n_primes_init(iter);
#if 0 /* just to check q=300 & q=600 values */
q = 2;
while (q++ < 1000)
#else
while ((q = n_primes_next(iter)) < 10000)
#endif
{
slong s;
ulong k;
acb_dirichlet_group_t G;
acb_dirichlet_conrey_t x;
acb_dirichlet_group_init(G, q);
acb_dirichlet_conrey_init(x, G);
acb_ptr theta, z;
arb_t z1;
arb_ptr t;
theta = _acb_vec_init(G->phi_q);
z = _acb_vec_init(G->phi_q);
arb_init(z1);
arb_const_pi(z1, prec);
arb_div_ui(z1, z1, q, prec);
arb_neg(z1, z1);
arb_exp(z1, z1, prec);
t = _arb_vec_init(q);
acb_dirichlet_arb_quadratic_powers(t, q, z1, prec);
for (s = 0; s <= 1; s++)
{
k = 0;
acb_dirichlet_conrey_one(x, G);
do {
acb_set_arb(z + k, t + x->n);
k++;
} while (acb_dirichlet_conrey_next(x, G) >= 0);
acb_dirichlet_dft_conrey(theta, z, G, prec);
for (k = 0; k < G->phi_q; k++)
{
if (acb_contains_zero(theta + k))
{
slong i;
acb_dirichlet_conrey_one(x, G);
for (i = 0; i < k; i++)
acb_dirichlet_conrey_next(x, G);
if (acb_dirichlet_conrey_conductor(G,x) < q || acb_dirichlet_conrey_parity(G, x) != s)
continue;
flint_printf("\ntheta null q = %wu, n = %wu\n",q, x->n);
acb_printd(theta + k, 10);
}
}
/* change t for odd characters */
for (k = 0; k < q; k++)
arb_mul_ui(t + k, t + k, k, prec);
}
_arb_vec_clear(t, q);
_acb_vec_clear(theta, G->phi_q);
_acb_vec_clear(z, G->phi_q);
arb_clear(z1);
acb_dirichlet_conrey_clear(x);
acb_dirichlet_group_clear(G);
}
flint_cleanup();
return EXIT_SUCCESS;
}

139
acb_dirichlet/profile/p-vec.c Normal file
View file

@ -0,0 +1,139 @@
/*
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 "acb_dirichlet.h"
#include "profiler.h"
typedef void (*dir_f) (ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv);
static void
dir_empty(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv)
{
return;
}
static void
vecloop(dir_f dir, ulong minq, ulong maxq, ulong * rand, ulong nr, ulong * v, ulong nv)
{
ulong q;
TIMEIT_ONCE_START
for (q = minq; q <= maxq; q++)
{
ulong r;
acb_dirichlet_group_t G;
acb_dirichlet_char_t chi;
acb_dirichlet_group_init(G, q);
acb_dirichlet_char_init(chi, G);
for (r = 0; r < nr; r++)
{
acb_dirichlet_char(chi, G, rand[r] % q);
dir(v, G, chi, nv);
}
acb_dirichlet_char_clear(chi);
acb_dirichlet_group_clear(G);
}
TIMEIT_ONCE_STOP
flint_printf("\n");
}
static void
acb_dirichlet_ui_chi_vec_sieve(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv)
{
slong k, p, pmax;
n_primes_t iter;
n_primes_init(iter);
pmax = (nv < G->q) ? nv : G->q;
v[1] = 0;
while ((p = n_primes_next(iter)) < pmax)
{
if (G->q % p == 0)
{
for (k = p; k < nv; k += p)
v[k] = ACB_DIRICHLET_CHI_NULL;
}
else
{
long chip;
chip = acb_dirichlet_ui_chi(G, chi, p);
for (k = p; k < nv; k += p)
if (v[k] != -1)
v[k] = nmod_add(v[k], chip, chi->order);
}
}
n_primes_clear(iter);
}
int main()
{
slong iter, k, nv, nref, r, nr;
ulong minq, maxq;
ulong * rand;
int i, ni = 5;
ulong q[5] = { 2, 1000, 3000, 10000, 100000 };
ulong qq[5] = { 500, 2000, 5000, 12000, 100500 };
ulong * v;
flint_rand_t state;
nr = 20;
flint_randinit(state);
rand = flint_malloc(nr * sizeof(ulong));
v = flint_malloc(nv * sizeof(ulong));
for (r = 0; r < nr; r++)
rand[r] = n_randprime(state, 42, 0);
for (i = 0; i < ni; i++)
{
ulong minq = q[i], maxq = qq[i];
nv = 2000;
flint_printf("%wu * chi(rand, 1..%wu) for all %wu <= q <= %wu....\n", nr, nv, minq, maxq);
fflush(stdout);
flint_printf("character only.......... ");
fflush(stdout);
vecloop(dir_empty, minq, maxq, rand, nr, v, nv);
flint_printf("big loop................ ");
fflush(stdout);
vecloop(acb_dirichlet_ui_chi_vec_loop, minq, maxq, rand, nr, v, nv);
flint_printf("med loop................ ");
fflush(stdout);
vecloop(acb_dirichlet_ui_chi_vec_primeloop, minq, maxq, rand, nr, v, nv);
flint_printf("sieve................... ");
fflush(stdout);
vecloop(acb_dirichlet_ui_chi_vec_sieve, minq, maxq, rand, nr, v, nv);
flint_printf("generic................. ");
fflush(stdout);
vecloop(acb_dirichlet_ui_chi_vec, minq, maxq, rand, nr, v, nv);
}
flint_free(v);
flint_free(rand);
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

42
acb_dirichlet/qseries_eval_arb.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 "acb_dirichlet.h"
/* assume a[0] = 0 */
void
acb_dirichlet_qseries_eval_arb(acb_t res, acb_srcptr a, const arb_t x, slong len, slong prec)
{
slong k;
arb_t xk2, dx, x2;
arb_init(xk2);
arb_init(dx);
arb_init(x2);
arb_set(dx, x);
arb_set(xk2, dx);
arb_mul(x2, dx, dx, prec);
acb_mul_arb(res, a + 1, xk2, prec);
/* TODO: reduce prec */
for (k = 2; k < len; k++)
{
arb_mul(dx, dx, x2, prec);
arb_mul(xk2, xk2, dx, prec);
acb_addmul_arb(res, a + k, xk2, prec);
}
arb_clear(xk2);
arb_clear(x2);
arb_clear(dx);
}

62
acb_dirichlet/si_poly_evaluate.c Normal file
View file

@ -0,0 +1,62 @@
/*
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 "acb_dirichlet.h"
/* bsgs evaluation */
void
acb_dirichlet_si_poly_evaluate(acb_t res, slong * v, slong len, const acb_t z, slong prec)
{
slong k, r, m;
acb_t sq;
acb_ptr zk;
if (len < 3)
{
if (len == 0)
{
acb_zero(res);
}
else if (len == 1)
{
acb_set_si(res, v[0]);
}
else if (len == 2)
{
acb_mul_si(res, z, v[1], prec);
acb_add_si(res, res, v[0], prec);
}
return;
}
m = n_sqrt(len) + 1;
zk = _acb_vec_init(m + 1);
_acb_vec_set_powers(zk, z, m + 1, prec);
acb_init(sq);
acb_zero(res);
k = len - 1;
r = k % m;
for (; k >= 0; r = m - 1)
{
acb_zero(sq);
for (; r >= 0; r--, k--)
acb_addmul_si(sq, zk + r, v[k], prec);
acb_mul(res, res, zk + m, prec);
acb_add(res, res, sq, prec);
}
_acb_vec_clear(zk, m + 1);
acb_clear(sq);
}

151
acb_dirichlet/test/t-chars.c Normal file
View file

@ -0,0 +1,151 @@
/*
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 "acb_dirichlet.h"
int main()
{
slong iter, bits;
flint_rand_t state;
flint_printf("chars....");
fflush(stdout);
flint_randinit(state);
for (bits = 5; bits <= 30; bits += 5)
{
for (iter = 0; iter < 50; iter++)
{
acb_dirichlet_group_t G;
acb_dirichlet_conrey_t x;
acb_dirichlet_char_t chi, chi2;
ulong q, iter2;
q = 2 + n_randint(state, 1 << bits);
acb_dirichlet_group_init(G, q);
acb_dirichlet_conrey_init(x, G);
acb_dirichlet_char_init(chi, G);
acb_dirichlet_char_init(chi2, G);
acb_dirichlet_group_dlog_precompute(G, 50);
/* check number char properties */
for (iter2 = 0; iter2 < 100; iter2++)
{
int par;
ulong m, n;
ulong order, chim1, pairing, cond;
do
m = n_randint(state, q);
while (n_gcd(q, m) > 1);
acb_dirichlet_char(chi, G, m);
acb_dirichlet_conrey_log(x, G, m);
acb_dirichlet_char_conrey(chi2, G, x);
if (!acb_dirichlet_char_eq(G, chi, chi2))
{
flint_printf("FAIL: init char\n\n");
flint_printf("q = %wu\n\n", q);
flint_printf("m = %wu\n\n", m);
acb_dirichlet_char_print(G, chi);
flint_printf("\n");
acb_dirichlet_char_print(G, chi2);
flint_printf("\n");
abort();
}
order = acb_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 = acb_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 = acb_dirichlet_ui_parity(G, m);
chim1 = acb_dirichlet_ui_chi(G, chi, q - 1);
if (acb_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", acb_dirichlet_char_parity(chi));
flint_printf("parity_ui = %d", par);
acb_dirichlet_char_print(G, chi);
abort();
}
do
n = n_randint(state, q);
while (n_gcd(q, n) > 1);
acb_dirichlet_char(chi2, G, n);
pairing = acb_dirichlet_ui_pairing(G, m, n);
if (pairing != acb_dirichlet_ui_chi(G, chi, n) * (G->expo / chi->order.n)
|| pairing != acb_dirichlet_ui_chi(G, chi2, m) * (G->expo / chi2->order.n))
{
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);
abort();
}
acb_dirichlet_conrey_next(x, G);
acb_dirichlet_char_next(chi, G);
acb_dirichlet_char_conrey(chi2, G, x);
if (!acb_dirichlet_char_eq(G, chi, chi2))
{
flint_printf("FAIL: next char\n\n");
flint_printf("q = %wu\n\n", q);
flint_printf("m = %wu\n\n", m);
acb_dirichlet_char_print(G, chi);
flint_printf("\n");
acb_dirichlet_char_print(G, chi2);
flint_printf("\n");
abort();
}
}
acb_dirichlet_char_clear(chi);
acb_dirichlet_char_clear(chi2);
acb_dirichlet_conrey_clear(x);
acb_dirichlet_group_clear(G);
}
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

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

@ -24,12 +24,14 @@ int main()
{
acb_t zn1, zn2, zn1n2, zn1zn2;
acb_dirichlet_group_t G;
acb_dirichlet_char_t chi;
ulong q, m, n1, n2, iter2;
int res;
q = 1 + n_randint(state, 1000);
acb_dirichlet_group_init(G, q);
acb_dirichlet_char_init(chi, G);
acb_init(zn1);
acb_init(zn2);
acb_init(zn1n2);
@ -42,12 +44,33 @@ int main()
m = 1 + n_randint(state, q);
} while (n_gcd(q, m) != 1);
acb_dirichlet_char(chi, G, m);
n1 = n_randint(state, 1000);
n2 = n_randint(state, 1000);
acb_dirichlet_chi(zn1, G, m, n1, 53);
acb_dirichlet_chi(zn2, G, m, n2, 53);
acb_dirichlet_chi(zn1n2, G, m, n1 * n2, 53);
acb_dirichlet_chi(zn1, G, chi, n1, 53);
acb_dirichlet_pairing(zn2, G, m, n1, 53);
if (!acb_overlaps(zn1, zn2))
{
acb_dirichlet_conrey_t x;
flint_printf("FAIL: overlap\n\n");
flint_printf("q = %wu\n\n", q);
flint_printf("m = %wu\n\n", m);
flint_printf("n = %wu\n\n", n1);
flint_printf("char = "); acb_printd(zn1, 15); flint_printf("\n\n");
flint_printf("pairing = "); acb_printd(zn2, 15); flint_printf("\n\n");
acb_dirichlet_char_print(G, chi);
acb_dirichlet_conrey_init(x, G);
acb_dirichlet_conrey_log(x, G, m);
flint_printf("log(m) = "); acb_dirichlet_conrey_print(G, x);
acb_dirichlet_conrey_log(x, G, n1);
flint_printf("log(n1) = "); acb_dirichlet_conrey_print(G, x);
abort();
}
acb_dirichlet_pairing(zn2, G, m, n2, 53);
acb_dirichlet_pairing(zn1n2, G, m, n1 * n2, 53);
acb_mul(zn1zn2, zn1, zn2, 53);
if (!acb_overlaps(zn1n2, zn1zn2))
@ -74,7 +97,7 @@ int main()
{
if (n_gcd(q, m) == 1)
{
acb_dirichlet_chi(zn2, G, m, n1, 53);
acb_dirichlet_pairing(zn2, G, m, n1, 53);
acb_add(zn1, zn1, zn2, 53);
}
}
@ -108,4 +131,3 @@ int main()
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

164
acb_dirichlet/test/t-conrey.c Normal file
View file

@ -0,0 +1,164 @@
/*
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 "acb_dirichlet.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("conrey....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 3000; iter++)
{
acb_dirichlet_group_t G;
acb_dirichlet_conrey_t x, y;
ulong q, n, k, sum;
long ref;
/*int * bits;*/
q = 1 + n_randint(state, 1000 * (1 + iter / 100));
acb_dirichlet_group_init(G, q);
acb_dirichlet_conrey_init(x, G);
acb_dirichlet_conrey_init(y, G);
/* check group size and elements */
acb_dirichlet_conrey_one(x, G);
sum = 1;
#if 1
for (n = 1; acb_dirichlet_conrey_next(x, G) >= 0; n++)
sum += x->n * x->n;
#else
/* iteration much faster than gcd below */
n = 1;
for (k = 2; k < G->q; k++)
{
if (n_gcd(k, G->q) > 1)
continue;
n++;
sum += k * k;
}
#endif
/* 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)
continue;
acb_dirichlet_conrey_first_primitive(x, G);
for (n = 1; acb_dirichlet_conrey_next_primitive(x, G) >= 0; n++);
ref = acb_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));
acb_dirichlet_conrey_log(x, G, m);
if (m != acb_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 = ");
acb_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 = acb_dirichlet_conrey_exp(x, G);
acb_dirichlet_conrey_log(y, G, m);
if (!acb_dirichlet_conrey_eq(G, x, y))
{
flint_printf("FAIL: conrey log and exp\n\n");
flint_printf("q = %wu\n\n", q);
flint_printf("conrey = ");
acb_dirichlet_conrey_print(G, x);
flint_printf("m = %wu\n\n", m);
flint_printf("log = ");
acb_dirichlet_conrey_print(G, y);
flint_printf("\n\nnumber = %wu\n\n", y->n);
abort();
}
acb_dirichlet_conrey_next_primitive(x, G);
m = x->n;
if (m != acb_dirichlet_conrey_exp(x, G))
{
flint_printf("FAIL: conrey number next primitive\n\n");
flint_printf("q = %wu\n\n", q);
flint_printf("conrey = ");
acb_dirichlet_conrey_print(G, y);
flint_printf(", m = %wu\n\n", y->n);
flint_printf("next primitive = ");
acb_dirichlet_conrey_print(G, x);
flint_printf(", m = %wu\n\n", m);
flint_printf("exp = %wu\n\n", x->n);
abort();
}
}
acb_dirichlet_conrey_clear(x);
acb_dirichlet_conrey_clear(y);
acb_dirichlet_group_clear(G);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

150
acb_dirichlet/test/t-dft.c Normal file
View file

@ -0,0 +1,150 @@
/*
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 "acb_dirichlet.h"
int main()
{
slong k;
slong prec = 100, digits = 30;
slong nq = 13;
ulong q[13] = { 2, 3, 4, 5, 6, 23, 10, 15, 30, 59, 308, 335, 961};
flint_rand_t state;
flint_printf("dft....");
fflush(stdout);
flint_randinit(state);
/* cyclic dft */
for (k = 0; k < nq; k++)
{
slong i;
acb_ptr v, w1, w2;
v = _acb_vec_init(q[k]);
w1 = _acb_vec_init(q[k]);
w2 = _acb_vec_init(q[k]);
for (i = 0; i < q[k]; i++)
acb_set_si(v + i, i);
acb_dirichlet_dft_pol(w1, v, q[k], prec);
acb_dirichlet_dft_fast(w2, v, q[k], prec);
for (i = 0; i < q[k]; i++)
{
if (!acb_overlaps(w1 + i, w2 + i))
{
flint_printf("differ from index %ld / %ld \n\n",i,q[k]);
flint_printf("pol =\n");
acb_vec_printd(w1, q[k], digits);
flint_printf("fast =\n");
acb_vec_printd(w2, q[k], digits);
flint_printf("\n\n");
abort();
}
}
_acb_vec_clear(v, q[k]);
_acb_vec_clear(w1, q[k]);
_acb_vec_clear(w2, q[k]);
}
/* Dirichlet group DFT */
for (k = 0; k < nq; k++)
{
slong i, j, len;
acb_dirichlet_group_t G;
acb_dirichlet_conrey_t x, y;
acb_t chiy;
acb_ptr v, w1, w2;
acb_dirichlet_group_init(G, q[k]);
len = G->phi_q;
v = _acb_vec_init(len);
w1 = _acb_vec_init(len);
w2 = _acb_vec_init(len);
acb_dirichlet_conrey_init(x, G);
acb_dirichlet_conrey_one(x, G);
for (i = 0; i < len; i++)
acb_randtest_precise(v + i, state, prec, 0);
/* naive */
acb_init(chiy);
acb_dirichlet_conrey_init(y, G);
acb_dirichlet_conrey_one(x, G);
for (i = 0; i < len; i++)
{
acb_zero(w1 + i);
acb_dirichlet_conrey_one(y, G);
for (j = 0; j < len; j++)
{
acb_dirichlet_pairing_conrey(chiy, G, x, y, prec);
acb_addmul(w1 + i, chiy, v + j, prec);
acb_dirichlet_conrey_next(y, G);
}
acb_dirichlet_conrey_next(x, G);
}
acb_clear(chiy);
acb_dirichlet_conrey_clear(y);
acb_dirichlet_conrey_clear(x);
/* dft */
acb_dirichlet_dft_conrey(w2, v, G, prec);
for (i = 0; i < len; i++)
{
if (!acb_overlaps(w1 + i, w2 + i))
{
flint_printf("FAIL\n\n");
flint_printf("q = %wu\n", q[k]);
flint_printf("v [size %wu]\n", len);
acb_vec_printd(v, len, digits);
flint_printf("\nDFT differ from index %ld / %ld \n", i, len);
flint_printf("\nnaive =\n");
acb_vec_printd(w1, len, digits);
flint_printf("\nfast =\n");
acb_vec_printd(w2, len, digits);
flint_printf("\n\n");
abort();
}
else if (acb_rel_accuracy_bits(w1 + i) < 30
|| acb_rel_accuracy_bits(w2 + i) < 30)
{
flint_printf("FAIL\n\n");
flint_printf("q = %wu\n", q[k]);
flint_printf("\nDFT inaccurate from index %ld / %ld \n", i, len);
flint_printf("\nnaive =\n");
acb_printd(w1 + i, digits);
flint_printf("\nfast =\n");
acb_printd(w2 + i, digits);
flint_printf("\nerrors %ld & %ld [prec = %wu]\n",
acb_rel_accuracy_bits(w1 + i),
acb_rel_accuracy_bits(w2 + i), prec
);
abort();
}
}
_acb_vec_clear(v, len);
_acb_vec_clear(w1, len);
_acb_vec_clear(w2, len);
acb_dirichlet_group_clear(G);
}
flint_randclear(state);
flint_printf("PASS\n");
}

86
acb_dirichlet/test/t-gauss.c Normal file
View file

@ -0,0 +1,86 @@
/*
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 "acb_dirichlet.h"
int main()
{
slong prec = 128;
ulong q;
flint_printf("gauss....");
fflush(stdout);
/* check Gauss sums */
for (q = 3; q < 250; q ++)
{
acb_dirichlet_group_t G;
acb_dirichlet_conrey_t x;
acb_dirichlet_char_t chi;
acb_t s1, s2, s3, s4;
acb_dirichlet_group_init(G, q);
acb_dirichlet_conrey_init(x, G);
acb_dirichlet_char_init(chi, G);
acb_init(s1);
acb_init(s2);
acb_init(s3);
acb_init(s4);
acb_dirichlet_conrey_one(x, G);
while (1) {
acb_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)
acb_dirichlet_gauss_sum_theta(s4, G, chi, prec);
else
acb_set(s4, s1);
if (!acb_overlaps(s1, s2)
|| !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);
acb_printd(s1, 25);
flint_printf("\ndefault ", q, x->n);
acb_printd(s2, 25);
flint_printf("\nfactor ", q, x->n);
acb_printd(s3, 25);
flint_printf("\ntheta ", q, x->n);
acb_printd(s4, 25);
abort();
}
if (acb_dirichlet_conrey_next(x, G) < 0)
break;
}
acb_clear(s1);
acb_clear(s2);
acb_clear(s3);
acb_clear(s4);
acb_dirichlet_group_clear(G);
acb_dirichlet_char_clear(chi);
acb_dirichlet_conrey_clear(x);
}
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

87
acb_dirichlet/test/t-jacobi.c Normal file
View file

@ -0,0 +1,87 @@
/*
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 "acb_dirichlet.h"
int main()
{
slong prec = 128;
ulong q;
flint_printf("jacobi....");
fflush(stdout);
/* check Jacobi sums */
for (q = 29 * 29; q > 1; q = q%2 ? 3*q+1 : q/2)
{
slong m1, m2;
acb_dirichlet_group_t G;
acb_dirichlet_char_t chi1, chi2;
acb_t s1, s2;
acb_dirichlet_group_init(G, q);
acb_dirichlet_char_init(chi1, G);
acb_dirichlet_char_init(chi2, G);
acb_init(s1);
acb_init(s2);
acb_dirichlet_char_one(chi1, G);
for (m1 = 0; m1 < 50; m1++)
{
acb_dirichlet_char_one(chi2, G);
for (m2 = 0; m2 < 50; m2++)
{
acb_dirichlet_jacobi_sum_naive(s1, G, chi1, chi2, prec);
acb_dirichlet_jacobi_sum(s2, G, chi1, chi2, prec);
if (!acb_overlaps(s1, s2))
{
flint_printf("FAIL: J_%wu(%wu,%wu)",
q, chi1->x->n, chi2->x->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,
acb_dirichlet_ui_conductor(G, nmod_mul(chi1->x->n, chi2->x->n, G->mod))
);
abort();
}
if (acb_dirichlet_char_next(chi2, G) < 0)
break;
}
if (acb_dirichlet_char_next(chi1, G) < 0)
break;
}
acb_clear(s1);
acb_clear(s2);
acb_dirichlet_group_clear(G);
acb_dirichlet_char_clear(chi1);
acb_dirichlet_char_clear(chi2);
}
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

219
acb_dirichlet/test/t-l.c Normal file
View file

@ -0,0 +1,219 @@
/*
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 "acb_dirichlet.h"
#define nq 5
#define nx 3
void
test_dft(ulong q)
{
ulong i;
slong prec = 100;
acb_dirichlet_group_t G;
acb_dirichlet_conrey_t x;
acb_dirichlet_char_t chi;
acb_t s, z;
acb_ptr v;
acb_dirichlet_group_init(G, q);
acb_dirichlet_conrey_init(x, G);
acb_dirichlet_char_init(chi, G);
acb_init(s);
acb_one(s);
acb_div_si(s, s, 2, prec);
v = _acb_vec_init(G->phi_q);
/* all at once */
acb_dirichlet_l_vec_hurwitz(v, s, G, prec);
/* check with complete loop */
i = 0;
acb_init(z);
acb_dirichlet_conrey_one(x, G);
do {
acb_dirichlet_char_conrey(chi, G, x);
acb_dirichlet_l_hurwitz(z, s, G, chi, prec);
if (!acb_overlaps(z, v + i))
{
flint_printf("\n L value differ");
flint_printf("\nL(1/2, %wu) single = ", x->n);
acb_printd(z, 20);
flint_printf("\nL(1/2, %wu) multi = ", x->n);
acb_printd(v + i, 20);
flint_printf("\n\n");
acb_vec_printd(v, G->phi_q, 10);
flint_printf("\n\n");
}
else if (acb_rel_accuracy_bits(z) < prec - 8
|| acb_rel_accuracy_bits(v + i) < prec - 8)
{
flint_printf("FAIL\n\n");
flint_printf("q = %wu\n", q);
flint_printf("\nL(1/2,chi_%wu(%wu,)) inaccurate\n", q, x->n);
flint_printf("\nsingle =\n");
acb_printd(z, 30);
flint_printf("\ndft =\n");
acb_printd(v + i, 30);
flint_printf("\nerrors %ld & %ld [prec = %wu]\n",
acb_rel_accuracy_bits(z),
acb_rel_accuracy_bits(v + i), prec);
abort();
}
i++;
} while (acb_dirichlet_conrey_next(x, G) >= 0);
acb_clear(s);
_acb_vec_clear(v, G->phi_q);
acb_dirichlet_char_clear(chi);
acb_dirichlet_conrey_clear(x);
acb_dirichlet_group_clear(G);
}
int main()
{
slong i, j;
ulong q[nq] = {3, 5, 61, 91, 800};
ulong m[nq] = {2, 4, 11, 2, 3};
slong prec = 150;
acb_ptr x;
/* cannot test at s = 1 with hurwitz */
const char * x_r[nx] = { "1", "0.5", "0.5" };
const char * x_i[nx] = { "1", "0", "6" };
acb_t ref, res;
/*
default(realprecision, 54)
X = [ 1 + I, 1/2, 1/2 + 6 * I ]
C = [Mod(2,3),Mod(4,5),Mod(11,61),Mod(2,91),Mod(3,800)]
v = concat([ [lfun(c,x) | x<-X] | c<-C])
apply(z->printf("\"%s\",\n",real(z)),v)
apply(z->printf("\"%s\",\n",imag(z)),v)
*/
const char * ref_r[nq * nx] = {
"0.655527984002548033786648216345221087939439503905627469",
"0.480867557696828626181220063235589877776829730832371863",
"1.56831301727577320813799211138797101541772722814204757",
"0.521271244517346991221550773660594765406476858135844321",
"0.231750947504015755883383661760877226427886696409005898",
"0.275543455389521803395512886745330595086898302178508437",
"0.598221809458540554839300433683735304093606595684903281",
"0.489264190003695740292779374874163221990017067040417393",
"0.573331076412428980263984182365365715292560207445592018",
"0.510279695870740409778738767334484809708615155539404548",
"0.635626509594367380604827545000418331455019188562281349",
"0.129304857274642475564179442785425797926079767522671163",
"1.18088858810025653590356481638012816019876881487868657",
"2.17175778983760437737667738885959799183430688287297767",
"3.41568550810774629867945639900431994221065497147578087"
};
const char * ref_i[nq * nx] = {
"0.220206044893215842652155131408935133596486560067476474",
"0",
"-0.969458654385732175077973304161399773236957587792986099",
"0.354614573267731219242838516784605303288232150760467423",
"0",
"-0.995392028773643947872231871832838543767598301887892796",
"1.04370497561090171487193145841005574472705644411957863",
"-0.108902811943905225853677097712717212629591264759957602",
"-0.232114369998608907337769019848201890558327186146689311",
"-0.133300066189980774635445078240315148544665020358019145",
"0.0119464572932630291870372694406253796888930803905106876",
"-0.567660589679294457801153713636532209809112025502518666",
"-0.654079942571300523223917775358845549990877148918886474",
"0.970337207245832214408380510793679653538607483205616894",
"-1.43652482351673593824956935036654893593947145947637807"
};
flint_printf("l....");
fflush(stdout);
x = _acb_vec_init(nx);
for (j = 0; j < nx; j++)
{
if (arb_set_str(acb_realref(x + j), x_r[j], prec) ||
arb_set_str(acb_imagref(x + j), x_i[j], prec) )
{
flint_printf("error while setting x[%ld] <- %s+I*%s\n",
j, x_r[j], x_i[j]);
abort();
}
}
acb_init(ref);
acb_init(res);
for (i = 0; i < nq; i++)
{
acb_dirichlet_group_t G;
acb_dirichlet_char_t chi;
acb_dirichlet_group_init(G, q[i]);
acb_dirichlet_char_init(chi, G);
acb_dirichlet_char(chi, G, m[i]);
for (j = 0; j < nx; j++)
{
if (arb_set_str(acb_realref(ref), ref_r[i * nx + j], prec - 10) ||
arb_set_str(acb_imagref(ref), ref_i[i * nx + j], prec - 10) )
{
flint_printf("error while setting ref <- %s+I*%s\n",
ref_r[i * nx + j], ref_i[i * nx + j]);
abort();
}
acb_dirichlet_l_hurwitz(res, x + j, G, chi, prec + 10);
if (!acb_contains(ref, res))
{
flint_printf("FAIL:\n\n");
flint_printf("q = %wu\n", q[i]);
flint_printf("m = %wu\n", m[i]);
flint_printf("x = ");
acb_printd(x, 54);
flint_printf("\nref = ");
acb_printd(ref, 54);
flint_printf("\nl(chi,x) = ");
acb_printd(res, 54);
flint_printf("\n\n");
abort();
}
}
acb_dirichlet_char_clear(chi);
acb_dirichlet_group_clear(G);
/* test using dft */
test_dft(q[i]);
}
acb_clear(ref);
acb_clear(res);
_acb_vec_clear(x, nx);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

120
acb_dirichlet/test/t-thetanull.c Normal file
View file

@ -0,0 +1,120 @@
/*
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 "acb_dirichlet.h"
int main()
{
slong prec = 64;
ulong q;
flint_printf("thetanull....");
fflush(stdout);
/* check the only theta functions
* theta(chi) = sum chi(k)* k^odd * exp(-Pi * k^2 / q)
* vanishing at 1 correspond to two specific
* characters of moduli 300 and 600 + conjugates
*/
for (q = 3; q < 1000; q ++)
{
acb_dirichlet_group_t G;
acb_dirichlet_char_t chi;
ulong * v, nv, k;
acb_t sum;
acb_ptr z;
arb_t eq;
arb_ptr t, kt, tt;
if (q % 4 == 2)
/* no primitive character mod q */
continue;
acb_dirichlet_group_init(G, q);
acb_dirichlet_char_init(chi, G);
z = _acb_vec_init(G->expo);
acb_dirichlet_vec_nth_roots(z, G->expo, prec);
nv = acb_dirichlet_theta_length_d(q, 1, prec);
v = flint_malloc(nv * sizeof(ulong));
arb_init(eq);
arb_const_pi(eq, prec);
arb_div_ui(eq, eq, q, prec);
arb_neg(eq, eq);
arb_exp(eq, eq, prec);
t = _arb_vec_init(nv);
acb_dirichlet_arb_quadratic_powers(t, nv, eq, prec);
kt = _arb_vec_init(nv);
for (k = 1; k < nv; k++)
arb_mul_ui(kt + k, t + k, k, prec);
/* theta function on primitive characters */
acb_init(sum);
acb_dirichlet_char_first_primitive(chi, G);
do {
ulong m;
acb_zero(sum);
acb_dirichlet_ui_chi_vec(v, G, chi, nv);
m = G->expo / chi->order.n;
tt = acb_dirichlet_char_parity(chi) ? kt : t;
for (k = 1; k < nv; k++)
if (v[k] != ACB_DIRICHLET_CHI_NULL)
acb_addmul_arb(sum, z + (v[k] * m), 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 (!acb_contains_zero(sum))
{
flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->x->n);
acb_printd(sum, 10);
flint_printf("\n");
acb_dirichlet_char_print(G, chi);
flint_printf("\n");
abort();
}
}
else if (acb_contains_zero(sum))
{
flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->x->n);
acb_printd(sum, 10);
flint_printf("\n");
acb_dirichlet_char_print(G, chi);
flint_printf("\n");
abort();
}
} while (acb_dirichlet_char_next_primitive(chi, G) >= 0);
_acb_vec_clear(z, G->expo);
_arb_vec_clear(t, nv);
acb_clear(sum);
arb_clear(eq);
flint_free(v);
acb_dirichlet_group_clear(G);
acb_dirichlet_char_clear(chi);
}
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

74
acb_dirichlet/test/t-vec.c Normal file
View file

@ -0,0 +1,74 @@
/*
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 "acb_dirichlet.h"
static ulong
vec_diff(ulong * v, ulong * ref, ulong nv)
{
ulong k;
for (k = 1; k < nv; k++)
if (ref[k] != v[k])
return k;
return 0;
}
int main()
{
ulong q;
flint_printf("vec....");
fflush(stdout);
for (q = 2; q < 600; q ++)
{
acb_dirichlet_group_t G;
acb_dirichlet_conrey_t x;
acb_dirichlet_char_t chi;
ulong * v1, * v2, nv, k;
acb_dirichlet_group_init(G, q);
acb_dirichlet_conrey_init(x, G);
acb_dirichlet_char_init(chi, G);
nv = 100;
v1 = flint_malloc(nv * sizeof(ulong));
v2 = flint_malloc(nv * sizeof(ulong));
acb_dirichlet_conrey_one(x, G);
do {
acb_dirichlet_char_conrey(chi, G, x);
acb_dirichlet_ui_chi_vec_loop(v1, G, chi, nv);
acb_dirichlet_ui_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);
abort();
}
} while (acb_dirichlet_conrey_next(x, G) >= 0);
flint_free(v1);
flint_free(v2);
acb_dirichlet_group_clear(G);
acb_dirichlet_char_clear(chi);
acb_dirichlet_conrey_clear(x);
}
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

39
acb_dirichlet/theta_length.c Normal file
View file

@ -0,0 +1,39 @@
/*
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 "acb_dirichlet.h"
#include <math.h>
#define PI 3.14159265358
#define LOG2 0.69314718055
ulong
acb_dirichlet_theta_length_d(ulong q, double t, slong prec)
{
double a, la;
a = PI / (double)q * t * t;
la = (a < .3) ? -log(2*a*(1-a)) : .8;
la = ((double)prec * LOG2 + la) / a;
return ceil(sqrt(la)+.5);
}
ulong
acb_dirichlet_theta_length(ulong q, const arb_t t, slong prec)
{
double dt;
ulong len;
arf_t at;
arf_init(at);
arb_get_lbound_arf(at, t, 53);
dt = arf_get_d(at, ARF_RND_DOWN);
len = acb_dirichlet_theta_length_d(q, dt, prec);
arf_clear(at);
return len;
}

32
acb_dirichlet/ui_chi.c Normal file
View file

@ -0,0 +1,32 @@
/*
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 "acb_dirichlet.h"
ulong
acb_dirichlet_ui_chi(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n)
{
if (n_gcd(G->q, n) > 1)
{
return ACB_DIRICHLET_CHI_NULL;
}
else
{
ulong v;
acb_dirichlet_conrey_t x;
acb_dirichlet_conrey_init(x, G);
acb_dirichlet_conrey_log(x, G, n);
v = acb_dirichlet_ui_chi_conrey(G, chi, x);
return v;
}
}

24
acb_dirichlet/ui_chi_conrey.c Normal file
View file

@ -0,0 +1,24 @@
/*
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 "acb_dirichlet.h"
ulong
acb_dirichlet_ui_chi_conrey(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const acb_dirichlet_conrey_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
acb_dirichlet/ui_chi_vec.c Normal file
View file

@ -0,0 +1,21 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_ui_chi_vec(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv)
{
if (3 * nv > G->q)
acb_dirichlet_ui_chi_vec_loop(v, G, chi, nv);
else
acb_dirichlet_ui_chi_vec_primeloop(v, G, chi, nv);
}

48
acb_dirichlet/ui_chi_vec_loop.c Normal file
View file

@ -0,0 +1,48 @@
/*
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 "acb_dirichlet.h"
/* loop over whole group */
void
acb_dirichlet_ui_chi_vec_loop(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv)
{
int j;
ulong t;
slong k;
acb_dirichlet_conrey_t x;
acb_dirichlet_conrey_init(x, G);
acb_dirichlet_conrey_one(x, G);
for (k = 0; k < nv; k++)
v[k] = ACB_DIRICHLET_CHI_NULL;
t = v[1] = 0;
while ( (j = acb_dirichlet_conrey_next(x, G)) >= 0 )
{
/* exponents were modified up to j */
for (k = G->num - 1; k >= j; k--)
t = nmod_add(t, chi->expo[k], chi->order);
if (x->n < nv)
v[x->n] = t;
}
/* fix result outside primes */
/*acb_dirichlet_vec_set_null(v, G, nv);*/
/* copy outside modulus */
for (k = G->q; k < nv ; k++ )
v[k] = v[k - G->q];
acb_dirichlet_conrey_clear(x);
}

69
acb_dirichlet/ui_chi_vec_primeloop.c Normal file
View file

@ -0,0 +1,69 @@
/*
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 "acb_dirichlet.h"
static void
chi_vec_evenpart(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_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
acb_dirichlet_ui_chi_vec_primeloop(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv)
{
slong k, l;
for (k = 1; k < nv; k++)
v[k] = 0;
if (G->neven)
chi_vec_evenpart(v, G, chi, nv);
for (l = G->neven; l < G->num; l++)
{
acb_dirichlet_prime_group_struct P = G->P[l];
/* FIXME: there may be some precomputed dlog in P if needed */
dlog_vec_add(v, nv, P.g, chi->expo[l], P.pe, P.phi, chi->order);
}
acb_dirichlet_ui_vec_set_null(v, G, nv);
}

53
acb_dirichlet/ui_conductor.c Normal file
View file

@ -0,0 +1,53 @@
/*
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 "acb_dirichlet.h"
/* TODO: product of local conductors */
ulong
acb_dirichlet_ui_conductor(const acb_dirichlet_group_t G, ulong a)
{
slong k;
ulong ap, cond;
cond = 1;
for (k = (G->neven == 2); k < G->num; k++)
{
ulong p;
nmod_t pe;
p = G->P[k].p;
pe = G->P[k].pe;
ap = a % pe.n;
if (ap == 1)
continue;
if (p == 2)
{
cond = 4;
if (a % 4 == 3)
ap = pe.n - ap;
}
else
{
cond *= p;
ap = nmod_pow_ui(ap, p - 1, pe);
}
while (ap != 1)
{
cond *= p;
ap = nmod_pow_ui(ap, p, pe);
}
}
return cond;
}

41
acb_dirichlet/ui_order.c Normal file
View file

@ -0,0 +1,41 @@
/*
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 "acb_dirichlet.h"
/* order of an element knowing the factorization of a multiple */
ulong
nmod_order_precomp(ulong a, nmod_t mod, ulong expo, n_factor_t fac)
{
int k;
ulong pe, ap, order = 1;
for (k = 0; k < fac.num; k++)
{
pe = n_pow(fac.p[k], fac.exp[k]);
ap = nmod_pow_ui(a, expo / pe, mod);
while ( ap != 1)
{
ap = nmod_pow_ui(ap, fac.p[k], mod);
order *= fac.p[k];
}
}
return order;
}
ulong
acb_dirichlet_ui_order(const acb_dirichlet_group_t G, ulong a)
{
n_factor_t fac;
n_factor_init(&fac);
n_factor(&fac, G->expo, 1);
return nmod_order_precomp(a, G->mod, G->expo, fac);
}

34
acb_dirichlet/ui_pairing.c Normal file
View file

@ -0,0 +1,34 @@
/*
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 "acb_dirichlet.h"
ulong
acb_dirichlet_ui_pairing(const acb_dirichlet_group_t G, ulong m, ulong n)
{
ulong x;
acb_dirichlet_conrey_t a, b;
if (n_gcd(G->q, m) > 1 || n_gcd(G->q, n) > 1)
return ACB_DIRICHLET_CHI_NULL;
acb_dirichlet_conrey_init(a, G);
acb_dirichlet_conrey_init(b, G);
acb_dirichlet_conrey_log(a, G, m);
acb_dirichlet_conrey_log(b, G, n);
x = acb_dirichlet_ui_pairing_conrey(G, a, b);
acb_dirichlet_conrey_clear(a);
acb_dirichlet_conrey_clear(b);
return x;
}

28
acb_dirichlet/ui_pairing_conrey.c Normal file
View file

@ -0,0 +1,28 @@
/*
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 "acb_dirichlet.h"
/* todo: modular arithmetic */
ulong
acb_dirichlet_ui_pairing_conrey(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t a, const acb_dirichlet_conrey_t b)
{
ulong x, k;
x = 0;
for (k = 0; k < G->num; k++)
x = (x + G->PHI[k] * a->log[k] * b->log[k]) % G->expo;
return x;
}

28
acb_dirichlet/ui_parity.c Normal file
View file

@ -0,0 +1,28 @@
/*
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 "acb_dirichlet.h"
int
acb_dirichlet_ui_parity(const acb_dirichlet_group_t G, ulong a)
{
int par;
par = 0;
if (G->neven && a % 4 == 3)
par++;
if (n_jacobi_unsigned(a, G->rad_q) == -1)
par++;
return par % 2;
}

25
acb_dirichlet/ui_theta_arb.c Normal file
View file

@ -0,0 +1,25 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_ui_theta_arb(acb_t res, const acb_dirichlet_group_t G, ulong a, const arb_t t, slong prec)
{
acb_dirichlet_char_t chi;
acb_dirichlet_char_init(chi, G);
acb_dirichlet_char(chi, G, a);
acb_dirichlet_chi_theta_arb(res, G, chi, t, prec);
acb_dirichlet_char_clear(chi);
}

31
acb_dirichlet/ui_vec_set_null.c Normal file
View file

@ -0,0 +1,31 @@
/*
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 "acb_dirichlet.h"
void
acb_dirichlet_ui_vec_set_null(ulong *v, const acb_dirichlet_group_t G, slong nv)
{
slong k, l;
if (G->q_even > 1)
{
for (k = 2; k < nv; k += 2)
v[k] = -1;
}
for (l = G->neven; l < G->num; l++)
{
ulong p = G->P[l].p;
for (k = p; k < nv; k += p)
v[k] = ACB_DIRICHLET_CHI_NULL;
}
}

91
acb_dirichlet/vec_nth_root.c Normal file
View file

@ -0,0 +1,91 @@
/*
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 "acb_dirichlet.h"
/* assume xs has size >= len * step */
void
_acb_vec_set_powers_step(acb_ptr xs, slong n, slong len, slong step, slong prec)
{
slong i, j;
prec += n_clog(n, 2);
for (i = 0, j = 0; i < len; i++, j += step)
{
if (i == 0)
acb_one(xs + j);
else if (i == 1)
acb_dirichlet_nth_root(xs + j, n, prec);
else if (i % 2 == 0)
acb_mul(xs + j, xs + j / 2, xs + j / 2, prec);
else
acb_mul(xs + j, xs + j - step, xs + step, prec);
}
}
/* assume len = p^e and z has size >= len * step */
void
_acb_vec_roots_pe(acb_ptr z, slong p, slong e, slong len, slong step, slong prec)
{
if (e <= 1)
{
_acb_vec_set_powers_step(z, p, len, step, prec);
}
else
{
slong q, r;
_acb_vec_roots_pe(z, p, e - 1, len / p, step * p, prec);
_acb_vec_set_powers_step(z, n_pow(p, e), p, step, prec);
for (q = p; q < len; q += p)
for (r = 1; r < p; r++)
acb_mul(z + (q + r) * step, z + q * step, z + r * step, prec);
}
}
void
acb_dirichlet_vec_nth_roots(acb_ptr z, slong len, slong prec)
{
slong i, q;
n_factor_t fac;
acb_one(z + 0);
n_factor_init(&fac);
n_factor(&fac, len, 0);
q = 1;
for (i = 0; i < fac.num; i++)
{
slong p, e, pe, mp, mq;
p = fac.p[i];
e = fac.exp[i];
pe = n_pow(p, e);
mp = len / pe;
mq = len / q;
_acb_vec_roots_pe(z, p, e, pe, mp, prec);
if (i > 0)
{
slong k, l;
for (k = mp; k < len; k += mp)
{
for (l = mq; l < len - k; l += mq)
acb_mul(z + k + l, z + k, z + l, prec);
for (; l < len; l += mq)
acb_mul(z + k + l - len, z + k, z + l, prec);
}
}
q *= pe;
}
}

273
dlog.h Normal file
View file

@ -0,0 +1,273 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2016 Pascal Molin
******************************************************************************/
#ifndef DLOG_H
#define DLOG_H
#ifdef DLOG_INLINES_C
#define DLOG_INLINE
#else
#define DLOG_INLINE static __inline__
#endif
#include "flint/ulong_extras.h"
#include "flint/nmod_vec.h"
enum
{
DLOG_MODPE, DLOG_CRT, DLOG_POWER, DLOG_BSGS, DLOG_TABLE, DLOG_23
};
typedef struct dlog_precomp_struct dlog_precomp_struct;
typedef struct dlog_precomp_struct * dlog_precomp_ptr;
/* log in (1+pZ/p^eZ), e small: use recursion formulas
* could use padic log instead but exponent is small
* for ulongs */
typedef struct
{
ulong inv1p; /* 1 / (1 + p) */
ulong invloga1; /* 1 / log(a^(p-1),1+p) */
}
dlog_1modpe_struct;
typedef dlog_1modpe_struct dlog_1modpe_t[1];
/* log in (Z/p^eZ)^* */
typedef struct
{
ulong p;
ulong e;
ulong pe1; /* p^(e-1) */
ulong inva;
nmod_t pe;
dlog_precomp_struct * modp;
dlog_1modpe_t modpe;
}
dlog_modpe_struct;
typedef dlog_modpe_struct dlog_modpe_t[1];
/* all logs precomputed in (Z/modZ)^ast */
typedef struct
{
ulong mod;
ulong * table;
}
dlog_table_struct;
typedef dlog_table_struct dlog_table_t[1];
/* bsgs table, already in flint */
typedef struct apow {
ulong k;
ulong ak;
} apow_t;
int apow_cmp(const apow_t * x, const apow_t * y);
typedef struct {
nmod_t mod;
ulong m;
ulong am;
ulong g;
apow_t * table;
} dlog_bsgs_struct;
typedef dlog_bsgs_struct dlog_bsgs_t[1];
/* typedef bsgs_t dlog_bsgs_t; */
/* Pollard rho */
typedef struct {
ulong a;
nmod_t n;
nmod_t mod;
int nisprime;
} dlog_rho_struct;
typedef dlog_rho_struct dlog_rho_t[1];
/* CRT decomposition (Pohlig-Hellman) */
typedef struct
{
nmod_t mod;
nmod_t n;
ulong num;
ulong * expo;
ulong * crt_coeffs;
dlog_precomp_ptr pre;
}
dlog_crt_struct;
typedef dlog_crt_struct dlog_crt_t[1];
/* dlog when generator has prime power order */
typedef struct
{
nmod_t mod;
ulong p;
ulong e;
ulong * apk;
dlog_precomp_struct * pre;
}
dlog_power_struct;
typedef dlog_power_struct dlog_power_t[1];
typedef ulong dlog_order23_t[1];
/* generic decomposition */
/*typedef */
struct
dlog_precomp_struct
{
int type;
ulong cost;
union
{
dlog_table_t table;
dlog_bsgs_t bsgs;
dlog_crt_t crt;
dlog_power_t power;
dlog_modpe_t modpe;
dlog_order23_t order23;
} t;
};
typedef dlog_precomp_struct dlog_precomp_t[1];
void dlog_precomp_modpe_init(dlog_precomp_t pre, ulong a, ulong p, ulong e, ulong pe, ulong num);
void dlog_precomp_small_init(dlog_precomp_t pre, ulong a, ulong mod, ulong n, ulong num);
void dlog_precomp_n_init(dlog_precomp_t pre, ulong a, ulong mod, ulong n, ulong num);
void dlog_precomp_p_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong num);
void dlog_precomp_pe_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong e, ulong pe, ulong num);
void dlog_precomp_clear(dlog_precomp_t pre);
ulong dlog_precomp(const dlog_precomp_t pre, ulong b);
ulong dlog_order23_init(dlog_order23_t t, ulong a);
ulong dlog_table_init(dlog_table_t t, ulong a, ulong mod);
ulong dlog_crt_init(dlog_crt_t t, ulong a, ulong mod, ulong n, ulong num);
ulong dlog_power_init(dlog_power_t t, ulong a, ulong mod, ulong p, ulong e, ulong num);
ulong dlog_modpe_init(dlog_modpe_t t, ulong a, ulong p, ulong e, ulong pe, ulong num);
ulong dlog_bsgs_init(dlog_bsgs_t t, ulong a, ulong mod, ulong n, ulong m);
void dlog_1modpe_init(dlog_1modpe_t t, ulong a1, ulong p, ulong e, nmod_t pe);
void dlog_rho_init(dlog_rho_t t, ulong a, ulong mod, ulong n);
/*#define dlog_bsgs_init(t, a, n, m) bsgs_table_init(t, a, n, m)*/
DLOG_INLINE void
dlog_order23_clear(dlog_order23_t t)
{
return;
}
DLOG_INLINE void
dlog_table_clear(dlog_table_t t)
{
flint_free(t->table);
}
void dlog_crt_clear(dlog_crt_t t);
DLOG_INLINE void
dlog_power_clear(dlog_power_t t)
{
flint_free(t->apk);
dlog_precomp_clear(t->pre);
flint_free(t->pre);
}
DLOG_INLINE void
dlog_modpe_clear(dlog_modpe_t t)
{
dlog_precomp_clear(t->modp);
flint_free(t->modp);
}
DLOG_INLINE void
dlog_bsgs_clear(dlog_bsgs_t t)
{
flint_free(t->table);
}
DLOG_INLINE void
dlog_rho_clear(dlog_rho_t t)
{
return;
}
DLOG_INLINE ulong
dlog_bsgs_size(ulong n, ulong num)
{
if (2 * num < n)
return (1 + n_sqrt(n)) * (1 + n_sqrt(num));
else
return n;
}
/*#define dlog_bsgs_clear(t) bsgs_table_clear(t)*/
ulong dlog_order23(const dlog_order23_t t, ulong b);
ulong dlog_table(const dlog_table_t t, ulong b);
ulong dlog_crt(const dlog_crt_t t, ulong b);
ulong dlog_power(const dlog_power_t t, ulong b);
ulong dlog_modpe(const dlog_modpe_t t, ulong b);
ulong dlog_mod2e(const dlog_modpe_t t, ulong b);
ulong dlog_bsgs(const dlog_bsgs_t t, ulong b);
ulong dlog_rho(const dlog_rho_t t, ulong b);
ulong dlog_1modpe_mod1p(ulong b1, ulong p, ulong e, ulong inv1p, nmod_t pe);
ulong dlog_1modpe(const dlog_1modpe_t t, ulong b1, ulong p, ulong e, nmod_t pe);
/*#define dlog_bsgs(t, b) n_discrete_log_bsgs_table(t, b)*/
#define DLOG_SMALL_LIM 50
#define DLOG_TABLE_LIM 50
#define DLOG_TABLE_P_LIM 50
#define DLOG_TABLE_MODPE_LIM 50
#define DLOG_TABLE_PE_LIM 50
#define DLOG_TABLE_N_LIM 50
#define DLOG_BSGS_LIM 500
#define DLOG_LOOP_MAX_FACTOR 6
#define DLOG_G_SMALL 0
#define DLOG_G_BIG 1
void dlog_n_factor_group(n_factor_t * fac, ulong bound);
#define DLOG_NOT_FOUND UWORD_MAX
#define DLOG_NONE UWORD_MAX
ulong dlog_vec_pindex_factorgcd(ulong * v, ulong nv, ulong p, nmod_t mod, ulong a, ulong na, ulong loga, ulong logm1, nmod_t order, int maxtry);
void dlog_vec_fill(ulong * v, ulong nv, ulong x);
void dlog_vec_set_not_found(ulong *v, ulong nv, nmod_t mod);
void dlog_vec_loop(ulong * v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
void dlog_vec_loop_add(ulong * v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
void dlog_vec_eratos_add(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
void dlog_vec_eratos(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
void dlog_vec_sieve_add(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
void dlog_vec_sieve(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
void dlog_vec_add(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
void dlog_vec(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
#endif

26
dlog/1modpe.c Normal file
View file

@ -0,0 +1,26 @@
/*
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 "dlog.h"
ulong
dlog_1modpe(const dlog_1modpe_t t, ulong b1, ulong p, ulong e, nmod_t pe)
{
if (e == 1)
return 0;
else
{
ulong logb1;
logb1 = dlog_1modpe_mod1p(b1, p, e, t->inv1p, pe);
/* only need mod p^(e-1) */
return nmod_mul(logb1, t->invloga1, pe);
}
}

32
dlog/1modpe_init.c Normal file
View file

@ -0,0 +1,32 @@
/*
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 "dlog.h"
void
dlog_1modpe_init(dlog_1modpe_t t, ulong a1, ulong p, ulong e, nmod_t pe)
{
if (e == 1)
{
t->inv1p = 1;
t->invloga1 = 0;
}
else
{
ulong loga1;
if (a1 == 1)
abort();
t->inv1p = nmod_inv(1 + p, pe); /* 1 - p + p^2 - ... */
loga1 = dlog_1modpe_mod1p(a1, p, e, t->inv1p, pe);
/* only need inverse mod p^(e-1) but does not hurt */
t->invloga1 = nmod_inv(loga1, pe);
}
}

39
dlog/1modpe_mod1p.c Normal file
View file

@ -0,0 +1,39 @@
/*
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 "dlog.h"
/* for odd prime p, assume b1 = 1 mod p */
ulong
dlog_1modpe_mod1p(ulong b1, ulong p, ulong e, ulong inv1p, nmod_t pe)
{
int f;
ulong x, xf, pf, pf1;
pf1 = 1;
pf = p;
x = 0;
for (f = 1; f < e; f++)
{
if (b1 % pf != 1)
{
flint_printf("ERROR dlog_1modpe_1modp: %wu %% %wu != 1 mod %wu\n\n",
b1, pf, pe.n);
abort();
}
xf = (b1 - 1) / pf;
xf = (xf % p) * pf1;
x += xf;
b1 = nmod_mul(b1, nmod_pow_ui(inv1p, xf, pe), pe);
pf1 = pf;
pf *= p;
}
return x;
}

34
dlog/bsgs.c Normal file
View file

@ -0,0 +1,34 @@
/*
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 <stdlib.h>
#include "dlog.h"
ulong
dlog_bsgs(const dlog_bsgs_t t, ulong b)
{
ulong i;
apow_t c, * x;
c.ak = b;
for (i = 0; i < t->g; i++)
{
x = bsearch(&c, t->table, t->m, sizeof(apow_t),
(int(*)(const void*,const void*))apow_cmp);
if (x != NULL)
return i * t->m + x->k;
c.ak = nmod_mul(c.ak, t->am, t->mod);
}
flint_printf("Exception (dlog_bsgs). discrete log not found.\n");
flint_printf(" table size %wu, cosize %wu mod %wu. %wu not found (a^-m=%wu)\n",
t->m, t->g, t->mod.n, b, t->am);
abort();
}

43
dlog/bsgs_init.c Normal file
View file

@ -0,0 +1,43 @@
/*
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 <stdlib.h>
#include "dlog.h"
int
apow_cmp(const apow_t * x, const apow_t * y)
{
return (x->ak < y->ak) ? -1 : (x->ak > y->ak);
}
/* set size of table m=sqrt(nk) to compute k logs in a group of size n */
ulong
dlog_bsgs_init(dlog_bsgs_t t, ulong a, ulong mod, ulong n, ulong m)
{
ulong k, ak;
if (m >= n) m = n + 1;
t->table = (apow_t *)flint_malloc(m * sizeof(apow_t));
nmod_init(&t->mod, mod);
t->m = m;
t->g = n / m + 1;
for (k = 0, ak = 1; k < m; k++)
{
t->table[k].k = k;
t->table[k].ak = ak;
ak = nmod_mul(ak, a, t->mod);
}
t->am = nmod_inv(ak, t->mod);
qsort(t->table, m, sizeof(apow_t), (int(*)(const void*,const void*))apow_cmp);
return t->g;
}

31
dlog/crt.c Normal file
View file

@ -0,0 +1,31 @@
/*
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 "dlog.h"
ulong
dlog_crt(const dlog_crt_t t, ulong b)
{
int k;
ulong r = 0;
for (k = 0; k < t->num; k++)
{
ulong bk, rk;
bk = nmod_pow_ui(b, t->expo[k], t->mod);
rk = dlog_precomp(t->pre + k, bk);
#if 0
flint_printf("##[crt-%d]: log(%wu)=log(%wu^%wu) = %wu [size %wu mod %wu]\n",
k, bk, b, t->expo[k], rk, t->n/t->expo[k], t->mod);
#endif
r = nmod_add(r, nmod_mul(t->crt_coeffs[k], rk, t->n), t->n);
}
return r;
}

23
dlog/crt_clear.c Normal file
View file

@ -0,0 +1,23 @@
/*
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 "dlog.h"
void
dlog_crt_clear(dlog_crt_t t)
{
int k;
flint_free(t->expo);
flint_free(t->crt_coeffs);
for (k = 0; k < t->num; k++)
dlog_precomp_clear(t->pre + k);
flint_free(t->pre);
}

58
dlog/crt_init.c Normal file
View file

@ -0,0 +1,58 @@
/*
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 "dlog.h"
ulong
dlog_crt_init(dlog_crt_t t, ulong a, ulong mod, ulong n, ulong num)
{
int k;
n_factor_t fac;
ulong * M, * u;
ulong cost = 0;
n_factor_init(&fac);
n_factor(&fac, n, 1);
t->num = fac.num;
nmod_init(&t->mod,mod);
nmod_init(&t->n, n);
M = t->expo = flint_malloc(t->num * sizeof(ulong));
u = t->crt_coeffs = flint_malloc(t->num * sizeof(ulong));
t->pre = flint_malloc(t->num * sizeof(dlog_precomp_struct));
for (k = 0; k < t->num; k++)
{
ulong p, e, mk;
p = fac.p[k];
e = fac.exp[k];
if (0 && mod % p == 0)
{
flint_printf("dlog_crt_init: modulus must be prime to order.\n");
abort();
}
mk = n_pow(p, e);
M[k] = n / mk;
u[k] = nmod_mul(M[k], n_invmod(M[k] % mk, mk), t->n);
/* depends on the power */
#if 0
flint_printf("[sub-crt -- init for size %wu mod %wu]\n", mk, mod);
#endif
dlog_precomp_pe_init(t->pre + k, nmod_pow_ui(a, M[k], t->mod), mod, p, e, mk, num);
cost += t->pre[k].cost;
}
#if 0
if (cost > 500)
flint_printf("[crt init for size %wu mod %wu -> cost %wu]\n", n,mod,cost);
#endif
return cost;
}

53
dlog/factor_group.c Normal file
View file

@ -0,0 +1,53 @@
/*
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 "dlog.h"
/* group components up to bound and return cofactor */
void
dlog_n_factor_group(n_factor_t * fac, ulong bound)
{
int i, j, k;
ulong m[FLINT_MAX_FACTORS_IN_LIMB];
ulong c = 1;
i = 0;
for (k = fac->num - 1; k >= 0; k--)
{
ulong qe = n_pow(fac->p[k], fac->exp[k]);
if (qe > bound)
c *= qe;
else
{
/* try to insert somewhere in m */
for (j = 0; j < i && (m[j] * qe > bound); j++);
if (j == i)
m[i++] = qe;
else
m[j] *= qe;
}
}
for (j = 0; j < i; j++)
{
fac->p[j] = m[j];
fac->exp[j] = DLOG_G_SMALL;
}
if (c > 1)
{
fac->p[i] = c;
fac->exp[i] = DLOG_G_BIG;
i++;
}
fac->num = i;
}

43
dlog/mod2e.c Normal file
View file

@ -0,0 +1,43 @@
/*
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 "dlog.h"
ulong
dlog_mod2e(const dlog_modpe_t t, ulong b1)
{
if (t->e == 2)
return (b1 % 4) == 3;
else
{
slong f;
ulong pf1, pf, x, xf;
pf1 = 1;
pf = 4;
x = 0;
for (f = 2; f < t->e; f++)
{
if (b1 % pf != 1)
{
flint_printf("ERROR dlog_mod2e: %wu %% %wu != 1 mod %wu\n\n",
b1, pf, t->pe.n);
abort();
}
xf = (b1 - 1) / pf;
xf = (f == 2) ? xf % 4 : (xf % 2) * (pf1 / 2);
b1 = nmod_mul(b1, nmod_pow_ui(t->inva, xf, t->pe), t->pe);
x += xf;
pf1 = pf;
pf *= 2;
}
return x;
}
}

41
dlog/modpe.c Normal file
View file

@ -0,0 +1,41 @@
/*
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 "dlog.h"
ulong
dlog_modpe(const dlog_modpe_t t, ulong b)
{
ulong x;
if (t->p == 2)
return dlog_mod2e(t, b);
x = dlog_precomp(t->modp, b % t->p);
if (t->e > 1)
{
ulong b1, y;
#if 0
b1 = nmod_mul(b, nmod_pow_ui(t->inva, x, t->pe), t->pe);
y = dlog_1modpe(t->modpe.rec, b1, t->p, t->e, t->pe);
y = y % t->pe1;
x = x + (t->p - 1) * y;
#else
b1 = nmod_pow_ui(b, t->p - 1, t->pe);
y = dlog_1modpe(t->modpe, b1, t->p, t->e, t->pe);
y = y % t->pe1;
x = n_submod(x, y % (t->p - 1), t->p - 1);
x = y + t->pe1 * x;
#endif
}
return x;
}

44
dlog/modpe_init.c Normal file
View file

@ -0,0 +1,44 @@
/*
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 "dlog.h"
/* todo: recursion could be made quadratic (not very useful for ulongs) */
ulong
dlog_modpe_init(dlog_modpe_t t, ulong a, ulong p, ulong e, ulong pe, ulong num)
{
ulong a1;
t->p = p;
t->e = e;
nmod_init(&t->pe, pe);
t->inva = nmod_inv(a, t->pe);
if (p == 2)
{
t->modp = NULL;
t->pe1 = (e <= 2) ? 2 : pe / 4;
t->modpe->inv1p = t->inva;
t->modpe->invloga1 = 1;
return e - 2;
}
else
{
t->modp = flint_malloc(sizeof(dlog_precomp_struct));
t->pe1 = pe / p;
dlog_precomp_n_init(t->modp, a, p, p - 1, num);
a1 = nmod_pow_ui(a, p - 1, t->pe);
dlog_1modpe_init(t->modpe, a1, p, e, t->pe);
return t->modp->cost + e;
}
}

Some files were not shown because too many files have changed in this diff Show more