{ "cells": [ { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Odd part of n=1001\n", "Factorization is [13, 11, 7] \n", "\n" ] } ], "source": [ "'Implementation of Theorem 4.6 from arXiv:2002.02059,'\n", "'An algorithm for factoring positive integers'\n", "\n", "def tfactor(n): #this is the main function that factors the number n\n", " factors=[] # factors will be collected in a list\n", " while n%2 == 0: #get rid of powers of 2 first\n", " factors.append(2)\n", " n=n/2 \n", " print('Odd part of n='+str(n)) # now the fun starts\n", " # spreads=[] # lines of code involving this list are to study the multiple of the factor that arises in the gcd\n", " while n!=1: #if n=1, we're done, otherwise we begin the main loop\n", " seen={} #we will use a dict to keep track of the numbers mod n that we've passed, \n", " k=0 # and how many hops it took to get there\n", " m=n #m is the placeholder for our position in Z/nZ. We start at 0 mod n\n", " '''Main while loop'''\n", " while k<=(n-1)/2: #if n is prime, all (n+1)/2 values will be different so we may need to check all of them\n", " if m in seen: \n", " # print('First occurrence of repeat at k=',seen[m]) # to see how far one must go to find factor\n", " break # exits this while loop because we have found a factor\n", " else:\n", " seen.update({m%n:k}) # otherwise, add the value to the dictionary,\n", " k+=1 #increment k to keep track that were going to hop again\n", " m-=2*k # and drop to the next negative triangular number\n", " m=m%n # take its value mod n\n", " continue # and go back to the beginning of the loop to see if this m value is old or new\n", " '''n is prime case'''\n", " if k==(n+1)/2: # notice while only goes to (n-1)/2, but we add one in the else conditition to arrive here\n", " factors.append(n) # conclude n is prime\n", " n=1 # 1 remains after we divide out, exits while loop\n", " '''Next apply proposition 4.5 to extract factor by taking gcd'''\n", " else: # if we exit the while loop because of a collision, then we:\n", " dist=k-seen[m] # take the difference of positions of the repeated values\n", " fac=gcd(dist,n) # factor is the gcd of n and this distance\n", " factors.append(fac) # add the factor to our list\n", " # spreads.append(dist/fac) # for studying spreads\n", " n=n/fac\n", " print('Factorization is ',factors, '\\n')\n", " #print('Spreads are', spreads)\n", "\n", "'''A slick implementation of Euclid's algorithm for the gcd'''\n", "def gcd(x,y): \n", " while(y): # while y is not 0\n", " x, y = y, x % y # replace by divisor and remainder\n", " return x # returns last non-zero remainder\n", "\n", "semis=[4, 6, 9, 10, 14, 15, 21, 22, 25, 26, 33, 34, 35, 38, 39, 46, 49, 51, 55, 57, 58, 62, 65, 69, 74, 77, 82, 85, 86, 87, 91, 93, 94, 95, 106, 111, 115, 118, 119, 121, 122, 123, 129, 133, 134, 141, 142, 143, 145, 146, 155, 158, 159, 161, 166, 169, 177, 178, 183, 185, 187, 194, 201, 202, 203, 205, 206, 209, 213, 214, 215, 217, 218, 219, 221, 226, 235, 237, 247, 249, 253, 254, 259, 262, 265, 267, 274, 278, 287, 289, 291, 295, 298, 299, 301, 302, 303, 305, 309, 314, 319, 321, 323, 326, 327, 329, 334, 335, 339, 341, 346, 355, 358, 361, 362, 365, 371, 377, 381, 382, 386, 391, 393, 394, 395, 398, 403, 407, 411, 413, 415, 417, 422, 427, 437, 445, 446, 447, 451, 453, 454, 458, 466, 469, 471, 473, 478, 481, 482, 485, 489, 493, 497, 501, 502, 505, 511, 514, 515, 517, 519, 526, 527, 529, 533, 535, 537, 538, 542, 543, 545, 551, 553, 554, 559, 562, 565, 566, 573, 579, 581, 583, 586, 589, 591, 597, 611, 614, 622, 623, 626, 629, 633, 634, 635, 649, 655, 662, 667, 669, 671, 674, 679, 681, 685, 687, 689, 694, 695, 697, 698, 699, 703, 706, 707, 713, 717, 718, 721, 723, 731, 734, 737, 745, 746, 749, 753, 755, 758, 763, 766, 767, 771, 778, 779, 781, 785, 789, 791, 793, 794, 799, 802, 803, 807, 813, 815, 817, 818, 831, 835, 838, 841, 842, 843, 849, 851, 862, 865, 866, 869, 871, 878, 879, 886, 889, 893, 895, 898, 899, 901, 905, 913, 914, 917, 921, 922, 923, 926, 933, 934, 939, 943, 949, 951, 955, 958, 959, 961, 965, 973, 974, 979, 982, 985, 989, 993, 995, 998\n", "]\n", "#for n in semis:\n", "# tfactor(n)\n", "tfactor(1001)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "7 * 11 * 13" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "factor(1001)" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Odd part of n=1562338765544279\n", "Factorization is [2, 2, 2, 23434997, 1321, 109, 463] \n", "\n", "None\n", "646.2216827869415\n" ] } ], "source": [ "'''If you want to time test. Built in SageMath factorization (which is actually built-in PARI \n", "factorization) is still much faster, but this uses a hybrid of several algorithms including\n", "Pollard's rho, Shanks', elliptic curves, and quadratic sieve. '''\n", "from time import time\n", "start=time()\n", "print(tfactor(12498710124354232))\n", "end=time()\n", "print(end-start)" ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.1", "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.3" } }, "nbformat": 4, "nbformat_minor": 2 }