change dlog_mope to make recursion by hand

dlog.h

@@ -26,15 +26,14 @@
 #ifndef DLOG_H
 #define DLOG_H
 
-#ifdef ACB_INLINES_C
-#define ACB_INLINE
+#ifdef DLOG_INLINES_C
+#define DLOG_INLINE
 #else
-#define ACB_INLINE static __inline__
+#define DLOG_INLINE static __inline__
 #endif
 
-#include "ulong_extras.h"
-#include "nmod_vec.h"
-#include "padic.h"
+#include "flint/ulong_extras.h"
+#include "flint/nmod_vec.h"
 
 enum
 {
@@ -44,25 +43,17 @@ enum
 typedef struct dlog_precomp_struct dlog_precomp_struct;
 typedef struct dlog_precomp_struct * dlog_precomp_ptr;
 
-/* log in (1+pZ/p^eZ), e large: compute via p-adic log */
-typedef struct
-{
-    padic_ctx_t ctx;     /* padic context */
-    padic_t invlog;      /* 1 / log_p(a^(p-1)) */
-}
-dlog_1modpe_padic_struct;
-
-typedef dlog_1modpe_padic_struct dlog_1modpe_padic_t[1];
-
-/* log in (1+pZ/p^eZ), e small: use recursion formulas */
+/* log in (1+pZ/p^eZ), e small: use recursion formulas
+ * could use padic log instead but exponent is small
+ * for ulongs */
 typedef struct
 {
     ulong inv1p;         /* 1 / (1 + p) */
     ulong invloga1;      /* 1 / log(a^(p-1),1+p) */
 }
-dlog_1modpe_rec_struct;
+dlog_1modpe_struct;
 
-typedef dlog_1modpe_rec_struct dlog_1modpe_rec_t[1];
+typedef dlog_1modpe_struct dlog_1modpe_t[1];
 
 /* log in (Z/p^eZ)^* */
 typedef struct
@@ -73,10 +64,7 @@ typedef struct
     ulong inva;
     nmod_t pe;
     dlog_precomp_struct * modp;
-    union {
-      dlog_1modpe_rec_t rec;     /* if e <= 3 && e < p */
-      dlog_1modpe_padic_t padic; /* otherwise */
-    } modpe;
+    dlog_1modpe_t modpe;
 }
 dlog_modpe_struct;
 
@@ -114,10 +102,10 @@ typedef dlog_bsgs_struct dlog_bsgs_t[1];
 
 /* Pollard rho */
 typedef struct {
-  ulong a;
-  nmod_t n;
-  nmod_t mod;
-  int nisprime;
+    ulong a;
+    nmod_t n;
+    nmod_t mod;
+    int nisprime;
 } dlog_rho_struct;
 
 typedef dlog_rho_struct dlog_rho_t[1];
@@ -186,33 +174,25 @@ ulong dlog_crt_init(dlog_crt_t t, ulong a, ulong mod, ulong n, ulong num);
 ulong dlog_power_init(dlog_power_t t, ulong a, ulong mod, ulong p, ulong e, ulong num);
 ulong dlog_modpe_init(dlog_modpe_t t, ulong a, ulong p, ulong e, ulong pe, ulong num);
 ulong dlog_bsgs_init(dlog_bsgs_t t, ulong a, ulong mod, ulong n, ulong m);
-void dlog_1modpe_rec_init(dlog_1modpe_rec_t t, ulong a1, ulong p, ulong e, nmod_t pe);
-void dlog_1modpe_padic_init(dlog_1modpe_padic_t t, ulong a1, ulong p, ulong e);
+void dlog_1modpe_init(dlog_1modpe_t t, ulong a1, ulong p, ulong e, nmod_t pe);
 void dlog_rho_init(dlog_rho_t t, ulong a, ulong mod, ulong n);
 /*#define dlog_bsgs_init(t, a, n, m) bsgs_table_init(t, a, n, m)*/
 
-ACB_INLINE void
+DLOG_INLINE void
 dlog_order23_clear(dlog_order23_t t)
 {
     return;
 }
 
-ACB_INLINE void
+DLOG_INLINE void
 dlog_table_clear(dlog_table_t t)
 {
-  flint_free(t->table);
-}
-
-ACB_INLINE void
-dlog_1modpe_padic_clear(dlog_1modpe_padic_t t)
-{
-    padic_clear(t->invlog);
-    padic_ctx_clear(t->ctx);
+    flint_free(t->table);
 }
 
 void dlog_crt_clear(dlog_crt_t t);
 
-ACB_INLINE void
+DLOG_INLINE void
 dlog_power_clear(dlog_power_t t)
 {
     flint_free(t->apk);
@@ -220,25 +200,23 @@ dlog_power_clear(dlog_power_t t)
     flint_free(t->pre);
 }
 
-ACB_INLINE void
+DLOG_INLINE void
 dlog_modpe_clear(dlog_modpe_t t)
 {
-  dlog_precomp_clear(t->modp);
-  if (0 && t->e > 2)
-      dlog_1modpe_padic_clear(t->modpe.padic);
-  flint_free(t->modp);
+    dlog_precomp_clear(t->modp);
+    flint_free(t->modp);
 }
 
