/*
  Author:  Pate Williams (c) 1997

  Algorithm 1.4.4 (Primitive Root). See "A Course
  in Computational Algebraic Number Theory" by
  Henri Cohen page 25. A primitive root is a
  number, a, such that a ^ p mod p = a.
*/

#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include "lip.h"

#define BOUND 10000000l

struct factor {int expon; verylong prime;};

void trial_division(long B, verylong zN, struct factor *f,
                    long *count)
{
  int e, one = 0, pri = 0;
  long p;
  verylong za = 0, zb = 0;

  *count = 0;
  zcopy(zN, &za);
  zpstart2();
  do {
    p = zpnext();
    if (zsmod(za, p) == 0) {
      e = 0;
      do {
        zsdiv(za, p, &zb);
        zcopy(zb, &za);
        e++;
      } while (zsmod(za, p) == 0);
      f[*count].expon = e;
      zintoz(p, &f[*count].prime);
      *count = *count + 1;
      one = zscompare(za, 1l) == 0;
      if (!one) pri = zprobprime(za, 5l);
    }
  } while (!one && !pri && p <= B);
  if (pri) {
    f[*count].expon = 1;
    zcopy(za, &f[*count].prime);
    *count = *count + 1;
  }
  if (p > B) {
    fprintf(stderr, "fatal error\ncan't factor number");
    exit(1);
  }
  zfree(&za);
  zfree(&zb);
}

void primitive_root(verylong zp, verylong *za)
{
  long i, k;
  struct factor f[32];
  verylong zb = 0, ze = 0, zq = 0, zr = 0, zp1 = 0;

  for (i = 0; i < 32; i++) f[i].prime = 0;
  zsadd(zp, - 1l, &zp1);
  trial_division(BOUND, zp1, f, &k);
  zone(za);
  L2:
  i = 0;
  zsadd(*za, 1l, &zb);
  zcopy(zb, za);
  L3:
  zdiv(zp1, f[i].prime, &zq, &zr);
  zexpmod(*za, zq, zp, &ze);
  if (zscompare(ze, 1l) == 0) goto L2;
  if (++i < k) goto L3;
  zfree(&zb);
  zfree(&ze);
  zfree(&zq);
  zfree(&zr);
  zfree(&zp1);
  for (i = 0; i < 32; i++) zfree(&f[i].prime);
}

int main(void)
{
  char answer[256];
  int prime;
  verylong za = 0, zp = 0;

  do {
    do {
      printf("p = "); zread(&zp);
      prime = zprobprime(zp, 5l);
      if (!prime) printf("*error*\np must be prime\n");
    } while (!prime);
    primitive_root(zp, &za);
    printf("a = "); zwriteln(za);
    printf("another root (n or y)? ");
    scanf("%s", answer);
  } while (tolower(answer[0]) == 'y');
  zfree(&za);
  zfree(&zp);
  return 0;
}