diff --git a/CLAUDE.md b/CLAUDE.md index 5b668d2..4d7b388 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -63,7 +63,8 @@ There is no test suite in this repo. ### `bruteforce/p2pk_bruteforce_gpu.cu` (GPU, CUDA, single file) - Same algorithmic strategy as the CPU version (Jacobian batch add + Montgomery batch inversion), but GMP doesn't run on-device, so 256-bit field arithmetic mod the secp256k1 prime is hand-written (`u256` = 4×uint64 limbs, schoolbook multiply via `unsigned __int128`, fast reduction using `2^256 ≡ 2^32+977 (mod p)`). Verified bit-for-bit against `secp256k1_ec_pubkey_create` before trusting it (see conversation/test methodology — no separate test file is checked in). - Each CUDA thread is an independent search lane (own random starting privkey, own local batch state), same as a CPU thread — just thousands of them instead of ~11. Per-thread batch size is `GPU_EC_BATCH_SIZE` (128), smaller than the CPU's 256, to keep per-thread local-memory footprint (Jacobian batch arrays) bounded across tens of thousands of concurrent threads. -- Scalar multiplication on-device is plain double-and-add (no wNAF/windowing/precomputed tables like libsecp256k1 has on CPU) — the main further-optimization opportunity if more speed is needed. +- Scalar multiplication (`scalar_mult_basic`, plain double-and-add) runs **once per kernel launch**, not once per batch: the base point P0 is then advanced by `GPU_EC_BATCH_SIZE·G` each round via a cheap affine EC add folded into the batch inversion (the `NEXT` slot in the Zj/prefix arrays), so the expensive scalar mult is amortized over all `GPU_OUTER_ITERS_PER_LAUNCH × GPU_EC_BATCH_SIZE` keys. There is exactly one field inversion per round (Montgomery batch trick), done with the libsecp256k1 addition-chain `modinv` (255 sqr + 15 mul), not naive square-and-multiply. +- Remaining optimization headroom: double-and-add scalar mult (no wNAF/windowing), schoolbook `mulmod`, and occupancy tuning (grid size / batch size vs per-thread local memory). - Host side still uses libsecp256k1 (CPU) only for one-time setup: precomputing G-multiples and validating/loading target keys into a host-sorted array + Bloom filter, both uploaded once to device memory. Device-side matching = Bloom filter check + binary search over the sorted target array. - `NVCC_ARCH` in the Makefile defaults to `sm_61` (Pascal) — must match the actual GPU (`make gpu-info` shows compute capability via `nvidia-smi`, then pick the matching `sm_XX`). - On WSL2, the NVIDIA driver comes from the Windows host (visible as `nvidia-smi` working out of the box) — only the CUDA Toolkit needs installing inside WSL, never a driver. diff --git a/bruteforce/README.md b/bruteforce/README.md index 5215d26..715309a 100644 --- a/bruteforce/README.md +++ b/bruteforce/README.md @@ -68,11 +68,12 @@ make gpu # compila (default: sm_61 / Pascal — sovrascrivi con NVCC_ARC ./p2pk_bruteforce_gpu ``` -Su una Quadro P4000 (Pascal, 2016) si osservano ~15M keys/sec, contro i -~300K-2M keys/sec della CPU — ma la moltiplicazione scalare su GPU non è -ottimizzata quanto quella di libsecp256k1 (niente wNAF/tabelle precalcolate), -quindi il margine di ulteriore miglioramento è ancora ampio su GPU più -recenti. +Su una Quadro P4000 (Pascal, 2016) si osservano ~30M keys/sec, contro i +~300K-2M keys/sec della CPU. La moltiplicazione scalare (costosa) viene +eseguita una sola volta per kernel launch e poi il punto base viene fatto +avanzare con semplici addizioni EC ripiegate nella batch inversion, quindi il +grosso del tempo è ora nell'aritmetica di campo; su GPU più recenti il +margine di miglioramento è ancora ampio. > **Nota (WSL2)**: il driver NVIDIA arriva dal sistema Windows host (basta > che `nvidia-smi` funzioni già) — dentro WSL serve installare solo il CUDA diff --git a/bruteforce/p2pk_bruteforce_gpu.cu b/bruteforce/p2pk_bruteforce_gpu.cu index 5baec18..e397853 100644 --- a/bruteforce/p2pk_bruteforce_gpu.cu +++ b/bruteforce/p2pk_bruteforce_gpu.cu @@ -54,9 +54,6 @@ struct u256 { uint64_t d[4]; }; __constant__ u256 FIELD_P = {{ 0xFFFFFFFEFFFFFC2FULL, 0xFFFFFFFFFFFFFFFFULL, 0xFFFFFFFFFFFFFFFFULL, 0xFFFFFFFFFFFFFFFFULL }}; -// p - 2, usato per l'inversione modulare via piccolo teorema di Fermat (a^(p-2) mod p) -__constant__ u256 FIELD_P_MINUS_2 = {{ 0xFFFFFFFEFFFFFC2DULL, 0xFFFFFFFFFFFFFFFFULL, - 0xFFFFFFFFFFFFFFFFULL, 0xFFFFFFFFFFFFFFFFULL }}; // 2^256 mod p = 2^32 + 977, costante per la riduzione veloce (forma speciale del primo) #define FIELD_K 0x1000003D1ULL @@ -188,23 +185,35 @@ __device__ inline u256 mulmod(const u256& a, const u256& b) { return reduce512(t); } -// Inversione modulare via piccolo teorema di Fermat: a^(p-2) mod p. -// Costosa (~256 quadrati + fino a 256 moltiplicazioni), usata una sola volta -// per batch (per invertire il prodotto totale nel trucco di Montgomery), -// quindi il costo è ammortizzato su GPU_EC_BATCH_SIZE chiavi. -__device__ inline u256 modinv_fermat(const u256& a) { - u256 result = { { 1, 0, 0, 0 } }; - u256 base = a; - for (int limb = 3; limb >= 0; limb--) { - uint64_t e = FIELD_P_MINUS_2.d[limb]; - for (int bit = 63; bit >= 0; bit--) { - result = mulmod(result, result); - if ((e >> bit) & 1ULL) { - result = mulmod(result, base); - } - } - } - return result; +__device__ inline u256 sqrmod(const u256& a) { return mulmod(a, a); } + +// Inversione modulare a^(p-2) mod p tramite la catena di addizione di +// libsecp256k1: 255 quadrati + 15 moltiplicazioni, contro i ~256 quadrati + +// ~128 moltiplicazioni del quadrato-e-moltiplica ingenuo (piccolo teorema di +// Fermat generico). Sfrutta la struttura di p-2 = 2^256 - 2^32 - 979, i cui +// bit sono blocchi di uni; si costruiscono i valori a^(2^n - 1) e si combinano. +// Usata una sola volta per round, quindi ammortizzata su GPU_EC_BATCH_SIZE chiavi. +__device__ inline u256 modinv(const u256& a) { + u256 x2, x3, x6, x9, x11, x22, x44, x88, x176, x220, x223, t1; + int j; + + x2 = sqrmod(a); x2 = mulmod(x2, a); + x3 = sqrmod(x2); x3 = mulmod(x3, a); + x6 = x3; for (j = 0; j < 3; j++) x6 = sqrmod(x6); x6 = mulmod(x6, x3); + x9 = x6; for (j = 0; j < 3; j++) x9 = sqrmod(x9); x9 = mulmod(x9, x3); + x11 = x9; for (j = 0; j < 2; j++) x11 = sqrmod(x11); x11 = mulmod(x11, x2); + x22 = x11; for (j = 0; j < 11; j++) x22 = sqrmod(x22); x22 = mulmod(x22, x11); + x44 = x22; for (j = 0; j < 22; j++) x44 = sqrmod(x44); x44 = mulmod(x44, x22); + x88 = x44; for (j = 0; j < 44; j++) x88 = sqrmod(x88); x88 = mulmod(x88, x44); + x176 = x88; for (j = 0; j < 88; j++) x176 = sqrmod(x176); x176 = mulmod(x176, x88); + x220 = x176; for (j = 0; j < 44; j++) x220 = sqrmod(x220); x220 = mulmod(x220, x44); + x223 = x220; for (j = 0; j < 3; j++) x223 = sqrmod(x223); x223 = mulmod(x223, x3); + + t1 = x223; for (j = 0; j < 23; j++) t1 = sqrmod(t1); t1 = mulmod(t1, x22); + for (j = 0; j < 5; j++) t1 = sqrmod(t1); t1 = mulmod(t1, a); + for (j = 0; j < 3; j++) t1 = sqrmod(t1); t1 = mulmod(t1, x2); + for (j = 0; j < 2; j++) t1 = sqrmod(t1); t1 = mulmod(t1, a); + return t1; } // ============================================================================ @@ -393,8 +402,10 @@ __device__ inline int check_match(const uint64_t* bloom, uint64_t bloom_mask, __constant__ u256 D_GX; __constant__ u256 D_GY; -__constant__ u256 D_PRECOMP_GX[GPU_EC_BATCH_MULT]; +__constant__ u256 D_PRECOMP_GX[GPU_EC_BATCH_MULT]; // 1G, 2G, ..., (GPU_EC_BATCH_MULT)G __constant__ u256 D_PRECOMP_GY[GPU_EC_BATCH_MULT]; +__constant__ u256 D_STEP_GX; // (GPU_EC_BATCH_SIZE)G, passo per avanzare il punto base P0 +__constant__ u256 D_STEP_GY; __global__ void search_kernel(const uint64_t* bloom_bits, uint64_t bloom_mask, const TargetRecord* sorted_targets, int n_targets, @@ -403,80 +414,106 @@ __global__ void search_kernel(const uint64_t* bloom_bits, uint64_t bloom_mask, int* d_found_flag, uint8_t* d_found_privkey, int* d_found_target_idx) { uint32_t tid = blockIdx.x * blockDim.x + threadIdx.x; - u256 privkey = random_start_privkey(global_seed, launch_id, tid); + // privkey di base del punto P0 corrente (avanza di GPU_EC_BATCH_SIZE per round) + u256 base_priv = random_start_privkey(global_seed, launch_id, tid); unsigned long long local_attempts = 0; - u256 Xj[GPU_EC_BATCH_MULT], Yj[GPU_EC_BATCH_MULT], Zj[GPU_EC_BATCH_MULT]; - u256 invZ[GPU_EC_BATCH_MULT], prefix[GPU_EC_BATCH_MULT]; - bool valid[GPU_EC_BATCH_MULT]; + // Il punto base P0 = base_priv * G viene calcolato UNA sola volta con la + // moltiplicazione scalare (costosa), poi avanzato di GPU_EC_BATCH_SIZE*G a + // ogni round tramite una semplice addizione EC ripiegata nella batch + // inversion — così l'unica scalar mult è ammortizzata sull'intero launch. + Jacobian P0j = scalar_mult_basic(base_priv, D_GX, D_GY); + u256 zinv = modinv(P0j.Z); + u256 zinv2 = mulmod(zinv, zinv); + u256 x0 = mulmod(P0j.X, zinv2); + u256 y0 = mulmod(P0j.Y, mulmod(zinv2, zinv)); + + // Indici 0..GPU_EC_BATCH_MULT-1: punti del batch P0+1G..P0+(MULT)G. + // Indice GPU_EC_BATCH_MULT: il prossimo P0 = P0 + GPU_EC_BATCH_SIZE*G, + // incluso nella stessa batch inversion per non pagare una seconda inversione. + u256 Xj[GPU_EC_BATCH_SIZE], Yj[GPU_EC_BATCH_SIZE], Zj[GPU_EC_BATCH_SIZE]; + u256 prefix[GPU_EC_BATCH_SIZE]; + bool valid[GPU_EC_BATCH_SIZE]; + const int NEXT = GPU_EC_BATCH_MULT; // ultimo slot = prossimo punto base for (int outer = 0; outer < GPU_OUTER_ITERS_PER_LAUNCH; outer++) { if (*d_found_flag) break; - Jacobian P0j = scalar_mult_basic(privkey, D_GX, D_GY); - u256 zinv = modinv_fermat(P0j.Z); - u256 zinv2 = mulmod(zinv, zinv); - u256 zinv3 = mulmod(zinv2, zinv); - u256 x0 = mulmod(P0j.X, zinv2); - u256 y0 = mulmod(P0j.Y, zinv3); - + // Chiave 0 del round: P0 stesso (privkey = base_priv) uint8_t key0[64]; u256_to_be32(x0, key0); u256_to_be32(y0, key0 + 32); int m = check_match(bloom_bits, bloom_mask, sorted_targets, n_targets, key0); if (m >= 0 && atomicCAS(d_found_flag, 0, 1) == 0) { - u256_to_be32(privkey, d_found_privkey); + u256_to_be32(base_priv, d_found_privkey); *d_found_target_idx = m; } - int valid_count = 0; + // Batch: P0 + iG (affine+affine con Z1=1) + il prossimo P0 nell'ultimo slot for (int i = 0; i < GPU_EC_BATCH_MULT; i++) { jacobian_add_affine_z1(x0, y0, D_PRECOMP_GX[i], D_PRECOMP_GY[i], Xj[i], Yj[i], Zj[i]); valid[i] = !(Zj[i].d[0] == 0 && Zj[i].d[1] == 0 && Zj[i].d[2] == 0 && Zj[i].d[3] == 0); - if (valid[i]) valid_count++; + } + jacobian_add_affine_z1(x0, y0, D_STEP_GX, D_STEP_GY, Xj[NEXT], Yj[NEXT], Zj[NEXT]); + valid[NEXT] = !(Zj[NEXT].d[0] == 0 && Zj[NEXT].d[1] == 0 && Zj[NEXT].d[2] == 0 && Zj[NEXT].d[3] == 0); + + // Prodotto prefisso dei soli Z validi + int last = -1; + for (int i = 0; i < GPU_EC_BATCH_SIZE; i++) { + if (!valid[i]) continue; + prefix[i] = (last < 0) ? Zj[i] : mulmod(prefix[last], Zj[i]); + last = i; } - if (valid_count > 0) { - int last = -1; - for (int i = 0; i < GPU_EC_BATCH_MULT; i++) { - if (!valid[i]) continue; - prefix[i] = (last < 0) ? Zj[i] : mulmod(prefix[last], Zj[i]); - last = i; + if (last < 0) { + // Caso impossibile in pratica (tutti Z=0): re-inizializza il base point + base_priv = random_start_privkey(global_seed, launch_id ^ 0x5A5A5A5A, tid + outer + 1); + P0j = scalar_mult_basic(base_priv, D_GX, D_GY); + zinv = modinv(P0j.Z); zinv2 = mulmod(zinv, zinv); + x0 = mulmod(P0j.X, zinv2); y0 = mulmod(P0j.Y, mulmod(zinv2, zinv)); + local_attempts += GPU_EC_BATCH_SIZE; + continue; + } + + // Un'unica inversione per l'intero round (trucco di Montgomery) + u256 inv_tmp = modinv(prefix[last]); + + // Passaggio all'indietro: calcola invZ[i] e consumalo subito + // (conversione affine + match), senza array invZ[] separato. + u256 next_x0 = x0, next_y0 = y0; // fallback se lo slot NEXT non fosse valido + for (int i = last; i >= 0; i--) { + if (!valid[i]) continue; + int prev = -1; + for (int j = i - 1; j >= 0; j--) { if (valid[j]) { prev = j; break; } } + u256 invZ = (prev >= 0) ? mulmod(inv_tmp, prefix[prev]) : inv_tmp; + inv_tmp = mulmod(inv_tmp, Zj[i]); + + u256 iv2 = mulmod(invZ, invZ); + u256 ax = mulmod(Xj[i], iv2); + u256 ay = mulmod(Yj[i], mulmod(iv2, invZ)); + + if (i == NEXT) { + // Non è una chiave da controllare: è il punto base del prossimo round + next_x0 = ax; next_y0 = ay; + continue; } - u256 inv_tmp = modinv_fermat(prefix[last]); - - for (int i = last; i >= 0; i--) { - if (!valid[i]) continue; - int prev = -1; - for (int j = i - 1; j >= 0; j--) { if (valid[j]) { prev = j; break; } } - invZ[i] = (prev >= 0) ? mulmod(inv_tmp, prefix[prev]) : inv_tmp; - inv_tmp = mulmod(inv_tmp, Zj[i]); - } - - for (int i = 0; i < GPU_EC_BATCH_MULT; i++) { - if (!valid[i]) continue; - u256 iv2 = mulmod(invZ[i], invZ[i]); - u256 iv3 = mulmod(iv2, invZ[i]); - u256 ax = mulmod(Xj[i], iv2); - u256 ay = mulmod(Yj[i], iv3); - - uint8_t key[64]; - u256_to_be32(ax, key); - u256_to_be32(ay, key + 32); - - int mi = check_match(bloom_bits, bloom_mask, sorted_targets, n_targets, key); - if (mi >= 0 && atomicCAS(d_found_flag, 0, 1) == 0) { - u256 found = privkey; - add_small_u256(found, (uint32_t)(i + 1)); - u256_to_be32(found, d_found_privkey); - *d_found_target_idx = mi; - } + uint8_t key[64]; + u256_to_be32(ax, key); + u256_to_be32(ay, key + 32); + int mi = check_match(bloom_bits, bloom_mask, sorted_targets, n_targets, key); + if (mi >= 0 && atomicCAS(d_found_flag, 0, 1) == 0) { + u256 found = base_priv; + add_small_u256(found, (uint32_t)(i + 1)); + u256_to_be32(found, d_found_privkey); + *d_found_target_idx = mi; } } local_attempts += GPU_EC_BATCH_SIZE; - add_small_u256(privkey, GPU_EC_BATCH_SIZE); + add_small_u256(base_priv, GPU_EC_BATCH_SIZE); + x0 = next_x0; + y0 = next_y0; } atomicAdd(d_total_attempts, local_attempts); @@ -571,6 +608,21 @@ int main(int argc, char** argv) { CUDA_CHECK(cudaMemcpyToSymbol(D_PRECOMP_GY, h_pgy.data(), sizeof(u256) * GPU_EC_BATCH_MULT)); } + // Passo di avanzamento del punto base: GPU_EC_BATCH_SIZE * G + { + uint8_t pk[32] = {0}; + uint32_t v = (uint32_t)GPU_EC_BATCH_SIZE; + pk[31] = v & 0xFF; pk[30] = (v >> 8) & 0xFF; pk[29] = (v >> 16) & 0xFF; + secp256k1_pubkey pub; + if (!secp256k1_ec_pubkey_create(ctx, &pub, pk)) { fprintf(stderr, "[ERROR] precompute step\n"); return 1; } + unsigned char buf[65]; size_t outlen = 65; + secp256k1_ec_pubkey_serialize(ctx, buf, &outlen, &pub, SECP256K1_EC_UNCOMPRESSED); + u256 sx = be32_to_u256(buf + 1); + u256 sy = be32_to_u256(buf + 33); + CUDA_CHECK(cudaMemcpyToSymbol(D_STEP_GX, &sx, sizeof(u256))); + CUDA_CHECK(cudaMemcpyToSymbol(D_STEP_GY, &sy, sizeof(u256))); + } + // Caricamento target keys (stesso formato/euristica della versione CPU) printf("[+] Caricamento chiavi target da %s...\n", target_file); std::vector h_sorted;