#### Aurifeuillian Factorisation

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),
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),
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),
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
```