GNFS_v1.0.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The function being used is: x^3 + x^2 + 3002023016*x + 3199785469 \n",
      "\n",
      "Rational Factor Base: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 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, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 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, 1361, 1367, 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, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997]\n",
      "Algebraic Factor Base: [(1, 7), (6, 7), (11, 17), (1, 19), (6, 19), (11, 19), (5, 23), (8, 23), (9, 23), (13, 29), (6, 31), (18, 31), (16, 37), (21, 47), (10, 59), (11, 59), (37, 59), (29, 61), (34, 61), (58, 61), (22, 67), (49, 71), (42, 73), (48, 79), (12, 97), (3, 109), (9, 113), (27, 113), (76, 113), (101, 127), (55, 131), (87, 149), (92, 149), (118, 149), (96, 157), (84, 163), (146, 167), (143, 173), (175, 179), (43, 191), (177, 197), (126, 211), (125, 229), (4, 233), (98, 241), (177, 241), (206, 241), (39, 251), (44, 251), (167, 251), (32, 263), (61, 263), (169, 263), (34, 269), (241, 271), (36, 277), (254, 283), (200, 313), (53, 317), (77, 317), (186, 317), (175, 331), (96, 337), (343, 347), (297, 349), (137, 373), (304, 373), (30, 379), (144, 379), (204, 379), (155, 389), (24, 397), (134, 397), (238, 397), (168, 401), (86, 409), (245, 419), (321, 421), (380, 431), (337, 439), (166, 443), (256, 457), (326, 457), (331, 457), (439, 461), (308, 467), (378, 487), (142, 491), (39, 499), (41, 509), (256, 521), (30, 523), (75, 541), (141, 541), (324, 541), (258, 547), (362, 571), (396, 577), (429, 587), (203, 601), (45, 607), (71, 619), (548, 619), (618, 619), (102, 631), (135, 641), (182, 643), (214, 643), (246, 643), (123, 647), (269, 653), (135, 659), (257, 659), (266, 659), (393, 673), (13, 677), (229, 683), (406, 691), (428, 691), (547, 691), (110, 701), (351, 709), (696, 719), (378, 727), (158, 733), (184, 739), (2, 743), (90, 751), (690, 751), (721, 751), (524, 757), (738, 761), (559, 769), (428, 787), (172, 797), (372, 809), (794, 811), (102, 823), (350, 823), (370, 823), (117, 829), (161, 829), (550, 829), (436, 839), (309, 853), (700, 857), (356, 859), (615, 859), (746, 859), (610, 863), (319, 877), (10, 883), (40, 883), (832, 883), (21, 887), (36, 911), (754, 919), (496, 929), (176, 937), (282, 937), (478, 937), (349, 947), (711, 953), (698, 967), (935, 977), (44, 983), (237, 983), (701, 983), (298, 997), (26, 1021), (144, 1021), (850, 1021), (218, 1031), (851, 1031), (992, 1031), (452, 1033), (396, 1039), (109, 1049), (309, 1061), (322, 1061), (429, 1061), (194, 1063), (161, 1093), (272, 1093), (659, 1093), (156, 1103), (192, 1103), (754, 1103), (290, 1109), (826, 1109), (1101, 1109), (600, 1117), (209, 1123), (488, 1129), (251, 1151), (1019, 1151), (1031, 1151), (381, 1153), (896, 1153), (1028, 1153), (1013, 1163), (846, 1171), (265, 1187), (975, 1187), (1133, 1187), (321, 1201), (272, 1213), (1203, 1223), (489, 1231), (777, 1237), (500, 1249), (1253, 1259), (267, 1279), (292, 1279), (719, 1279), (1063, 1283), (871, 1289), (857, 1291), (181, 1301), (399, 1301), (720, 1301), (1035, 1303), (365, 1307), (308, 1321), (6, 1327), (538, 1327), (782, 1327), (46, 1367), (268, 1373), (447, 1381), (1087, 1381), (1227, 1381), (887, 1423), (1275, 1427), (1111, 1429), (7, 1433), (97, 1433), (1328, 1433), (654, 1439), (427, 1451), (813, 1453), (1188, 1483), (1212, 1487), (243, 1493), (1322, 1543), (131, 1549), (1463, 1549), (1503, 1549), (1268, 1553), (355, 1559), (178, 1567), (76, 1579), (580, 1583), (488, 1597), (793, 1607), (706, 1609), (1262, 1619), (596, 1621), (1434, 1637), (62, 1657), (779, 1657), (815, 1657), (475, 1663), (561, 1663), (626, 1663), (482, 1669), (759, 1693), (167, 1697), (340, 1697), (1189, 1697), (906, 1699), (1080, 1699), (1411, 1699), (319, 1709), (958, 1721), (1284, 1723), (280, 1733), (821, 1741), (974, 1741), (1686, 1741), (711, 1747), (737, 1753), (520, 1759), (48, 1777), (1316, 1787), (316, 1789), (127, 1801), (1701, 1801), (1773, 1801), (994, 1811), (623, 1823), (1367, 1823), (1655, 1823), (539, 1831), (607, 1867), (1484, 1867), (1642, 1867), (1822, 1871), (1082, 1877), (1358, 1889), (1798, 1901), (1033, 1913), (211, 1933), (66, 1949), (758, 1949), (1124, 1949), (837, 1979), (2, 1997), (604, 1997), (1390, 1997), (301, 2003), (497, 2003), (1204, 2003), (202, 2011), (1077, 2017), (1780, 2027), (1938, 2039), (1199, 2063), (1858, 2069), (1311, 2081), (1369, 2081), (1481, 2081), (413, 2083), (1799, 2083), (1953, 2083), (639, 2089), (1764, 2089), (1774, 2089), (549, 2099), (760, 2099), (789, 2099), (1188, 2111), (1248, 2113), (750, 2131), (801, 2137), (1413, 2137), (2059, 2137), (214, 2141), (1268, 2143), (1438, 2143), (1579, 2143), (964, 2179), (789, 2213), (806, 2221), (1298, 2239), (1468, 2243), (699, 2251), (19, 2273), (409, 2273), (1844, 2273), (211, 2281), (1243, 2293), (1331, 2293), (2011, 2293), (100, 2297), (2203, 2333), (114, 2341), (496, 2347), (1203, 2371), (1006, 2377), (449, 2381), (2010, 2381), (2302, 2381), (993, 2383), (149, 2389), (2261, 2393), (323, 2399), (493, 2411), (945, 2423), (1639, 2423), (2261, 2423), (1012, 2437), (796, 2447), (1289, 2459), (2437, 2467), (1452, 2473), (1184, 2521), (1550, 2539), (704, 2549), (118, 2551), (965, 2551), (1467, 2551), (524, 2591), (1199, 2609), (1823, 2609), (2195, 2609), (1121, 2621), (1736, 2621), (2384, 2621), (679, 2633), (232, 2647), (1673, 2663), (2602, 2671), (1366, 2677), (2525, 2687), (2312, 2689), (856, 2693), (138, 2699), (235, 2699), (2325, 2699), (1975, 2707), (2342, 2713), (1438, 2719), (2266, 2731), (2428, 2741), (359, 2753), (882, 2753), (1511, 2753), (2154, 2767), (1803, 2777), (1568, 2789), (908, 2791), (2081, 2797), (685, 2801), (2616, 2833), (1264, 2837), (598, 2843), (535, 2857), (432, 2861), (459, 2861), (1969, 2861), (2129, 2879), (2138, 2887), (973, 2897), (1872, 2903), (2812, 2909), (2889, 2917), (1523, 2927), (435, 2939), (586, 2939), (1917, 2939), (4, 2957), (1240, 2957), (1712, 2957), (637, 2963), (2723, 2969), (2238, 2971), (1237, 2999), (2371, 3001), (462, 3037), (1016, 3041), (305, 3049), (2007, 3061), (1271, 3067), (1624, 3079), (729, 3083), (2618, 3083), (2818, 3083), (620, 3109), (1089, 3109), (1399, 3109), (2620, 3163), (1386, 3169), (880, 3181), (40, 3191), (59, 3221), (590, 3221), (2571, 3221), (1585, 3229), (2117, 3229), (2755, 3229), (189, 3253), (235, 3259), (734, 3259), (2289, 3259), (1133, 3299), (662, 3301), (1716, 3313), (2553, 3323), (1290, 3329), (2504, 3329), (2863, 3329), (1285, 3347), (191, 3359), (892, 3359), (2275, 3359), (2737, 3371), (964, 3389), (2984, 3391), (2950, 3407), (3322, 3457), (404, 3463), (3188, 3463), (3333, 3463), (325, 3467), (915, 3467), (2226, 3467), (108, 3469), (846, 3469), (2514, 3469), (3429, 3491), (1510, 3499), (2036, 3499), (3451, 3499), (1035, 3511), (1213, 3517), (2749, 3527), (1636, 3539), (2882, 3557), (108, 3581), (1421, 3593), (436, 3613), (2249, 3617), (2640, 3631), (26, 3637), (980, 3643), (1827, 3697), (819, 3701), (3444, 3709), (200, 3767), (1002, 3793), (1602, 3797), (1186, 3821), (2992, 3823), (2904, 3833), (2068, 3851), (2752, 3851), (2881, 3851), (3220, 3853), (3786, 3863), (6, 3877), (146, 3881), (1320, 3889), (1504, 3907), (1989, 3919), (2815, 3919), (3033, 3919), (3346, 3923), (164, 3929), (2256, 3943), (824, 3947), (1689, 3967), (2882, 3967), (3362, 3967), (3315, 3989), (2037, 4003), (2541, 4003), (3427, 4003), (2445, 4007), (2763, 4007), (2805, 4007), (3670, 4013), (2872, 4019), (3527, 4021), (166, 4049), (711, 4051), (2604, 4057), (1156, 4073), (1851, 4093), (2613, 4093), (3721, 4093), (3285, 4099), (3879, 4111), (2681, 4129), (2374, 4139), (2466, 4139), (3437, 4139), (978, 4153), (2511, 4177), (2693, 4177), (3149, 4177), (739, 4201), (2729, 4217), (1741, 4219), (3081, 4229), (146, 4231), (618, 4243), (1480, 4243), (2144, 4243), (2498, 4259), (1980, 4271), (2219, 4273), (3887, 4283), (345, 4327), (2916, 4339), (506, 4349), (495, 4357), (1221, 4357), (2640, 4357), (3461, 4363), (3440, 4373), (2656, 4391), (4297, 4409), (2890, 4421), (325, 4441), (694, 4447), (2396, 4457), (803, 4481), (1515, 4481), (2162, 4481), (201, 4493), (1313, 4507), (1633, 4517), (1850, 4519), (2173, 4547), (2508, 4547), (4412, 4547), (3763, 4561), (66, 4567), (955, 4567), (3545, 4567), (3182, 4583), (912, 4603), (2919, 4621), (2805, 4637), (1501, 4639), (3658, 4639), (4118, 4639), (1585, 4643), (272, 4649), (4398, 4649), (4627, 4649), (2195, 4651), (2557, 4651), (4549, 4651), (2395, 4657), (2786, 4663), (4374, 4673), (2049, 4679), (1075, 4691), (424, 4721), (2134, 4721), (2162, 4721), (1452, 4751), (1468, 4759), (1000, 4783), (1722, 4783), (2060, 4783), (323, 4813), (1760, 4813), (2729, 4813), (1879, 4817), (3143, 4817), (4611, 4817), (3090, 4861), (2231, 4877), (3596, 4889), (900, 4919), (1847, 4919), (2171, 4919), (3931, 4933), (2138, 4943), (900, 4951), (2663, 4987), (2812, 4987), (4498, 4987), (2774, 4999), (193, 5003), (1160, 5003), (3649, 5003), (1564, 5009), (2420, 5011), (2926, 5011), (4675, 5011), (222, 5023), (1805, 5023), (2995, 5023), (112, 5039), (461, 5059), (3957, 5077), (2914, 5081), (3887, 5099), (1955, 5101), (1404, 5107), (3893, 5107), (4916, 5107), (1117, 5113), (2735, 5119), (4386, 5153), (3987, 5167), (1091, 5171), (558, 5197), (4131, 5209), (1181, 5227), (3990, 5231), (3085, 5233), (3106, 5233), (4274, 5233), (1640, 5237), (2723, 5261), (3689, 5261), (4109, 5261), (4546, 5279), (2842, 5281), (1871, 5309), (3472, 5323), (2604, 5333), (3172, 5347), (3489, 5347), (4032, 5347), (3653, 5351), (4165, 5387), (1500, 5399), (2377, 5417), (4460, 5431), (2745, 5437), (3966, 5441), (5110, 5443), (2719, 5501), (3540, 5501), (4742, 5501), (2186, 5503), (170, 5519), (4124, 5521), (384, 5527), (2407, 5531), (515, 5557), (2128, 5557), (2913, 5557), (2943, 5563), (2697, 5569), (3649, 5569), (4791, 5569), (279, 5591), (2556, 5623), (212, 5639), (5520, 5639), (5545, 5639), (277, 5647), (2378, 5651), (1249, 5653), (3496, 5659), (1884, 5669), (461, 5683), (2355, 5683), (2866, 5683), (4603, 5693), (918, 5717), (4852, 5717), (5663, 5717), (3383, 5741), (3820, 5743), (1, 5749), (340, 5749), (5407, 5749), (5736, 5783), (3238, 5801), (4097, 5801), (4266, 5801), (844, 5813), (3028, 5827), (4546, 5839), (1129, 5851), (608, 5857), (2579, 5857), (2669, 5857), (4433, 5861), (2243, 5867), (2011, 5869), (3884, 5879), (3902, 5879), (3971, 5879), (4657, 5881), (3051, 5897), (4831, 5903), (4429, 5927), (2464, 6007), (4228, 6011), (524, 6043), (4838, 6089), (4940, 6091), (2146, 6101), (2589, 6113), (1000, 6121), (3350, 6131), (4770, 6133), (2622, 6151), (2681, 6173), (3421, 6199), (2, 6203), (2508, 6211), (3813, 6211), (6100, 6211), (1761, 6217), (4496, 6217), (6176, 6217), (4401, 6221), (5277, 6229), (5679, 6247), (585, 6263), (1828, 6263), (3849, 6263), (4195, 6271), (3633, 6329), (6028, 6337), (5454, 6343), (3377, 6353), (4451, 6353), (4877, 6353), (4019, 6359), (2668, 6367), (3609, 6421), (3376, 6427), (4127, 6427), (5350, 6427), (4067, 6449), (1155, 6451), (1878, 6451), (3417, 6451), (2224, 6473), (2505, 6481), (4780, 6481), (5676, 6481), (5376, 6491), (1380, 6521), (1501, 6521), (3639, 6521), (1494, 6529), (5847, 6547), (106, 6551), (2476, 6553), (5914, 6563), (2707, 6569), (5162, 6571), (5434, 6581), (6238, 6619), (5470, 6637), (5568, 6653), (3397, 6659), (551, 6661), (913, 6661), (5196, 6661), (1690, 6673), (5673, 6689), (1413, 6701), (282, 6709), (1621, 6709), (4805, 6709), (4569, 6719), (2423, 6733), (649, 6737), (1793, 6761), (5163, 6761), (6565, 6761), (1264, 6763), (2282, 6781), (5153, 6781), (6126, 6781), (4152, 6791), (4392, 6793), (3141, 6803), (2025, 6827), (2313, 6827), (2488, 6827), (5216, 6829), (867, 6841), (698, 6863), (505, 6907), (2256, 6907), (4145, 6907), (51, 6911), (528, 6911), (6331, 6911), (2143, 6959), (6561, 6961), (2911, 6967), (239, 6971), (332, 6971), (6399, 6971), (2281, 6983), (4713, 6983), (6971, 6983), (387, 6991), (1966, 6991), (4637, 6991), (3814, 6997), (5233, 7001), (3422, 7027), (6056, 7043), (131, 7057), (2361, 7057), (4564, 7057), (2198, 7069), (5826, 7069), (6113, 7069), (7056, 7079), (663, 7103), (918, 7103), (5521, 7103), (5063, 7129), (6294, 7151), (367, 7177), (6871, 7177), (7115, 7177), (2996, 7187), (5185, 7187), (6192, 7187), (7183, 7193), (3121, 7207), (1204, 7219), (6059, 7219), (7174, 7219), (6828, 7229), (997, 7237), (3995, 7243), (3105, 7247), (6651, 7283), (1049, 7297), (1968, 7297), (4279, 7297), (2119, 7307), (1133, 7309), (1542, 7309), (4633, 7309), (3875, 7321), (4795, 7349), (953, 7369), (6038, 7393), (4317, 7417), (5155, 7451), (5240, 7457), (516, 7487), (1003, 7487), (5967, 7487), (889, 7489), (1006, 7489), (5593, 7489), (815, 7507), (7195, 7529), (1619, 7541), (4581, 7561), (7196, 7589), (521, 7603), (2119, 7603), (4962, 7603), (6998, 7621), (5228, 7639), (1271, 7643), (3081, 7643), (3290, 7643), (1983, 7649), (847, 7669), (2331, 7669), (4490, 7669), (3456, 7673), (458, 7691), (473, 7691), (6759, 7691), (2565, 7699), (7349, 7703), (4723, 7717), (3560, 7723), (6035, 7727), (3634, 7741), (4704, 7741), (7143, 7741), (3811, 7759), (7426, 7789), (4304, 7817), (5087, 7817), (6242, 7817), (4034, 7823), (4843, 7823), (6768, 7823), (3898, 7853), (5407, 7853), (6400, 7853), (4225, 7867), (2328, 7877), (5819, 7877), (7606, 7877), (1393, 7883), (1650, 7907), (610, 7919), (5913, 7927), (478, 7933), (2333, 7933), (5121, 7933), (2005, 7937), (214, 7949), (2192, 7949), (5542, 7949), (6224, 7951), (5740, 7963), (2469, 8017), (4037, 8039), (4865, 8053), (4887, 8059), (4053, 8081), (4958, 8081), (7150, 8081), (6694, 8093), (4647, 8101), (5555, 8101), (5999, 8101), (1, 8111), (2610, 8111), (5499, 8111), (5756, 8123), (4735, 8147), (3308, 8161), (4701, 8167), (4623, 8171), (314, 8191), (1569, 8209), (7085, 8209), (7763, 8209), (6774, 8219), (3351, 8231), (3614, 8237), (593, 8243), (1740, 8243), (5909, 8243), (2207, 8263), (697, 8273), (3543, 8273), (4032, 8273), (1423, 8287), (215, 8291), (149, 8293), (1646, 8297), (3330, 8311), (4635, 8317), (4787, 8329), (7086, 8363), (5656, 8369), (6837, 8377), (2174, 8387), (508, 8389), (1344, 8389), (6536, 8389), (2584, 8419), (6200, 8423), (3767, 8429), (2595, 8443), (7996, 8447), (2000, 8467), (4736, 8501), (4401, 8521), (5001, 8521), (7639, 8521), (6608, 8527), (87, 8537), (450, 8539), (1898, 8539), (6190, 8539), (1146, 8563), (959, 8573), (1533, 8573), (6080, 8573), (4009, 8597), (5848, 8597), (7336, 8597), (7344, 8599), (6033, 8609), (4589, 8623), (1485, 8627), (7626, 8629), (3616, 8663), (6066, 8681), (8201, 8689), (8566, 8693), (3686, 8699), (3803, 8707), (4281, 8713), (4768, 8719), (3368, 8737), (8153, 8761), (2147, 8803), (5175, 8819), (4837, 8821), (7276, 8837), (5661, 8849), (343, 8861), (2094, 8861), (6423, 8861), (8323, 8863), (1575, 8867), (7428, 8867), (8730, 8867), (1608, 8887), (2865, 8893), (8730, 8923), (2036, 8929), (2983, 8929), (3909, 8929), (3703, 8951), (6853, 8951), (7345, 8951), (4584, 8971), (36, 8999), (5332, 9001), (5531, 9001), (7138, 9001), (3188, 9007), (632, 9011), (1831, 9011), (6547, 9011), (3920, 9029), (8649, 9041), (7610, 9043), (5397, 9049), (5784, 9103), (4507, 9109), (3022, 9133), (2056, 9157), (905, 9161), (2061, 9161), (6194, 9161), (6116, 9173)]\n",
      "Quadratic Character Base: [(6212, 9187), (758, 9199), (1316, 9203), (6870, 9209), (2928, 9221), (74, 9257), (3228, 9257), (5954, 9257), (1357, 9277), (1519, 9283), (2114, 9293)] \n",
      "\n",
      "Sieving to find more than 2012 smooth pairs.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "You're going to need to increase the sieve size to get about 1978 more rows.\n",
      "34 total smooth pairs found.\n",
      "1 total bases found.\n",
      "\n",
      "\n",
      "Square Root 1: 285072678265828550720 Square Root 2: 2449833115033385834889081654\n",
      "Factors are: 141597143119199 and 353223804178417\n",
      "Non-trivial factorizations: 1\n",
      "Elapsed seconds: 131.74826645851135\n",
      "Elapsed CPU time: 131.515625\n"
     ]
    }
   ],
   "source": [
    "import ast\n",
    "import copy\n",
    "import time\n",
    "\n",
    "n = 1415960099 #ChangeME\n",
    "#m = 31 #ChangeME\n",
    "\n",
    "f = None\n",
    "d = None\n",
    "depth = None\n",
    "\n",
    "a_lb = None\n",
    "a_ub = None\n",
    "b_lb = None\n",
    "b_ub = None\n",
    "\n",
    "configs = True #Set various variables below using config file config.txt in same directory\n",
    "find_factor_bases = True #Automatically retrieve factor bases / quadratic character base based on your upper bound size limit\n",
    "estimate_good_m = True #Get an m that is close to n^(1/d) (you  will need to decide your degree d beforehand)\n",
    "sieve_size_adjust = True #If not enough bases (solutions) found, incrementally increase sieve size to get more smooth pairs\n",
    "write_to_output_file = True #Write program output to file\n",
    "use_depthier_alg_factors = True #FIXME\n",
    "record_execution_time = True #Record how long it takes for the main chunk of the program to execute (printed at the end)\n",
    "verbose = False # -------------\n",
    "\n",
    "rat_factor_base = None\n",
    "alg_factor_base = None\n",
    "quad_character_base = None\n",
    "\n",
    "ub_rat_factor_base = None\n",
    "ub_alg_factor_base = None\n",
    "ub_quad_character_base = None\n",
    "\n",
    "if (configs == True):\n",
    "    with open(\"config.txt\") as file:\n",
    "        lines = file.readlines()\n",
    "        \n",
    "        for line in lines:\n",
    "            variable = line.split(\"=\")[0]\n",
    "            \n",
    "            if (variable[0] == \"#\"):\n",
    "                continue\n",
    "                \n",
    "            value = line.split(\"=\")[1]\n",
    "            \n",
    "            if (variable == (\"n\").lower()):\n",
    "                n = sage.rings.integer.Integer(value)\n",
    "            elif (variable == (\"m\").lower()):\n",
    "                m = sage.rings.integer.Integer(value)\n",
    "            elif (variable == (\"f\").lower()):\n",
    "                f = ast.literal_eval(value)\n",
    "            elif (variable == (\"d\").lower()):\n",
    "                d = sage.rings.integer.Integer(value)\n",
    "            elif (variable == (\"depth\").lower()):\n",
    "                depth = sage.rings.integer.Integer(value)\n",
    "            elif (variable == (\"r_factorbase\").lower()):\n",
    "                rat_factor_base = ast.literal_eval(value)\n",
    "            elif (variable == (\"a_factorbase\").lower()):\n",
    "                alg_factor_base = ast.literal_eval(value)\n",
    "            elif (variable == (\"q_characterbase\").lower()):\n",
    "                quad_character_base = ast.literal_eval(value)\n",
    "            elif (variable == (\"a_lb\")):\n",
    "                a_lb = sage.rings.integer.Integer(value)\n",
    "            elif (variable == (\"a_ub\")):\n",
    "                a_ub = sage.rings.integer.Integer(value)\n",
    "            elif (variable == (\"b_lb\")):\n",
    "                b_lb = sage.rings.integer.Integer(value)\n",
    "            elif (variable == (\"b_ub\")):\n",
    "                b_ub = sage.rings.integer.Integer(value)\n",
    "            elif (variable == (\"ub_rat_factor_base\")):\n",
    "                ub_rat_factor_base = sage.rings.integer.Integer(value)\n",
    "            elif (variable == (\"ub_alg_factor_base\")):\n",
    "                ub_alg_factor_base = sage.rings.integer.Integer(value)\n",
    "            elif (variable == (\"ub_quad_character_base\")):\n",
    "                ub_quad_character_base = sage.rings.integer.Integer(value)\n",
    "                \n",
    "        file.close()\n",
    "        \n",
    "\n",
    "def BaseMExpansion(n,m):\n",
    "#Pseudocode of this function at: https://sites.pitt.edu/~bonidie/cs441/Chapter4-2.pdf\n",
    "    result = []\n",
    "    q = n\n",
    "    while (q != 0):\n",
    "        a = q % m\n",
    "        q = q // m\n",
    "        result.append(a)\n",
    "    return result[::-1] \n",
    "\n",
    "def useFunct(f,x,n):\n",
    "    list2 = f[::-1]\n",
    "    result = 0\n",
    "    for i in range(len(list2)):\n",
    "        if (i == 0):\n",
    "            result += list2[0]\n",
    "        else:\n",
    "            result += list2[i] * power_mod(x,i,n)\n",
    "    return result\n",
    "\n",
    "def GetGoodM(n,d):\n",
    "    \n",
    "    result = int(n^(1/d))\n",
    "    \n",
    "    while ( (n // (result^d)) != 1):\n",
    "        result-=1\n",
    "        \n",
    "    return result\n",
    "\n",
    "def GetRationalFactorBase(lenRatFactorBase):\n",
    "    result = []\n",
    "    currentPrime = 2\n",
    "\n",
    "    while (len(result) < lenRatFactorBase):\n",
    "        \n",
    "        result.append(currentPrime)\n",
    "        \n",
    "        if (currentPrime == 2):\n",
    "            currentPrime = 3\n",
    "            \n",
    "        else:\n",
    "            currentPrime = next_prime(currentPrime + 1)\n",
    "        \n",
    "    return result\n",
    "\n",
    "def GetAlgebraicFactorBase(f,lenAlgFactorBase,d):\n",
    "    result = []\n",
    "    currentPrime = 2\n",
    "    while (len(result) < lenAlgFactorBase):\n",
    "        countedPrimes = 0\n",
    "        for r in range(currentPrime + 1):\n",
    "            if ( (useFunct(f,r,currentPrime) % currentPrime) == 0 and countedPrimes <= d):\n",
    "                result.append((r,currentPrime))\n",
    "                countedPrimes+=1\n",
    "                if (countedPrimes == d):\n",
    "                    break\n",
    "        if (currentPrime == 2):\n",
    "            currentPrime = 3\n",
    "        else:\n",
    "            currentPrime = next_prime(currentPrime + 1)\n",
    "    return result\n",
    "\n",
    "            \n",
    "def GetQuadraticCharacterBase(f,lenQuadCharBase,d,alg_factor_base):\n",
    "    usedPrimes = set()\n",
    "    result = []\n",
    "    \n",
    "    currentPrime = 2\n",
    "    \n",
    "    for i in range(len(alg_factor_base)):\n",
    "        usedPrimes.add(alg_factor_base[i][1])\n",
    "    \n",
    "    while (len(result) < lenQuadCharBase):\n",
    "        countedPrimes = 0\n",
    "        \n",
    "        for s in range(currentPrime + 1):\n",
    "            if (currentPrime in usedPrimes):\n",
    "                break\n",
    "                \n",
    "            if ( (useFunct(f,s,currentPrime) % currentPrime) == 0 and countedPrimes <= d):\n",
    "                result.append((s,currentPrime))\n",
    "                countedPrimes+=1\n",
    "                if (countedPrimes == d):\n",
    "                    break\n",
    "        if (currentPrime == 2):\n",
    "            currentPrime = 3\n",
    "        else:\n",
    "            currentPrime = next_prime(currentPrime + 1)\n",
    "    return result\n",
    "\n",
    "\n",
    "\n",
    "def compute_integer_product_of_pairs(polynomial_ring, f, m, N, integer_pairs, vector):\n",
    "        prod  = 1;\n",
    "\n",
    "        for j in range(len(vector)):\n",
    "            if (1 == vector[j]):\n",
    "                prod  = prod*(integer_pairs[j][0] + m*integer_pairs[j][1]) % N;\n",
    "\n",
    "        return prod;\n",
    "\n",
    "\n",
    "def find_square(f, beta_squared, polynomial_ring):\n",
    "    # computes the square root of beta_squared\n",
    "    # (a polynomial over ZZ of degree less than d, that represents the element in the order ZZ[theta])\n",
    "    # and returns the polynomial representing the square root.\n",
    "\n",
    "    K.<a> = NumberField(f);\n",
    "    beta2 = beta_squared(a);\n",
    "    if(not beta2.is_square()):\n",
    "        if (verbose):\n",
    "            print('The supposedly found square is not square!')\n",
    "        return False;\n",
    "    return polynomial_ring(beta2.sqrt().list());\n",
    "\n",
    "\n",
    "#Modifications to this function from https://web.archive.org/web/20220724010535/https://groups.google.com/g/sage-support/c/LTQjSUDoT8Q\n",
    "def compute_numberfield_product_of_pairs(polynomial_ring, f, integer_pairs, vector):\n",
    "\n",
    "    #I = polynomial_ring.ideal(f); #FIXED\n",
    "\n",
    "    product_polynomial = polynomial_ring([1]);\n",
    "\n",
    "    for j in range(len(vector)):\n",
    "\n",
    "        if (1 == vector[j]):\n",
    "\n",
    "            linear_poly = polynomial_ring([integer_pairs[j][0], integer_pairs[j][1]]);\n",
    "            product_polynomial = product_polynomial * linear_poly;\n",
    "            #product_polynomial.mod(I); #FIXED\n",
    "            product_polynomial = product_polynomial.mod(f) #FIXED\n",
    "\n",
    "    return product_polynomial;\n",
    "\n",
    "\n",
    "\n",
    "def compute_difference_of_squares(polynomial_ring, f, m, N, integer_pairs, vector):\n",
    "\n",
    "    found_squares = False;\n",
    "\n",
    "    u_plus_v = 0; \n",
    "    u_minus_v = 0;\n",
    "    u = None\n",
    "    v = None\n",
    "\n",
    "    ints_mod_N = ZZ.quo(ZZ.ideal(N));\n",
    "\n",
    "    vsquared = ints_mod_N(compute_integer_product_of_pairs(polynomial_ring, f, m, N, integer_pairs, vector));\n",
    "\n",
    "    if ( is_square(vsquared)):\n",
    "\n",
    "        beta_squared = compute_numberfield_product_of_pairs(polynomial_ring, f, integer_pairs, vector);\n",
    "        beta = find_square(f, beta_squared, polynomial_ring);\n",
    "\n",
    "        if (False != beta):\n",
    "\n",
    "            u = (beta(m)) % N;\n",
    "            v = vsquared.sqrt().lift();\n",
    "\n",
    "            found_squares = True;\n",
    "\n",
    "            return (found_squares,beta_squared,beta,u,v);\n",
    "\n",
    "        else:\n",
    "            if (verbose):\n",
    "                print('Failed to find a square root in number field.')\n",
    "    else:\n",
    "        if (verbose):\n",
    "            print('Integer was not square.')\n",
    "\n",
    "    return (found_squares,'Filler Word.');\n",
    "\n",
    "\n",
    "def runIt(n,m,f,d,a_lb,a_ub,b_lb,b_ub,depth,lengthRow,rat_factor_base,alg_factor_base,quad_character_base):\n",
    "    \n",
    "    r_mat = matrix(0,lengthRow)\n",
    "    tuples = []\n",
    "    \n",
    "    rational_remainingPrimes = set()\n",
    "    rational_remainingPrimes2 = set()\n",
    "    rational_remainingPrimes3 = set()\n",
    "    \n",
    "    def get_factors(numA,algebraic=False):\n",
    "        result = []\n",
    "        len_divisors = None\n",
    "        \n",
    "        if (algebraic):\n",
    "            len_divisors = len(alg_factor_base)\n",
    "        else:\n",
    "            len_divisors = len(rat_factor_base)\n",
    "            \n",
    "        for i in range(len_divisors):\n",
    "            if (numA == 1):\n",
    "                break\n",
    "                \n",
    "            if (algebraic):\n",
    "                thisFactor = alg_factor_base[i][1]\n",
    "            else:\n",
    "                thisFactor = rat_factor_base[i]\n",
    "                \n",
    "            exponent_tracker = 0\n",
    "            \n",
    "            while(thisFactor.divides(numA)):\n",
    "                    \n",
    "                exponent_tracker += 1\n",
    "                numA = numA // thisFactor\n",
    "                \n",
    "            if (exponent_tracker > 0):\n",
    "                result.append((thisFactor,exponent_tracker))\n",
    "                \n",
    "        \n",
    "        if (numA != 1 and is_prime(numA) and (not algebraic)):\n",
    "            if (numA in rational_remainingPrimes):\n",
    "                if (numA in rational_remainingPrimes2):\n",
    "                    rational_remainingPrimes3.add(numA)\n",
    "                else:\n",
    "                    rational_remainingPrimes2.add(numA)\n",
    "            else:\n",
    "                rational_remainingPrimes.add(numA)\n",
    "                \n",
    "        CompletelyFactored = False\n",
    "        \n",
    "        if (numA == 1):\n",
    "            CompletelyFactored = True\n",
    "        \n",
    "        return result,CompletelyFactored\n",
    "    \n",
    "    \n",
    "        \n",
    "    for a in range(a_lb,a_ub):\n",
    "        for b in range(b_lb,b_ub):\n",
    "                \n",
    "            r = (a + (b*m))\n",
    "            r_alg = (power_mod((-1*b),d,n) * (useFunct(f,(a*power_mod(b,-1,n)*-1),n))) % n\n",
    "            r_alg_2 = abs(r_alg-n)\n",
    "            \n",
    "            if (r == 0 or r_alg == 0):\n",
    "                continue\n",
    "                \n",
    "            depth_additions = []\n",
    "            \n",
    "            for i in range(depth):\n",
    "                depth_additions.append(r_alg_2 + (n*(i)))\n",
    "            \n",
    "            factor_r = get_factors(int(r))\n",
    "            r_factors = factor_r[0]\n",
    "            rat_fact_base_match = factor_r[1] #Completely factored over the rational factor base\n",
    "            if (not rat_fact_base_match): #If a + bm factors are not in the rational factor base\n",
    "                continue\n",
    "            \n",
    "            factor_r_alg = get_factors(int(r_alg),True)\n",
    "            r_alg_factors = factor_r_alg[0]\n",
    "            alg_fact_base_match = factor_r_alg[1] #If a + bθ factors completely over the algebraic factor base\n",
    "        \n",
    "            depth_alg_fact_base_match = False #All r_alg factors are in one depthier algebraic factor base\n",
    "            \n",
    "            if (not alg_fact_base_match):\n",
    "                for i in range(depth):\n",
    "                    get_r_alg_depthier_factors = get_factors(int(depth_additions[i])%n,True)\n",
    "                    r_alg_depth_factors = get_r_alg_depthier_factors[0]\n",
    "                    depth_alg_fact_base_match = get_r_alg_depthier_factors[1]\n",
    "                    \n",
    "                    #r_alg_factors_depthier.append(factor(int(depth_additions[i])%n)) #Modulo n, fixme?\n",
    "                    #r_alg_factors_depthier.append(r_alg_depth_factors)\n",
    "                    \n",
    "                    if (depth_alg_fact_base_match):\n",
    "                        r_alg_factors = r_alg_depth_factors\n",
    "                        break\n",
    "                \n",
    "            if (not (alg_fact_base_match or depth_alg_fact_base_match)): #If a + bθ is not factorable over any algebraic factor base\n",
    "                continue\n",
    "            \n",
    "\n",
    "            #if (r_alg_factors == None):\n",
    "            #    print(r_alg_factors)\n",
    "            #    print('R_alg_factors is none, uh oh!','\\n','a:',a,'b:',b,'r_alg:',r_alg)\n",
    "            #    continue\n",
    "            assert(r_alg_factors != None) #FIXME \n",
    "\n",
    "            new_row_r = [0 for j in range(lengthRow)]\n",
    "\n",
    "            if (r >= 0):\n",
    "                new_row_r[0] = 0\n",
    "            else:\n",
    "                new_row_r[0] = 1\n",
    "\n",
    "            for i in range(len(rat_factor_base)):\n",
    "                for j in range(len(r_factors)):\n",
    "                    if (rat_factor_base[i] == r_factors[j][0]):\n",
    "                        new_row_r[i+1] = r_factors[j][1]%2\n",
    "                        break\n",
    "\n",
    "\n",
    "            usedPrimes = set()\n",
    "            thisPrime = None\n",
    "            for i in range(len(alg_factor_base)):\n",
    "\n",
    "                thisR = alg_factor_base[i][0]\n",
    "                thisPrime = alg_factor_base[i][1]\n",
    "\n",
    "                #\"A first degree prime ideal represented by the pair (r, p) divides <a + bθ> if and only if p\n",
    "                #....divides N(a + bθ), which occurs if and only if a ≡ −br (mod p)\"\n",
    "                #.......\n",
    "                #\"Only one such pair can have a ≡ −br (mod p). Such an (r, p) pair is the one\n",
    "                #...that will be “responsible” for counting the number of times p divides N(a + bθ)\"\n",
    "                #...--Matthew E. Briggs\n",
    "                for j in range(len(r_alg_factors)):\n",
    "                    if (thisPrime == r_alg_factors[j][0]):\n",
    "                        if (thisPrime in usedPrimes):\n",
    "                            if ( (a % thisPrime) == (-1*b*thisR)%thisPrime):\n",
    "                                new_row_r[i+1+len(rat_factor_base)] = r_alg_factors[j][1]%2\n",
    "                                #Unsure if I should comment this line out based on the \"only one\" comment from Briggs above      \n",
    "                        else:\n",
    "                            if ((a % thisPrime) == (-1*b*thisR)%thisPrime):\n",
    "                                new_row_r[i+1+len(rat_factor_base)] = r_alg_factors[j][1]%2\n",
    "                                usedPrimes.add(thisPrime)\n",
    "\n",
    "            #counter1 = 0\n",
    "            for i in range(len(quad_character_base)):\n",
    "                s = quad_character_base[i][0]\n",
    "                q = quad_character_base[i][1]\n",
    "\n",
    "                if (kronecker((a + (b*s)),q) == 1):\n",
    "                    new_row_r[i+1+len(rat_factor_base)+len(alg_factor_base)] = 0\n",
    "                    #counter1 += 1\n",
    "                elif (kronecker((a + (b*s)),q) == 0):\n",
    "                    new_row_r[i+1+len(rat_factor_base)+len(alg_factor_base)] = 1\n",
    "                    #counter1 += 1\n",
    "                else:\n",
    "                    new_row_r[i+1+len(rat_factor_base)+len(alg_factor_base)] = 1\n",
    "                    #counter1 += 1\n",
    "\n",
    "            #if (counter1 >= SOME_THRESHOLD):\n",
    "            #    print('Higher probability.')\n",
    "\n",
    "            r_mat = r_mat.insert_row(r_mat.nrows(), new_row_r)\n",
    "            tuples.append((a,b))\n",
    "                \n",
    "                \n",
    "    if (r_mat.nrows() < lengthRow):\n",
    "        print('You\\'re going to need to increase the sieve size to get about',lengthRow-r_mat.nrows(),'more rows.')\n",
    "        return r_mat,tuples,False,lengthRow-r_mat.nrows(),rational_remainingPrimes2,rational_remainingPrimes3\n",
    "    return r_mat,tuples,True,rational_remainingPrimes2,rational_remainingPrimes3\n",
    "\n",
    "\n",
    "if (record_execution_time == True):\n",
    "    theTime = (time.time(), time.process_time())\n",
    "    \n",
    "if (estimate_good_m):\n",
    "    assert(d != None) #FIXME!!\n",
    "    m = GetGoodM(n,d)\n",
    "\n",
    "R.<x> = ZZ[];\n",
    "polynomial_ring = R        \n",
    "    \n",
    "    \n",
    "if (f == None):\n",
    "    f = BaseMExpansion(n,m)\n",
    "    f_check = polynomial_ring(f[::-1])\n",
    "    print('The function being used is:',f_check,'\\n')\n",
    "    assert (f_check.is_monic() and f_check.is_irreducible()) #FIXME!!\n",
    "\n",
    "if (d == None):\n",
    "    d = len(f)-1\n",
    "\n",
    "if (ub_rat_factor_base == None or ub_alg_factor_base == None or ub_quad_character_base == None):\n",
    "    print('Using some default factor base sizes.')\n",
    "if (ub_rat_factor_base == None):\n",
    "    ub_rat_factor_base = 15\n",
    "if (ub_alg_factor_base == None):\n",
    "    ub_alg_factor_base = 30\n",
    "if (ub_quad_character_base == None):\n",
    "    ub_quad_character_base = 9\n",
    "\n",
    "if (find_factor_bases or (alg_factor_base == None or quad_character_base == None)):\n",
    "    rat_factor_base = GetRationalFactorBase(ub_rat_factor_base)\n",
    "    alg_factor_base = GetAlgebraicFactorBase(f,ub_alg_factor_base,d)\n",
    "    quad_character_base = GetQuadraticCharacterBase(f,ub_quad_character_base,d,alg_factor_base)\n",
    "\n",
    "print('Rational Factor Base:',rat_factor_base)\n",
    "print('Algebraic Factor Base:',alg_factor_base)\n",
    "print('Quadratic Character Base:',quad_character_base,'\\n')\n",
    "    \n",
    "lengthRow = len(rat_factor_base)+len(alg_factor_base)+len(quad_character_base)+1\n",
    "\n",
    "error = False\n",
    "\n",
    "if (a_lb == None or a_ub == None or b_lb == None or b_ub == None):\n",
    "    print('Using some default bounds for sieving interval.')\n",
    "if (a_lb == None):\n",
    "    a_lb = -1000\n",
    "if (a_ub == None):\n",
    "    a_ub = 1000\n",
    "if (b_lb == None):\n",
    "    b_lb = 1\n",
    "if (b_ub == None):\n",
    "    b_ub = 60\n",
    "\n",
    "print('Sieving to find more than',lengthRow,'smooth pairs.')\n",
    "result = runIt(n,m,f,d,a_lb,a_ub,b_lb,b_ub,depth,lengthRow,rat_factor_base,alg_factor_base,quad_character_base)\n",
    "r_mat = result[0]\n",
    "tuples = result[1]\n",
    "remainingPrimes3 = result[-1]\n",
    "remainingPrimes2 = result[-2]\n",
    "\n",
    "print(len(tuples),'total smooth pairs found.')\n",
    "\n",
    "M = matrix(GF(2), len(tuples), r_mat)\n",
    "solutions = M.kernel().basis_matrix().rows();\n",
    "\n",
    "print(len(solutions),'total bases found.')\n",
    "print('\\n')\n",
    "\n",
    "nonTrivialFactorizations = 0\n",
    "\n",
    "if (result[2] == True or len(solutions) > 0):\n",
    "    \n",
    "    f = polynomial_ring(f[::-1])\n",
    "    \n",
    "    for i in range(len(solutions)):\n",
    "\n",
    "        tempRow = solutions[i]\n",
    "\n",
    "        integer_pairs = tuples\n",
    "        vector = tempRow\n",
    "\n",
    "        result = compute_difference_of_squares(R,f,m,n,integer_pairs,vector)\n",
    "\n",
    "        if (result[0] == True):\n",
    "            u = result[3]\n",
    "            v = result[4]\n",
    "\n",
    "            if (u != v):\n",
    "                if (gcd(n,u+v) != 1 and gcd(n,u+v) != n):\n",
    "                    if (verbose):\n",
    "                        print(result)\n",
    "                    else:\n",
    "                        print('Square Root 1:',result[3],'Square Root 2:',result[4])\n",
    "                    print('Factors are:',gcd(n,u-v),'and',gcd(n,u+v))\n",
    "                    nonTrivialFactorizations += 1\n",
    "                else:\n",
    "                    if (verbose):\n",
    "                        print('Trivial factorization found: u=',u,'v=',v,'gcd(n,u-v) =',1,'and','gcd(n,u+v) =',n)\n",
    "            else:\n",
    "                if (verbose):\n",
    "                    print('Trivial factorization found:',u,',',v)\n",
    "                    \n",
    "\n",
    "f = open(\"more_rational_factorbase_primes.txt\", \"w\")\n",
    "#print('More primes you could\\'ve added to the rational factor base:')\n",
    "f.write('More primes you could\\'ve added to the rational factor base:')\n",
    "#print('Strong candidates:')\n",
    "f.write('\\n')\n",
    "f.write('Strong candidates:')\n",
    "#print('\\n')\n",
    "f.write('\\n')\n",
    "for x in remainingPrimes3:\n",
    "    f.write(str(x)+'\\n')\n",
    "    #print(x,end=\" \")\n",
    "#print('Mediocre candidates:')\n",
    "for x in remainingPrimes2:\n",
    "    f.write(str(x)+'\\n')\n",
    "#    print(x,end=\" \")\n",
    "f.close()\n",
    "\n",
    "print('Non-trivial factorizations:',nonTrivialFactorizations)\n",
    "\n",
    "if (record_execution_time == True):\n",
    "    theTime2 = (time.time(), time.process_time())\n",
    "    print('Elapsed seconds:',theTime2[0]-theTime[0])\n",
    "    print('Elapsed CPU time:',theTime2[1]-theTime[1])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.3",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}