math/mpreduce.h: Missing include files.

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <mLib/alloc.h>
#include <mLib/darray.h>
#include <mLib/dstr.h>
#include <mLib/quis.h>
#include <mLib/report.h>

#include "gf.h"
#include "gfreduce.h"
#include "mp.h"

#include "factor.h"

static int primitivep(fact_v *f, mp *p)
{
  gfreduce r;
  unsigned i;
  mp *x = MP_NEW;
  int rc = 1;

  if (!gf_irreduciblep(p))
    return (0);
#ifdef DEBUG
  MP_PRINTX("p", p);
#endif
  gfreduce_create(&r, p);
  for (i = 0; i < DA_LEN(f); i++) {
    x = gfreduce_exp(&r, x, MP_TWO, DA(f)[i].r);
#ifdef DEBUG
    MP_PRINT("	r", DA(f)[i].r);
    MP_PRINTX("	 x", x);
#endif
    if (MP_EQ(x, MP_ONE)) {
      rc = 0;
      break;
    }
  }
  gfreduce_destroy(&r);
  MP_DROP(x);
  return (rc);
}

static int dofip(fact_v *f, 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(f, m, p, n - 1, i) : primitivep(f, *p))
      return (1);
    *p = mp_clearbit(*p, *p, i);
  }
  return (0);
}

static mp *fip(fact_v *f, unsigned m)
{
  mp *p = MP_ONE;
  unsigned n;

  p = mp_setbit(p, p, m);
  n = 0;
  while (!dofip(f, m, &p, n, m))
    n += 2;
  return (p);
}

static void findprim(unsigned m)
{
  fact_v f = DA_INIT;
  mp *q = mp_lsl(MP_NEW, MP_ONE, m);
  mp *p;
  unsigned i;

  q = mp_sub(q, q, MP_ONE);
  factor(q, &f);

#ifdef DEBUG
  {
    size_t i;
    for (i = 0; i < DA_LEN(&f); i++) {
      mp_writefile(DA(&f)[i].p, stdout, 10);
      printf("^%u = ", DA(&f)[i].e);
      mp_writefile(DA(&f)[i].n, stdout, 10);
      fputs(" (", stdout);
      mp_writefile(DA(&f)[i].r, stdout, 10);
      fputs(")\n", stdout);
    }
  }
#endif

  p = fip(&f, m);
  MP_PRINTX("p", p);
  for (i = m + 1; i--; ) {
    if (mp_testbit(p, i))
      printf("%u ", i);
  }
  putchar('\n');
  mp_drop(p);
  mp_drop(q);
  freefactors(&f);
}

int main(int argc, char *argv[])
{
  ego(argv[0]);
  if (argc != 2) {
    fprintf(stderr, "Usage: %s M\n", QUIS);
    exit(1);
  }
  findprim(atoi(argv[1]));
  return (0);
}