+#include <stdio.h>
+#include <stdlib.h>
+
+#include <mLib/alloc.h>
+
+#include "field.h"
+#include "gf.h"
+#include "mptext.h"
+#include "fibrand.h"
+#include "mprand.h"
+#include "gfn.h"
+
+/*----- Polynomials over finite fields ------------------------------------*/
+
+typedef struct poly {
+ unsigned sz, len;
+ mp **v;
+} poly;
+
+#define POLY_INIT { 0, 0, 0 }
+
+#define POLY_ZEROP(p) (!(p)->len)
+#define POLY_CONSTANTP(p) ((p)->len == 1)
+#define POLY_DEGREE(p) ((p)->len - 1)
+
+void poly_free(poly *p)
+{
+ unsigned i;
+
+ for (i = 0; i < p->sz; i++)
+ MP_DROP(p->v[i]);
+ xfree(p->v); p->v = 0;
+ p->sz = p->len = 0;
+}
+
+void poly_ensure(field *f, poly *p, unsigned sz)
+{
+ mp **x;
+ unsigned i;
+ unsigned n;
+
+ if (p->sz >= sz)
+ return;
+ n = p->sz; if (!n) n = 2; do n <<= 1; while (n < sz);
+ x = xmalloc(n * sizeof(mp *));
+ for (i = 0; i < p->sz; i++)
+ x[i] = p->v[i];
+ for (; i < n; i++)
+ x[i] = MP_COPY(f->zero);
+ xfree(p->v);
+ p->v = x;
+ p->sz = n;
+}
+
+void poly_fix(field *f, poly *p, unsigned n)
+{
+ while (n && F_ZEROP(f, p->v[n - 1]))
+ n--;
+ p->len = n;
+}
+
+#define MAX(a, b) ((a) > (b) ? (a) : (b))
+
+static void putmphex(field *f, const char *name, mp *x)
+{
+ mp *z;
+
+ printf("%s = 0x", name);
+ z = F_OUT(f, MP_NEW, x);
+ mp_writefile(z, stdout, 16);
+ MP_DROP(z);
+ printf("\n");
+}
+
+void poly_dump(field *f, const char *name, poly *p)
+{
+ unsigned i;
+ mp *z = MP_NEW;
+
+ printf("%s = (", name);
+ for (i = p->len; i--; ) {
+ printf("0x");
+ z = F_OUT(f, z, p->v[i]);
+ mp_writefile(z, stdout, 16);
+ if (i) printf(", ");
+ }
+ mp_drop(z);
+ printf(")\n");
+}
+
+void poly_add(field *f, poly *d, poly *x, poly *y)
+{
+ unsigned i;
+ unsigned n = MAX(x->len, y->len);
+
+/* poly_dump(f, "add x", x); */
+/* poly_dump(f, "add y", y); */
+ poly_ensure(f, d, n);
+ for (i = 0; i < n; i++) {
+ mp *a = (i < x->len) ? x->v[i] : f->zero;
+ mp *b = (i < y->len) ? y->v[i] : f->zero;
+ d->v[i] = F_ADD(f, d->v[i], a, b);
+ }
+ poly_fix(f, d, n);
+/* poly_dump(f, "add d", d); */
+}
+
+void poly_sub(field *f, poly *d, poly *x, poly *y)
+{
+ unsigned i;
+ unsigned n = MAX(x->len, y->len);
+
+/* poly_dump(f, "sub x", x); */
+/* poly_dump(f, "sub y", y); */
+ poly_ensure(f, d, n);
+ for (i = 0; i < n; i++) {
+ mp *a = (i < x->len) ? x->v[i] : f->zero;
+ mp *b = (i < y->len) ? y->v[i] : f->zero;
+ d->v[i] = F_SUB(f, d->v[i], a, b);
+ }
+ poly_fix(f, d, n);
+/* poly_dump(f, "sub d", d); */
+}
+
+static void omuladd(field *f, poly *x, poly *y, mp *z, unsigned o)
+{
+ unsigned i;
+ unsigned n = MAX(x->len, y->len + o);
+ mp *t = MP_NEW;
+
+/* poly_dump(f, "omuladd x", x); */
+/* poly_dump(f, "omuladd y", y); */
+/* putmphex(f, "omuladd z", z); */
+/* printf("omuladd o = %u\n", o); */
+ poly_ensure(f, x, n);
+ for (i = o; i < n; i++) {
+ if (i < y->len + o) {
+ mp *a = (i < x->len) ? x->v[i] : f->zero;
+ t = F_MUL(f, t, y->v[i - o], z);
+ x->v[i] = F_ADD(f, x->v[i], a, t);
+ }
+ }
+ mp_drop(t);
+ poly_fix(f, x, n);
+/* poly_dump(f, "omuladd d", x); */
+}
+
+static void omulsub(field *f, poly *x, poly *y, mp *z, unsigned o)
+{
+ unsigned i;
+ unsigned n = MAX(x->len, y->len + o);
+ mp *t = MP_NEW;
+
+/* poly_dump(f, "omulsub x", x); */
+/* poly_dump(f, "omulsub y", y); */
+/* putmphex(f, "omulsub z", z); */
+/* printf("omulsub o = %u\n", o); */
+ poly_ensure(f, x, n);
+ for (i = o; i < n; i++) {
+ if (i < y->len + o) {
+ mp *a = (i < x->len) ? x->v[i] : f->zero;
+ t = F_MUL(f, t, y->v[i - o], z);
+ x->v[i] = F_SUB(f, x->v[i], a, t);
+ }
+ }
+ mp_drop(t);
+ poly_fix(f, x, n);
+/* poly_dump(f, "omulsub d", x); */
+}
+
+void poly_mul(field *f, poly *d, poly *x, poly *y)
+{
+ poly dd = POLY_INIT;
+ unsigned i;
+
+/* poly_dump(f, "mul x", x); */
+/* poly_dump(f, "mul y", y); */
+ poly_ensure(f, &dd, x->len + y->len);
+ for (i = 0; i < y->len; i++) {
+ if (!F_ZEROP(f, y->v[i]))
+ omuladd(f, &dd, x, y->v[i], i);
+ }
+ poly_free(d);
+ *d = dd;
+/* poly_dump(f, "mul d", d); */
+}
+
+void poly_copy(field *f, poly *d, poly *p)
+{
+ unsigned i;
+
+ poly_ensure(f, d, p->len);
+ for (i = 0; i < p->len; i++) {
+ MP_DROP(d->v[i]);
+ d->v[i] = MP_COPY(p->v[i]);
+ }
+ d->len = p->len;
+}
+
+void poly_div(field *f, poly *q, poly *r, poly *x, poly *y)
+{
+ poly qq = POLY_INIT, rr = POLY_INIT;
+ mp *i = MP_NEW;
+ mp *z = MP_NEW;
+ unsigned j;
+ unsigned o, oo = 0;
+
+/* poly_dump(f, "div x", x); */
+/* poly_dump(f, "div y", y); */
+ assert(y->len);
+ i = F_INV(f, MP_NEW, y->v[y->len - 1]);
+ poly_copy(f, &rr, x);
+ if (rr.len >= y->len) {
+ o = oo = rr.len - y->len + 1;
+ j = rr.len;
+ if (q) poly_ensure(f, &qq, o);
+ while (o) {
+ o--; j--;
+ if (!F_ZEROP(f, rr.v[j])) {
+ z = F_MUL(f, z, rr.v[j], i);
+ if (q) qq.v[o] = MP_COPY(z);
+ omulsub(f, &rr, y, z, o);
+ }
+ }
+ }
+ if (q) {
+ poly_fix(f, &qq, oo);
+ poly_free(q);
+ *q = qq;
+/* poly_dump(f, "div q", q); */
+ }
+/* poly_dump(f, "div r", &rr); */
+ if (!r)
+ poly_free(&rr);
+ else {
+ poly_free(r);
+ *r = rr;
+ }
+ mp_drop(i);
+ mp_drop(z);
+}
+
+void poly_monic(field *f, poly *d, poly *p)
+{
+ mp *z;
+ unsigned i, n;
+
+ assert(p->len);
+ n = p->len;
+/* poly_dump(f, "monic p", p); */
+ poly_ensure(f, d, n);
+ z = F_INV(f, MP_NEW, p->v[n - 1]);
+ for (i = 0; i < n; i++)
+ d->v[i] = F_MUL(f, d->v[i], p->v[i], z);
+ poly_fix(f, d, n);
+/* poly_dump(f, "monic d", d); */
+ mp_drop(z);
+}
+
+void poly_gcd(field *f, poly *g, poly *x, poly *y)
+{
+ poly uu = POLY_INIT, vv = POLY_INIT, rr = POLY_INIT;
+ poly *u = &uu, *v = &vv, *r = &rr;
+ poly *t;
+
+/* poly_dump(f, "gcd x", x); */
+/* poly_dump(f, "gcd y", y); */
+ poly_copy(f, u, x); poly_copy(f, v, y);
+ if (u->len < v->len) { t = u; u = v; v = t; }
+ while (!POLY_ZEROP(v)) {
+ poly_div(f, 0, r, u, v);
+ t = u; u = v; v = r; r = t;
+ }
+ poly_monic(f, g, u);
+ poly_free(u);
+ poly_free(v);
+ poly_free(r);
+/* poly_dump(f, "gcd g", g); */
+}
+
+mp *poly_solve(field *f, mp *d, mp *p, grand *r)
+{
+ poly g = POLY_INIT;
+ poly ut = POLY_INIT;
+ poly c = POLY_INIT;
+ poly h = POLY_INIT;
+ mp *z;
+ unsigned m = f->nbits;
+ unsigned i;
+
+ poly_ensure(f, &g, m + 1);
+ for (i = 0; i <= m; i++)
+ g.v[i] = mp_testbit(p, i) ? MP_COPY(f->one) : MP_COPY(f->zero);
+ poly_fix(f, &g, m + 1);
+ poly_ensure(f, &ut, 2);
+ while (POLY_DEGREE(&g) > 1) {
+ ut.v[1] = mprand(ut.v[1], m, r, 0);
+ poly_fix(f, &ut, 2);
+ poly_copy(f, &c, &ut);
+/* poly_dump(f, "c-in", &c); */
+
+ for (i = 1; i < m; i++) {
+ poly_mul(f, &c, &c, &c);
+ poly_add(f, &c, &c, &ut);
+ poly_div(f, 0, &c, &c, &g);
+/* putchar('.'); fflush(stdout); */
+ }
+/* poly_dump(f, "c-out", &c); */
+ poly_gcd(f, &h, &c, &g);
+/* poly_dump(f, "h", &h); */
+ if (POLY_CONSTANTP(&h) || POLY_DEGREE(&h) == POLY_DEGREE(&g))
+ continue;
+ if (2 * POLY_DEGREE(&h) > POLY_DEGREE(&g))
+ poly_div(f, &g, 0, &g, &h);
+ else
+ poly_copy(f, &g, &h);
+ }
+ z = MP_COPY(g.v[0]);
+ poly_free(&g);
+ poly_free(&ut);
+ poly_free(&c);
+ poly_free(&h);
+ mp_drop(d);
+ return (z);
+}
+
+/*----- Other stuff -------------------------------------------------------*/
+
+static int primep(unsigned x)
+{
+ unsigned i;
+
+ for (i = 2; i*i < x; i++) {
+ if (x%i == 0)
+ return (0);
+ }
+ return (1);
+}
+
+static unsigned gcd(unsigned u, unsigned v)
+{
+ unsigned r;
+
+ if (u < v) { r = u; u = v; v = r; }
+ for (;;) {
+ r = u%v;
+ if (!r) return (v);
+ u = v; v = r;
+ }
+}
+
+static int onbtype(unsigned m)
+{
+ unsigned t;
+ unsigned p, h;
+ unsigned k, x, d;
+
+ if (m%8 == 0)
+ return (0);
+ for (t = 1; t <= 2; t++) {
+ p = t*m + 1;
+ if (!primep(p))
+ continue;
+ for (x = 2, k = 1; x != 1; x = (2*x)%p, k++);
+ h = t*m/k;
+ d = gcd(h, m);
+ if (d == 1)
+ return (t);
+ }
+ return (0);
+}
+
+static mp *fieldpoly(unsigned m, int t)
+{
+ mp *p, *q, *r, *z;
+ unsigned i;
+
+ switch (t) {
+ case 1:
+ p = MP_ZERO;
+ for (i = 0; i <= m; i++)
+ p = mp_setbit(p, p, i);
+ break;
+ case 2:
+ q = MP_ONE;
+ p = MP_THREE;
+ r = MP_NEW;
+ for (i = 1; i < m; i++) {
+ r = mp_lsl(r, p, 1);
+ r = mp_xor(r, r, q);
+ z = r; r = q; q = p; p = z;
+ }
+ mp_drop(q);
+ mp_drop(r);
+ break;
+ default:
+ abort();
+ }
+ return (p);
+}
+
+static int dofip(unsigned m, mp **p, unsigned n, unsigned x)
+{
+ unsigned i;
+
+ for (i = n + 1; i < x; i++) {
+ *p = mp_setbit(*p, *p, i);
+ if (n ? dofip(m, p, n - 1, i) : gf_irreduciblep(*p))
+ return (1);
+ *p = mp_clearbit(*p, *p, i);
+ }
+ return (0);
+}
+
+static mp *fip(unsigned m)
+{
+ mp *p = MP_ONE;
+ unsigned n;
+
+ p = mp_setbit(p, p, m);
+ n = 0;
+ while (!dofip(m, &p, n, m))
+ n += 2;
+ return (p);
+}
+
+static mp *fnb(mp *p)
+{
+ unsigned t;
+ field *f;
+ grand *r = fibrand_create(0);
+ mp *q, *x;
+ unsigned m = mp_bits(p) - 1;
+
+ if ((t = onbtype(m)) == 0)
+ return (0);
+ f = field_binpoly(p);
+ q = fieldpoly(m, t);
+ x = poly_solve(f, MP_NEW, q, r);
+ MP_DROP(q);
+ F_DESTROY(f);
+ return (x);
+}
+
+int main(int argc, char *argv[])
+{
+ int m = atoi(argv[1]);
+ int i;
+ mp *p;
+ mp *beta;
+
+ if (!argv[2])
+ p = fip(m);
+ else {
+ p = MP_ONE;
+ p = mp_setbit(p, p, m);
+ for (i = 2; i < argc; i++)
+ p = mp_setbit(p, p, atoi(argv[i]));
+ }
+ if (!gf_irreduciblep(p)) {
+ printf("not irreducible\n");
+ return (1);
+ }
+
+ MP_PRINTX("p", p);
+ for (i = m + 1; i--; ) {
+ if (mp_testbit(p, i))
+ printf("%u ", i);
+ }
+ putchar('\n');
+
+ beta = fnb(p);
+ if (!beta)
+ printf("no onb\n");
+ else
+ MP_PRINTX("beta", beta);
+
+ mp_drop(p);
+ mp_drop(beta);
+ assert(mparena_count(MPARENA_GLOBAL) == 0);
+ return (0);
+}
projects / u / mdw / catacomb / commitdiff