Increase speed by 20-35% by optimizing mp_mul_mod_word_sub().

.gitignore

@@ -1,5 +1,7 @@
 # Makefile artifacts
 *.o
 *.so
+*.x64
 cache-opencl.*
 bin
+__pycache__/

Makefile

@@ -7,7 +7,7 @@ EXECUTABLE=profanity2.x64
 UNAME_S := $(shell uname -s)
 ifeq ($(UNAME_S),Darwin)
 	LDFLAGS=-framework OpenCL
-	CFLAGS=-c -std=c++11 -Wall -mmmx -O2
+	CFLAGS=-c -std=c++11 -Wall -O2
 else
 	LDFLAGS=-s -lOpenCL -mcmodel=large
 	CFLAGS=-c -std=c++11 -Wall -mmmx -O2 -mcmodel=large 

profanity.cl

@@ -277,31 +277,51 @@ mp_word mp_mul_word_add_extra(mp_number * const r, const mp_number * const a, co
 	return *extra < cM ? 1 : (*extra == cM ? cA : 0);
 }
 
-// Multiplies a number with a word, potentially adds modhigher to it, and then subtracts it from en existing number, no extra words, no overflow
-// This is a special function only used for modular multiplication
+// Multiplies a number with a word, potentially adds modhigher to it, and then subtracts it from
+// an existing number, no extra words, no overflow.
+//
+// This is a special function only used for modular multiplication.
+//
+// Optimization (secp256k1 fast reduction) contributed by Rodrigo Madera (@madera).
+//
+// The secp256k1 prime has the special form:
+//
+//   p = 2^256 - pmod, where pmod = 2^32 + 977 = 0x1000003D1
+//
+// Therefore, working modulo 2^256:
+//
+//   q * p = q * (2^256 - pmod) = q * 2^256 - q * pmod == -q * pmod  (mod 2^256)
+//
+//   (r - q * p) mod 2^256 == (r + q * pmod) mod 2^256
+//
+// So instead of multiplying q by the full 256-bit p and subtracting, we multiply q by the
+// 33-bit pmod and add. This reduces the amount of bits used, giving us 20-35% speed improvements.
+//
 void mp_mul_mod_word_sub(mp_number * const r, const mp_word w, const bool withModHigher) {
-	// Having these numbers declared here instead of using the global values in __constant address space seems to lead
-	// to better optimizations by the compiler on my GTX 1070.
-	mp_number mod = { { 0xfffffc2f, 0xfffffffe, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff} };
-	mp_number modhigher = { {0x00000000, 0xfffffc2f, 0xfffffffe, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff} };
-
-	mp_word cM = 0; // Carry for multiplication
-	mp_word cS = 0; // Carry for subtraction
-	mp_word tS = 0; // Temporary storage for subtraction
-	mp_word tM = 0; // Temporary storage for multiplication
-	mp_word cA = 0; // Carry for addition of modhigher
-
-	for (mp_word i = 0; i < MP_WORDS; ++i) {
-		tM = (mod.d[i] * w + cM);
-		cM = mul_hi(mod.d[i], w) + (tM < cM);
-
-		tM += (withModHigher ? modhigher.d[i] : 0) + cA;
-		cA = tM < (withModHigher ? modhigher.d[i] : 0) ? 1 : (tM == (withModHigher ? modhigher.d[i] : 0) ? cA : 0);
-
-		tS = r->d[i] - tM - cS;
-		cS = tS > r->d[i] ? 1 : (tS == r->d[i] ? cS : 0);
-
-		r->d[i] = tS;
+	const mp_word lo977 = 977u * w;
+	const mp_word hi977 = mul_hi(977u, w);
+
+	const mp_word p0 = lo977;
+	const ulong p1_full = (ulong)w + hi977 + (withModHigher ? 0x000003D1u : 0u);
+	const mp_word p1 = (mp_word)p1_full;
+	const mp_word p2 = (mp_word)(p1_full >> 32) + (withModHigher ? 1u : 0u);
+
+	ulong s = (ulong)r->d[0] + p0;
+	r->d[0] = (mp_word)s;
+	mp_word c = (mp_word)(s >> 32);
+
+	s = (ulong)r->d[1] + p1 + c;
+	r->d[1] = (mp_word)s;
+	c = (mp_word)(s >> 32);
+
+	s = (ulong)r->d[2] + p2 + c;
+	r->d[2] = (mp_word)s;
+	c = (mp_word)(s >> 32);
+
+	for (mp_word i = 3; i < MP_WORDS; ++i) {
+		s = (ulong)r->d[i] + c;
+		r->d[i] = (mp_word)s;
+		c = (mp_word)(s >> 32);
 	}
 }
 

profanity.cpp

@@ -33,24 +33,46 @@ std::string readFile(const char * const szFilename)
 	return contents.str();
 }
 
