add discrete log precomputations

acb_dirichlet/dlog.h

@@ -0,0 +1,178 @@
+/*=============================================================================
+
+    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
+
+******************************************************************************/
+
+#ifndef DLOG_H
+#define DLOG_H
+
+#include "padic.h"
+#include "ulong_extras.h"
+
+enum
+{
+    DLOG_MODPE, DLOG_CRT, DLOG_POWER, DLOG_BSGS, DLOG_TABLE
+};
+
+typedef struct dlog_precomp_struct dlog_precomp_struct;
+
+/* log in (1+pZ/p^eZ): compute via p-adic log */
+typedef struct
+{
+    ulong p;
+    ulong e;
+    padic_ctx_t ctx;     /* padic context */
+    padic_t invlog;
+    padic_t x;
+    fmpz_t r;
+}
+dlog_1modpe_struct;
+
+typedef dlog_1modpe_struct dlog_1modpe_t[1];
+
+/* log in (Z/p^eZ)^* */
+typedef struct
+{
+    ulong p;
+    ulong e;
+    ulong pe;
+    /*dlog_precomp_t modp;*/
+    dlog_precomp_struct * modp;
+    dlog_1modpe_t modpe;
+}
+dlog_modpe_struct;
+
+typedef dlog_modpe_struct dlog_modpe_t[1];
+
+/* all logs precomputed in (Z/modZ)^ast */
+typedef struct
+{
+    ulong mod;
+    ulong * table;
+}
+dlog_table_struct;
+
+typedef dlog_table_struct dlog_table_t[1];
+
+/* bsgs table, already in flint */
+
+typedef struct apow {
+    ulong k;
+    ulong ak;
+} apow_t;
+
+typedef struct {
+    ulong n;
+    double ninv;
+    ulong m;
+    ulong am;
+    apow_t * table;
+} dlog_bsgs_struct;
+
+typedef dlog_bsgs_struct dlog_bsgs_t[1];
+/* typedef bsgs_t dlog_bsgs_t; */
+
+/* CRT decomposition (Pohlig-Hellman) */
+typedef struct
+{
+    ulong mod;
+    ulong n;
+    ulong num;
+    ulong * expo;
+    ulong * crt_coeffs;
+    dlog_precomp_struct ** pre;
+    /*
+       void * pre;
+       dlog_precomp_t * pre;
+       */
+}
+dlog_crt_struct;
+
+typedef dlog_crt_struct dlog_crt_t[1];
+
+/* dlog when generator has prime power order */
+typedef struct
+{
+    ulong mod;
+    ulong p;
+    ulong e;
+    ulong * apk;
+    dlog_precomp_struct * pre;
+    /*
+       void * pre;
+       dlog_precomp_t * pre;
+       */
+}
+dlog_power_struct;
+
+typedef dlog_power_struct dlog_power_t[1];
+
+/* generic decomposition */
+/*typedef */
+struct
+dlog_precomp_struct
+{
+    int type;
+    union
+    {
+        dlog_table_t table;
+        dlog_bsgs_t bsgs;
+        dlog_crt_t crt;
+        dlog_power_t power;
+        dlog_modpe_t modpe;
+    } t;
+};
+
+typedef dlog_precomp_struct dlog_precomp_t[1];
+
+void dlog_table_init(dlog_table_t t, ulong a, ulong mod);
+void dlog_1modpe_init(dlog_1modpe_t t, ulong a1, ulong p, ulong e);
+void dlog_crt_init(dlog_crt_t t, ulong a, ulong mod, ulong n, ulong num);
+void dlog_power_init(dlog_power_t t, ulong a, ulong mod, ulong p, ulong e, ulong num);
+void dlog_modpe_init(dlog_modpe_t t, ulong a, ulong p, ulong e, ulong pe, ulong num);
+void dlog_bsgs_init(dlog_bsgs_t t, ulong a, ulong n, ulong m);
+/*#define dlog_bsgs_init(t, a, n, m) bsgs_table_init(t, a, n, m)*/
+
+void dlog_table_clear(dlog_table_t t);
+void dlog_1modpe_clear(dlog_1modpe_t t);
+void dlog_crt_clear(dlog_crt_t t);
+void dlog_power_clear(dlog_power_t t);
+void dlog_modpe_clear(dlog_modpe_t t);
+void dlog_bsgs_clear(dlog_bsgs_t t);
+/*#define dlog_bsgs_clear(t) bsgs_table_clear(t)*/
+
+ulong dlog_table(const dlog_table_t t, ulong b);
+ulong dlog_crt(const dlog_crt_t t, ulong b);
+ulong dlog_power(const dlog_power_t t, ulong b);
+ulong dlog_modpe(const dlog_modpe_t t, ulong b);
+ulong dlog_1modpe(const dlog_1modpe_t t, ulong b);
+ulong dlog_bsgs(const dlog_bsgs_t t, ulong b);
+/*#define dlog_bsgs(t, b) n_discrete_log_bsgs_table(t, b)*/
+
+void dlog_precomp_n_init(dlog_precomp_t pre, ulong a, ulong mod, ulong n, ulong num);
+void dlog_precomp_p_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong num);
+void dlog_precomp_pe_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong e, ulong pe, ulong num);
+void dlog_precomp_clear(dlog_precomp_t pre);
+ulong dlog_precomp(const dlog_precomp_t pre, ulong b);
+
+#endif

