Algorithm for prime number factorization (according to Knuth)

diff --git a/src/primefactors.cpp b/src/primefactors.cpp
new file mode 100644 (file)
index 0000000..20c77bd
--- /dev/null
+++ b/src/primefactors.cpp
@@ -0,0 +1,118 @@
+/*
+ * primefactors.cpp
+ *
+ *  Created on: 25.03.2017
+ *      Author: christoph
+ */
+
+/* Build Log:
+ *
+make all
+Building file: ../src/primefactors.cpp
+Invoking: GCC C++ Compiler
+g++ -O0 -g3 -Wall -c -fmessage-length=0 -MMD -MP -MF"src/primefactors.d" -MT"src/primefactors.o" -o "src/primefactors.o" "../src/primefactors.cpp"
+Finished building: ../src/primefactors.cpp
+
+Building target: PrimeFactors
+Invoking: GCC C++ Linker
+g++  -o "PrimeFactors"  ./src/primefactors.o   -lgmpxx -lgmp
+Finished building target: PrimeFactors
+ *
+ */
+
+#include <iostream>
+#include <vector>
+#include <cmath>
+#include <gmpxx.h>
+
+using namespace std;
+
+/* Algorithm A from Knuth's Art of Computer Programm 2: Seminumerical Algorithms */
+/* page 380 */
+void factor_int(mpz_class N, mpz_class* d, unsigned long& nooffactors, mpz_class *factors)
+{
+       /* Step A1: initialize t = 0, k = 0, n = N */
+       unsigned long t=0;
+       unsigned long k=0;
+       mpz_class n = N;
+
+    mpz_class* p  = new mpz_class[n.get_ui()];
+
+       mpz_class q, r;
+
+       /* Step A2: n=1? */
+    while(n != 1)
+    {
+           /* Step A3: q = floor(n/dk), r = n mod dk */
+               q = n/d[k];
+               r = n%d[k];
+
+               /* Step A4: r != 0 */
+               if(r != 0)
+               {
+                       /* Step A6: Low quotient */
+                       if(q > d[k])
+                       {
+                               k++;
+                       }
+                       else
+                       {
+                               p[t] = n;
+                               n = 1;
+                       }
+               }
+               else
+               {
+                       /* Step A5: Factor found */
+                       p[t] = d[k];
+                       t++;
+                       n = q;
+               }
+       }
+
+    for(size_t i=0;i<=t;i++)
+       factors[i] = p[i];
+    nooffactors = t;
+
+}
+
+int main()
+{
+   mpz_class N;
+   mpz_class *trialdivisors, *factors;
+   unsigned long max;
+   unsigned long nooffactors;
+
+   cout << "Enter Integer to be factored.... ";
+   cin >> N;
+   cout << endl;
+
+   max = (unsigned long)(sqrt(N.get_ui()));
+
+   trialdivisors = new mpz_class[max];
+   factors = new mpz_class[max];
+
+   trialdivisors[0] = 2;
+   trialdivisors[1] = 3;
+   trialdivisors[2] = 5;
+   for(unsigned int i=3;i<max;i++)
+   {
+      if(i%2 == 1)
+      {
+         trialdivisors[i] = trialdivisors[i-1]+2;
+      }
+      else
+      {
+         trialdivisors[i] = trialdivisors[i-1]+4;
+      }
+   }
+
+   factor_int(N, trialdivisors, nooffactors, factors);
+
+   for(unsigned long i=0;i<=nooffactors;i++)
+          cout << factors[i] << " ";
+   cout << endl;
+
+   return (0);
+}
+