-std::vector<cl_device_id> getAllDevices(cl_device_type deviceType = CL_DEVICE_TYPE_GPU)
-{
+std::vector<cl_device_id> getAllDevices(cl_device_type deviceType = CL_DEVICE_TYPE_GPU) {
 	std::vector<cl_device_id> vDevices;
 
 	cl_uint platformIdCount = 0;
-	clGetPlatformIDs (0, NULL, &platformIdCount);
+	cl_int err = clGetPlatformIDs(0, NULL, &platformIdCount);
+	if (err != CL_SUCCESS || platformIdCount == 0) {
+		std::cerr << "warning: no OpenCL platforms found, err = " << err << std::endl;
+		return vDevices;
+	}
+
+	std::vector<cl_platform_id> platformIds(platformIdCount);
+	err = clGetPlatformIDs(platformIdCount, platformIds.data(), NULL);
+	if (err != CL_SUCCESS) {
+		std::cerr << "warning: failed to enumerate OpenCL platforms, err = " << err << std::endl;
+		return vDevices;
+	}
 
-	std::vector<cl_platform_id> platformIds (platformIdCount);
-	clGetPlatformIDs (platformIdCount, platformIds.data (), NULL);
+	for (auto it = platformIds.cbegin(); it != platformIds.cend(); ++it) {
+		cl_uint countDevice = 0;
 
-	for( auto it = platformIds.cbegin(); it != platformIds.cend(); ++it ) {
-		cl_uint countDevice;
-		clGetDeviceIDs(*it, deviceType, 0, NULL, &countDevice);
+		err = clGetDeviceIDs(*it, deviceType, 0, NULL, &countDevice);
+		if (err != CL_SUCCESS || countDevice == 0) {
+			char platformName[256] = {0};
+			clGetPlatformInfo(*it, CL_PLATFORM_NAME, sizeof(platformName), platformName, NULL);
+			std::cerr << "warning: skipping OpenCL platform without usable GPU devices: "
+			          << platformName << ", err = " << err << std::endl;
+			continue;
+		}
 
 		std::vector<cl_device_id> deviceIds(countDevice);
-		clGetDeviceIDs(*it, deviceType, countDevice, deviceIds.data(), &countDevice);
+		err = clGetDeviceIDs(*it, deviceType, countDevice, deviceIds.data(), &countDevice);
+		if (err != CL_SUCCESS) {
+			char platformName[256] = {0};
+			clGetPlatformInfo(*it, CL_PLATFORM_NAME, sizeof(platformName), platformName, NULL);
+			std::cerr << "warning: failed to get GPU devices from platform: "
+			          << platformName << ", err = " << err << std::endl;
+			continue;
+		}
 
-		std::copy( deviceIds.begin(), deviceIds.end(), std::back_inserter(vDevices) );
+		std::copy(deviceIds.begin(), deviceIds.end(), std::back_inserter(vDevices));
 	}
 
 	return vDevices;

tests/Makefile

@@ -0,0 +1,32 @@
+# Builds the mp-math correctness tests and benchmarks. Mirrors the style of the
+# main Makefile. Requires an OpenCL runtime; on machines without a GPU, PoCL
+# (http://portablecl.org, e.g. `apt install pocl-opencl-icd`) works fine.
+#
+# Usage:
+#   make
+#   ./test_correctness_pr49.x64 [num_random_cases]
+#   ./bench_mod_mul_pr49.x64 [global_size] [iterations] [repetitions]
+
+CC=g++
+CDEFINES=
+EXECUTABLES=test_correctness_pr49.x64 bench_mod_mul_pr49.x64
+
+UNAME_S := $(shell uname -s)
+ifeq ($(UNAME_S),Darwin)
+	LDFLAGS=-framework OpenCL
+	CFLAGS=-std=c++11 -Wall -O2
+else
+	LDFLAGS=-lOpenCL
+	CFLAGS=-std=c++11 -Wall -O2
+endif
+
+all: $(EXECUTABLES)
+
+test_correctness_pr49.x64: test_correctness_pr49.cpp testutil.hpp ../types.hpp
+	$(CC) $(CFLAGS) $(CDEFINES) test_correctness_pr49.cpp $(LDFLAGS) -o $@
+
+bench_mod_mul_pr49.x64: bench_mod_mul_pr49.cpp testutil.hpp ../types.hpp
+	$(CC) $(CFLAGS) $(CDEFINES) bench_mod_mul_pr49.cpp $(LDFLAGS) -o $@
+
+clean:
+	rm -rf $(EXECUTABLES)

tests/bench_mod_mul_pr49.cpp

@@ -0,0 +1,118 @@
+/* bench_mod_mul_pr49.cpp
+ * ======================
+ * Benchmark: mp_mod_mul with the old (pre-PR#49) vs new (optimized) mp_mul_mod_word_sub.
+ *
+ * Each work item performs a dependency-chained sequence of modular multiplications
+ * x = x * y (mod p) using the real kernels compiled from profanity.cl + tests/harness.cl.
+ * The chained data dependency prevents the compiler from eliminating work.
+ *
+ * Both variants are also cross-checked: starting from identical inputs, their final
+ * outputs must be bit-exact identical.
+ *
+ * Note: absolute numbers on a CPU OpenCL implementation (e.g. PoCL) are not
+ * representative of GPU throughput, but the relative old-vs-new comparison exercises
+ * exactly the code path the PR optimizes.
+ *
+ * Build & run (see tests/Makefile):
+ *   cd tests && make && ./bench_mod_mul_pr49.x64 [global_size] [iterations] [repetitions]
+ */
+
+#include <chrono>
+#include <cstdio>
+#include <random>
+#include <vector>
+
+#include "testutil.hpp"
+
+struct BenchResult {
+	double bestSeconds;
+	std::vector<mp_number> output;
+};
+
+static BenchResult benchVariant(const ClSetup & s, const char * const kernelName,
+	const std::vector<mp_number> & x, const std::vector<mp_number> & y,
+	const uint32_t iterations, const int reps)
+{
+	const size_t count = x.size();
+	cl_int err;
+	cl_kernel kernel = clCreateKernel(s.program, kernelName, &err);
+	clCheck(err, "clCreateKernel");
+
+	cl_mem yBuf = clCreateBuffer(s.context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, count * sizeof(mp_number), (void *)y.data(), &err);
+	clCheck(err, "clCreateBuffer(y)");
+
+	BenchResult result;
+	result.bestSeconds = -1;
+	result.output.resize(count);
+
+	// reps timed runs + one warmup run (which also produces the output used
+	// for the old-vs-new cross-check)
+	for (int rep = 0; rep <= reps; ++rep) {
+		cl_mem xBuf = clCreateBuffer(s.context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, count * sizeof(mp_number), (void *)x.data(), &err);
+		clCheck(err, "clCreateBuffer(x)");
+		clCheck(clSetKernelArg(kernel, 0, sizeof(cl_mem), &xBuf), "clSetKernelArg(0)");
+		clCheck(clSetKernelArg(kernel, 1, sizeof(cl_mem), &yBuf), "clSetKernelArg(1)");
+		clCheck(clSetKernelArg(kernel, 2, sizeof(uint32_t), &iterations), "clSetKernelArg(2)");
+		clCheck(clFinish(s.queue), "clFinish(pre)");
+
+		const auto t0 = std::chrono::steady_clock::now();
+		clCheck(clEnqueueNDRangeKernel(s.queue, kernel, 1, NULL, &count, NULL, 0, NULL, NULL), "clEnqueueNDRangeKernel");
+		clCheck(clFinish(s.queue), "clFinish");
+		const double seconds = std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count();
+
+		if (rep == 0) {
+			clCheck(clEnqueueReadBuffer(s.queue, xBuf, CL_TRUE, 0, count * sizeof(mp_number), result.output.data(), 0, NULL, NULL), "clEnqueueReadBuffer");
+			clCheck(clFinish(s.queue), "clFinish(read)");
+		} else if (result.bestSeconds < 0 || seconds < result.bestSeconds) {
+			result.bestSeconds = seconds;
+		}
+
+		clReleaseMemObject(xBuf);
+	}
+
+	clReleaseMemObject(yBuf);
+	clReleaseKernel(kernel);
+	return result;
+}
+
+int main(int argc, char * * argv) {
+	const size_t globalSize = argc > 1 ? (size_t)std::atol(argv[1]) : 4096;
+	const uint32_t iterations = argc > 2 ? (uint32_t)std::atol(argv[2]) : 2000;
+	const int reps = argc > 3 ? std::atoi(argv[3]) : 5;
+
+	std::mt19937_64 rng(0xBEEF);
+	std::vector<mp_number> x(globalSize), y(globalSize);
+	for (size_t i = 0; i < globalSize; ++i) {
+		for (int j = 0; j < MP_NWORDS; ++j) {
+			x[i].d[j] = (uint32_t)rng();
+			y[i].d[j] = (uint32_t)rng();
+		}
+	}
+
+	const ClSetup s = clSetup();
+	std::printf("Work items: %zu, chained mod-muls per item: %u, repetitions: %d (best time taken)\n",
+		globalSize, iterations, reps);
+
+	const BenchResult oldResult = benchVariant(s, "bench_mod_mul_old", x, y, iterations, reps);
+	const BenchResult newResult = benchVariant(s, "bench_mod_mul_new", x, y, iterations, reps);
+
+	size_t mismatches = 0;
+	for (size_t i = 0; i < globalSize; ++i) {
+		if (!mpEqual(oldResult.output[i], newResult.output[i])) {
+			++mismatches;
+		}
+	}
+	if (mismatches) {
+		std::printf("CROSS-CHECK FAILED: outputs differ in %zu/%zu work items\n", mismatches, globalSize);
+		return 1;
+	}
+	std::printf("Cross-check: outputs of old and new variants are bit-exact identical (%zu work items x %u chained mod-muls).\n\n",
+		globalSize, iterations);
+
+	const double totalMuls = (double)globalSize * iterations;
+	std::printf("  old (master):  %9.2f ms  (%8.2f Mmul/s)\n", oldResult.bestSeconds * 1e3, totalMuls / oldResult.bestSeconds / 1e6);
+	std::printf("  new (PR #49):  %9.2f ms  (%8.2f Mmul/s)\n", newResult.bestSeconds * 1e3, totalMuls / newResult.bestSeconds / 1e6);
+	std::printf("  speedup: x%.3f (%+.1f%%)\n", oldResult.bestSeconds / newResult.bestSeconds,
+		(oldResult.bestSeconds / newResult.bestSeconds - 1) * 100);
+	return 0;
+}