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

This commit is contained in:
Pascal 2016-09-14 15:29:04 +02:00
commit f3a3dfaf94
17 changed files with 64 additions and 73 deletions

View file

@ -1,31 +1,15 @@
/*=============================================================================
*
*
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
@ -203,8 +187,6 @@ void acb_dirichlet_char_denormalize(acb_dirichlet_char_t chi, const acb_dirichle
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);

View file

@ -29,7 +29,7 @@ acb_dirichlet_chi_vec(acb_ptr v, const acb_dirichlet_group_t G, const acb_dirich
if (a[k] != ACB_DIRICHLET_CHI_NULL)
acb_dirichlet_power(v + k, t, a[k], prec);
else
*(v + k) = *(v + 0);
acb_zero(v + k);
}
acb_dirichlet_powers_clear(t);

View file

@ -16,6 +16,6 @@ acb_dirichlet_conrey_pow(acb_dirichlet_conrey_t c, const acb_dirichlet_group_t G
{
ulong k;
for (k = 0; k < G->num ; k++)
c->log[k] = (a->log[k] * n) % G->P[k].phi;
c->log[k] = n_mulmod2(a->log[k], n, G->P[k].phi);
c->n = nmod_pow_ui(a->n, n, G->mod);
}

30
dlog.h
View file

@ -1,27 +1,13 @@
/*=============================================================================
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
******************************************************************************/
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 DLOG_H
#define DLOG_H

View file

@ -37,7 +37,7 @@ dlog_precomp_clear(dlog_precomp_t pre)
dlog_order23_clear(pre->t.order23);
break;
default:
printf("THE TYPE IS %d\n", pre->type);
flint_printf("dlog_precomp_clear: unknown type %d\n", pre->type);
abort();
break;
}

View file

