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).