some documentation + syntax

acb_dirichlet.h

@@ -1,15 +1,31 @@
-/*
+/*=============================================================================
+ *
+ *
+
+    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
@@ -108,6 +124,7 @@ int acb_dirichlet_conrey_eq(const acb_dirichlet_group_t G, const acb_dirichlet_c
 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);
@@ -196,15 +213,16 @@ void acb_dirichlet_ui_chi_vec(ulong *v, const acb_dirichlet_group_t G, const acb
 typedef struct
 {
     ulong order;
-    acb_ptr z;
-    ulong m;
-    ulong M;
-    acb_ptr Z;
+    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);
@@ -215,12 +233,15 @@ void acb_dirichlet_chi(acb_t res, const acb_dirichlet_group_t G, const acb_diric
 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 x, 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);
@@ -232,6 +253,28 @@ void acb_dirichlet_si_poly_evaluate(acb_t res, slong * v, slong len, const acb_t
 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

acb_dirichlet/chi_theta_arb.c

@@ -26,16 +26,16 @@
 #include "acb_dirichlet.h"
 #include "acb_poly.h"
 
-/* x = Pi / q * t^2 */
+/* q(t) = Pi / q * t^2 */
 static void
-acb_dirichlet_arb_theta_argt(arb_t x, ulong q, const arb_t t, slong prec)
+acb_dirichlet_arb_theta_xt(arb_t xt, ulong q, const arb_t t, slong prec)
 {
-    arb_const_pi(x, prec);
-    arb_div_ui(x, x, q, prec);
-    arb_mul(x, x, t, prec);
-    arb_mul(x, x, t, prec);
-    arb_neg(x, x);
-    arb_exp(x, x, 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
@@ -43,7 +43,7 @@ acb_dirichlet_chi_theta_arb(acb_t res, const acb_dirichlet_group_t G, const acb_
 {
     slong len;
     ulong * a;
-    arb_t x;
+    arb_t xt;
     acb_dirichlet_powers_t z;
 
     len = acb_dirichlet_theta_length(G->q, t, prec);
@@ -52,11 +52,12 @@ acb_dirichlet_chi_theta_arb(acb_t res, const acb_dirichlet_group_t G, const acb_
     acb_dirichlet_ui_chi_vec(a, G, chi, len);
     acb_dirichlet_powers_init(z, chi->order.n, len, prec);
 
-    arb_init(x);
-    acb_dirichlet_arb_theta_argt(x, G->q, t, prec);
-    acb_dirichlet_arb_theta_naive(res, x, chi->parity, a, z, len, prec);
+    arb_init(xt);
+    acb_dirichlet_arb_theta_xt(xt, G->q, t, prec);
+    /* TODO: switch to theta smallorder at some point */
+    acb_dirichlet_arb_theta_naive(res, xt, chi->parity, a, z, len, prec);
 
-    arb_clear(x);
+    arb_clear(xt);
     flint_free(a);
     acb_dirichlet_powers_clear(z);
 }

acb_dirichlet/power.c

@@ -26,15 +26,14 @@
 #include "acb_dirichlet.h"
 #include "acb_poly.h"
 
-/* assume won't be modified */
 void
 acb_dirichlet_power(acb_t z, const acb_dirichlet_powers_t t, ulong n, slong prec)
 {
     if (n < t->m)
     {
-        /* FIXME: I do not want to copy: is this all right? */
-        *z = *(t->z + n);
-        /* acb_set(z, t->z + n); */
+        /* TODO: could avoid copy ? line below does not work
+           z = *(t->z + n); */
+        acb_set(z, t->z + n);
     }
     else
     {

acb_dirichlet/powers_init.c

@@ -32,6 +32,7 @@ acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong num, slon
     ulong m;
     acb_t zeta;
     t->order = order;
+
     m = (num == 1) ? 1 : num * (prec / 64 + 1);
     if (m > order)
         m = order;
@@ -54,4 +55,5 @@ acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong num, slon
         t->M = 0;
         t->Z = NULL;
     }
+    acb_clear(zeta);
 }

acb_dirichlet/si_poly_evaluate.c

@@ -25,6 +25,7 @@
 
 #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)
 {

acb_dirichlet/theta_length.c

@@ -29,25 +29,25 @@
 #define LOG2 0.69314718055
 
 ulong
-acb_dirichlet_theta_length_d(ulong q, double x, slong prec)
+acb_dirichlet_theta_length_d(ulong q, double t, slong prec)
 {
     double a, la;
-    a = PI / (double)q * x * x;
+    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 x, slong prec)
+acb_dirichlet_theta_length(ulong q, const arb_t t, slong prec)
 {
-    double dx;
+    double dt;
     ulong len;
-    arf_t ax;
-    arf_init(ax);
-    arb_get_lbound_arf(ax, x, 53);
-    dx = arf_get_d(ax, ARF_RND_DOWN);
-    len = acb_dirichlet_theta_length_d(q, dx, prec);
-    arf_clear(ax);
+    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;
 }

acb_dirichlet/ui_theta_arb.c

@@ -0,0 +1,39 @@
+/*=============================================================================
+
+    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
+
+******************************************************************************/
+
+#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);
+}

doc/source/acb_dirichlet.rst

@@ -227,6 +227,42 @@ unity.
 
     There are no restrictions on *n*.
 
+Roots of unity
+-------------------------------------------------------------------------------
+
+.. function:: void acb_dirichlet_nth_root(acb_t res, ulong order, slong prec)
+
+   sets *res* to `\exp(\frac{2i\pi}{\mathrm{order}})` to precision *prec*.
+
+.. type:: acb_dirichlet_powers_struct
+
+.. type:: acb_dirichlet_powers_t
+
+   this structure allows to compute *n* powers of a fixed root of unity of order *m*
+   using precomputations. Extremal cases are
+
+   - all powers are stored: `O(m)` initialization + storage, `O(n)` eval
+
+   - nothing stored: `O(1)` initialization + storage, `O(\log(m)n)` eval
+
+   - `k` step decomposition: `O(k m^{\frac1k})` init + storage, `O((k-1)n)` eval.
+
+   Currently, only baby-step giant-step decomposition (i.e. `k=2`)
+   is implemented, allowing to obtain each power using one multiplication.
+
+.. function:: void acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong num, slong prec)
+
+   initialize the powers structure for *num* evaluations of powers of the root of unity
+   of order *order*.
+
+.. function:: void acb_dirichlet_powers_clear(acb_dirichlet_powers_t t)
+
+   clears *t*.
+
+.. function:: void acb_dirichlet_power(acb_t z, const acb_dirichlet_powers_t t, ulong n, slong prec)
+
+   sets *z* to `x^n` where *t* contains precomputed powers of `x`.
+
 Vector evaluation
 -------------------------------------------------------------------------------
 
@@ -251,6 +287,54 @@ Operations
 
    take the power of some character
 
+Theta sums
+-------------------------------------------------------------------------------
+
+We call Theta series of a Dirichlet character the quadratic series
+
+.. math::
+
+   \Theta_q(a) = \sum_{n\geq 0} \chi_q(a, n) n^p x^{n^2}
+
+where `p` is the parity of the character `\chi_q(a,\cdot)`.
+
+For `\Re(t)>0` we write `x(t)=\exp(-\frac{\pi}{N}t^2)` and define
+
+.. math::
+
+   \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_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)
+   
+   should be used, which is not done automatically (to avoid recomputing the
+   Gauss sum).
+
+.. function:: ulong acb_dirichlet_theta_length(ulong q, const arb_t t, slong prec)
+
+   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_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)
+
+   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*.
+
+
 Gauss and Jacobi sums
 -------------------------------------------------------------------------------