Based on "On computing factors of cyclotomic polynomials" by Richard P. Brent.
For square-free n>1 the cyclotomic polynomial Pn(x) satisfies the
identity of Aurifeuille, Le Lasseur and Lucas which states that
if n is odd then Pn((-1)(n-1)/2x)=Cn2-nxDn2, or
if n is even then, ±Pn/2(-x2)=Cn2-nxDn2.
The first few polynomials C(x) and D(x) can be found here.
Factors of the first few values of nn±1 can be found here.
The following code requires the Multi-Precision Integer Arithmetic package.
/* factors1(N, K, F, G) is true if the strings F and G are two of the */ /* factors, not necessarily prime, of N^((2K+1)N)-1 if N=1 mod 4, or of */ /* N^((2K+1)N)+1 otherwise, where N and K are integers, and N is */ /* square-free. */ /* e.g. factors1(37, 0, F, G) gives factors F and G of 37^37-1: */ /* F=`46959719470144429555105032871`, G=`6243610407478181159725577611` */ /* factors1(38, 0, F, G) gives factors F and G of 38^38+1: */ /* F=`73515350937824503550370503117`, G=`10128165505710094110937686497` */ factors1(N, K, F, G):- lucas(N, Cs, Ds), power([N], [K], M), multiply(M, M, MM), multiply_short(MM, N, X), evaluate(Cs, X, Csx), evaluate(Ds, X, Dsx), multiply(Dsx, M, T1), multiply_short(T1, N, T2), add(Csx, T2, F0), subtract(Csx, T2, G0), bigint_string(F0, F), bigint_string(G0, G). /* factors2(N, M, F, G) is true if the strings F and G are two of the */ /* factors, not necessarily prime, of (M^2 N)^N-1 if N=1 mod 4, or of */ /* (M^2 N)^N+1 otherwise, where N and M are integers, and N is square-free.*/ /* e.g. factors2(21, 22, F, G) gives factors F and G of 10164^21-1: */ /* F=`1153738090869575800135471`, G=`1053479076989729217964291` */ /* factors2(23, 3, F, G) gives factors F and G of 207^23+1: */ /* F=`41606227894419123631183627`, G=`21384360489997575328139509` */ factors2(N, M, F, G):- lucas(N, Cs, Ds), multiply_short([M], M, MM), multiply_short(MM, N, X), evaluate(Cs, X, Csx), evaluate(Ds, X, Dsx), multiply_short(Dsx, M, T1), multiply_short(T1, N, T2), add(Csx, T2, F0), subtract(Csx, T2, G0), bigint_string(F0, F), bigint_string(G0, G). % For square-free n>1 the cyclotomic polynomial CP_n(x) satisfies the % identity of Aurifeuille, Le Lasseur and Lucas: % For odd n, CP_n((-1)^((n-1)/2)x)=C_n^2-nxD_n^2 % For even n, +/-CP_(n/2)(-x^2)=C_n^2-nxD_n^2 % e.g. lucas(22, Cs, Ds) gives % Cs=[1,11,27,33,21,11,21,33,27,11,1], Ds=[1,4,7,6,3,3,6,7,4,1] lucas(N, [1|Cs], [1|Ds]):- N > 1, N mod 4 > 0, N mod 9 > 0, N mod 25 > 0, N mod 49 > 0, N mod 121 > 0, n_prime(N, N1), naive_totient(N1, TwoE), E is TwoE // 2, lucas_1(2, TwoE, N, N1, Qs), % Start at 2 to avoid jacobi(N, 1, Q) lucas_2(1, E, N, [1|Qs], [1], [1], Cs, Ds). % Calculates all the elements, except the first, of Qs, with an additional % zero at the end required for the last element of Ds, which is not used. % e.g. lucas_1(2, 20, 33, 33, Qs) gives % Qs=[1,0,1,-1,-2,-1,1,0,1,0,-2,-1,1,0,1,1,-2,-1,1,0] lucas_1(K, E, _, _, [0]):- K > E, !. lucas_1(K, E, N, N1, [Q|Qs]):- 1 =:= K mod 2, !, % K is odd jacobi(N, K, Q), % q_k=jacobi(n,k) K1 is K + 1, lucas_1(K1, E, N, N1, Qs). lucas_1(K, E, N, N1, [Q|Qs]):- % K is even int_gcd(N1, K, GK), NGK is N1 // GK, mobius(NGK, Mobius), % mobius(n'/gcd(n',k)) naive_totient(GK, Totient), % totient(gcd(n',k)) cosine(N, K, Cosine), % cos((n-1)k pi/4) Q is Mobius * Totient * Cosine, K1 is K + 1, lucas_1(K1, E, N, N1, Qs). % Calculates all the elements of Cs and Ds lucas_2(E, E, N, Qs, Cs0, Ds0, [C], []):-!, lucas_3(Cs0, Ds0, Qs, N, 0, C0, 0, _), C is C0 // (E + E). lucas_2(K, E, N, Qs, Cs0, Ds0, [C|Cs], [D|Ds]):- lucas_3(Cs0, Ds0, Qs, N, 0, C0, 0, D0), C is C0 // (K + K), D is (D0 + C) // (K + K + 1), K1 is K + 1, lucas_2(K1, E, N, Qs, [C|Cs0], [D|Ds0], Cs, Ds). % Calculates the k-th element of Cs and Ds lucas_3([], [], _, _, C, C, D, D). lucas_3([CC|Cs], [DD|Ds], [Q1,Q2,Q3|Qs], N, C0, C, D0, D):- C1 is C0 + N*Q1*DD - Q2*CC, % Add nq_(2k-2j-1)d_j - q_(2k-2j)c_j D1 is D0 + Q3*CC - Q2*DD, % Add q_(2k-2j+1)c_j - q_(2k-2j)d_j lucas_3(Cs, Ds, [Q3|Qs], N, C1, C, D1, D). % Returns n if n = 1 mod 4, 2n otherwise n_prime(N, N):-1 =:= N mod 4, !. n_prime(N, N1):-N1 is N + N. % Returns cos((n-1)*k*pi/4) for even k cosine(N, K, 1):-0 =:= ((N-1)*K // 2) mod 4, !. cosine(N, K, -1):-2 =:= ((N-1)*K // 2) mod 4, !. cosine(_, _, 0). /* evaluate(Polynomial, X, Y) is true if the bigint Y is the value of the */ /* integer-valued Polynomial evaluated for the bigint value X. */ evaluate([A|As], X, Y):- int_bigint(A, Ab), evaluate_1(As, X, Ab, Y). evaluate_1([], _, Y, Y). evaluate_1([A0|As], X, Y0, Y):- A0 > 0, !, int_bigint(A0, A), multiply(X, Y0, XY0), add(XY0, A, Y1), evaluate_1(As, X, Y1, Y). evaluate_1([A0|As], X, Y0, Y):- A1 is -A0, int_bigint(A1, A), multiply(X, Y0, XY0), subtract(XY0, A, Y1), evaluate_1(As, X, Y1, Y). /* naive_totient(N, Totient) is true if Totient is Euler's totient function, */ /* defined as the number of positive integers R (1 <= R < N) that are */ /* relatively prime to N. */ /* e.g. naive_totient(42, X) gives X=12. */ naive_totient(N, Totient):- naive_totient_1(0, N, 0, Totient). naive_totient_1(N, N, Totient, Totient):-!. naive_totient_1(R, N, Totient0, Totient):- int_gcd(R, N, 1), !, Totient1 is Totient0 + 1, R1 is R + 1, naive_totient_1(R1, N, Totient1, Totient). naive_totient_1(R, N, Totient0, Totient):- R1 is R + 1, naive_totient_1(R1, N, Totient0, Totient). /* int_gcd(X, Y, Z) is true if Z is the greatest common divisor of X and Y. */ int_gcd(X, 0, X):-!. int_gcd(X, Y, Z):-X1 is abs(X), Y1 is abs(Y), int_gcd_1(X1, Y1, Z). int_gcd_1(X, 0, Z):-!, Z=X. int_gcd_1(X, Y, Z):-R is X mod Y, int_gcd_1(Y, R, Z). /* jacobi(A, N, J) is true if J is the Jacobi symbol for A and N. */ /* e.g. jacobi(12079, 17, 1), jacobi(12079, 19, -1). */ /* Thus 12079 is a quadratic residue modulo 17, but a quadratic */ /* non-residue modulo 19. */ jacobi(A, N, J):- A >= 0, N > 1, jacobi_1(A, N, J). jacobi_1(A, _, J):- A =< 1, !, J = A. jacobi_1(A, N, J):- A mod 2 =:= 1, A mod 4 =:= 3, N mod 4 =:= 3, !, N_A is N mod A, jacobi(N_A, A, J0), J is -J0. jacobi_1(A, N, J):- A mod 2 =:= 1, !, N_A is N mod A, jacobi(N_A, A, J). jacobi_1(A, N, J):- N mod 8 =:= 1, !, A_2 is A // 2, jacobi(A_2, N, J). jacobi_1(A, N, J):- N mod 8 =:= 7, !, A_2 is A // 2, jacobi(A_2, N, J). jacobi_1(A, N, J):- A_2 is A // 2, jacobi(A_2, N, J0), J is -J0. /* 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