From: Christoph Fuerst <ch.fuerst@gmx.at>
Date: Sat, 25 Mar 2017 08:58:48 +0000 (+0100)
Subject: Algorithm for prime number factorization (according to Knuth)
X-Git-Url: http://git.risc.jku.at/gitweb/?a=commitdiff_plain;h=0af05f4277c5711fa4ccfbf3984f283bbd574355;p=cfuerst%2Fformal-numbers.git

Algorithm for prime number factorization (according to Knuth)
---

diff --git a/src/primefactors.cpp b/src/primefactors.cpp
new file mode 100644
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);
+}
+