+#include <math.h>
#include <stdio.h>
#include <stdlib.h>
poly_mul(f, &c, &c, &c);
poly_add(f, &c, &c, &ut);
poly_div(f, 0, &c, &c, &g);
-/* putchar('.'); fflush(stdout); */
+/* putchar('.'); fflush(stdout); */
}
-/* poly_dump(f, "c-out", &c); */
+/* poly_dump(f, "c-out", &c); */
poly_gcd(f, &h, &c, &g);
-/* poly_dump(f, "h", &h); */
+/* poly_dump(f, "h", &h); */
if (POLY_CONSTANTP(&h) || POLY_DEGREE(&h) == POLY_DEGREE(&g))
continue;
if (2 * POLY_DEGREE(&h) > POLY_DEGREE(&g))
}
}
+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)
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;
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);
}
n = 0;
while (!dofip(m, &p, n, m))
n += 2;
- return (p);
+ return (p);
}
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);