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

split gauss and jacobi source files

Browse source
This commit is contained in:
Pascal 2016-09-16 16:53:56 +02:00
parent 71250d6d28
commit 3a9dc880e5
6 changed files with 146 additions and 110 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
acb_dirichlet.h
View file

@ -259,6 +259,7 @@ void acb_dirichlet_root_number(acb_t res, const acb_dirichlet_group_t G, const a
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);
ulong jacobi_one_prime(ulong p, ulong e, ulong pe, ulong cond);
void acb_dirichlet_jacobi_sum_factor(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_gauss(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);

15
acb_dirichlet/gauss_sum.c
View file

@ -74,21 +74,6 @@ gauss_sum_non_primitive(acb_t res, const acb_dirichlet_group_t G, const acb_diri
}
}
void
acb_dirichlet_gauss_sum_order2(acb_t res, const acb_dirichlet_char_t chi, slong prec)
{
if (chi->parity)
{
arb_zero(acb_realref(res));
arb_sqrt_ui(acb_imagref(res), chi->q, prec);
}
else
{
arb_zero(acb_imagref(res));
arb_sqrt_ui(acb_realref(res), chi->q, prec);
}
}
void
acb_dirichlet_gauss_sum(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
{

27
acb_dirichlet/gauss_sum_order2.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"
void
acb_dirichlet_gauss_sum_order2(acb_t res, const acb_dirichlet_char_t chi, slong prec)
{
if (chi->parity)
{
arb_zero(acb_realref(res));
arb_sqrt_ui(acb_imagref(res), chi->q, prec);
}
else
{
arb_zero(acb_imagref(res));
arb_sqrt_ui(acb_realref(res), chi->q, prec);
}
}

96
acb_dirichlet/jacobi_sum.c
View file

@ -12,7 +12,7 @@
#include "acb_dirichlet.h"
/* J_N(1,a) = sum on x = 1 mod some p | q */
static ulong
ulong
jacobi_one_prime(ulong p, ulong e, ulong pe, ulong cond)
{
if (e > 1 && cond % (p*p) == 0)
@ -40,103 +40,9 @@ jacobi_one(const acb_dirichlet_group_t G, ulong cond)
return r;
}
/* should use only for prime power modulus */
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);
}
void
acb_dirichlet_jacobi_sum_factor(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);

79
acb_dirichlet/jacobi_sum_factor.c Normal file
View file

@ -0,0 +1,79 @@
/*
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_factor(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);
}

38
acb_dirichlet/jacobi_sum_gauss.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"
/* should use only for prime power modulus */
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);
}