Cyclotomic Polynomials

Calculates the coefficients of the N-th cyclotomic polynomial. These procedures require the Polynomial Operations package.
The first few cyclotomic polynomials can be found here.
A polynomial is represented by a list of its coefficients, high-order first. e.g. x3 + 3x2 - 2x + 4 is represented by [1,3,-2,4].

/* cyclotomic(N, CP) is true if the list CP contains the coefficients of the */
/*   N-th cyclotomic polynomial.                                             */
/* e.g. cyclotomic(105, CP) gives                                            */
/*   CP=[1,1,1,0,0,-1,-1,-2,-1,-1,0,0,1,1,1,1,1,1,0,0,-1,0,-1,0,-1,0,-1,0,   */
/*       -1,0,0,1,1,1,1,1,1,0,0,-1,-1,-2,-1,-1,0,0,1,1,1]                    */
/* This version is significantly faster and more compact than cyclo_poly/2.  */
cyclotomic(1, [1,-1]):-!.
cyclotomic(N, CP):-
  N > 1,
  x_n_1(N, Num),                % x^N-1
  cyclotomic_1(1, N, [1], Den),
  p_div(Num, Den, CP, [0]).

cyclotomic_1(D, N, Den, Den):-
  D + D > N, !.                 % No more divisors less than N
cyclotomic_1(D, N, Den0, Den):-
  0 =:= N mod D, !,             % D is a divisor of N
  cyclotomic(D, CPD),
  p_mult(Den0, CPD, Den1),
  D1 is D + 1,
  cyclotomic_1(D1, N, Den1, Den).
cyclotomic_1(D, N, Den0, Den):- % D is not a divisor of N
  D1 is D + 1,
  cyclotomic_1(D1, N, Den0, Den).
  
/* cyclo_poly(N, CP) is true if the list CP contains the coefficients of the */
/*   N-th cyclotomic polynomial.                                             */
/* e.g. cyclo_poly(105, CP) gives                                            */
/*   CP=[1,1,1,0,0,-1,-1,-2,-1,-1,0,0,1,1,1,1,1,1,0,0,-1,0,-1,0,-1,0,-1,0,   */
/*       -1,0,0,1,1,1,1,1,1,0,0,-1,-1,-2,-1,-1,0,0,1,1,1]                    */
/* This version is significantly slower than cyclotomic/2.                   */
cyclo_poly(1, [1,-1]):-!.
cyclo_poly(N, CP):-
  N > 1,
  cyclo_poly_1(1, N, [1], [1], CP).

cyclo_poly_1(D, N, Num, Den, CP):-
  D > N, !,         % No more divisors less than or equal to N
  p_div(Num, Den, CP, [0]).
cyclo_poly_1(D, N, Num, Den, CP):-
  0 =:= N mod D, !, % D is a divisor of N
  NdivD is N // D,
  mobius(NdivD, Mu),
  cyclo_poly_2(Mu, D, Num, Num1, Den, Den1),
  D1 is D + 1,
  cyclo_poly_1(D1, N, Num1, Den1, CP).
cyclo_poly_1(D, N, Num, Den, CP):-
  D1 is D + 1,      % D is not a divisor of N
  cyclo_poly_1(D1, N, Num, Den, CP).

cyclo_poly_2(1, D, Num, Num1, Den, Den):-!,
  x_n_1(D, XD1), % x^D-1
  p_mult(Num, XD1, Num1).
cyclo_poly_2(-1, D, Num, Num, Den, Den1):-!,
  x_n_1(D, XD1), % x^D-1
  p_mult(Den, XD1, Den1).
cyclo_poly_2(0, _, Num, Num, Den, Den).

/* mobius(N, Mu) is true if Mu is the Mobius function of N defined by:       */
/*   mobius(1)=1; mobius(p1p2...pr)=(-1)^r if p1,p2,...,pr are distinct      */
/*   primes; mobius(n)=0 if n is divisible by the square of a prime.         */
mobius(N, Mu):-
  N > 0,
  mobius_1(2, N, 1, Mu).

mobius_1(P, N, Mu, Mu):-
  P > N, !.               % N=1 or all prime factors of N are distinct
mobius_1(P, N, Mu0, Mu):-
  0 =\= N mod P, !,       % P is not a factor of N
  P1 is P + 1,
  mobius_1(P1, N, Mu0, Mu).
mobius_1(P, N, Mu0, Mu):- % P is a factor of N
  N1 is N // P,
  0 =\= N1 mod P, !,      % P^2 is not a factor of N
  P1 is P + 1,
  Mu1 is -Mu0,
  mobius_1(P1, N1, Mu1, Mu).
mobius_1(_, _, _, 0).     % P^2 is a factor of N

/* x_n_1(N, Poly) is true if Poly is the polynomial x^N-1.                   */
/* e.g. x_n_1(3, [1,0,0,-1]).                                                */
x_n_1(N, [1|Poly]):-x_n_1a(N, Poly).
 
x_n_1a(1, [-1]):-!.
x_n_1a(N, [0|Poly]):-
  N1 is N - 1,
  x_n_1a(N1, Poly).

/*
 * Applications of cyclotomic polynomials
 */
  
/* facts_m(N, Polys) is true if the list of polynomials in Polys constitute  */
/*   the irreducible factors of the polynomial x^N-1.                        */
/* e.g. facts_m(15, [[1,-1],[1,1,1],[1,1,1,1,1],[1,-1,0,1,-1,1,0,-1,1]]).    */
/*   i.e. x^15-1=(x-1)*(x^2+x+1)*(x^4+x^3+x^2+x+1)*(x^8-x^7+x^5-x^4+x^3-x+1) */
facts_m(N, Polys):-
  N > 0,
  facts_m_1(1, N, Polys).

facts_m_1(D, N, []):-
  D > N, !.              % No more divisors of N
facts_m_1(D, N, [Poly|Polys]):-
  0 =:= N mod D, !,      % D is a divisor of N
  cyclotomic(D, Poly),
  D1 is D + 1,
  facts_m_1(D1, N, Polys).
facts_m_1(D, N, Polys):- % D is not a divisor of N
  D1 is D + 1,
  facts_m_1(D1, N, Polys).

/* facts_p(N, Polys) is true if the list of polynomials in Polys constitute  */
/*   the irreducible factors of the polynomial x^N+1.                        */
/* e.g. facts_p(14, [[1,0,1],[1,0,-1,0,1,0,-1,0,1,0,-1,0,1]]).               */
/*   i.e. x^14+1=(x^2+1)*(x^12-x^10+x^8-x^6+x^4-x^2+1)                       */
facts_p(1, [[1,1]]):-!.
facts_p(N, Polys):-  
  N > 1,
  TwoN is 2 * N,
  facts_m(TwoN, Fs),         % Factors of x^(2*N)-1
  facts_m(N, Gs),            % Factors of x^N-1
  delete_all(Gs, Fs, Polys). % Factors of x^N+1
  
/* delete_all(Xs, Ys, Zs) is true if Zs is the result of removing the first  */
/*   occurrence of each element of the list Xs from the list Ys.             */
delete_all([], Ys, Ys).
delete_all([X|Xs], Ys, Zs):-efface(Ys, X, Ws), delete_all(Xs, Ws, Zs).

/* efface(Xs, Y, Zs) is true if Zs is the result of deleting the first       */
/*   occurrence of the element Y from the list Xs.                           */
efface([], _, []).
efface([Y|Xs], Y, Xs):-!.
efface([X|Xs], Y, [X|Zs]):-efface(Xs, Y, Zs).

LPA Index     Home Page

Valid HTML 4.0!