ptab: Replace the Catacomb groups.

[u/mdw/catacomb] / utils / fnb.c

diff --git a/utils/fnb.c b/utils/fnb.c
index 8b7c1ab..98941c9 100644 (file)
--- a/utils/fnb.c
+++ b/utils/fnb.c
@@ -1,3 +1,4 @@
+#include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 
@@ -349,19 +350,28 @@ static unsigned gcd(unsigned u, unsigned v)
   }
 }
 
+static unsigned order(unsigned x, unsigned p)
+{
+  unsigned y, k;
+
+  if (!x || x == 1) return (0);
+  for (y = x, k = 1; y != 1; y = (y*x)%p, k++);
+  return (k);
+}
+
 static int onbtype(unsigned m)
 {
   unsigned t;
   unsigned p, h;
-  unsigned k, x, d;
+  unsigned k, d;
 
   if (m%8 == 0)
     return (0);
-  for (t = 1; t <= 2; t++) {
+  for (t = 1;; t++) {
     p = t*m + 1;
     if (!primep(p))
       continue;
-    for (x = 2, k = 1; x != 1; x = (2*x)%p, k++);
+    k = order(2, p);
     h = t*m/k;
     d = gcd(h, m);
     if (d == 1)
@@ -370,7 +380,9 @@ static int onbtype(unsigned m)
   return (0);
 }
 
-static mp *fieldpoly(unsigned m, int t)
+#define PI 3.1415926535897932384626433832795028841971693993751
+
+static mp *fieldpoly(unsigned m, int t, grand *rr)
 {
   mp *p, *q, *r, *z;
   unsigned i;
@@ -393,8 +405,48 @@ static mp *fieldpoly(unsigned m, int t)
       mp_drop(q);
       mp_drop(r);
       break;
-    default:
-      abort();
+    default: {
+#ifdef notdef
+      unsigned pp = t*m + 1;
+      unsigned uu;
+      unsigned j;
+      struct cplx { double r, i; };
+      struct cplx e, n;
+      struct cplx *f;
+
+      do uu = GR_RANGE(rr, pp); while (order(uu, pp) != t);
+      f = xmalloc(sizeof(struct cplx) * (m + 1));
+      for (i = 0; i <= m; i++) f[i].r = f[i].i = 0;
+      f[0].r = 1;
+      printf("poly init; type = %u\n", t);
+      for (i = m + 1; i--; )
+       printf("%3u: %g + %g i\n", i, f[i].r, f[i].i);
+      putchar('\n');
+      for (i = 1; i <= m; i++) {
+       e.r = e.i = 0;
+       for (j = 0; j < t; j++) {
+         double z = (pow(2, i) * pow(uu, j) * PI)/pp;
+         e.r -= cos(z); e.i -= sin(z);
+       }
+       printf("new factor: %g + %g i\n", e.r, e.i);
+       for (j = i; j--; ) {
+         f[j + 1].r += f[j].r;
+         f[j + 1].i += f[j].i;
+         n.r = f[j].r * e.r - f[j].i * e.i;
+         n.i = f[j].r * e.i + f[j].i * e.r;
+         f[j] = n;
+       }
+       printf("poly after %u iters\n", i);
+       for (j = m + 1; j--; )
+         printf("%3u: %g + %g i\n", j, f[j].r, f[j].i);
+       putchar('\n');
+      }
+      xfree(f);
+      p = MP_ZERO;
+#else
+    abort();
+#endif
+    } break;      
   }
   return (p);
 }
@@ -432,10 +484,10 @@ static mp *fnb(mp *p)
   mp *q, *x;
   unsigned m = mp_bits(p) - 1;
 
-  if ((t = onbtype(m)) == 0)
+  if ((t = onbtype(m)) == 0 || t > 2)
     return (0);
   f = field_binpoly(p);
-  q = fieldpoly(m, t);
+  q = fieldpoly(m, t, r);
   x = poly_solve(f, MP_NEW, q, r);
   MP_DROP(q);
   F_DESTROY(f);