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

acb_dirichlet.h

@@ -255,11 +255,15 @@ void _acb_dirichlet_theta_argument_at_arb(arb_t xt, ulong q, const arb_t t, slon
 void acb_dirichlet_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_order2(acb_t res, const acb_dirichlet_char_t chi, 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_root_number(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec);
+void acb_dirichlet_root_number(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);

acb_dirichlet/dft_fast.c

@@ -161,7 +161,7 @@ _acb_dft_cyc(acb_ptr w, acb_srcptr v, dft_cyc_step * cyc, slong num, slong prec)
         /* reorder */
         for (i = 0; i < m; i++)
             for (j = 0; j < M; j++)
-                w[j + M * i] = t[i + m * j];
+                acb_set(w + j + M * i, t + i + m * j);
         _acb_vec_clear(t, m * M);
     }
 }
@@ -216,7 +216,7 @@ _acb_dft_prod(acb_ptr w, acb_srcptr v, dft_prod_step * cyc, slong num, slong pre
         /* reorder */
         for (i = 0; i < m; i++)
             for (j = 0; j < M; j++)
-                w[j + M * i] = t[i + m * j];
+                acb_set(w + j + M * i, t + i + m * j);
         _acb_vec_clear(t, m * M);
     }
 }

acb_dirichlet/gauss_sum.c

@@ -74,6 +74,20 @@ 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)
 {
@@ -84,16 +98,7 @@ acb_dirichlet_gauss_sum(acb_t res, const acb_dirichlet_group_t G, const acb_diri
     }
     else if (chi->order.n <= 2)
     {
-        if (chi->parity)
+        acb_dirichlet_gauss_sum_order2(res, chi, prec);
-        {
-            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_theta.c

@@ -14,38 +14,21 @@
 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 */
+        flint_printf("gauss_sum_theta: non available for non primitive character"
-        acb_dirichlet_gauss_sum_naive(res, G, chi, prec);
+                "or exceptional characters chi_300(71,.), chi_300(131,.), "
-        return;
+                "chi_600(11,.) and chi_600(491,.)\n");
-    }
+        abort();
-
-    arb_one(x);
-    acb_dirichlet_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));
+        acb_t iq;
-        arb_sqrt_ui(acb_realref(res), G->q, prec);
+        acb_init(iq);
+        acb_dirichlet_gauss_sum_order2(iq, chi, prec);
+        _acb_dirichlet_root_number(res, G, chi, prec);
+        acb_mul(res, res, iq, prec);
+        acb_clear(iq);
     }
-
-    acb_mul(res, res, eps, prec);
-
-    arb_clear(x);
-    acb_clear(eps);
 }

acb_dirichlet/root_number.c

@@ -0,0 +1,52 @@
+/*
+    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_root_number(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);
+    arb_one(x);
+    acb_dirichlet_theta_arb(res, G, chi, x, prec);
+    acb_init(eps);
+    acb_conj(eps, res);
+    acb_div(res, res, eps, prec);
+
+    arb_clear(x);
+    acb_clear(eps);
+}
+
+void
+acb_dirichlet_root_number(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
+{
+    if (chi->conductor < G->q)
+    {
+        flint_printf("root number: need primitive character\n");
+        abort();
+    }
+    else if (G->num > 1)
+    {
+        acb_t iq;
+        acb_init(iq);
+        acb_dirichlet_gauss_sum_order2(iq, chi, prec);
+        acb_dirichlet_gauss_sum(res, G, chi, prec);
+        acb_div(res, res, iq, prec);
+        acb_clear(iq);
+    }
+    else
+    {
+        _acb_dirichlet_root_number(res, G, chi, prec);
+    }
+}