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; }