Add support for Christensen in CPU aberth

6 files changed, 66 insertions, 4 deletions

diff --git a/aberth/compute_host.cpp b/aberth/compute_host.cpp
index 343318f..0af20df 100644
--- a/aberth/compute_host.cpp
+++ b/aberth/compute_host.cpp
@@ -14,6 +14,8 @@ vector<int> computeHost(int W, int H, Com bottomLeft, Com topRight) {
 	constexpr const int numThreads = 4;
 	static_assert(ispow2(numThreads));
 
+	cerr << "Finding roots for " << PolyGen::numPolys() << " polynomials" << endl;
+
 	vector<int> counts(W * H);
 	mutex countsMutex;
 
diff --git a/aberth/defs.h b/aberth/defs.h
index 41d1f89..1f0ab83 100644
--- a/aberth/defs.h
+++ b/aberth/defs.h
@@ -6,6 +6,7 @@
 using namespace std;
 
 
-constexpr const int N = 18;
+constexpr const int N = 18;  // e.g. for Derbyshire
+// constexpr const int N = 5;  // e.g. for Christensen
 
 using Com = complex<double>;
diff --git a/aberth/main.cpp b/aberth/main.cpp
index 152640f..3f87aed 100644
--- a/aberth/main.cpp
+++ b/aberth/main.cpp
@@ -74,6 +74,7 @@ int main(int argc, char **argv) {
 
 	if (argc <= 1) {
 		W = H = 900;
+		// Use 2.5 everywhere for Christensen
 		const Com bottomLeft = Com(-1.5, -1.5);
 		const Com topRight = Com(1.5, 1.5);
 
diff --git a/aberth/polygen.cpp b/aberth/polygen.cpp
index 6df89db..e606dc8 100644
--- a/aberth/polygen.cpp
+++ b/aberth/polygen.cpp
@@ -28,9 +28,43 @@ namespace PolyGen::Derbyshire {
         return poly;
     }
 
+    int numPolys() {
+        return 1 << N;
+    }
+}
+
+namespace PolyGen::Christensen {
+    // Returns whether we just looped around
+    bool next(Poly &poly) {
+        for (int i = 0; i < (int)poly.size(); i++) {
+            if (poly[i] < kCoeffBound) {
+                poly[i]++;
+                return false;
+            }
+            poly[i] = -kCoeffBound;
+        }
+        return true;
+    }
+
+    Poly atIndex(int index) {
+        Poly poly;
+        for (int i = 0; i < (int)poly.size(); i++) {
+            poly[i] = index % (2 * kCoeffBound + 1) - kCoeffBound;
+            index /= 2 * kCoeffBound + 1;
+        }
+        assert(index == 0);
+        return poly;
+    }
+
+    int numPolys() {
+        return ipow(2 * kCoeffBound + 1, N + 1);
+    }
+}
+
+namespace PolyGen {
     vector<Job> genJobs(int targetJobs) {
-        int njobs = min(1 << N, ceil2(targetJobs));
-        int jobsize = (1 << N) / njobs;
+        int njobs = min(numPolys(), ceil2(targetJobs));
+        int jobsize = numPolys() / njobs;
 
         vector<Job> jobs(njobs);
         for (int i = 0; i < njobs; i++) {
diff --git a/aberth/polygen.h b/aberth/polygen.h
index 1381444..b0022f1 100644
--- a/aberth/polygen.h
+++ b/aberth/polygen.h
@@ -18,8 +18,21 @@ namespace PolyGen {
 
         Poly atIndex(int index);
 
-        vector<Job> genJobs(int targetJobs);
+        int numPolys();
+    }
+
+    namespace Christensen {
+        constexpr int kCoeffBound = 4;
+
+        // Returns whether we just looped around
+        bool next(Poly &poly);
+
+        Poly atIndex(int index);
+
+        int numPolys();
     }
 
     using namespace Derbyshire;
+
+    vector<Job> genJobs(int targetJobs);
 }
diff --git a/aberth/util.h b/aberth/util.h
index a276e2a..49195b5 100644
--- a/aberth/util.h
+++ b/aberth/util.h
@@ -27,3 +27,14 @@ constexpr static T ceil2(T value) {
 		if (value2 == 0) return value << 1;
 	}
 }
+
+template <typename T>
+constexpr static T ipow(T base, T ex) {
+	T x = 1, y = base;
+	while (ex > 0) {
+		if ((ex & 1) == 1) x *= y;
+		y *= y;
+		ex >>= 1;
+	}
+	return x;
+}