acb_dirichlet/dlog_1modpe.c

@@ -0,0 +1,82 @@
+/*=============================================================================
+
+    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"
+#include "padic.h"
+
+typedef struct
+{
+    ulong p;
+    padic_ctx_t ctx;     /* padic context */
+    padic_t invlog;
+}
+dlog_1modpe_struct;
+
+typedef dlog_1modpe_struct dlog_1modpe_t[1];
+
+void
+dlog_1modpe_init(dlog_1modpe_t t, ulong a1, ulong p, ulong e)
+{
+    fmpz_t tmp;
+    t->p = p;
+    fmpz_init(tmp);
+    padic_init(t->invlog);
+
+    fmpz_set_ui(tmp, p);
+    padic_ctx_init(t->ctx , tmp , 0 , e, PADIC_TERSE);
+
+    padic_set_ui(t->invlog, a1, t->ctx);
+    padic_inv(t->invlog, t->invlog, t->ctx);
+
+    fmpz_clear(tmp);
+}
+
+void
+dlog_1modpe_clear(dlog_1modpe_t t)
+{
+    padic_clear(t->invlog);
+    padic_ctx_clear(t->ctx);
+}
+
+/* assume b = 1 mod p, not checked */
+ulong
+dlog_1modpe(const dlog_1modpe_t t, ulong b)
+{
+    padic_t px;
+    fmpz_t ix;
+    ulong ux;
+    padic_init(px);
+    fmpz_init(ix);
+
+    padic_set_ui(px, b, t->ctx);
+    padic_log(px, px, t->ctx);
+    padic_mul(px, px, t->invlog, t->ctx);
+    padic_get_fmpz(ix, px, t->ctx);
+    ux =  fmpz_get_ui(ix);
+
+    padic_clear(px);
+    fmpz_clear(ix);
+    return ux;
+}

acb_dirichlet/dlog_bsgs.c

@@ -0,0 +1,84 @@
+/*=============================================================================
+
+    This file is part of FLINT.
+
+    FLINT 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.
+
+    FLINT 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 FLINT; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2013 Mike Hansen
+    Copyright (C) 2016 Pascal Molin
+
+******************************************************************************/
+
+#include <stdlib.h>
+
+#include "dlog.h"
+#include "ulong_extras.h"
+
+static int
+apow_cmp(const apow_t * x, const apow_t * y)
+{
+    return (x->ak < y->ak) ? -1 : (x->ak > y->ak);
+}
+
+/* set size of table m=sqrt(nk) to compute k logs in a group of size n */
+void
+dlog_bsgs_init(dlog_bsgs_t t, ulong a, ulong n, ulong m)
+{
+    ulong k, ak;
+    t->table = (apow_t *)flint_malloc(m * sizeof(apow_t));
+
+    t->n = n;
+    t->ninv = n_precompute_inverse(n);
+    t->m = m;
+
+    for (k = 0, ak = 1; k < m; k++)
+    {
+        t->table[k].k = k;
+        t->table[k].ak = ak;
+        ak = n_mulmod_precomp(ak, a, n, t->ninv);
+    }
+
+    t->am = n_invmod(ak, n);
+    qsort(t->table, m, sizeof(apow_t), (int(*)(const void*,const void*))apow_cmp);
+}
+
+void
+dlog_bsgs_clear(dlog_bsgs_t t)
+{
+    flint_free(t->table);
+}
+
+ulong
+dlog_bsgs(const dlog_bsgs_t t, ulong b)
+{
+    ulong i;
+    apow_t c, * x;
+
+    c.k = 0;
+    c.ak = b;
+    for (i = 0; i < t->m; i++)
+    {
+        x = bsearch(&c, t->table, t->m, sizeof(apow_t),
+            (int(*)(const void*,const void*))apow_cmp);
+        if (x != NULL)
+            return i * t->m + x->k;
+        c.ak = n_mulmod_precomp(c.ak, t->am, t->n, t->ninv);
+    }
+    flint_printf("Exception (n_discrete_log_bsgs).  discrete log not found.\n");
+    abort();
+}

