delaybench.cpp raw
1 // delaybench.cpp - Limenka fork delay-function benchmark
2 //
3 // Candidates for the per-block sequential delay (D = 60s, CPU-bound,
4 // low hardware variance, non-shortcuttable, verified by recomputation):
5 //
6 // a1. M=1 schoolbook long division, hardware 128/64 DIV per step
7 // a2. M=1 schoolbook long division, Barrett reciprocal estimate
8 // (division-free: mul/add/sub/compare only - the uniform variant)
9 // b1. M=4 Knuth D, hardware-DIV qhat
10 // b2. M=4 Knuth D, Barrett reciprocal qhat (division-free)
11 // c. quadratic map x = x*x + c mod 2^128 (record only; shortcut
12 // status unresolved)
13 // d. SHA256 single chain (hardware-variance control baseline)
14 //
15 // The dividend stream comes from xorshift64* seeded from the block
16 // input: the work binds to the block, cannot be precomputed, and the
17 // working set is a few words regardless of stream length.
18 //
19 // Build: g++ -O2 -march=native -o delaybench delaybench.cpp
20 #include <cstdint>
21 #include <cstdio>
22 #include <cstring>
23 #include <chrono>
24 #include <cstdlib>
25
26 using u64 = uint64_t;
27 using u128 = unsigned __int128;
28 using u32 = uint32_t;
29 using u8 = uint8_t;
30
31 static u64 prng_state;
32
33 static inline u64 prng_next() {
34 u64 x = prng_state;
35 x ^= x >> 12;
36 x ^= x << 25;
37 x ^= x >> 27;
38 prng_state = x;
39 return x * 0x2545F4914F6CDD1Dull;
40 }
41
42 // ---------------- a1. M=1, hardware 128/64 DIV ----------------
43 static u64 run_div1(u64 d, u64 words, bool verbose) {
44 u64 r = 0;
45 for (u64 i = 0; i < words; i++) {
46 u64 w = prng_next();
47 r = (u64)((((u128)r << 64) | w) % d);
48 if (verbose)
49 printf(" div1 step %llu: r=%016llx\n", (unsigned long long)i, (unsigned long long)r);
50 }
51 return r;
52 }
53
54 // high 128 bits of (u1:u0) * (v1:v0) - GCC truncates u128*u128 to the
55 // low 128 bits, so the schoolbook is done explicitly
56 static inline u64 mulhi128(u64 u1, u64 u0, u64 v1, u64 v0) {
57 u64 hi0 = (u64)((u128)u0 * v0 >> 64);
58 u128 t1 = (u128)u0 * v1 + hi0;
59 u64 lo1 = (u64)t1, hi1 = (u64)(t1 >> 64);
60 u128 t2 = (u128)u1 * v0 + lo1;
61 u64 lo2 = (u64)t2, hi2 = (u64)(t2 >> 64);
62 (void)lo2;
63 return (u64)((u128)u1 * v1 + hi1 + hi2);
64 }
65
66 // ---------------- a2. M=1, Barrett reciprocal (division-free) ----------------
67 // v = floor((2^128 - 1) / d), precomputed once per block.
68 // qhat = (u * v) >> 128 for u = (r << 64 | w) < d*2^64; estimate is
69 // exact or 1 too small; refine with one conditional subtract.
70 static inline u64 barrett_step(u64 r, u64 w, u64 d, u128 v) {
71 u128 u = ((u128)r << 64) | w;
72 u64 q = mulhi128(r, w, (u64)(v >> 64), (u64)v);
73 u128 rem = u - (u128)q * d;
74 if (rem >= d) {
75 rem -= d;
76 q++;
77 }
78 (void)q;
79 return (u64)rem;
80 }
81
82 static u64 run_div1_barrett(u64 d, u64 words, bool verbose) {
83 u128 v = (((u128)0 - 1) << 64 | (u128)0xffffffffffffffffULL) / d; // floor((2^128-1)/d)
84 u64 r = 0;
85 for (u64 i = 0; i < words; i++) {
86 u64 w = prng_next();
87 r = barrett_step(r, w, d, v);
88 if (verbose)
89 printf(" bar1 step %llu: r=%016llx\n", (unsigned long long)i, (unsigned long long)r);
90 }
91 return r;
92 }
93
94 // ---------------- b1. M=4 Knuth D, hardware-DIV qhat ----------------
95 static inline void knuth_step_div(u64* u, const u64* d, u64 w) {
96 u64 q = (u64)((((u128)u[3] << 64) | u[2]) / d[3]);
97 u64 r[5] = {w, u[0], u[1], u[2], u[3]};
98 bool neg = false;
99 for (int i = 0; i < 4; i++) {
100 u128 prod = (u128)q * d[i];
101 u64 lo = (u64)prod, hi = (u64)(prod >> 64);
102 u64 t = r[i] - lo;
103 u64 b0 = (r[i] < lo);
104 r[i] = t;
105 // r[i+1] -= hi + b0, split to avoid the hi+b0 wrap when
106 // hi == 2^64-1 and b0 == 1
107 t = r[i + 1] - b0;
108 u64 b1 = (r[i + 1] < b0);
109 r[i + 1] = t;
110 t = r[i + 1] - hi;
111 b1 |= (r[i + 1] < hi);
112 r[i + 1] = t;
113 for (int j = i + 2; j < 5 && b1; j++) {
114 t = r[j] - 1;
115 b1 = (r[j] < 1);
116 r[j] = t;
117 }
118 if (b1) neg = true; // borrow escaped the 5-word value
119 }
120 // negative test: top bit of the 5-word value (two's complement).
121 // The add-back carry must propagate into r[4].
122 while ((r[4] >> 63) != 0) {
123 u64 carry = 0;
124 for (int i = 0; i < 4; i++) {
125 u64 s = r[i] + d[i];
126 u64 c1 = (s < r[i]);
127 u64 s2 = s + carry;
128 u64 c2 = (s2 < s);
129 r[i] = s2;
130 carry = c1 | c2;
131 }
132 r[4] += carry;
133 }
134 for (int rep = 0; rep < 4; rep++) {
135 // r >= d? (equality also subtracts once)
136 int ge = 1;
137 for (int i = 3; i >= 0; i--) {
138 if (r[i] > d[i]) break;
139 if (r[i] < d[i]) { ge = 0; break; }
140 }
141 if (!ge) break;
142 u64 borrow = 0;
143 for (int i = 0; i < 4; i++) {
144 u64 t = r[i] - d[i];
145 u64 b1 = (r[i] < d[i]);
146 u64 t2 = t - borrow;
147 u64 b2 = (t < borrow);
148 r[i] = t2;
149 borrow = b1 | b2;
150 }
151 }
152 u[0] = r[0]; u[1] = r[1]; u[2] = r[2]; u[3] = r[3];
153 }
154
155 static u64 run_knuth4_div(const u64* d, u64 words, bool verbose) {
156 u64 u[4] = {0, 0, 0, 0};
157 for (u64 i = 0; i < words; i++) {
158 u64 w = prng_next();
159 knuth_step_div(u, d, w);
160 if (verbose)
161 printf(" knuthdiv step %llu: u=%016llx %016llx %016llx %016llx\n",
162 (unsigned long long)i,
163 (unsigned long long)u[0], (unsigned long long)u[1],
164 (unsigned long long)u[2], (unsigned long long)u[3]);
165 }
166 return u[0] ^ u[1] ^ u[2] ^ u[3];
167 }
168
169 // ---------------- b2. M=4 Knuth D, Barrett reciprocal qhat ----------------
170 // v3 = floor((2^128 - 1) / d3), precomputed per block. qhat =
171 // ((u3:u2) * v3) >> 128, exact or off by 1; refined by two compares.
172 static u64 run_knuth4_barrett(const u64* d, u64 words, bool verbose) {
173 u64 u[4] = {0, 0, 0, 0};
174 u128 v3 = (((u128)0 - 1) << 64 | (u128)0xffffffffffffffffULL) / d[3];
175 for (u64 i = 0; i < words; i++) {
176 u64 w = prng_next();
177 u64 q = mulhi128(u[3], u[2], (u64)(v3 >> 64), (u64)v3);
178 // refine: while (u3:u2) < q*d3, q--; while >= (q+1)*d3, q++
179 {
180 u128 top = (u128)u[3] << 64 | u[2];
181 if ((u128)q * d[3] > top) q--;
182 else if (top >= ((u128)q + 1) * d[3]) q++;
183 }
184 u64 r[5] = {w, u[0], u[1], u[2], u[3]};
185 bool neg = false;
186 for (int i = 0; i < 4; i++) {
187 u128 prod = (u128)q * d[i];
188 u64 lo = (u64)prod, hi = (u64)(prod >> 64);
189 u64 t = r[i] - lo;
190 u64 b0 = (r[i] < lo);
191 r[i] = t;
192 // r[i+1] -= hi + b0, split to avoid the hi+b0 wrap
193 t = r[i + 1] - b0;
194 u64 b1 = (r[i + 1] < b0);
195 r[i + 1] = t;
196 t = r[i + 1] - hi;
197 b1 |= (r[i + 1] < hi);
198 r[i + 1] = t;
199 for (int j = i + 2; j < 5 && b1; j++) {
200 t = r[j] - 1;
201 b1 = (r[j] < 1);
202 r[j] = t;
203 }
204 if (b1) neg = true;
205 }
206 // negative test: top bit of the 5-word value (two's
207 // complement). The add-back carry must propagate into r[4].
208 while ((r[4] >> 63) != 0) {
209 u64 carry = 0;
210 for (int i = 0; i < 4; i++) {
211 u64 s = r[i] + d[i];
212 u64 c1 = (s < r[i]);
213 u64 s2 = s + carry;
214 u64 c2 = (s2 < s);
215 r[i] = s2;
216 carry = c1 | c2;
217 }
218 r[4] += carry;
219 }
220 for (int rep = 0; rep < 4; rep++) {
221 int ge = 1;
222 for (int i = 3; i >= 0; i--) {
223 if (r[i] > d[i]) break;
224 if (r[i] < d[i]) { ge = 0; break; }
225 }
226 if (!ge) break;
227 u64 borrow = 0;
228 for (int i = 0; i < 4; i++) {
229 u64 t = r[i] - d[i];
230 u64 b1 = (r[i] < d[i]);
231 u64 t2 = t - borrow;
232 u64 b2 = (t < borrow);
233 r[i] = t2;
234 borrow = b1 | b2;
235 }
236 }
237 u[0] = r[0]; u[1] = r[1]; u[2] = r[2]; u[3] = r[3];
238 if (verbose)
239 printf(" knuthbar step %llu: u=%016llx %016llx %016llx %016llx\n",
240 (unsigned long long)i,
241 (unsigned long long)u[0], (unsigned long long)u[1],
242 (unsigned long long)u[2], (unsigned long long)u[3]);
243 }
244 return u[0] ^ u[1] ^ u[2] ^ u[3];
245 }
246
247 // ---------------- c. quadratic map x = x*x + c mod 2^128 ----------------
248 static u64 run_quad128(u64 c, u64 words) {
249 u128 x = c;
250 for (u64 i = 0; i < words; i++) {
251 x = x * x + c;
252 }
253 return (u64)(x >> 64) ^ (u64)x;
254 }
255
256 // ---------------- d. SHA256 chain (scalar reference) ----------------
257 static inline u32 rotr32(u32 x, int n) { return (x >> n) | (x << (32 - n)); }
258
259 static void sha256_compress(u32* h, const u8* block) {
260 static const u32 K[64] = {
261 0x428a2f98,0x71374491,0xb5c0fbcf,0xe9b5dba5,0x3956c25b,0x59f111f1,0x923f82a4,0xab1c5ed5,
262 0xd807aa98,0x12835b01,0x243185be,0x550c7dc3,0x72be5d74,0x80deb1fe,0x9bdc06a7,0xc19bf174,
263 0xe49b69c1,0xefbe4786,0x0fc19dc6,0x240ca1cc,0x2de92c6f,0x4a7484aa,0x5cb0a9dc,0x76f988da,
264 0x983e5152,0xa831c66d,0xb00327c8,0xbf597fc7,0xc6e00bf3,0xd5a79147,0x06ca6351,0x14292967,
265 0x27b70a85,0x2e1b2138,0x4d2c6dfc,0x53380d13,0x650a7354,0x766a0abb,0x81c2c92e,0x92722c85,
266 0xa2bfe8a1,0xa81a664b,0xc24b8b70,0xc76c51a3,0xd192e819,0xd6990624,0xf40e3585,0x106aa070,
267 0x19a4c116,0x1e376c08,0x2748774c,0x34b0bcb5,0x391c0cb3,0x4ed8aa4a,0x5b9cca4f,0x682e6ff3,
268 0x748f82ee,0x78a5636f,0x84c87814,0x8cc70208,0x90befffa,0xa4506ceb,0xbef9a3f7,0xc67178f2,
269 };
270 u32 w[64];
271 for (int i = 0; i < 16; i++)
272 w[i] = (block[i * 4] << 24) | (block[i * 4 + 1] << 16) | (block[i * 4 + 2] << 8) | block[i * 4 + 3];
273 for (int i = 16; i < 64; i++) {
274 u32 s0 = rotr32(w[i - 15], 7) ^ rotr32(w[i - 15], 18) ^ (w[i - 15] >> 3);
275 u32 s1 = rotr32(w[i - 2], 17) ^ rotr32(w[i - 2], 19) ^ (w[i - 2] >> 10);
276 w[i] = w[i - 16] + s0 + w[i - 7] + s1;
277 }
278 u32 a = h[0], b = h[1], c2 = h[2], d2 = h[3], e = h[4], f = h[5], g = h[6], hh = h[7];
279 for (int i = 0; i < 64; i++) {
280 u32 S1 = rotr32(e, 6) ^ rotr32(e, 11) ^ rotr32(e, 25);
281 u32 ch = (e & f) ^ (~e & g);
282 u32 t1 = hh + S1 + ch + K[i] + w[i];
283 u32 S0 = rotr32(a, 2) ^ rotr32(a, 13) ^ rotr32(a, 22);
284 u32 maj = (a & b) ^ (a & c2) ^ (b & c2);
285 u32 t2 = S0 + maj;
286 hh = g; g = f; f = e; e = d2 + t1;
287 d2 = c2; c2 = b; b = a; a = t1 + t2;
288 }
289 h[0] += a; h[1] += b; h[2] += c2; h[3] += d2;
290 h[4] += e; h[5] += f; h[6] += g; h[7] += hh;
291 }
292
293 static u64 run_sha256chain(u64 words) {
294 u32 h[8] = {0x6a09e667,0xbb67ae85,0x3c6ef372,0xa54ff53a,0x510e527f,0x9b05688c,0x1f83d9ab,0x5be0cd19};
295 u8 block[64] = {0};
296 for (u64 i = 0; i < words; i++) {
297 u64 w0 = prng_next();
298 memcpy(block, &w0, 8);
299 sha256_compress(h, block);
300 }
301 return h[0] ^ h[7];
302 }
303
304 int main(int argc, char** argv) {
305 u64 words = argc > 1 ? strtoull(argv[1], nullptr, 10) : (1ull << 24);
306 bool verbose = argc > 2 && strcmp(argv[2], "--verbose") == 0;
307 prng_state = 0x9e3779b97f4a7c15ull;
308
309 u64 d1 = 0x9e3779b97f4a7c15ull | 1;
310 u64 d4[4] = {0x243f6a8885a308d3 | 1, 0x13198a2e03707344, 0xa4093822299f31d0, 0x082efa98ec4e6c89 | (1ull << 63)};
311
312 double secs[6];
313 u64 res[6];
314 const char* names[6] = {
315 "M=1 div hw ", "M=1 barrett ", "M=4 knuth div ",
316 "M=4 knuth barr", "x^2+c mod 2^128", "sha256 chain ",
317 };
318 const u64 seed = 0x9e3779b97f4a7c15ull;
319 auto run = [&](int idx, u64 r) { res[idx] = r; };
320
321 prng_state = seed;
322 auto t0 = std::chrono::steady_clock::now();
323 run(0, run_div1(d1, words, verbose));
324 auto t1 = std::chrono::steady_clock::now();
325 prng_state = seed;
326 run(1, run_div1_barrett(d1, words, verbose));
327 auto t2 = std::chrono::steady_clock::now();
328 prng_state = seed;
329 run(2, run_knuth4_div(d4, words, verbose));
330 auto t3 = std::chrono::steady_clock::now();
331 prng_state = seed;
332 run(3, run_knuth4_barrett(d4, words, verbose));
333 auto t4 = std::chrono::steady_clock::now();
334 prng_state = seed;
335 run(4, run_quad128(d1, words));
336 auto t5 = std::chrono::steady_clock::now();
337 prng_state = seed;
338 run(5, run_sha256chain(words));
339 auto t6 = std::chrono::steady_clock::now();
340
341 secs[0] = std::chrono::duration<double>(t1 - t0).count();
342 secs[1] = std::chrono::duration<double>(t2 - t1).count();
343 secs[2] = std::chrono::duration<double>(t3 - t2).count();
344 secs[3] = std::chrono::duration<double>(t4 - t3).count();
345 secs[4] = std::chrono::duration<double>(t5 - t4).count();
346 secs[5] = std::chrono::duration<double>(t6 - t5).count();
347
348 printf("words %llu\n", (unsigned long long)words);
349 for (int i = 0; i < 6; i++)
350 printf("%s: %8.3f s %6.2f ns/step (%llx)\n", names[i], secs[i],
351 secs[i] / words * 1e9, (unsigned long long)res[i]);
352 return 0;
353 }
354