This code requires the Multi-Precision Integer Arithmetic package.
% p_1(Integer, Factor) is true if Factor is a factor of the Integer,
% where Integer and Factor are expressed as strings. If Factor = `1`,
% no factor has been found.
% e.g. p_1(`11111111111111111`, X) gives X=`20711723`
% e.g. power([2],[139],A),subtract(A,[1],B),bigint_string(B,C),p_1(C,D)
% factorizes the Mersenne number 2^139-1 given enough time
p_1(N0, G0):-
string_bigint(N0, N),
p_1_1(N, G0).
p_1_1(N, `3`):-
divide(N, [3], _, Mod),
Mod = [0], !.
p_1_1(N, G0):-
% Note: Ss contains a list of primes: each integer between 2 and 2000
% which is the power of a prime has been replaced by that prime.
Ss = [2,3,2,5,7,2,3,11,13,2,17,19,23,5,3,29,31,2,37,41,43,47,7,53,59,61,
2,67,71,73,79,3,83,89,97,101,103,107,109,113,11,5,127,2,131,137,
139,149,151,157,163,167,13,173,179,181,191,193,197,199,211,223,227,
229,233,239,241,3,251,2,257,263,269,271,277,281,283,17,293,307,311,
313,317,331,337,7,347,349,353,359,19,367,373,379,383,389,397,401,
409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,
503,509,2,521,523,23,541,547,557,563,569,571,577,587,593,599,601,
607,613,617,619,5,631,641,643,647,653,659,661,673,677,683,691,701,
709,719,727,3,733,739,743,751,757,761,769,773,787,797,809,811,821,
823,827,829,839,29,853,857,859,863,877,881,883,887,907,911,919,929,
937,941,947,953,31,967,971,977,983,991,997,1009,1013,1019,1021,2,
1031,1033,1039,1049,1051,1061,1063,1069,1087,1091,1093,1097,1103,
1109,1117,1123,1129,1151,1153,1163,1171,1181,1187,1193,1201,1213,
1217,1223,1229,1231,1237,1249,1259,1277,1279,1283,1289,1291,1297,
1301,1303,1307,1319,1321,1327,11,1361,1367,37,1373,1381,1399,1409,
1423,1427,1429,1433,1439,1447,1451,1453,1459,1471,1481,1483,1487,
1489,1493,1499,1511,1523,1531,1543,1549,1553,1559,1567,1571,1579,
1583,1597,1601,1607,1609,1613,1619,1621,1627,1637,1657,1663,1667,
1669,41,1693,1697,1699,1709,1721,1723,1733,1741,1747,1753,1759,
1777,1783,1787,1789,1801,1811,1823,1831,1847,43,1861,1867,1871,
1873,1877,1879,1889,1901,1907,1913,1931,1933,1949,1951,1973,1979,
1987,1993,1997,1999],
B = [3],
subtract(B, [1], Q),
p_1_2(Ss, N, B, Q, G),
bigint_string(G, G0).
p_1_2([], _, _, _, [1]):-!. % No factor found
p_1_2(_, _, [1], [0], [1]):-!. % No factor found
p_1_2(_, N, _, Q, G):-
% G = gcd(Q, N),
gcd(Q, N, G),
\+ equal([1], G), !. % Factor found
p_1_2([S|Ss], N, B, Q, G):-
% B1 = (B ^ S) mod N,
powermod(B, [S], N, B1),
% Q1 = Q * (B1 - 1) mod N,
subtract(B1, [1], B1_1),
multiply(Q, B1_1, Q_B1_1),
divide(Q_B1_1, N, _, Q1),
p_1_2(Ss, N, B1, Q1, G).