-ACB_INLINE void
+DLOG_INLINE void
 dlog_bsgs_clear(dlog_bsgs_t t)
 {
     flint_free(t->table);
 }
 
-ACB_INLINE void
+DLOG_INLINE void
 dlog_rho_clear(dlog_rho_t t)
 {
-  return;
+    return;
 }
 /*#define dlog_bsgs_clear(t) bsgs_table_clear(t)*/
 
@@ -247,11 +225,11 @@ 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_mod2e(const dlog_modpe_t t, ulong b);
 ulong dlog_bsgs(const dlog_bsgs_t t, ulong b);
 ulong dlog_rho(const dlog_rho_t t, ulong b);
 ulong dlog_1modpe_mod1p(ulong b1, ulong p, ulong e, ulong inv1p, nmod_t pe);
-ulong dlog_1modpe_padic(const dlog_1modpe_padic_t t, ulong b1);
-ulong dlog_1modpe_rec(const dlog_1modpe_rec_t t, ulong b1, ulong p, ulong e, nmod_t pe);
+ulong dlog_1modpe(const dlog_1modpe_t t, ulong b1, ulong p, ulong e, nmod_t pe);
 /*#define dlog_bsgs(t, b) n_discrete_log_bsgs_table(t, b)*/
 
 #define DLOG_SMALL_LIM 50

dlog/1modpe.c

@@ -26,7 +26,7 @@
 #include "dlog.h"
 
 ulong
