chiark - git - mdw - catacomb/blob - utils/prim.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4
   5 #include <mLib/alloc.h>
   6 #include <mLib/darray.h>
   7 #include <mLib/dstr.h>
   8 #include <mLib/quis.h>
   9 #include <mLib/report.h>
  10
  11 #include "gf.h"
  12 #include "gfreduce.h"
  13 #include "mp.h"
  14
  15 #include "factor.h"
  16
  17 static int primitivep(fact_v *f, mp *p)
  18 {
  19   gfreduce r;
  20   unsigned i;
  21   mp *x = MP_NEW;
  22   int rc = 1;
  23
  24   if (!gf_irreduciblep(p))
  25     return (0);
  26 #ifdef DEBUG
  27   MP_PRINTX("p", p);
  28 #endif
  29   gfreduce_create(&r, p);
  30   for (i = 0; i < DA_LEN(f); i++) {
  31     x = gfreduce_exp(&r, x, MP_TWO, DA(f)[i].r);
  32 #ifdef DEBUG
  33     MP_PRINT("  r", DA(f)[i].r);
  34     MP_PRINTX("  x", x);
  35 #endif
  36     if (MP_EQ(x, MP_ONE)) {
  37       rc = 0;
  38       break;
  39     }
  40   }
  41   gfreduce_destroy(&r);
  42   MP_DROP(x);
  43   return (rc);
  44 }
  45
  46 static int dofip(fact_v *f, unsigned m, mp **p, unsigned n, unsigned x)
  47 {
  48   unsigned i;
  49
  50   for (i = n + 1; i < x; i++) {
  51     *p = mp_setbit(*p, *p, i);
  52     if (n ? dofip(f, m, p, n - 1, i) : primitivep(f, *p))
  53       return (1);
  54     *p = mp_clearbit(*p, *p, i);
  55   }
  56   return (0);
  57 }
  58
  59 static mp *fip(fact_v *f, unsigned m)
  60 {
  61   mp *p = MP_ONE;
  62   unsigned n;
  63
  64   p = mp_setbit(p, p, m);
  65   n = 0;
  66   while (!dofip(f, m, &p, n, m))
  67     n += 2;
  68   return (p);
  69 }
  70
  71 static void findprim(unsigned m)
  72 {
  73   fact_v f = DA_INIT;
  74   mp *q = mp_lsl(MP_NEW, MP_ONE, m);
  75   mp *p;
  76   unsigned i;
  77
  78   q = mp_sub(q, q, MP_ONE);
  79   factor(q, &f);
  80
  81 #ifdef DEBUG
  82   {
  83     size_t i;
  84     for (i = 0; i < DA_LEN(&f); i++) {
  85       mp_writefile(DA(&f)[i].p, stdout, 10);
  86       printf("^%u = ", DA(&f)[i].e);
  87       mp_writefile(DA(&f)[i].n, stdout, 10);
  88       fputs(" (", stdout);
  89       mp_writefile(DA(&f)[i].r, stdout, 10);
  90       fputs(")\n", stdout);
  91     }
  92   }
  93 #endif
  94
  95   p = fip(&f, m);
  96   MP_PRINTX("p", p);
  97   for (i = m + 1; i--; ) {
  98     if (mp_testbit(p, i))
  99       printf("%u ", i);
 100   }
 101   putchar('\n');
 102   mp_drop(p);
 103   mp_drop(q);
 104   freefactors(&f);
 105 }
 106
 107 int main(int argc, char *argv[])
 108 {
 109   ego(argv[0]);
 110   if (argc != 2) {
 111     fprintf(stderr, "Usage: %s M\n", QUIS);
 112     exit(1);
 113   }
 114   findprim(atoi(argv[1]));
 115   return (0);
 116 }