cleanup: Big pile of whitespace fixes, all at once.

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

#include <math.h>
#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 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, d;

  if (m%8 == 0)
    return (0);
  for (t = 1;; t++) {
    p = t*m + 1;
    if (!primep(p))
      continue;
    k = order(2, p);
    h = t*m/k;
    d = gcd(h, m);
    if (d == 1)
      return (t);
  }
  return (0);
}

#define PI 3.1415926535897932384626433832795028841971693993751

static mp *fieldpoly(unsigned m, int t, grand *rr)
{
  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: {
#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);
}

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 || t > 2)
    return (0);
  f = field_binpoly(p);
  q = fieldpoly(m, t, r);
  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);
}