-dlog_1modpe_rec(const dlog_1modpe_rec_t t, ulong b1, ulong p, ulong e, nmod_t pe)
+dlog_1modpe(const dlog_1modpe_t t, ulong b1, ulong p, ulong e, nmod_t pe)
 {
     if (e == 1)
         return 0;

dlog/1modpe_init.c

@@ -26,7 +26,7 @@
 #include "dlog.h"
 
 void
-dlog_1modpe_rec_init(dlog_1modpe_rec_t t, ulong a1, ulong p, ulong e, nmod_t pe)
+dlog_1modpe_init(dlog_1modpe_t t, ulong a1, ulong p, ulong e, nmod_t pe)
 {
     if (e == 1)
     {

dlog/1modpe_mod1p.c

@@ -25,6 +25,7 @@
 
 #include "dlog.h"
 
+/* for odd prime p, assume b1 = 1 mod p */
 ulong
 dlog_1modpe_mod1p(ulong b1, ulong p, ulong e, ulong inv1p, nmod_t pe)
 {
@@ -36,7 +37,11 @@ dlog_1modpe_mod1p(ulong b1, ulong p, ulong e, ulong inv1p, nmod_t pe)
     for (f = 1; f < e; f++)
     {      
         if (b1 % pf != 1)
+        {
+            flint_printf("ERROR dlog_1modpe_1modp: %wu %% %wu != 1 mod %wu\n\n",
+                    b1, pf, pe.n);
             abort();
+        }
         xf = (b1 - 1) / pf;
         xf = (xf % p) * pf1;
         x += xf;

dlog/1modpe_padic.c

@@ -1,63 +0,0 @@
-/*=============================================================================
-
-    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"
-
-/* assume b = 1 mod p, not checked */
-ulong
-dlog_1modpe_padic(const dlog_1modpe_padic_t t, ulong b1)
-{
-    padic_t px;
-    fmpz_t ix;
-    ulong ux;
-
-    if (b1 == 1)
-        return 0;
-
-    padic_init(px);
-    fmpz_init(ix);
-
-    padic_set_ui(px, b1, t->ctx);
-    flint_printf("set %wu -> ", b1);
-    padic_print(px, t->ctx);
-
-    padic_log(px, px, t->ctx);
-    flint_printf("\n\n compute log -> ");
-    padic_print(px, t->ctx);
-
-    flint_printf("\n\n 1/log(a^(p-1)) -> ");
-    padic_print(t->invlog, t->ctx);
-
-    padic_mul(px, px, t->invlog, t->ctx);
-    flint_printf("\n\n divide by log(a^(p-1)) -> ");
-    padic_print(px, t->ctx);
-
-    padic_get_fmpz(ix, px, t->ctx);
-    ux =  fmpz_get_ui(ix);
-    flint_printf("\n\nlog_p(%wu)/log_p(a) = %wu\n", b1, ux);
-    padic_clear(px);
-    fmpz_clear(ix);
-    return ux;
-}

dlog/mod2e.c

@@ -25,23 +25,33 @@
 
 #include "dlog.h"
 
-void
-dlog_1modpe_padic_init(dlog_1modpe_padic_t t, ulong a1, ulong p, ulong e)
+ulong
+dlog_mod2e(const dlog_modpe_t t, ulong b1)
 {
-    fmpz_t tmp;
-
-    fmpz_init(tmp);
-    padic_init(t->invlog);
-
-    fmpz_set_ui(tmp, p);
-    padic_ctx_init(t->ctx , tmp , 0 , e, PADIC_SERIES);
-
-    padic_set_ui(t->invlog, a1, t->ctx);
-    flint_printf("set %wu -> ", a1);
-    flint_printf("\n\n compute log -> ");
-    padic_log(t->invlog, t->invlog, t->ctx);
-    padic_print(t->invlog, t->ctx);
-    padic_inv(t->invlog, t->invlog, t->ctx);
-
-    fmpz_clear(tmp);
+    if (t->e == 2)
+        return (b1 % 4) == 3;
+    else
+    {
+        slong f;
+        ulong pf1, pf, x, xf;
+        pf1 = 1;
+        pf = 4;
+        x = 0;
+        for (f = 2; f < t->e; f++)
+        {
+            if (b1 % pf != 1)
+            {
+                flint_printf("ERROR dlog_mod2e: %wu %% %wu != 1 mod %wu\n\n",
+                        b1, pf, t->pe.n);
+                abort();
+            }
+            xf = (b1 - 1) / pf;
+            xf = (f == 2) ? xf % 4 : (xf % 2) * (pf1 / 2);
+            b1 = nmod_mul(b1, nmod_pow_ui(t->inva, xf, t->pe), t->pe);
+            x += xf;
+            pf1 = pf;
+            pf *= 2;
+        }
+        return x;
+    }
 }

dlog/modpe.c

@@ -29,23 +29,21 @@ ulong
 dlog_modpe(const dlog_modpe_t t, ulong b)
 {
     ulong x;
+    if (t->p == 2)
+        return dlog_mod2e(t, b);
     x = dlog_precomp(t->modp, b % t->p);
     if (t->e > 1)
     {
         ulong b1, y;
 #if 0
         b1 = nmod_mul(b, nmod_pow_ui(t->inva, x, t->pe), t->pe);
-#else
-        b1 = nmod_pow_ui(b, t->p - 1, t->pe);
-#endif
-        if (1 || t->e <= 2)
-            y = dlog_1modpe_rec(t->modpe.rec, b1, t->p, t->e, t->pe);
-        else
-            y = dlog_1modpe_padic(t->modpe.padic, b1);
+        y = dlog_1modpe(t->modpe.rec, b1, t->p, t->e, t->pe);
         y = y % t->pe1;
-#if 0
         x = x + (t->p - 1) * y;
 #else
+        b1 = nmod_pow_ui(b, t->p - 1, t->pe);
+        y = dlog_1modpe(t->modpe, b1, t->p, t->e, t->pe);
+        y = y % t->pe1;
         x = n_submod(x, y % (t->p - 1), t->p - 1);
         x = y + t->pe1 * x;
 #endif

dlog/modpe_init.c

@@ -32,18 +32,26 @@ dlog_modpe_init(dlog_modpe_t t, ulong a, ulong p, ulong e, ulong pe, ulong num)
 
     t->p = p;
     t->e = e;
-    t->pe1 = pe / p;
     nmod_init(&t->pe, pe);
     t->inva = nmod_inv(a, t->pe);
 
-    t->modp = flint_malloc(sizeof(dlog_precomp_struct));
-    dlog_precomp_n_init(t->modp, a, p, p - 1, num);
-
-    a1 = nmod_pow_ui(a, p - 1, t->pe);
-    if (1 || e <= 2)
-        dlog_1modpe_rec_init(t->modpe.rec, a1, p, e, t->pe);
+    if (p == 2)
+    {
+        t->modp = NULL;
+        t->pe1 = (e <= 2) ? 2 : pe / 4;
+        t->modpe->inv1p = t->inva;
+        t->modpe->invloga1 = 1;
+        return e - 2;
+    }
     else
-        dlog_1modpe_padic_init(t->modpe.padic, a1, p, e);
+    {
+        t->modp = flint_malloc(sizeof(dlog_precomp_struct));
+        t->pe1 = pe / p;
+        dlog_precomp_n_init(t->modp, a, p, p - 1, num);
 
-    return t->modp->cost + e;
+        a1 = nmod_pow_ui(a, p - 1, t->pe);
+        dlog_1modpe_init(t->modpe, a1, p, e, t->pe);
+
+        return t->modp->cost + e;
+    }
 }

dlog/precomp_clear.c

@@ -28,6 +28,8 @@
 void
 dlog_precomp_clear(dlog_precomp_t pre)
 {
+    if (pre == NULL)
+        return;
     switch (pre->type)
     {
         case DLOG_MODPE:

dlog/precomp_n_init.c

@@ -30,7 +30,7 @@
 void
 dlog_precomp_n_init(dlog_precomp_t pre, ulong a, ulong mod, ulong n, ulong num)
 {
-    if (n%2 && n_is_probabprime(n))
+    if (n % 2 && n_is_probabprime(n))
         dlog_precomp_p_init(pre, a, mod, n, num);
     else {
         if (n < DLOG_TABLE_N_LIM)