Factoring by Trial Division

We implement the algorithm for finding the prime factors  of an integer by trial division by primes.  We assume that  a table of prime up to 10^5 is available, either using the Sieve of Eratosthenes or using the built-in Prime command of Mathematica.   This table  is loaded  before we define the function.  In the first program we divide  until all the factors are found or the primes are exhausted.  If the number is not fully factored the method does not check if the remaining factor is prime.  The variable i  is used to index the current prime. The first While statement stops if  i> t, the number of primes or if m=1, i.e. the number is factored.  Inside the loop, the second while statement determines the  exponent of the  current prime divisor. If non--zero this is added to the list of factors. Once we exit the main loop, we must still check if the number is fully factored. If  m \ne 1 then we add this number to the list of factors. At this point, one can insert a  primality test (if a fast test is available) to see if the complete factorization is obtained.  Note that we use the variable  m in the program as we cannot assign to the variable n  which is the argument of the function.

Trial Division by Primes

P=Table[Prime[i], {i,1,100000}]
trialdivide[n_]:=Module[ {factors={}, t,e,m,i},
           t=Length[P];e=0; m=n;
           i=1;
           While[ i<=t && m!=1,
             e=0;
             While[ Mod[m,P[[i]]]==0,e++;
               m=Quotient[m,P[[i]]]];
             If[ e!=0, AppendTo[factors,{P[[i]],e}]];
             i++];
           If[m!=1,   AppendTo[factors, {m,1}]];
           Return[factors]]

trialdivide[1001]

Let us test this on a larger number.

trialdivide[10^17-1]

Here   two small prime factors are found and then a much larger number is returned. This number
 has no prime factor less than $10^5$, but is not necessarily prime as  it is greater than $10^{10}$.  We can use the built--in {\bf PrimeQ} to check that the number is not prime.

PrimeQ[(10^17-1)/9]

To test the speed of the process, we can time it using the {\bf Timing} command. This can be placed after the command using two  slashes. The output is a list with the time taken as the first argument and the result of the computation as the second argument.

trialdivide[ 911 911 Prime[10000]]//Timing

The  program given above does not check if the  current prime $p$ exceeds $\sqrt{m}$. In this case we know that the remaining integer is prime.  This change can be incorporated into the program as follows.  As computing the square root is a time consuming operation, we should do this only when $m$ changes, i.e. when a factor has been found. For this purpose, we introduce an extra variable $x$.  Now the algorithm can terminate for three reasons. The number of primes is exhausted, all the factors are found or the remaining factor is prime.  We test these conditions by adding a Boolean variable {\bf done} which becomes  {\bf True} when the second two conditions are satisfied.  We check this variable and  $m$ after we exit the main loop to  test  under which condition the algorithm terminated. Note that incorporating this extra variable if not necessary, and one can have three conditions in the first statement of the While loop. The procedure prints a message when an incomplete factorization is found.

newtrialdivide[n_]:= 
            Module[ {i=1,factorlist={},x,t,m,done=False},
                 m=n; x=N[Sqrt[m]];t=Length[P];
                 While[!done && i<=t,
                        p=P[[i]];
                        e=0;
                        While[Mod[m,p]==0, e=e+1;
                            m=m/p];
                        If[e!=0, AppendTo[factorlist,{p,e}];
                                  x=N[Sqrt[m]]];
                        If[m==1 || p>x , 
                              done=True];
                        i=i+1];
                 If[m==1,Return[factorlist],
                    AppendTo[factorlist, {m,1}];
                 If[done, Return[factorlist],
                           Print["Incomplete factorization"];
                           Return[factorlist]]]]

b=Prime[6000]

a b

newtrialdivide[4855506841]//Timing

Compare our program with the built--in FactorInteger.

FactorInteger[ 4855506841]//Timing

The predefined function is much faster as the algorithm is written in machine language.  For routine factorization programs it is  reasonable to use the built--in function.  Unlike our function, which will terminate after a fixed number of divisions, {\bf FactorInteger} does not terminate.   The second function actually turns out to be slower, though it performs fewer divisions because of the additional variables and the square root computation. This happens in {\it Mathematica} because it is an interpreted language, not compiled like C or Pascal.

trialdivide[ 4855506841]//Timing

newtrialdivide[22456197671 ]//Timing

FactorInteger[ 22456197671]//Timing

trialdivide[ 22456197671]//Timing

When using {\bf FactorInteger} for large numbers it is better to stop when a factor has been found.  This can be accomplished by setting the flag {\bf FactorComplete} to be {\bf False}.

n= Prime[6000] 115751 100000000063

??FactorInteger

newtrialdivide[n]//Timing

trialdivide[n]//Timing

FactorInteger[n, FactorComplete->False]//Timing

The complete factorization can also be obtained in this case using {\bf FactorInteger} as $n$ has two  prime factors less than $10^7$ and the square root of the largest prime factor is less than  $10^6$.

FactorInteger[n]//Timing

To see the limitations of the procedure we can try  a number which is the product of two primes each of which has more than $10$ digits.  The process did not stop after about 8 hours and needed to be aborted.

FactorInteger[ 10000000283 8748934379]//Timing