acb_dirichlet/dlog_crt.c

@@ -0,0 +1,92 @@
+/*=============================================================================
+
+    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 "dlog.h"
+
+void
+dlog_crt_init(dlog_crt_t t, ulong a, ulong mod, ulong n, ulong num)
+{
+    int k;
+    n_factor_t fac;
+    ulong * M, * u;
+
+    n_factor_init(&fac);
+    n_factor(&fac, n, 1);
+
+    t->num = fac.num;
+    t->mod = mod;
+    t->n = n;
+
+    M = t->expo = flint_malloc(t->num * sizeof(ulong));
+    u = t->crt_coeffs = flint_malloc(t->num * sizeof(ulong));
+    t->pre = flint_malloc(t->num * sizeof(dlog_precomp_t));
+    for (k = 0; k < t->num; k++)
+        t->pre[k] = flint_malloc(sizeof(dlog_precomp_struct));
+
+    for (k = 0; k < t->num; k++)
+    {
+        ulong p, e, mk;
+        p = fac.p[k];
+        e = fac.exp[k];
+        mk = n_pow(p, e);
+        M[k] = n / mk;
+        u[k] = M[k] * n_invmod(M[k] % mk, mk) % n;
+        /* depends on the power */
+        dlog_precomp_pe_init(t->pre[k], n_powmod(a, M[k], mod), mod, p, e, mk, num);
+    }
+}
+
+void
+dlog_crt_clear(dlog_crt_t t)
+{
+    int k;
+    flint_free(t->expo);
+    flint_free(t->crt_coeffs);
+    for (k = 0; k < t->num; k++)
+    {
+        dlog_precomp_clear(t->pre[k]);
+        flint_free(t->pre[k]);
+    }
+    flint_free(t->pre);
+}
+
+ulong
+dlog_crt(const dlog_crt_t t, ulong b)
+{
+    int k;
+    ulong r = 0;
+    for (k = 0; k < t->num; k++)
+    {
+        ulong bk, rk;
+        bk = n_powmod(b, t->expo[k], t->mod);
+        rk = dlog_precomp(t->pre[k], bk);
+#if 0
+        flint_printf("##[crt-%d]: log(%wu)=log(%wu^%wu) = %wu [size %wu mod %wu]\n",
+                k, bk, b, t->expo[k], rk, t->n/t->expo[k], t->mod);
+#endif
+        r = (r + t->crt_coeffs[k] * rk) % t->n;
+    }
+    return r;
+}

acb_dirichlet/dlog_modpe.c

@@ -0,0 +1,55 @@
+/*=============================================================================
+
+    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 "dlog.h"
+#include "padic.h"
+
+void
+dlog_modpe_init(dlog_modpe_t t, ulong a, ulong p, ulong e, ulong pe, ulong num)
+{
+    t->p = p;
+    t->e = e;
+    t->pe = pe;
+    dlog_precomp_n_init(t->modp, a, p, p-1, num);
+    dlog_1modpe_init(t->modpe, n_powmod(a, p - 1, pe), p, e);
+}
+
+void
+dlog_modpe_clear(dlog_modpe_t t)
+{
+  dlog_precomp_clear(t->modp);
+  dlog_1modpe_clear(t->modpe);
+}
+
+ulong
+dlog_modpe(const dlog_modpe_t t, ulong b)
+{
+    ulong x;
+    x = dlog_precomp(t->modp, b % t->p);
+    /*b = b * n_powmod(t->a, -x, t->pe);*/
+    b = n_powmod(b, t->p - 1, t->pe);
+    x = x + (t->p-1) * dlog_1modpe(t->modpe, b);
+    return x;
+}

acb_dirichlet/dlog_power.c