@ -15,12 +15,13 @@
void
dlog_precomp_modpe_init(dlog_precomp_t pre, ulong a, ulong p, ulong e, ulong pe, ulong num)
{
if ( pe < DLOG_TABLE_MODPE_LIM )
if (pe < DLOG_TABLE_MODPE_LIM)
{
dlog_precomp_small_init(pre, a, pe, pe - pe / p, num);
return;
}
else {
else
{
if (e > 1)
{
pre->type = DLOG_MODPE;

View file

@ -18,7 +18,8 @@ dlog_precomp_n_init(dlog_precomp_t pre, ulong a, ulong mod, ulong n, ulong num)
{
if (n % 2 && n_is_probabprime(n))
dlog_precomp_p_init(pre, a, mod, n, num);
else {
else
{
if (n < DLOG_TABLE_N_LIM)
{
dlog_precomp_small_init(pre, a, mod, n, num);

View file

@ -10,13 +10,12 @@
*/
#include "dlog.h"
#include "math.h"
/* we known the order is prime */
void
dlog_precomp_p_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong num)
{
if ( p < DLOG_TABLE_P_LIM )
if (p < DLOG_TABLE_P_LIM)
{
dlog_precomp_small_init(pre, a, mod, p, num);
}

View file

@ -14,13 +14,13 @@
void
dlog_precomp_pe_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong e, ulong pe, ulong num)
{
if ( pe < DLOG_TABLE_PE_LIM )
if (pe < DLOG_TABLE_PE_LIM)
{
dlog_precomp_small_init(pre, a, mod, pe, num);
}
else
{
if ( e == 1)
if (e == 1)
{
dlog_precomp_p_init(pre, a, mod, p, num);
}

View file

@ -21,7 +21,7 @@ dlog_precomp_small_init(dlog_precomp_t pre, ulong a, ulong mod, ulong n, ulong n
}
else
{
if ( mod < DLOG_TABLE_LIM )
if (mod < DLOG_TABLE_LIM)
{
pre->type = DLOG_TABLE;
pre->cost = dlog_table_init(pre->t.table, a, mod);

View file

@ -10,7 +10,6 @@
*/
#include "dlog.h"
#include <math.h>
void
dlog_rho_init(dlog_rho_t t, ulong a, ulong mod, ulong n)

View file

@ -73,6 +73,10 @@ int main()
}
dlog_modpe_clear(modpe);
/* multiplication can overflow on 32-bit */
if ((double) pe * p > LIM)
break;
}
}

View file

@ -51,7 +51,7 @@ int main()
fflush(stdout);
flint_randinit(state);
for (bits = 10; bits <= 35; bits += 5)
for (bits = 10; bits <= FLINT_MIN(35, FLINT_BITS); bits += 5)
{
for (nv = 10; nv <= 10000; nv *= 10)

View file

@ -26,6 +26,6 @@ dlog_vec_loop(ulong * v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod
vx = nmod_add(vx, va, order);
}
while (x != 1);
for(x = mod.n + 1; x < nv; x++)
for (x = mod.n + 1; x < nv; x++)
v[x] = v[x - mod.n];
}

View file

@ -168,9 +168,12 @@ Character type
.. function:: void acb_dirichlet_char_init(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
Initializes *chi* to an element of the group *G* and sets its value
to the principal character.
.. function:: void acb_dirichlet_char_clear(acb_dirichlet_char_t chi)
Initializes and clear *chi*.
Clears *chi*.
.. function:: void acb_dirichlet_char(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, ulong n)
@ -189,6 +192,22 @@ Character type
Return 1 if *chi* is the principal character mod *q*.
.. function:: void acb_dirichlet_char_one(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
Sets *chi* to the principal character.
.. function:: int acb_dirichlet_char_next(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
Sets *x* to the next character in *G* with lexicographic Conrey ordering
(see :func:`acb_dirichlet_conrey_next`). The return value
is the index of the last updated exponent of *x*, or *-1* if the last
element has been reached.
.. function:: int acb_dirichlet_char_next_primitive(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
Like :func:`acb_dirichlet_char_next`, but only generates primitive
characters.
Character properties
-------------------------------------------------------------------------------
@ -226,10 +245,10 @@ No discrete log computation is performed.
.. function:: ulong acb_dirichlet_char_order(const acb_dirichlet_char_t chi)
Return the order of `\chi_q(a,\cdot)` which is the order of `a\mod q`.
Return the order of `\chi_q(a,\cdot)` which is the order of `a\bmod q`.
This number is precomputed for the *char* type.
.. function:: acb_dirichlet_char_is_real(const acb_dirichlet_char_t chi)
.. function:: int acb_dirichlet_char_is_real(const acb_dirichlet_char_t chi)
Return 1 if *chi* is a real character (iff it has order `\leq 2`).
@ -328,7 +347,7 @@ Gauss and Jacobi sums
.. math::
G_q(a) = \sum_{x \mod q} \chi_q(a, x)e^{\frac{2i\pi x}q}
G_q(a) = \sum_{x \bmod q} \chi_q(a, x)e^{\frac{2i\pi x}q}
.. function:: 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)
@ -336,7 +355,7 @@ Gauss and Jacobi sums
.. math::
J_q(a,b) = \sum_{x \mod q} \chi_q(a, x)\chi_q(b, 1-x)
J_q(a,b) = \sum_{x \bmod q} \chi_q(a, x)\chi_q(b, 1-x)
Theta sums
-------------------------------------------------------------------------------
@ -355,16 +374,16 @@ For `\Re(t)>0` we write `x(t)=\exp(-\frac{\pi}{N}t^2)` and define
\Theta_q(a,t) = \sum_{n\geq 0} \chi_q(a, n) x(t)^{n^2}.
.. function:: 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);
.. function:: 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)
.. function:: void acb_dirichlet_ui_theta_arb(acb_t res, const acb_dirichlet_group_t G, ulong a, const arb_t t, slong prec);
.. function:: void acb_dirichlet_ui_theta_arb(acb_t res, const acb_dirichlet_group_t G, ulong a, const arb_t t, slong prec)
Compute the theta series `\Theta_q(a,t)` for real argument `t>0`.
Beware that if `t<1` the functional equation
.. math::
t \theta(a,t) = \epsilon(\chi) \theta(\frac1a, \frac1t)
t \theta(a,t) = \epsilon(\chi) \theta\left(\frac1a, \frac1t\right)
should be used, which is not done automatically (to avoid recomputing the
Gauss sum).

View file

@ -141,7 +141,7 @@ Basic properties and manipulation
.. function:: void acb_poly_truncate(acb_poly_t poly, slong n)
Truncates *poly* to have length at most *n*, i.e. degree
strictly smaller than *n*.
strictly smaller than *n*. We require that *n* is nonnegative.
.. function:: slong acb_poly_valuation(const acb_poly_t poly)
@ -364,7 +364,7 @@ Arithmetic
.. function:: void _acb_poly_divrem(acb_ptr Q, acb_ptr R, acb_srcptr A, slong lenA, acb_srcptr B, slong lenB, slong prec)
.. function:: void acb_poly_divrem(acb_poly_t Q, acb_poly_t R, const acb_poly_t A, const acb_poly_t B, slong prec)
.. function:: int acb_poly_divrem(acb_poly_t Q, acb_poly_t R, const acb_poly_t A, const acb_poly_t B, slong prec)
Performs polynomial division with remainder, computing a quotient `Q` and
a remainder `R` such that `A = BQ + R`. The implementation reverses the

View file

@ -135,7 +135,7 @@ Basic manipulation
.. function:: void arb_poly_truncate(arb_poly_t poly, slong n)
Truncates *poly* to have length at most *n*, i.e. degree
strictly smaller than *n*.
strictly smaller than *n*. We require that *n* is nonnegative.
.. function:: slong arb_poly_valuation(const arb_poly_t poly)