Ottimizza il kernel GPU: ~2x throughput sulla Quadro P4000

Tre modifiche, tutte a parità di risultato (verificate end-to-end):

- Ammortamento della moltiplicazione scalare: prima scalar_mult_basic
  (l'operazione di gran lunga più costosa) veniva rieseguita a ogni
  batch, cioè 40 volte per kernel launch. Ora viene eseguita UNA sola
  volta per launch; il punto base P0 viene poi fatto avanzare di
  GPU_EC_BATCH_SIZE*G a ogni round con una semplice addizione EC
  affine, ripiegata nella batch inversion già presente (nuovo slot
  NEXT negli array Zj/prefix). La scalar mult è così ammortizzata su
  ~5120 chiavi invece di 128.

- Una sola inversione di campo per round invece di due: la
  normalizzazione del nuovo P0 condivide l'inversione di Montgomery
  del batch.

- Inversione modulare tramite la catena di addizione di libsecp256k1
  (255 quadrati + 15 moltiplicazioni) al posto del quadrato-e-moltiplica
  ingenuo (~256 + ~128). Rimosso il vecchio modinv_fermat e la costante
  FIELD_P_MINUS_2 non più usata.

- Rimosso l'array locale invZ[] (passaggio all'indietro fuso con la
  conversione affine + match), riducendo la memoria locale per-thread e
  migliorando l'occupancy.

Misurato nello stesso ambiente/scheda: 15.0M -> ~30M keys/sec.
Correttezza verificata end-to-end con privkey nota a offset 0, 50, 127,
128 (avanzamento del punto base), 200 e 5119 (40 avanzamenti
consecutivi): tutte trovate con la chiave privata corretta. Nessun
falso positivo su chiave pubblica valida casuale.

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.

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
+Su una Quadro P4000 (Pascal, 2016) si osservano ~30M keys/sec, contro i
-~300K-2M keys/sec della CPU — ma la moltiplicazione scalare su GPU non è
+~300K-2M keys/sec della CPU. La moltiplicazione scalare (costosa) viene
-ottimizzata quanto quella di libsecp256k1 (niente wNAF/tabelle precalcolate),
+eseguita una sola volta per kernel launch e poi il punto base viene fatto
-quindi il margine di ulteriore miglioramento è ancora ampio su GPU più
+avanzare con semplici addizioni EC ripiegate nella batch inversion, quindi il
-recenti.
+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

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.
+__device__ inline u256 sqrmod(const u256& a) { return mulmod(a, a); }
-// Costosa (~256 quadrati + fino a 256 moltiplicazioni), usata una sola volta
+
-// per batch (per invertire il prodotto totale nel trucco di Montgomery),
+// Inversione modulare a^(p-2) mod p tramite la catena di addizione di
-// quindi il costo è ammortizzato su GPU_EC_BATCH_SIZE chiavi.
+// libsecp256k1: 255 quadrati + 15 moltiplicazioni, contro i ~256 quadrati +
-__device__ inline u256 modinv_fermat(const u256& a) {
+// ~128 moltiplicazioni del quadrato-e-moltiplica ingenuo (piccolo teorema di
-    u256 result = { { 1, 0, 0, 0 } };
+// Fermat generico). Sfrutta la struttura di p-2 = 2^256 - 2^32 - 979, i cui
-    u256 base = a;
+// bit sono blocchi di uni; si costruiscono i valori a^(2^n - 1) e si combinano.
-    for (int limb = 3; limb >= 0; limb--) {
+// Usata una sola volta per round, quindi ammortizzata su GPU_EC_BATCH_SIZE chiavi.
-        uint64_t e = FIELD_P_MINUS_2.d[limb];
+__device__ inline u256 modinv(const u256& a) {
-        for (int bit = 63; bit >= 0; bit--) {
+    u256 x2, x3, x6, x9, x11, x22, x44, x88, x176, x220, x223, t1;
-            result = mulmod(result, result);
+    int j;
-            if ((e >> bit) & 1ULL) {
+
-                result = mulmod(result, base);
+    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);
-    return result;
+    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];
+    // Il punto base P0 = base_priv * G viene calcolato UNA sola volta con la
-    u256 invZ[GPU_EC_BATCH_MULT], prefix[GPU_EC_BATCH_MULT];
+    // moltiplicazione scalare (costosa), poi avanzato di GPU_EC_BATCH_SIZE*G a
-    bool valid[GPU_EC_BATCH_MULT];
+    // 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);
+        // Chiave 0 del round: P0 stesso (privkey = base_priv)
-        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);
-
         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);
 
-        if (valid_count > 0) {
+        // Prodotto prefisso dei soli Z validi
         int last = -1;
-            for (int i = 0; i < GPU_EC_BATCH_MULT; i++) {
+        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;
         }
 
-            u256 inv_tmp = modinv_fermat(prefix[last]);
+        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; } }
-                invZ[i] = (prev >= 0) ? mulmod(inv_tmp, prefix[prev]) : inv_tmp;
+            u256 invZ = (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++) {
+            u256 iv2 = mulmod(invZ, invZ);
-                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);
+            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;
+            }
 
             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;
+                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<TargetRecord> h_sorted;