pub const POINTWISE_CU: &str = "// Stage 4: H = A*B - C, elementwise, and the two encodings it is written in.\n//\n// Twin of crates/metal/src/shaders/pointwise.metal. Depends on bn254_fr.cuh being\n// concatenated ahead of it. No #include: NVRTC has no filesystem.\n//\n// ORDERING NOTE. In the Metal backend this file is concatenated BEFORE ntt.metal, because\n// the fused NTT epilogue calls g16_store_h and MSL needs the definition in scope. Here\n// snarkrs_gpu_kernels::unit_stages fixes the order as gather / ntt / pointwise, so ntt.cu\n// carries a forward declaration of g16_store_h and the definition stays here, next to the\n// stage it belongs to. Same translation unit, so the function still inlines into the NTT\n// epilogue. If you move this file\'s contents, move the declaration in ntt.cu with it.\n//\n// ============================================================================\n// NO DIVISION BY Z(coset). THIS IS NOT AN OMISSION.\n// ============================================================================\n//\n// snarkjs\' joinABC computes exactly a*b - c and feeds it straight to the H multiexp. The\n// division by the vanishing polynomial is folded into the section 9 bases at setup: they\n// are the odd Lagrange polynomials of the 2n domain, and P = A*B - C vanishes on the even\n// points (those are the constraint rows), so sum_i P(inc^(2i+1)) * hExps[i] already\n// equals [P(tau)]_1 = [H(tau) * Z(tau)]_1. Dividing here as well double-counts Z and\n// produces a proof that fails verification with nothing else to go on. The CPU backend\n// (snarkrs_groth16::cpu::CpuCircuit::compute_h) and the Metal backend make the same choice and\n// the same argument; this has been checked against real proofs, it is not a guess.\n//\n// ============================================================================\n// TWO OUTPUT BUFFERS, AND WHY BOTH\n// ============================================================================\n//\n// H leaves this stage in two encodings written in the same pass:\n//\n// h_mont Montgomery limbs (layout::PackedFr). This is what any further field\n// arithmetic wants, and it is what a host readback compares against the CPU\n// backend, since PackedFr::to_fr is the identity on arkworks\' internal form.\n// h_std standard form, the integer in [0, r) (layout::PackedScalar). This is what\n// stage 9\'s Pippenger MSM wants, because a window digit of a Montgomery\n// representative is a digit of a*R mod r, which is a different number. Getting\n// that backwards yields a proof wrong by a factor of R.\n//\n// The alternative is one buffer plus a conversion launch later, which costs an extra\n// full read and write of the domain (2n * 32 bytes) plus a launch, against n * 32 bytes\n// of extra writes here. On a stage that is already memory bound the second write is the\n// cheaper of the two, and it removes a coordination point between this stage and the MSM\n// stage entirely. On a discrete card there is a second reason: h_std is the only buffer\n// the MSM stage needs, so writing it here is what lets H stay device resident\n// (snarkrs_groth16::HPoly::Device) and never cross PCIe at all.\n\n#ifndef G16_POINTWISE_CU\n#define G16_POINTWISE_CU\n\n// Writes one H coefficient in both encodings. `v` is Montgomery form.\n//\n// fr_from_mont is a Montgomery multiply by the integer 1, i.e. a bare Montgomery\n// reduction, so h_std costs one multiply per element and no branch.\n//\n// Declared in ntt.cu, which is concatenated ahead of this file; the signature there must\n// stay character-for-character the same as this one.\n__device__ __forceinline__ void g16_store_h(Fr* h_mont, Fr* h_std, u32 i, Fr v) {\n h_mont[i] = v;\n h_std[i] = fr_from_mont(v);\n}\n\n// Standalone stage 4. The production path does not launch this: the same arithmetic is\n// fused into the store epilogue of the final NTT batch (see G16_STORE_JOIN in ntt.cu),\n// which saves a whole read of A, B and C and a whole write of C. This kernel is kept\n// because it is the unfused reference the fused path is checked against, and because a\n// caller that wants the three coset vectors materialised (a debugging dump, say) needs a\n// way to finish without them being consumed in flight.\n//\n// LAUNCH: grid_dim.x = ceil(n / block_dim.x), block_dim.x free (256). One thread does one\n// element; a block is 256 independent elements with no sharing, no barrier and no shared\n// memory.\nextern \"C\" __global__ void g16_h_join(\n const Fr* a,\n const Fr* b,\n const Fr* c,\n Fr* h_mont,\n Fr* h_std,\n const u32 n)\n{\n u32 gid = blockIdx.x * blockDim.x + threadIdx.x;\n // Load bearing on CUDA, unlike on Metal: the grid is a whole number of blocks, so the\n // tail block runs threads past n and without this they would write past h_mont.\n if (gid >= n) {\n return;\n }\n g16_store_h(h_mont, h_std, gid, fr_sub(fr_mul(a[gid], b[gid]), c[gid]));\n}\n\n// Montgomery limbs to standard limbs, element-wise. Stage 0.5 of the witness path: the\n// witness is uploaded once, in Montgomery form for the gather, and this converts it in\n// place on the device into the standard-form copy the MSM digit decomposition needs, so\n// the MSM stage never packs or uploads the witness a second time (that used to be a\n// second full PCIe transfer plus a host-side Montgomery reduction per element, per\n// proof). One thread per element, no sharing, tail-guarded like everything else.\n// Byte-identical twin of fr_mont_to_std in msm.cu, and the duplication is forced by the unit\n// split: unit_stages does not carry msm.cu. A change to fr_from_mont\'s callers is made in\n// both.\nextern \"C\" __global__ void g16_mont_to_std(\n const Fr* in,\n Fr* out,\n const u32 n)\n{\n u32 gid = blockIdx.x * blockDim.x + threadIdx.x;\n if (gid >= n) {\n return;\n }\n out[gid] = fr_from_mont(in[gid]);\n}\n\n#endif // G16_POINTWISE_CU\n";Expand description
Stage 4: H = A*B - C, fused into the tail of the last transform where possible.