From 9441c452ece07597605637beca585147641db887 Mon Sep 17 00:00:00 2001
From: Pascal <molin.maths@gmail.com>
Date: Thu, 4 Aug 2016 16:32:31 +0200
Subject: [PATCH] fix dft on group

---
 acb_dirichlet/dft_fast.c   | 16 ++++++--
 acb_dirichlet/test/t-dft.c | 75 +++++++++++++++++++++++++++++++++++++-
 2 files changed, 87 insertions(+), 4 deletions(-)

diff --git a/acb_dirichlet/dft_fast.c b/acb_dirichlet/dft_fast.c
index 59df21df..fd92acf1 100644
--- a/acb_dirichlet/dft_fast.c
+++ b/acb_dirichlet/dft_fast.c
@@ -202,7 +202,10 @@ _acb_dft_prod(acb_ptr w, acb_srcptr v, dft_prod_step * cyc, slong num, slong pre
 {
     dft_prod_step c;
     if (num == 0)
+    {
+        flint_printf("error: reached num = 0 in dft_prod\n");
         abort(); /* or just copy v to w */
+    }
     c = cyc[0];
     if (num == 1)
     {
@@ -221,7 +224,7 @@ _acb_dft_prod(acb_ptr w, acb_srcptr v, dft_prod_step * cyc, slong num, slong pre
         {
             _acb_dft_prod(wi, vi, cyc + 1, num - 1, prec);
             wi += M;
-            vi += dv;
+            vi += dv; /* here = M */
         }
         /* after first pass */
         for (j = 0; j < M; j++)
@@ -237,7 +240,14 @@ _acb_dft_prod(acb_ptr w, acb_srcptr v, dft_prod_step * cyc, slong num, slong pre
 void
 acb_dirichlet_dft_prod_precomp(acb_ptr w, acb_srcptr v, acb_dft_prod_t t, slong prec)
 {
-    _acb_dft_prod(w, v, t->cyc, t->num, prec);
+    if (t->num == 0)
+    {
+        acb_set(w + 0, v + 0);
+    }
+    else
+    {
+        _acb_dft_prod(w, v, t->cyc, t->num, prec);
+    }
 }
 
 void
@@ -258,7 +268,7 @@ acb_dirichlet_dft_prod_init(acb_dft_prod_t t, slong * cyc, slong num, slong prec
         len /= m;
         t->cyc[i].m = m;
         t->cyc[i].M = len;
-        t->cyc[i].dv = dv;
+        t->cyc[i].dv = len;
         _acb_dirichlet_dft_cyc_init(t->cyc[i].pre, m, len, prec);
         dv *= m;
     }
diff --git a/acb_dirichlet/test/t-dft.c b/acb_dirichlet/test/t-dft.c
index c8c13c09..90278fd8 100644
--- a/acb_dirichlet/test/t-dft.c
+++ b/acb_dirichlet/test/t-dft.c
@@ -36,6 +36,7 @@ int main()
     flint_printf("dft....");
     fflush(stdout);
 
+    /* cyclic dft */
     for (k = 0; k < nq; k++)
     {
         slong i;
@@ -64,11 +65,83 @@ int main()
                 abort();
             }
         }
-        
+
         _acb_vec_clear(v, q[k]);
         _acb_vec_clear(w1, q[k]);
         _acb_vec_clear(w2, q[k]);
     }
 
+    /* Dirichlet group DFT */
+    for (k = 0; k < nq; k++)
+    {
+        slong i, j, len;
+        acb_dirichlet_group_t G;
+        acb_dirichlet_conrey_t x, y;
+        acb_t chiy;
+        acb_ptr v, w1, w2;
+
+        acb_dirichlet_group_init(G, q[k]);
+
+        len = G->phi_q;
+        v = _acb_vec_init(len);
+        w1 = _acb_vec_init(len);
+        w2 = _acb_vec_init(len);
+
+        acb_dirichlet_conrey_init(x, G);
+        acb_dirichlet_conrey_one(x, G);
+        for (i = 0; i < len; i++)
+        {
+            acb_set_si(v + i, x->n);
+            acb_dirichlet_conrey_next(x, G);
+        }
+
+        /* naive */
+        acb_init(chiy);
+        acb_dirichlet_conrey_init(y, G);
+        acb_dirichlet_conrey_one(x, G);
+        for (i = 0;  i < len; i++)
+        {
+            acb_zero(w1 + i);
+            acb_dirichlet_conrey_one(y, G);
+            for (j = 0; j < len; j++)
+            {
+                acb_dirichlet_pairing_conrey(chiy, G, x, y, prec);
+                acb_addmul(w1 + i, chiy, v + j, prec);
+                acb_dirichlet_conrey_next(y, G);
+            }
+            acb_dirichlet_conrey_next(x, G);
+        }
+        acb_clear(chiy);
+        acb_dirichlet_conrey_clear(y);
+        acb_dirichlet_conrey_clear(x);
+
+        /* dft */
+        acb_dirichlet_dft_conrey(w2, v, G, prec);
+
+        for (i = 0; i < len; i++)
+        {
+            if (!acb_overlaps(w1 + i, w2 + i))
+            {
+                flint_printf("FAIL\n\n");
+                flint_printf("q = %wu\n", q[k]);
+                flint_printf("v [size %wu]\n", len);
+                acb_vec_printd(v, len, 10);
+                flint_printf("\nDFT differ from index %ld / %ld \n", i, len);
+                flint_printf("\nnaive =\n");
+                acb_vec_printd(w1, len, 10);
+                flint_printf("\nfast =\n");
+                acb_vec_printd(w2, len, 10);
+                flint_printf("\n\n");
+                abort();
+            }
+        }
+
+        _acb_vec_clear(v, len);
+        _acb_vec_clear(w1, len);
+        _acb_vec_clear(w2, len);
+
+        acb_dirichlet_group_clear(G);
+    }
+
     flint_printf("PASS\n");
 }