How to compute cyclotomic numbers in GMP

The following program is for computing cyclotomic numbers in GMP written by A. Kruppa.

/* Compute a cyclotomic number, Phi_N(x). First command line parameter is */ 
/* N, second is x. */
/* Written in 2001 by A. Kruppa. Public Domain. */
/* Thanks to Yves Gallot, for his paper "Cyclotomic Polynomials and Prime
   Numbers" <http://perso.wanadoo.fr/yves.gallot/papers/cyclotomic.html>,
   and his advice */

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

/* Computes M"obius function */

int moeb(unsigned int n) {
  unsigned int i;
  int r=1;
  
  for (i=2; i<=n; i++) {
    if (n%i == 0) {
      n/=i;
      if (n%i == 0) return 0; /* Not squarefree */
      r=-r;
    }
  }
  return r;
}

/* Tests whether n is a perfect power of a prime. If so, returns that prime */
/* Returns 0 otherwise */

unsigned int isprimepower(unsigned int n) {
  unsigned int i;
  
  for (i=2; i*i <= n; i++) {
    if (n%i == 0) {
      while (n%i == 0) n/=i;
      if (n != 1) return 0;
      return i;
    }
  }
  /* n is prime */
  return n;
}


/* Computes Phi_N(x), where Phi_N is the Nth cyclotomic polynomial */

int main(int argc, char *argv[]) {
  int tN, x;
  unsigned int i, N;
  mpz_t Rp, Rq,t; /* R = Rp/Rq */

  if (argc < 3) {
    printf("%s <N> <x>\nComputes the value of the <N>th cyclotomic polynomial "
    "at point <x>.\nN must be a positive integer, x must be an integer\n", 
    argv[0]);
    return 0;
  }
  tN = atoi(argv[1]);
  x = atoi(argv[2]);

  if (tN < 1) {
    printf("N must be a positive integer\n");
    exit(EXIT_FAILURE);
  }
  N = (unsigned int)tN;

  if (x == -1) {
    if (N & 1) N *= 2;
    else N /= 2;
    x = 1;
  }

  if ((x == 1) && (N > 1)) {
    if ((i = isprimepower(N))) 
      printf("%d\n", i);
    else
      printf("1\n");
    return 0;
  }

  mpz_init_set_ui(Rp, 1);
  mpz_init_set_ui(Rq, 1);
  mpz_init(t);

  for (i=1; i<=N; i++) {
    if (N%i == 0) {
      mpz_set_si(t, x);
      mpz_pow_ui(t, t, i);
      mpz_sub_ui(t, t, 1);
      switch (moeb(N/i)) {
        case  1: {mpz_mul(Rp, Rp, t); break;}
        case -1: {mpz_mul(Rq, Rq, t); break;}
        case  0:
      }
    }
  }

  mpz_tdiv_qr(Rp, t, Rp, Rq);
  if (mpz_cmp_ui(t, 0) != 0) {
    printf("Program error, result is not an integer\n");
    mpz_clear(t);
    mpz_clear(Rq);
    mpz_clear(Rp);
    exit(EXIT_FAILURE);
  }
  mpz_out_str(stdout, 10, Rp);
  printf("\n");

  mpz_clear(t);
  mpz_clear(Rq);
  mpz_clear(Rp);

  return 0;
}

index