merge

acb_dirichlet.h

@@ -120,7 +120,6 @@ acb_dirichlet_conrey_copy(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t
 }
 
 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);
 ulong acb_dirichlet_conrey_order(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x);
@@ -184,6 +183,16 @@ 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);
+ACB_DIRICHLET_INLINE int
+acb_dirichlet_char_is_principal(const acb_dirichlet_char_t chi)
+{
+    return (chi->x->n == 1);
+}
+ACB_DIRICHLET_INLINE int
+acb_dirichlet_char_is_real(const acb_dirichlet_char_t chi)
+{
+    return (chi->order.n <= 2);
+}
 
 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);
@@ -236,7 +245,7 @@ void acb_dirichlet_chi_vec(acb_ptr v, const acb_dirichlet_group_t G, const acb_d
 void acb_dirichlet_arb_quadratic_powers(arb_ptr v, slong nv, const arb_t x, slong prec);
 void acb_dirichlet_qseries_arb(acb_t res, acb_srcptr a, const arb_t x, slong len, slong prec);
 void acb_dirichlet_qseries_arb_powers_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_qseries_arb_powers(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_qseries_arb_powers_smallorder(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, 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);

acb_dirichlet/char_init.c

@@ -14,6 +14,13 @@
 void
 acb_dirichlet_char_init(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
 {
+    slong k;
     acb_dirichlet_conrey_init(chi->x, G);
     chi->expo = flint_malloc(G->num * sizeof(ulong));
+    chi->q = G->q;
+    chi->conductor = 1;
+    chi->parity = 0;
+    for (k = 0; k < G->num; k++)
+        chi->expo[k] = 0;
+    nmod_init(&(chi->order), 1);
 }

acb_dirichlet/char_one.c

@@ -14,10 +14,12 @@
 void
 acb_dirichlet_char_one(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
 {
+    slong k;
     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);
+    for (k = 0; k < G->num; k++)
+        chi->expo[k] = 0;
+    nmod_init(&(chi->order), 1);
 }

acb_dirichlet/conrey.c

@@ -14,6 +14,7 @@
 void
 acb_dirichlet_conrey_init(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G) {
     x->log = flint_malloc(G->num * sizeof(ulong));
+    acb_dirichlet_conrey_one(x, G);
 }
 
 void

acb_dirichlet/qseries_arb_powers.c

@@ -51,7 +51,7 @@ acb_dirichlet_qseries_arb_powers_naive(acb_t res, const arb_t x, int parity, con
 
 /* small order, multiply by chi at the end */
 void
-acb_dirichlet_qseries_arb_powers(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec)
+acb_dirichlet_qseries_arb_powers_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;

acb_dirichlet/test/t-thetanull.c

@@ -106,6 +106,7 @@ int main()
         } while (acb_dirichlet_char_next_primitive(chi, G) >= 0);
 
         _acb_vec_clear(z, G->expo);
+        _arb_vec_clear(kt, nv);
         _arb_vec_clear(t, nv);
         acb_clear(sum);
         arb_clear(eq);

acb_dirichlet/theta_arb.c

@@ -22,7 +22,7 @@ _acb_dirichlet_theta_arb_smallorder(acb_t res, const acb_dirichlet_group_t G, co
     acb_dirichlet_ui_chi_vec(a, G, chi, len);
 
     _acb_dirichlet_powers_init(z, chi->order.n, 0, 0, prec);
-    acb_dirichlet_qseries_arb_powers(res, xt, chi->parity, a, z, len, prec);
+    acb_dirichlet_qseries_arb_powers_smallorder(res, xt, chi->parity, a, z, len, prec);
     acb_dirichlet_powers_clear(z);
 
     flint_free(a);

acb_dirichlet/theta_length.c

@@ -45,21 +45,24 @@ mag_tail_kexpk2_arb(mag_t res, const arb_t a, ulong n)
     mag_t m;
     mag_init(m);
     arb_get_mag_lower(m, a);
-    /* a < 1/4 ? */
+    /* a < 1/4 */
     if (mag_cmp_2exp_si(m, -2) <= 0)
     {
         mag_t c;
         mag_init(c);
-        mag_mul_ui(c, m, 2);
-        mag_addmul(c, c, c);
-        mag_mul_ui(res, m, n*n-n+1);
+        mag_mul_ui_lower(res, m, n*n-n+1);
         mag_expinv(res, res);
+        /* c = 2a(1+2a) */
+        mag_mul_ui(m, m, 2);
+        mag_set_ui(c, 1);
+        mag_add_lower(c, m, c);
+        mag_mul_lower(c, m, c);
         mag_div(res, res, c);
         mag_clear(c);
     }
     else
     {
-        mag_mul_ui(res, m, n*n-n-1);
+        mag_mul_ui_lower(res, m, n*n-n-1);
         mag_expinv(res, res);
         mag_mul_ui(res, res, 2);
     }

acb_dirichlet/ui_chi_vec_primeloop.c

@@ -51,7 +51,7 @@ acb_dirichlet_ui_chi_vec_primeloop(ulong *v, const acb_dirichlet_group_t G, cons
 {
     slong k, l;
 
-    for (k = 1; k < nv; k++)
+    for (k = 0; k < nv; k++)
         v[k] = 0;
 
     if (G->neven)

dlog/vec_sieve.c

@@ -91,6 +91,7 @@ dlog_vec_sieve(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod
                 smooth, limcount, mod.n, logcost, logcount, sievecount, missed);
 #endif
     n_primes_clear(iter);
+    dlog_precomp_clear(pre);
     for (k = mod.n + 1; k < nv; k++)
         v[k] = v[k - mod.n];
 }

doc/source/acb_dirichlet.rst

@@ -109,11 +109,27 @@ Conrey elements
 
 .. function:: int acb_dirichlet_conrey_next(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G)
 
-    This function allows to iterate on the elements of *G* looping on the *index*.
-    It produces elements in seemingly random *number* order. The return value
-    is the index of the last updated exponent of *x*, or *G->num* if the last
+    Sets *x* to the next conrey index in *G* with lexicographic ordering.
+    The return value
+    is the index of the last updated exponent of *x*, or *-1* if the last
     element has been reached.
 
+    This function allows to iterate on the elements of *G* looping on their *index*.
+    Note that it produces elements in seemingly random *number* order.
+
+    The following template can be used to loop over all elements *x* in *G*::
+
+        acb_conrey_one(x, G);
+        do {
+            /* use Conrey element x */
+        } while (acb_dirichlet_conrey_next(x, G) >= 0);
+
+
+.. function:: int acb_dirichlet_conrey_eq(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x, const acb_dirichlet_conrey_t y)
+
+   Return 1 if *x* equals *y*. This function checks both *number* and *index*,
+   writing ``(x->n == y->n)`` gives a faster check.
+
 Dirichlet characters
 -------------------------------------------------------------------------------
 
@@ -130,12 +146,12 @@ the group *G* is isomorphic to its dual.
 
 .. function:: ulong acb_dirichlet_ui_pairing_conrey(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t a, const acb_dirichlet_conrey_t b)
 
-    Compute the value of the Dirichlet pairing on numbers *m* and *n*, as
-    exponent modulo *G->expo*.
-    The second form takes the index *a* and *b*, and does not take discrete
-    logarithms.
+   Compute the value of the Dirichlet pairing on numbers *m* and *n*, as
+   exponent modulo *G->expo*.
+   The second form takes the Conrey index *a* and *b*, and does not take discrete
+   logarithms.
 
-    The returned value is the numerator of the actual value exponent mod the group exponent *G->expo*.
+   The returned value is the numerator of the actual value exponent mod the group exponent *G->expo*.
 
 Character type
 -------------------------------------------------------------------------------
@@ -165,6 +181,14 @@ Character type
 
     Sets *chi* to the Dirichlet character of Conrey index *x*.
 
+.. function:: int acb_dirichlet_char_eq(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2)
+
+   Return 1 if *chi1* equals *chi2*.
+
+.. function:: acb_dirichlet_char_is_principal(const acb_dirichlet_char_t chi)
+
+   Return 1 if *chi* is the principal character mod *q*.
+
 Character properties
 -------------------------------------------------------------------------------
 
@@ -205,6 +229,10 @@ No discrete log computation is performed.
    Return the order of `\chi_q(a,\cdot)` which is the order of `a\mod q`.
    This number is precomputed for the *char* type.
 
+.. function:: 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`).
+
 Character evaluation
 -------------------------------------------------------------------------------
 
@@ -333,11 +361,11 @@ For `\Re(t)>0` we write `x(t)=\exp(-\frac{\pi}{N}t^2)` and define
 
    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)
-   
+
    should be used, which is not done automatically (to avoid recomputing the
    Gauss sum).
 
@@ -346,13 +374,13 @@ For `\Re(t)>0` we write `x(t)=\exp(-\frac{\pi}{N}t^2)` and define
    Compute the number of terms to be summed in the theta series of argument *t*
    so that the tail is less than `2^{-\mathrm{prec}}`.
 
-.. function:: 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)
+.. function:: void acb_dirichlet_qseries_powers_naive(acb_t res, const arb_t x, int p, const ulong * a, const acb_dirichlet_powers_t z, slong len, slong prec)
 
-.. function:: 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)
+.. function:: void acb_dirichlet_qseries_powers_smallorder(acb_t res, const arb_t x, int p, const ulong * a, const acb_dirichlet_powers_t z, slong len, slong prec)
 
    Compute the series `\sum n^p z^{a_n} x^{n^2}` for exponent list *a*,
    precomputed powers *z* and parity *p* (being 0 or 1).
-   
+
    The *naive* version sums the series as defined, while the *smallorder*
    variant evaluates the series on the quotient ring by a cyclotomic polynomial
    before evaluating at the root of unity, ignoring its argument *z*.
@@ -385,7 +413,7 @@ the Fourier transform on Conrey labels as
    Compute the DFT of *v* using Conrey indices.
    This function assumes *v* and *w* are vectors
    of size *G->phi_q*, whose values correspond to a lexicographic ordering
-   of Conrey indices.
+   of Conrey indices (as obtained using :func:`acb_dirichlet_conrey_next`).
 
    For example, if `q=15`, the Conrey elements are stored in following
    order

doc/source/dlog.rst

@@ -71,7 +71,7 @@ Evaluation
 
 .. function:: ulong dlog_precomp(const dlog_precomp_t pre, ulong b)
 
-   Returns `\log(b)` for the group described in *pre*
+   Returns `\log(b)` for the group described in *pre*.
 
 Vector evaluations
 -------------------------------------------------------------------------------
@@ -100,14 +100,14 @@ Algorithms
 
 Several discrete logarithms strategies are implemented:
 
-- Complete lookup table for small groups
+- Complete lookup table for small groups.
 
-- Baby-step giant-step table
+- Baby-step giant-step table.
 
-combined with mathematical reductions
+combined with mathematical reductions:
 
 - Pohlig-Hellman decomposition (Chinese remainder decomposition on the
-  order of the group and base `p` decomposition for primepower order)
+  order of the group and base `p` decomposition for primepower order).
 
 - p-adic log for primepower modulus `p^e`.
 
@@ -115,7 +115,7 @@ For *dlog_vec* functions which compute the vector of discrete logarithms
 of successive integers `1\dots n`:
 
 - A simple loop on group elements avoiding all logarithms is done when
-  the group size is comparable with the number of elements requested
+  the group size is comparable with the number of elements requested.
 
 - Otherwise the logarithms are computed on primes and propagated by
   Eratosthene-like sieving on composite numbers.