aboutsummaryrefslogtreecommitdiff
path: root/aberth
diff options
context:
space:
mode:
authortomsmeding <tom.smeding@gmail.com>2019-07-16 21:57:21 +0200
committertomsmeding <tom.smeding@gmail.com>2019-07-16 21:57:29 +0200
commit3b1767a3f4fd2f0ff2b03d7daa9307fce115b5ea (patch)
tree6aeff844fb24693880ee9ae3b48cd9d3d7b9b6f1 /aberth
parent5e23e6b4986bb31e15e37309ee6b55224a4a630a (diff)
Add support for Christensen in CPU aberth
Diffstat (limited to 'aberth')
-rw-r--r--aberth/compute_host.cpp2
-rw-r--r--aberth/defs.h3
-rw-r--r--aberth/main.cpp1
-rw-r--r--aberth/polygen.cpp38
-rw-r--r--aberth/polygen.h15
-rw-r--r--aberth/util.h11
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;
+}