Make factorize() iterative instead of recursive

1 files changed, 50 insertions, 22 deletions

diff --git a/doc/examples/factor.cpp b/doc/examples/factor.cpp
index 24335ceb4..baa3e9a57 100644
--- a/doc/examples/factor.cpp
+++ b/doc/examples/factor.cpp
@@ -41,43 +41,71 @@ BigInt rho(const BigInt& n)
    return 0;
    }
 
-void concat(std::vector<BigInt>& x, const std::vector<BigInt>& y)
-   {
-   for(u32bit j = 0; j != y.size(); j++)
-      x.push_back(y[j]);
-   }
-
-std::vector<BigInt> factorize(const BigInt& n)
+std::vector<BigInt> remove_small_factors(BigInt& n)
    {
    std::vector<BigInt> factors;
 
-   if(n <= 1) // no prime factors at all
-      return factors;
+   while(n.is_even())
+      {
+      factors.push_back(2);
+      n /= 2;
+      }
 
-   if(is_prime(n)) // just n itself
+   for(u32bit j = 0; j != PRIME_TABLE_SIZE; j++)
       {
-      factors.push_back(n);
-      return factors;
+      if(n < PRIMES[j])
+         break;
+
+      BigInt x = gcd(n, PRIMES[j]);
+
+      if(x != 1)
+         {
+         n /= x;
+
+         u32bit occurs = 0;
+         while(x != 1)
+            {
+            x /= PRIMES[j];
+            occurs++;
+            }
+
+         for(u32bit k = 0; k != occurs; k++)
+            factors.push_back(PRIMES[j]);
+         }
       }
 
-   if(n.is_even())
+   return factors;
+   }
+
+std::vector<BigInt> factorize(const BigInt& n_in)
+   {
+   BigInt n = n_in;
+   std::vector<BigInt> factors = remove_small_factors(n);
+
+   if(is_prime(n))
       {
-      factors.push_back(2);
-      concat(factors, factorize(n / 2));
+      factors.push_back(n);
       return factors;
       }
 
-   BigInt factor = 0;
-   while(factor == 0)
-      factor = rho(n);
+   while(n != 1)
+      {
+      if(is_prime(n))
+         {
+         factors.push_back(n);
+         break;
+         }
 
-   concat(factors, factorize(factor));
-   concat(factors, factorize(n / factor));
+      BigInt a_factor = 0;
+      while(a_factor == 0)
+         a_factor = rho(n);
 
+      factors.push_back(a_factor);
+      n /= a_factor;
+      }
    return factors;
    }
 
-
 int main(int argc, char* argv[])
    {
    if(argc != 2)
@@ -91,8 +119,8 @@ int main(int argc, char* argv[])
       LibraryInitializer init;
 
       BigInt n(argv[1]);
-      std::vector<BigInt> factors = factorize(n);
 
+      std::vector<BigInt> factors = factorize(n);
       std::sort(factors.begin(), factors.end());
 
       std::cout << n << ": ";