From 284ce56948aae60134a1073c6499c7fb97ffb326 Mon Sep 17 00:00:00 2001 From: tomsmeding Date: Sat, 20 Apr 2019 21:15:46 +0200 Subject: More reorganisation --- aberth/Makefile | 6 +- aberth/compute_host.cpp | 58 +++++++++++++ aberth/compute_host.h | 9 ++ aberth/defs.h | 2 - aberth/host_aberth.cpp | 96 +++++++++++++++++++++ aberth/host_aberth.h | 12 +++ aberth/kernel.h | 1 + aberth/main.cpp | 224 ++---------------------------------------------- aberth/poly.cpp | 28 ++++++ aberth/poly.h | 42 +++++++++ aberth/polygen.cpp | 43 ++++++++++ aberth/polygen.h | 25 ++++++ aberth/util.cpp | 27 ------ aberth/util.h | 2 - 14 files changed, 323 insertions(+), 252 deletions(-) create mode 100644 aberth/compute_host.cpp create mode 100644 aberth/compute_host.h create mode 100644 aberth/host_aberth.cpp create mode 100644 aberth/host_aberth.h create mode 100644 aberth/poly.cpp create mode 100644 aberth/poly.h create mode 100644 aberth/polygen.cpp create mode 100644 aberth/polygen.h delete mode 100644 aberth/util.cpp diff --git a/aberth/Makefile b/aberth/Makefile index 99f0e16..cd59309 100644 --- a/aberth/Makefile +++ b/aberth/Makefile @@ -1,6 +1,6 @@ CXX = g++ -CXXFLAGS = -Wall -Wextra -std=c++17 -fwrapv $(OPTFLAGS) -OPTFLAGS = -O3 -ffast-math -march=native -mtune=native +CXXFLAGS = -Wall -Wextra -std=c++17 $(OPTFLAGS) +OPTFLAGS = -O3 -ffast-math -march=native -mtune=native -flto LDFLAGS = -pthread ifeq ($(shell uname),Darwin) LDFLAGS += -framework OpenCL @@ -24,7 +24,7 @@ remake: clean aberth: $(OBJECTS) aberth_kernel.o ../lodepng.o - $(CXX) -o $@ $^ $(LDFLAGS) + $(CXX) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) aberth_kernel.h: aberth_kernel.fut futhark opencl --library $^ diff --git a/aberth/compute_host.cpp b/aberth/compute_host.cpp new file mode 100644 index 0000000..38bf43d --- /dev/null +++ b/aberth/compute_host.cpp @@ -0,0 +1,58 @@ +#include +#include +#include +#include +#include "compute_host.h" +#include "host_aberth.h" +#include "polygen.h" +#include "util.h" + +using namespace std; + + +vector computeHost(int W, int H, Com bottomLeft, Com topRight) { + constexpr const int numThreads = 4; + static_assert(ispow2(numThreads)); + + vector counts(W * H); + mutex countsMutex; + + vector jobs = PolyGen::genJobs(numThreads); + assert(jobs.size() == numThreads); + + vector threads(jobs.size()); + for (int i = 0; i < (int)jobs.size(); i++) { + threads[i] = thread([W, H, &counts, &countsMutex, job = jobs[i], bottomLeft, topRight]() { + auto calcIndex = [](double value, double left, double right, int steps) -> int { + return (value - left) / (right - left) * (steps - 1) + 0.5; + }; + auto calcPos = [W, H, bottomLeft, topRight, &calcIndex](Com z) -> pair { + return make_pair( + calcIndex(z.real(), bottomLeft.real(), topRight.real(), W), + calcIndex(z.imag(), bottomLeft.imag(), topRight.imag(), H) + ); + }; + + vector localCounts(W * H); + + Poly poly = job.init; + for (int i = 0; i < job.numItems; i++) { + for (Com z : aberth(poly)) { + int x, y; + tie(x, y) = calcPos(z); + if (0 <= x && x < W && 0 <= y && y < H) { + localCounts[W * y + x]++; + } + } + PolyGen::next(poly); + } + + lock_guard guard(countsMutex); + for (int i = 0; i < W * H; i++) counts[i] += localCounts[i]; + }); + } + + for (thread &th : threads) th.join(); + + return counts; +} diff --git a/aberth/compute_host.h b/aberth/compute_host.h new file mode 100644 index 0000000..7a2ad36 --- /dev/null +++ b/aberth/compute_host.h @@ -0,0 +1,9 @@ +#pragma once + +#include +#include "defs.h" + +using namespace std; + + +vector computeHost(int W, int H, Com bottomLeft, Com topRight); diff --git a/aberth/defs.h b/aberth/defs.h index 40eacf0..41d1f89 100644 --- a/aberth/defs.h +++ b/aberth/defs.h @@ -9,5 +9,3 @@ using namespace std; constexpr const int N = 18; using Com = complex; -using Poly = array; -using AApprox = array; diff --git a/aberth/host_aberth.cpp b/aberth/host_aberth.cpp new file mode 100644 index 0000000..c567f30 --- /dev/null +++ b/aberth/host_aberth.cpp @@ -0,0 +1,96 @@ +#include +#include +#include "host_aberth.h" + +using namespace std; + + +struct AberthState { + // boundPoly is 's' in the stop condition formulated at p.189-190 of + // https://link.springer.com/article/10.1007%2FBF02207694 + + const Poly &poly; + Poly deriv; + Poly boundPoly; + AApprox approx; + double radius; + + AberthState(const Poly &poly); + void regenerate(); + bool step(); + void iterate(); +}; + +static thread_local minstd_rand randgenerator = minstd_rand(random_device()()); + +AberthState::AberthState(const Poly &poly) + : poly(poly), deriv(derivative(poly)), radius(maxRootNorm(poly)) { + + regenerate(); + for (int i = 0; i <= N; i++) { + boundPoly[i] = abs(poly[i]) * (4 * i + 1); + } +} + +void AberthState::regenerate() { + auto genCoord = [this]() { + return uniform_real_distribution(-radius, radius)(randgenerator); + }; + + for (int i = 0; i < N; i++) { + approx[i] = Com(genCoord(), genCoord()); + } +} + +// Lagrange-style step where the new elements are computed in parallel from the previous values +bool AberthState::step() { + array pairs; + for (int i = 0; i < N - 1; i++) { + for (int j = i + 1; j < N; j++) { + pairs[N * i + j] = 1.0 / (approx[i] - approx[j]); + } + } + + bool allConverged = true; + + AApprox newapprox; + AApprox offsets; + for (int i = 0; i < N; i++) { + Com pval = eval(poly, approx[i]); + Com derivval = eval(deriv, poly.size() - 1, approx[i]); + Com quo = pval / derivval; + Com sum = 0; + for (int j = 0; j < i; j++) sum -= pairs[N * j + i]; + for (int j = i + 1; j < N; j++) sum += pairs[N * i + j]; + offsets[i] = quo / (1.0 - quo * sum); + + newapprox[i] = approx[i] - offsets[i]; + + double sval = eval(boundPoly, abs(newapprox[i])); + if (abs(pval) > 1e-5 * sval) allConverged = false; + } + + approx = newapprox; + + return allConverged; +} + +void AberthState::iterate() { + int tries = 1, stepIdx = 1; + + while (!step()) { + stepIdx++; + + if (stepIdx > tries * 100) { + regenerate(); + stepIdx = 0; + tries++; + } + } +} + +AApprox aberth(const Poly &poly) { + AberthState state(poly); + state.iterate(); + return state.approx; +} diff --git a/aberth/host_aberth.h b/aberth/host_aberth.h new file mode 100644 index 0000000..566a84e --- /dev/null +++ b/aberth/host_aberth.h @@ -0,0 +1,12 @@ +#pragma once + +#include +#include "poly.h" + +using namespace std; + + +using AApprox = array; + + +AApprox aberth(const Poly &poly); diff --git a/aberth/kernel.h b/aberth/kernel.h index bd8ea32..3ab2772 100644 --- a/aberth/kernel.h +++ b/aberth/kernel.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include "defs.h" extern "C" { diff --git a/aberth/main.cpp b/aberth/main.cpp index 1cdea3a..152640f 100644 --- a/aberth/main.cpp +++ b/aberth/main.cpp @@ -1,232 +1,20 @@ #include #include -#include -#include -#include #include -#include -#include -#include -#include +#include #include -#include -#include -#include +#include #include +#include #include #include "../lodepng.h" +#include "compute_host.h" #include "defs.h" #include "kernel.h" -#include "util.h" using namespace std; -template -static T eval(const Poly &p, int nterms, T pt) { - T value = p[nterms - 1]; - for (int i = nterms - 2; i >= 0; i--) { - value = pt * value + (double)p[i]; - } - return value; -} - -template -static T eval(const Poly &p, T pt) { - return eval(p, p.size(), pt); -} - -static Poly derivative(const Poly &p) { - Poly res; - for (int i = res.size() - 2; i >= 0; i--) { - res[i] = (i+1) * p[i+1]; - } - return res; -} - -static double maxRootNorm(const Poly &poly) { - // Cauchy's bound: https://en.wikipedia.org/wiki/Geometrical_properties_of_polynomial_roots#Lagrange's_and_Cauchy's_bounds - - double value = 0; - double last = (double)poly.back(); - for (int i = 0; i < (int)poly.size() - 1; i++) { - value = max(value, abs(poly[i] / last)); - } - return 1 + value; -} - -static thread_local minstd_rand randgenerator = minstd_rand(random_device()()); - -struct AberthState { - const Poly &poly; - Poly deriv; - Poly boundPoly; - AApprox approx; - double radius; - - void regenerate() { - auto genCoord = [this]() { - return uniform_real_distribution(-radius, radius)(randgenerator); - }; - for (int i = 0; i < N; i++) { - approx[i] = Com(genCoord(), genCoord()); - } - } - - // boundPoly is 's' in the stop condition formulated at p.189-190 of - // https://link.springer.com/article/10.1007%2FBF02207694 - - AberthState(const Poly &poly) - : poly(poly), deriv(derivative(poly)), radius(maxRootNorm(poly)) { - - regenerate(); - for (int i = 0; i <= N; i++) { - boundPoly[i] = abs(poly[i]) * (4 * i + 1); - } - } - - // Lagrange-style step where the new elements are computed in parallel from the previous values - bool step() { - array pairs; - for (int i = 0; i < N - 1; i++) { - for (int j = i + 1; j < N; j++) { - pairs[N * i + j] = 1.0 / (approx[i] - approx[j]); - } - } - - bool allConverged = true; - - AApprox newapprox; - AApprox offsets; - for (int i = 0; i < N; i++) { - Com pval = eval(poly, approx[i]); - Com derivval = eval(deriv, poly.size() - 1, approx[i]); - Com quo = pval / derivval; - Com sum = 0; - for (int j = 0; j < i; j++) sum -= pairs[N * j + i]; - for (int j = i + 1; j < N; j++) sum += pairs[N * i + j]; - offsets[i] = quo / (1.0 - quo * sum); - - // approx[i] -= offsets[i]; - newapprox[i] = approx[i] - offsets[i]; - - double sval = eval(boundPoly, abs(newapprox[i])); - if (abs(pval) > 1e-5 * sval) allConverged = false; - } - - approx = newapprox; - - return allConverged; - } - - void iterate() { - int tries = 1, stepIdx = 1; - while (!step()) { - stepIdx++; - - if (stepIdx > tries * 100) { - regenerate(); - stepIdx = 0; - tries++; - } - } - } -}; - -static AApprox aberth(const Poly &poly) { - AberthState state(poly); - state.iterate(); - return state.approx; -} - -// Returns whether we just looped around -static bool nextDerbyshire(Poly &poly) { - for (int i = 1; i < (int)poly.size(); i++) { - if (poly[i] == -1) { - poly[i] = 1; - return false; - } - poly[i] = -1; - } - return true; -} - -static Poly derbyshireAtIndex(int index) { - Poly poly; - poly[0] = 1; - for (int i = 1; i <= N; i++) { - poly[i] = index & 1 ? 1 : -1; - index >>= 1; - } - assert(index == 0); - return poly; -} - -struct Job { - Poly init; - int numItems; -}; - -static vector derbyshireJobs(int targetJobs) { - int njobs = min(1 << N, ceil2(targetJobs)); - int jobsize = (1 << N) / njobs; - - vector jobs(njobs); - for (int i = 0; i < njobs; i++) { - jobs[i].init = derbyshireAtIndex(i * jobsize); - jobs[i].numItems = jobsize; - } - - return jobs; -} - -static vector computeCounts(int W, int H, Com bottomLeft, Com topRight) { - constexpr const int numThreads = 4; - static_assert(ispow2(numThreads)); - - vector counts(W * H); - mutex countsMutex; - - vector jobs = derbyshireJobs(numThreads); - assert(jobs.size() == numThreads); - - vector threads(jobs.size()); - for (int i = 0; i < (int)jobs.size(); i++) { - threads[i] = thread([W, H, &counts, &countsMutex, job = jobs[i], bottomLeft, topRight]() { - auto calcIndex = [](double value, double left, double right, int steps) -> int { - return (value - left) / (right - left) * (steps - 1) + 0.5; - }; - auto calcPos = [W, H, bottomLeft, topRight, &calcIndex](Com z) -> pair { - return make_pair( - calcIndex(z.real(), bottomLeft.real(), topRight.real(), W), - calcIndex(z.imag(), bottomLeft.imag(), topRight.imag(), H) - ); - }; - - vector localCounts(W * H); - - Poly poly = job.init; - for (int i = 0; i < job.numItems; i++) { - for (Com z : aberth(poly)) { - int x, y; - tie(x, y) = calcPos(z); - if (0 <= x && x < W && 0 <= y && y < H) { - localCounts[W * y + x]++; - } - } - nextDerbyshire(poly); - } - - lock_guard guard(countsMutex); - for (int i = 0; i < W * H; i++) counts[i] += localCounts[i]; - }); - } - - for (thread &th : threads) th.join(); - - return counts; -} - static void writeCounts(int W, int H, const vector &counts, const char *fname) { ofstream f(fname); f << W << ' ' << H << '\n'; @@ -289,12 +77,12 @@ int main(int argc, char **argv) { const Com bottomLeft = Com(-1.5, -1.5); const Com topRight = Com(1.5, 1.5); - // counts = computeCounts(W, H, bottomLeft, topRight); + // counts = computeHost(W, H, bottomLeft, topRight); Kernel().run_chunked(counts, W, H, bottomLeft, topRight, 42, 1 << 14); // Kernel().run_all(counts, W, H, bottomLeft, topRight, 42); writeCounts(W, H, counts, "out.txt"); - } else if (argc == 2) { + } else if (argc == 2 && strcmp(argv[1], "-h") != 0) { tie(W, H, counts) = readCounts(argv[1]); } else { cerr << "Usage: " << argv[0] << " -- compute and draw" << endl; diff --git a/aberth/poly.cpp b/aberth/poly.cpp new file mode 100644 index 0000000..cf19816 --- /dev/null +++ b/aberth/poly.cpp @@ -0,0 +1,28 @@ +#include +#include "poly.h" + + +ostream& operator<<(ostream &os, const Poly &p) { + static const char *supers[10] = { + "⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹" + }; + + os << p[0]; + + for (int i = 1; i < (int)p.size(); i++) { + if (p[i] < 0) os << " - " << -p[i]; + else if (p[i] > 0) os << " + " << p[i]; + else continue; + + os << "x"; + + if (i == 1) continue; + + ostringstream ss; + ss << i; + string s = ss.str(); + for (char c : s) os << supers[c - '0']; + } + + return os; +} diff --git a/aberth/poly.h b/aberth/poly.h new file mode 100644 index 0000000..82f04a2 --- /dev/null +++ b/aberth/poly.h @@ -0,0 +1,42 @@ +#pragma once + +#include "defs.h" + + +using Poly = array; + + +template +T eval(const Poly &p, int nterms, T pt) { + T value = p[nterms - 1]; + for (int i = nterms - 2; i >= 0; i--) { + value = pt * value + (double)p[i]; + } + return value; +} + +template +T eval(const Poly &p, T pt) { + return eval(p, p.size(), pt); +} + +inline Poly derivative(const Poly &p) { + Poly res; + for (int i = res.size() - 2; i >= 0; i--) { + res[i] = (i+1) * p[i+1]; + } + return res; +} + +inline double maxRootNorm(const Poly &poly) { + // Cauchy's bound: https://en.wikipedia.org/wiki/Geometrical_properties_of_polynomial_roots#Lagrange's_and_Cauchy's_bounds + + double value = 0; + double last = (double)poly.back(); + for (int i = 0; i < (int)poly.size() - 1; i++) { + value = max(value, abs(poly[i] / last)); + } + return 1 + value; +} + +ostream& operator<<(ostream &os, const Poly &p); diff --git a/aberth/polygen.cpp b/aberth/polygen.cpp new file mode 100644 index 0000000..6df89db --- /dev/null +++ b/aberth/polygen.cpp @@ -0,0 +1,43 @@ +#include "polygen.h" +#include "util.h" + +using namespace std; + + +namespace PolyGen::Derbyshire { + // Returns whether we just looped around + bool next(Poly &poly) { + for (int i = 1; i < (int)poly.size(); i++) { + if (poly[i] == -1) { + poly[i] = 1; + return false; + } + poly[i] = -1; + } + return true; + } + + Poly atIndex(int index) { + Poly poly; + poly[0] = 1; + for (int i = 1; i <= N; i++) { + poly[i] = index & 1 ? 1 : -1; + index >>= 1; + } + assert(index == 0); + return poly; + } + + vector genJobs(int targetJobs) { + int njobs = min(1 << N, ceil2(targetJobs)); + int jobsize = (1 << N) / njobs; + + vector jobs(njobs); + for (int i = 0; i < njobs; i++) { + jobs[i].init = atIndex(i * jobsize); + jobs[i].numItems = jobsize; + } + + return jobs; + } +} diff --git a/aberth/polygen.h b/aberth/polygen.h new file mode 100644 index 0000000..1381444 --- /dev/null +++ b/aberth/polygen.h @@ -0,0 +1,25 @@ +#pragma once + +#include +#include "poly.h" + +using namespace std; + + +namespace PolyGen { + struct Job { + Poly init; + int numItems; + }; + + namespace Derbyshire { + // Returns whether we just looped around + bool next(Poly &poly); + + Poly atIndex(int index); + + vector genJobs(int targetJobs); + } + + using namespace Derbyshire; +} diff --git a/aberth/util.cpp b/aberth/util.cpp deleted file mode 100644 index b8f6edc..0000000 --- a/aberth/util.cpp +++ /dev/null @@ -1,27 +0,0 @@ -#include "util.h" - - -ostream& operator<<(ostream &os, const Poly &p) { - static const char *supers[10] = { - "⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹" - }; - - os << p[0]; - - for (int i = 1; i < (int)p.size(); i++) { - if (p[i] < 0) os << " - " << -p[i]; - else if (p[i] > 0) os << " + " << p[i]; - else continue; - - os << "x"; - - if (i == 1) continue; - - ostringstream ss; - ss << i; - string s = ss.str(); - for (char c : s) os << supers[c - '0']; - } - - return os; -} diff --git a/aberth/util.h b/aberth/util.h index c21a7bb..a276e2a 100644 --- a/aberth/util.h +++ b/aberth/util.h @@ -27,5 +27,3 @@ constexpr static T ceil2(T value) { if (value2 == 0) return value << 1; } } - -ostream& operator<<(ostream &os, const Poly &p); -- cgit v1.2.3