/*
  Author:  Pate Williams (c) 1997

  Algorithm 1.7.5 (Prime Power Test). See "A Course
  in Computational Algebraic Number Theory" by Henri
  Cohen page 42.
*/

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

int Miller_Rabin(long t, verylong zn, verylong *za)
{
  int value = 1;
  long i, j, s = 0;
  verylong zb = 0, zn1 = 0, zr = 0, zy = 0;

  if (zodd(zn)) {
    if (zscompare(zn, 3l) > 0) {
      zsadd(zn, - 1l, &zn1);
      zcopy(zn1, &zr);
      while (!zodd(zr)) {
        s++;
        zrshift(zr, 1l, &zb);
        zcopy(zb, &zr);
      }
      for (i = 1; value && i <= t; i++) {
        do zrandomb(zn1, za); while (zscompare(*za, 2l) < 0);
        zexpmod(*za, zr, zn, &zy);
        if (zscompare(zy, 1l) != 0 &&
            zcompare(zy, zn1) != 0) {
          j = 1;
          while (value && j <= s - 1 && zcompare(zy, zn1) != 0) {
            zmulmod(zy, zy, zn, &zb);
            zcopy(zb, &zy);
            if (zscompare(zy, 1l) == 0) value = 0;
            j++;
          }
          if (value && zcompare(zy, zn1) != 0) value = 0;
        }
      }
    }
    else value = 1;
  }
  else value = zscompare(zn, 2l) == 0;
  zfree(&zb);
  zfree(&zn1);
  zfree(&zr);
  zfree(&zy);
  return value;
}

int prime_power(verylong zn, verylong *zp)
{
  int mod, one, power;
  verylong za = 0, zb = 0, zc = 0, zd = 0, zq = 0;

  if (!zodd(zn)) {
    zintoz(2l, zp);
    goto L4;
  }
  zcopy(zn, &zq);
  L2:
  if (Miller_Rabin(5l, zq, &za)) {
    zcopy(zq, zp);
    goto L4;
  }
  zexpmod(za, zq, zq, &zb);
  zsubmod(zb, za, zq, &zc);
  zgcd(zc, zq, &zd);

  if (zscompare(zd, 1l) == 0 || zcompare(zd, zq) == 0) goto L5;
  zcopy(zd, &zq);
  goto L2;
  L4:
  if (!Miller_Rabin(5l, *zp, &za)) {
    zcopy(*zp, &zq);
    goto L2;
  }
  zcopy(zn, &za);
  do {
    zdiv(za, *zp, &zb, &zc);
    zcopy(zb, &za);
    zmod(za, *zp, &zb);
    mod = zscompare(zb, 0l) == 0;
    one = zscompare(za, 1l) == 0;
  } while (!one && mod);
  power = one;
  goto L6;
  L5:
  power = 0;
  L6:
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&zq);
  return power;
}

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

  do {
    printf("enter the number to be tested below:\n");
    zread(&zn);
    if (prime_power(zn, &zp)) {
      printf("number is a prime power of ");
      zwriteln(zp);
    }
    else printf("number is not a prime power\n");
    printf("another number (n or y)? ");
    scanf("%s", answer);
  } while (tolower(answer[0]) == 'y');
  zfree(&zn);
  zfree(&zp);
  return 0;
}