@@ -0,0 +1,81 @@
+/*=============================================================================
+
+    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 "dlog.h"
+
+void
+dlog_power_init(dlog_power_t t, ulong a, ulong mod, ulong p, ulong e, ulong num)
+{
+    int k;
+    t->mod = mod;
+    t->p = p;
+    t->e = e;
+
+    t->apk = flint_malloc(e * sizeof(ulong));
+    t->pre = flint_malloc(sizeof(dlog_precomp_struct));
+
+    t->apk[0] = n_invmod(a, mod);
+    for (k = 1; k < e; k++)
+        t->apk[k] = n_powmod(t->apk[k-1], p, mod);
+
+    dlog_precomp_p_init(t->pre, n_invmod(t->apk[e-1], mod), mod, p, num);
+}
+
+void
+dlog_power_clear(dlog_power_t t)
+{
+    flint_free(t->apk);
+    dlog_precomp_clear(t->pre);
+    flint_free(t->pre);
+}
+
+/* 3^30*2+1, 2^30*3+1 are primes */
+
+ulong
+dlog_power(const dlog_power_t t, ulong b)
+{
+    int k;
+#ifdef C99
+    ulong pk[t->e];
+#else
+    ulong pk[30];
+#endif
+    ulong x;
+
+    pk[0] = 1;
+    for (k = 1; k < t->e; k++)
+       pk[k] = pk[k-1]*t->p;
+    
+    x = 0;
+    for(k = 0; k < t->e; k++)
+    {
+      ulong bk, xk;
+      bk = n_powmod(b, pk[t->e-1-k], t->mod);
+      xk = dlog_precomp(t->pre, bk);
+      b = b * n_powmod(t->apk[k], xk, t->mod) % t->mod;
+      x += xk * pk[k]; /* cannot overflow */
+    }
+    return x;
+}

acb_dirichlet/dlog_precomp.c

@@ -0,0 +1,174 @@
+/*=============================================================================
+
+    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 "math.h"
+#include "dlog.h"
+
+#define TABLE_P_LIM 50
+#define TABLE_PE_LIM 50
+#define TABLE_N_LIM 50
+#define BSGS_LIM  500 
+
+
+/* log mod p^e */
+void
+dlog_precomp_modpe_init(dlog_precomp_t pre, ulong a, ulong p, ulong e, ulong pe, ulong num)
+{
+    if ( pe < TABLE_PE_LIM )
+    {
+        pre->type = DLOG_TABLE;
+        dlog_table_init(pre->t.table, a, pe);
+        return;
+    }
+    if (e > 1)
+    {
+        pre->type = DLOG_MODPE;
+        dlog_modpe_init(pre->t.modpe, a, p, e, pe, num);
+    }
+    else
+    {
+        dlog_precomp_n_init(pre, a, p, p - 1, num);
+    }
+}
+
+/* group of order n modulo mod, mod a prime and no information on n */
+void
+dlog_precomp_n_init(dlog_precomp_t pre, ulong a, ulong mod, ulong n, ulong num)
+{
+    if (n%2 && n_is_probabprime(n))
+        dlog_precomp_p_init(pre, a, mod, n, num);
+    else {
+        if ( mod < TABLE_N_LIM )
+        {
+           pre->type = DLOG_TABLE;
+           dlog_table_init(pre->t.table, a, mod);
+        } else {
+            if (n < BSGS_LIM)
+            {
+                ulong m;
+                m = (2 * num < n) ? ceil(sqrt((double) n * num)) : n;
+                pre->type = DLOG_BSGS;
+                dlog_bsgs_init(pre->t.bsgs, a, mod, m);
+            } else {
+                pre->type = DLOG_CRT;
+                dlog_crt_init(pre->t.crt, a, mod, n, num);
+            }
+        }
+    }
+}
+
+/* we known the order is prime */
+void
+dlog_precomp_p_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong num)
+{
+    if ( p < TABLE_P_LIM )
+    {
+        pre->type = DLOG_TABLE;
+        dlog_table_init(pre->t.table, a, mod);
+    }
+    else
+    {
+        ulong m;
+        m = (2 * num < p) ? ceil(sqrt((double) p * num)) : p;
+        pre->type = DLOG_BSGS;
+        dlog_bsgs_init(pre->t.bsgs, a, mod, m);
+    }
+}
+
+void
+dlog_precomp_pe_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong e, ulong pe, ulong num)
+{
+    if ( pe < TABLE_PE_LIM )
+    {
+        pre->type = DLOG_TABLE;
+        dlog_table_init(pre->t.table, a, mod);
+    }
+    else
+    {
+        if ( e == 1)
+        {
+            dlog_precomp_p_init(pre, a, mod, p, num);
+        }
+        else
+        {
+            pre->type = DLOG_POWER;
+            dlog_power_init(pre->t.power, a, mod, p, e, num);
+        }
+    }
+}
+
+void
+dlog_precomp_clear(dlog_precomp_t pre)
+{
+    switch (pre->type)
+    {
+      case DLOG_MODPE:
+        dlog_modpe_clear(pre->t.modpe);
+        break;
+      case DLOG_CRT:
+        dlog_crt_clear(pre->t.crt);
+        break;
+      case DLOG_POWER:
+        dlog_power_clear(pre->t.power);
+        break;
+      case DLOG_TABLE:
+        dlog_table_clear(pre->t.table);
+        break;
+      case DLOG_BSGS:
+        dlog_bsgs_clear(pre->t.bsgs);
+        break;
+      default:
+        abort();
+        break;
+    }
+}
+
+ulong
+dlog_precomp(const dlog_precomp_t pre, ulong b)
+{
+    if (b == 1)
+        return 0;
+    switch (pre->type)
+    {
+      case DLOG_MODPE:
+        return dlog_modpe(pre->t.modpe, b);
+        break;
+      case DLOG_CRT:
+        return dlog_crt(pre->t.crt, b);
+        break;
+      case DLOG_POWER:
+        return dlog_power(pre->t.power, b);
+        break;
+      case DLOG_TABLE:
+        return dlog_table(pre->t.table, b);
+        break;
+      case DLOG_BSGS:
+        return dlog_bsgs(pre->t.bsgs, b);
+        break;
+      default:
+        abort();
+        break;
+    }
+}

acb_dirichlet/dlog_table.c

@@ -0,0 +1,55 @@
+/*=============================================================================
+
+    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 "dlog.h"
+
+void
+dlog_table_init(dlog_table_t t, ulong a, ulong mod)
+{
+  int k;
+  ulong ak;
+  t->mod = mod;
+  t->table = flint_malloc(mod * sizeof(ulong));
+  ak = 1; k = 0;
+  /* warning: do not check a is invertible modulo mod */
+  do
+  {
+    t->table[ak] = k++;
+    ak = (ak * a) % mod;
+  }
+  while (ak != 1);
+}
+
+void
+dlog_table_clear(dlog_table_t t)
+{
+  flint_free(t->table);
+}
+
+ulong
+dlog_table(const dlog_table_t t, ulong b)
+{
+    return t->table[b % t->mod];
+}

acb_dirichlet/test/t-dlog.c

@@ -0,0 +1,77 @@
+/*=============================================================================
+
+    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/dlog.h"
+#include <math.h>
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("dlog mod p....");
+    fflush(stdout);
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000; iter++)
+    {
+        dlog_table_t table;
+        dlog_bsgs_t bsgs;
+        dlog_precomp_t pre1, pre100;
+        ulong p, a, k;
+
+        p = n_randprime(state, 15, 0);
+        a = n_primitive_root_prime(p);
+
+        dlog_table_init(table, a, p);
+        dlog_bsgs_init(bsgs, a, p, ceil(sqrt((double)p)));
+        dlog_precomp_n_init(pre1, a, p, p-1, 1);
+        dlog_precomp_n_init(pre100, a, p, p-1, 100);
+
+        for (k = 1; k < 100 && k < p; k++)
+        {
+            ulong l1, l2, l3, l4;
+            l1 = dlog_table(table, k);
+            l2 = dlog_bsgs(bsgs, k);
+            l3 = dlog_precomp(pre1, k);
+            l4 = dlog_precomp(pre100, k);
+            if (l1 != l2 || l3 != l4 || l1 != l3)
+            {
+                flint_printf("FAIL: log(%wu,%wu) mod %wu: [%wu, %wu, %wu, %wu]\n",
+                        k, a, p, l1, l2, l3, l4);
+                abort();
+            }
+        }
+        dlog_table_clear(table);
+        dlog_bsgs_clear(bsgs);
+        dlog_precomp_clear(pre1);
+        dlog_precomp_clear(pre100);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}