From a40fce0f4fff22f86b2a8e3dd33637245db72e18 Mon Sep 17 00:00:00 2001 From: Oleksandr Brezhniev Date: Fri, 10 Jul 2026 13:41:14 +0100 Subject: [PATCH 1/6] perf: scalar-size partitioned MSM + wall-clock-aware window selection Partition scalars by significant bits before the bucket method: zeros are dropped, ones are accumulated with plain mixed additions, scalars up to 64 bits run Pippenger over a single word (~4 windows instead of ~16), and only full-width scalars pay full window costs. Fast paths skip the gather when one class covers (nearly) all points. Window size is now chosen by a wall-clock cost model (chunks execute in waves of nThreads) instead of minimizing total additions, with a total-work term as tie-break, so large MSMs no longer leave cores idle when the chunk count is below the thread count. Also: chunk-major slicedScalars layout (unit-stride digit reads in the bucket fill loop), per-partition nBits with +2 carry headroom (fixes a latent dropped carry for scalars using all scalarSize*8 bits), and a getBucketIndex guard against shifting past the scalar. Measured on 20 cores, G1 n=2^20: uniform scalars -14%, all-binary scalars -87%, realistic iden3 witness -14%. G2 unchanged. End to end: sha256_test (2M binary wires) proof -23%, small identity circuits unchanged (zkey load dominated). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01JpMJeCdh75Shs3hBRFFLkq --- c/msm.cpp | 181 ++++++++++++++++++++++++++++++++++++++++++++++++++++-- c/msm.hpp | 63 ++++++++++++++++--- 2 files changed, 230 insertions(+), 14 deletions(-) diff --git a/c/msm.cpp b/c/msm.cpp index 204e2cd..003ed73 100644 --- a/c/msm.cpp +++ b/c/msm.cpp @@ -1,4 +1,6 @@ #include +#include +#include #include "msm.hpp" #include "misc.hpp" @@ -9,6 +11,177 @@ void MSM::run(typename Curve::Point &r, uint64_t _scalarSize, uint64_t _n, uint64_t _nThreads) +{ + if (_n == 0) { + g.copy(r, g.zero()); + return; + } + if (_n == 1) { + g.mulByScalar(r, _bases[0], _scalars, _scalarSize); + return; + } + if (_scalarSize < 8) { + runPartition(r, _bases, _scalars, _scalarSize, _scalarSize*8, _n); + return; + } + + ThreadPool &threadPool = ThreadPool::defaultPool(); + + scalars = _scalars; + scalarSize = _scalarSize; + + const uint64_t nBlocks = std::min(threadPool.getThreadCount()*4, _n); + const uint64_t blockSize = (_n + nBlocks - 1) / nBlocks; + + enum ScalarClass : uint8_t { CLS_ZERO = 0, CLS_ONE = 1, CLS_SMALL = 2, CLS_BIG = 3 }; + + std::unique_ptr classes(new uint8_t[_n]); + std::unique_ptr blockCounts(new uint64_t[nBlocks*3]); + std::unique_ptr blockMaxBits(new uint64_t[nBlocks*2]); + + threadPool.parallelFor(0, nBlocks, [&] (int begin, int end, int numThread) { + for (int b = begin; b < end; b++) { + const uint64_t i0 = (uint64_t)b*blockSize; + const uint64_t i1 = std::min(i0 + blockSize, _n); + uint64_t nSmall = 0, nBig = 0, nOnes = 0; + uint64_t maxSmall = 0, maxBig = 0; + + for (uint64_t i = i0; i < i1; i++) { + const uint64_t bits = significantBits(_scalars + i*_scalarSize); + uint8_t cls; + + if (bits == 0) { + cls = CLS_ZERO; + } else if (bits == 1) { + cls = CLS_ONE; + nOnes++; + } else if (bits <= SMALL_SCALAR_BITS) { + cls = CLS_SMALL; + nSmall++; + if (bits > maxSmall) maxSmall = bits; + } else { + cls = CLS_BIG; + nBig++; + if (bits > maxBig) maxBig = bits; + } + classes[i] = cls; + } + blockCounts[b*3] = nSmall; + blockCounts[b*3+1] = nBig; + blockCounts[b*3+2] = nOnes; + blockMaxBits[b*2] = maxSmall; + blockMaxBits[b*2+1] = maxBig; + } + }); + + uint64_t nSmall = 0, nBig = 0, nOnes = 0; + uint64_t maxSmallBits = 0, maxBigBits = 0; + + std::unique_ptr blockOffsets(new uint64_t[nBlocks*2]); + + for (uint64_t b = 0; b < nBlocks; b++) { + blockOffsets[b*2] = nSmall; + blockOffsets[b*2+1] = nBig; + nSmall += blockCounts[b*3]; + nBig += blockCounts[b*3+1]; + nOnes += blockCounts[b*3+2]; + if (blockMaxBits[b*2] > maxSmallBits) maxSmallBits = blockMaxBits[b*2]; + if (blockMaxBits[b*2+1] > maxBigBits) maxBigBits = blockMaxBits[b*2+1]; + } + + const uint64_t overallMaxBits = std::max(maxBigBits, std::max(maxSmallBits, (uint64_t)(nOnes ? 1 : 0))); + + // When almost every scalar is full width (e.g. the H MSM, whose scalars + // are uniform field elements) partitioning saves nothing: run the whole + // input in place instead of paying the gather. + if (nBig >= _n - _n/16) { + runPartition(r, _bases, _scalars, _scalarSize, overallMaxBits + 2, _n); + return; + } + + // All scalars fit in 64 bits: gather only the scalars, bases stay in place. + if (nSmall == _n) { + std::unique_ptr smallScalars(new uint64_t[_n]); + + threadPool.parallelFor(0, _n, [&] (int begin, int end, int numThread) { + for (int i = begin; i < end; i++) { + std::memcpy(&smallScalars[i], _scalars + (uint64_t)i*_scalarSize, sizeof(uint64_t)); + } + }); + + runPartition(r, _bases, (uint8_t *)smallScalars.get(), sizeof(uint64_t), maxSmallBits + 2, _n); + return; + } + + std::unique_ptr smallScalars(nSmall ? new uint64_t[nSmall] : nullptr); + std::unique_ptr smallBases(nSmall ? new typename Curve::PointAffine[nSmall] : nullptr); + std::unique_ptr bigScalars(nBig ? new uint8_t[nBig*_scalarSize] : nullptr); + std::unique_ptr bigBases(nBig ? new typename Curve::PointAffine[nBig] : nullptr); + std::unique_ptr onesAcc(new typename Curve::Point[nBlocks]); + + threadPool.parallelFor(0, nBlocks, [&] (int begin, int end, int numThread) { + for (int b = begin; b < end; b++) { + const uint64_t i0 = (uint64_t)b*blockSize; + const uint64_t i1 = std::min(i0 + blockSize, _n); + uint64_t smallCur = blockOffsets[b*2]; + uint64_t bigCur = blockOffsets[b*2+1]; + + g.copy(onesAcc[b], g.zero()); + + for (uint64_t i = i0; i < i1; i++) { + switch (classes[i]) { + case CLS_ONE: + g.add(onesAcc[b], onesAcc[b], _bases[i]); + break; + case CLS_SMALL: + std::memcpy(&smallScalars[smallCur], _scalars + i*_scalarSize, sizeof(uint64_t)); + smallBases[smallCur] = _bases[i]; + smallCur++; + break; + case CLS_BIG: + std::memcpy(&bigScalars[bigCur*_scalarSize], _scalars + i*_scalarSize, _scalarSize); + bigBases[bigCur] = _bases[i]; + bigCur++; + break; + default: + break; + } + } + } + }); + + typename Curve::Point acc; + + g.copy(acc, onesAcc[0]); + for (uint64_t b = 1; b < nBlocks; b++) { + g.add(acc, acc, onesAcc[b]); + } + + if (nSmall > 0) { + typename Curve::Point rSmall; + + runPartition(rSmall, smallBases.get(), (uint8_t *)smallScalars.get(), + sizeof(uint64_t), maxSmallBits + 2, nSmall); + g.add(acc, acc, rSmall); + } + if (nBig > 0) { + typename Curve::Point rBig; + + runPartition(rBig, bigBases.get(), bigScalars.get(), + _scalarSize, maxBigBits + 2, nBig); + g.add(acc, acc, rBig); + } + + g.copy(r, acc); +} + +template +void MSM::runPartition(typename Curve::Point &r, + typename Curve::PointAffine *_bases, + uint8_t* _scalars, + uint64_t _scalarSize, + uint64_t _nBits, + uint64_t _n) { ThreadPool &threadPool = ThreadPool::defaultPool(); @@ -21,7 +194,7 @@ void MSM::run(typename Curve::Point &r, #ifdef MSM_BITS_PER_CHUNK bitsPerChunk = MSM_BITS_PER_CHUNK; #else - bitsPerChunk = calcBitsPerChunk(nPoints, scalarSize); + bitsPerChunk = calcBitsPerChunk(nPoints, _nBits, nThreads); #endif if (nPoints == 0) { @@ -33,7 +206,7 @@ void MSM::run(typename Curve::Point &r, return; } - const uint64_t nChunks = calcChunkCount(scalarSize, bitsPerChunk); + const uint64_t nChunks = calcChunkCount(_nBits, bitsPerChunk); const uint64_t nBuckets = calcBucketCount(bitsPerChunk); const uint64_t matrixSize = nThreads * nBuckets; const uint64_t nSlices = nChunks*nPoints; @@ -57,7 +230,7 @@ void MSM::run(typename Curve::Point &r, carry = 0; } - slicedScalars[i*nChunks + j] = bucketIndex; + slicedScalars[j*nPoints + i] = bucketIndex; } } }); @@ -73,7 +246,7 @@ void MSM::run(typename Curve::Point &r, } for (int i = 0; i < nPoints; i++) { - const int bucketIndex = slicedScalars[i*nChunks + j]; + const int bucketIndex = slicedScalars[j*nPoints + i]; if (bucketIndex > 0) { g.add(buckets[bucketIndex-1], buckets[bucketIndex-1], _bases[i]); diff --git a/c/msm.hpp b/c/msm.hpp index bdfdb69..f5ffd19 100644 --- a/c/msm.hpp +++ b/c/msm.hpp @@ -8,34 +8,48 @@ class MSM { const uint64_t MIN_CHUNK_SIZE_BITS = 3; const uint64_t MAX_CHUNK_SIZE_BITS = 16; + // Scalars of at most this many significant bits go to the small + // partition, which runs the bucket method over a single 64-bit word + // and so pays ~4 windows instead of ~16. Scalars 0 and 1 are cheaper + // still: they need no scalar multiplication at all. + static const uint64_t SMALL_SCALAR_BITS = 64; + Curve &g; uint8_t *scalars; uint64_t scalarSize; uint64_t bitsPerChunk; private: - uint64_t calcAddsCount(uint64_t nPoints, uint64_t scalarSize, uint64_t bitsPerChunk) const { - return calcChunkCount(scalarSize, bitsPerChunk) - * (nPoints + ((uint64_t)1 << bitsPerChunk) + bitsPerChunk + 1); + // Estimated wall-clock cost in point additions. Chunks run in parallel, + // so the elapsed time is the per-chunk cost times the number of waves of + // nThreads chunks; the second term charges 1/8 of the total work so that + // among near-equal wall costs the one burning fewer total additions wins. + uint64_t calcCost(uint64_t nPoints, uint64_t nBits, uint64_t bitsPerChunk, + uint64_t nThreads) const { + const uint64_t chunkCost = nPoints + ((uint64_t)1 << bitsPerChunk) + bitsPerChunk + 1; + const uint64_t nChunks = calcChunkCount(nBits, bitsPerChunk); + const uint64_t waves = (nChunks + nThreads - 1) / nThreads; + + return waves*chunkCost + nChunks*chunkCost/(8*nThreads); } - uint64_t calcBitsPerChunk(uint64_t n, uint64_t scalarSize) const { + uint64_t calcBitsPerChunk(uint64_t n, uint64_t nBits, uint64_t nThreads) const { uint64_t bitsPerChunk = MIN_CHUNK_SIZE_BITS; - uint64_t minAdds = calcAddsCount(n, scalarSize, bitsPerChunk); + uint64_t minCost = calcCost(n, nBits, bitsPerChunk, nThreads); for (uint64_t k = MIN_CHUNK_SIZE_BITS + 1; k <= MAX_CHUNK_SIZE_BITS; k++) { - const uint64_t curAdds = calcAddsCount(n, scalarSize, k); + const uint64_t curCost = calcCost(n, nBits, k, nThreads); - if (curAdds < minAdds) { - minAdds = curAdds; + if (curCost < minCost) { + minCost = curCost; bitsPerChunk = k; } } return bitsPerChunk; } - uint64_t calcChunkCount(uint64_t scalarSize, uint64_t bitsPerChunk) const { - return ((scalarSize * 8 - 1 ) / bitsPerChunk) + 1; + uint64_t calcChunkCount(uint64_t nBits, uint64_t bitsPerChunk) const { + return ((nBits - 1) / bitsPerChunk) + 1; } uint64_t calcBucketCount(uint64_t bitsPerChunk) const { @@ -44,6 +58,11 @@ class MSM { uint64_t getBucketIndex(uint64_t scalarIdx, uint64_t chunkIdx) const { uint64_t bitStart = chunkIdx*bitsPerChunk; + + // Chunks past the scalar bytes exist only to absorb the signed-digit + // carry; their digit is zero. + if (bitStart >= scalarSize*8) return 0; + uint64_t byteStart = bitStart/8; uint64_t efectiveBitsPerChunk = bitsPerChunk; @@ -59,15 +78,39 @@ class MSM { return uint64_t(v); } + // Number of significant bits of a little-endian scalar; 0 for a zero scalar. + uint64_t significantBits(const uint8_t *scalar) const { + for (int64_t k = (int64_t)scalarSize - 1; k >= 0; k--) { + if (scalar[k]) { + return (uint64_t)k*8 + (32 - __builtin_clz((uint32_t)scalar[k])); + } + } + return 0; + } + public: MSM(Curve &_g): g(_g) {} + // Partitions the scalars by significant bits (0, 1, up to 64 bits, wider) + // and runs the bucket method separately per partition, so the mostly-0/1 + // scalars of a circom witness don't pay full-width window costs. void run(typename Curve::Point &r, typename Curve::PointAffine *_bases, uint8_t* _scalars, uint64_t _scalarSize, uint64_t _n, uint64_t _nThreads=0); + + // One bucket-method pass over all n points using only the lowest nBits + // bits of each scalar. nBits must exceed the largest scalar's significant + // bit count by at least 2, so a signed-digit carry can never propagate + // out of the top chunk. + void runPartition(typename Curve::Point &r, + typename Curve::PointAffine *bases, + uint8_t* scalars, + uint64_t scalarSize, + uint64_t nBits, + uint64_t n); }; #include "msm.cpp" From cb8e99c92470f7edc403260265b74a92402d96be Mon Sep 17 00:00:00 2001 From: Oleksandr Brezhniev Date: Fri, 10 Jul 2026 19:35:09 +0100 Subject: [PATCH 2/6] perf: task-based MSM with point-splitting for many-core scaling Split the MSM into prepare (classify/gather/recode), collectTasks (one closure per partition x slice x window over a per-thread bucket arena) and finish (reduction), so several MSMs -- even over different curves -- can execute their bucket work in one parallel region. run() keeps the old single-call behavior on top. The cost model now picks window size and a point-split factor jointly: a partition short on windows splits its points into slices, gaining parallelism at the work-optimal window size instead of shrinking the window (which inflates total additions). Measured on 20 cores, interleaved: u64 scalars -46%, witness-mix -24%, G2 witness-mix -27%; uniform full-width scalars unchanged. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01JpMJeCdh75Shs3hBRFFLkq --- c/msm.cpp | 381 ++++++++++++++++++++++++++++++++++++------------------ c/msm.hpp | 170 +++++++++++++++++------- 2 files changed, 383 insertions(+), 168 deletions(-) diff --git a/c/msm.cpp b/c/msm.cpp index 003ed73..8118976 100644 --- a/c/msm.cpp +++ b/c/msm.cpp @@ -5,31 +5,98 @@ #include "misc.hpp" template -void MSM::run(typename Curve::Point &r, - typename Curve::PointAffine *_bases, - uint8_t* _scalars, - uint64_t _scalarSize, - uint64_t _n, - uint64_t _nThreads) +void MSM::preparePartition(Partition &p, uint64_t nThreads) { + ThreadPool &threadPool = ThreadPool::defaultPool(); + +#ifdef MSM_BITS_PER_CHUNK + p.bitsPerChunk = MSM_BITS_PER_CHUNK; + p.nSlices = 1; +#else + calcChunkConfig(p.n, p.nBits, nThreads, p.bitsPerChunk, p.nSlices); +#endif + + p.nChunks = calcChunkCount(p.nBits, p.bitsPerChunk); + p.nBuckets = calcBucketCount(p.bitsPerChunk); + p.digits.reset(new int32_t[p.nChunks * p.n]); + p.partials.reset(new typename Curve::Point[p.nSlices * p.nChunks]); + + // recode context for getBucketIndex + scalars = p.scalars; + scalarSize = p.scalarSize; + bitsPerChunk = p.bitsPerChunk; + + const uint64_t nChunks = p.nChunks; + const uint64_t nBuckets = p.nBuckets; + const uint64_t nPoints = p.n; + int32_t *digits = p.digits.get(); + + threadPool.parallelFor(0, nPoints, [&, nChunks, nBuckets, nPoints] (int begin, int end, int numThread) { + + for (int i = begin; i < end; i++) { + int carry = 0; + + for (uint64_t j = 0; j < nChunks; j++) { + int bucketIndex = getBucketIndex(i, j) + carry; + + if (bucketIndex >= (int)nBuckets) { + bucketIndex -= nBuckets*2; + carry = 1; + } else { + carry = 0; + } + + digits[j*nPoints + i] = bucketIndex; + } + } + }); +} + +template +void MSM::prepare(typename Curve::PointAffine *_bases, + uint8_t *_scalars, + uint64_t _scalarSize, + uint64_t _n, + uint64_t parallelismShare) +{ + ThreadPool &threadPool = ThreadPool::defaultPool(); + + const uint64_t nThreads = parallelismShare ? parallelismShare + : threadPool.getThreadCount(); + + partitions.clear(); + partitions.reserve(2); + onesAcc.reset(); + nOnesBlocks = 0; + trivial = false; + prepared = true; + if (_n == 0) { - g.copy(r, g.zero()); + trivial = true; + g.copy(trivialResult, g.zero()); return; } if (_n == 1) { - g.mulByScalar(r, _bases[0], _scalars, _scalarSize); - return; - } - if (_scalarSize < 8) { - runPartition(r, _bases, _scalars, _scalarSize, _scalarSize*8, _n); + trivial = true; + g.mulByScalar(trivialResult, _bases[0], _scalars, _scalarSize); return; } - ThreadPool &threadPool = ThreadPool::defaultPool(); - scalars = _scalars; scalarSize = _scalarSize; + if (_scalarSize < 8) { + partitions.emplace_back(); + Partition &p = partitions.back(); + p.bases = _bases; + p.scalars = _scalars; + p.scalarSize = _scalarSize; + p.n = _n; + p.nBits = _scalarSize*8; + preparePartition(p, nThreads); + return; + } + const uint64_t nBlocks = std::min(threadPool.getThreadCount()*4, _n); const uint64_t blockSize = (_n + nBlocks - 1) / nBlocks; @@ -95,52 +162,93 @@ void MSM::run(typename Curve::Point &r, // are uniform field elements) partitioning saves nothing: run the whole // input in place instead of paying the gather. if (nBig >= _n - _n/16) { - runPartition(r, _bases, _scalars, _scalarSize, overallMaxBits + 2, _n); + partitions.emplace_back(); + Partition &p = partitions.back(); + p.bases = _bases; + p.scalars = _scalars; + p.scalarSize = _scalarSize; + p.n = _n; + p.nBits = overallMaxBits + 2; + preparePartition(p, nThreads); return; } // All scalars fit in 64 bits: gather only the scalars, bases stay in place. if (nSmall == _n) { - std::unique_ptr smallScalars(new uint64_t[_n]); + partitions.emplace_back(); + Partition &p = partitions.back(); + p.ownScalars64.reset(new uint64_t[_n]); + + uint64_t *s64 = p.ownScalars64.get(); - threadPool.parallelFor(0, _n, [&] (int begin, int end, int numThread) { + threadPool.parallelFor(0, _n, [&, s64] (int begin, int end, int numThread) { for (int i = begin; i < end; i++) { - std::memcpy(&smallScalars[i], _scalars + (uint64_t)i*_scalarSize, sizeof(uint64_t)); + std::memcpy(&s64[i], _scalars + (uint64_t)i*_scalarSize, sizeof(uint64_t)); } }); - runPartition(r, _bases, (uint8_t *)smallScalars.get(), sizeof(uint64_t), maxSmallBits + 2, _n); + p.bases = _bases; + p.scalars = (uint8_t *)s64; + p.scalarSize = sizeof(uint64_t); + p.n = _n; + p.nBits = maxSmallBits + 2; + preparePartition(p, nThreads); return; } - std::unique_ptr smallScalars(nSmall ? new uint64_t[nSmall] : nullptr); - std::unique_ptr smallBases(nSmall ? new typename Curve::PointAffine[nSmall] : nullptr); - std::unique_ptr bigScalars(nBig ? new uint8_t[nBig*_scalarSize] : nullptr); - std::unique_ptr bigBases(nBig ? new typename Curve::PointAffine[nBig] : nullptr); - std::unique_ptr onesAcc(new typename Curve::Point[nBlocks]); + Partition *small = NULL; + Partition *big = NULL; - threadPool.parallelFor(0, nBlocks, [&] (int begin, int end, int numThread) { + if (nSmall > 0) { + partitions.emplace_back(); + small = &partitions.back(); + small->ownScalars64.reset(new uint64_t[nSmall]); + small->ownBases.reset(new typename Curve::PointAffine[nSmall]); + small->bases = small->ownBases.get(); + small->scalars = (uint8_t *)small->ownScalars64.get(); + small->scalarSize = sizeof(uint64_t); + small->n = nSmall; + small->nBits = maxSmallBits + 2; + } + if (nBig > 0) { + partitions.emplace_back(); + big = &partitions.back(); + big->ownScalars.reset(new uint8_t[nBig*_scalarSize]); + big->ownBases.reset(new typename Curve::PointAffine[nBig]); + big->bases = big->ownBases.get(); + big->scalars = big->ownScalars.get(); + big->scalarSize = _scalarSize; + big->n = nBig; + big->nBits = maxBigBits + 2; + } + + nOnesBlocks = nBlocks; + onesAcc.reset(new typename Curve::Point[nBlocks]); + + typename Curve::Point *ones = onesAcc.get(); + + threadPool.parallelFor(0, nBlocks, [&, ones] (int begin, int end, int numThread) { for (int b = begin; b < end; b++) { const uint64_t i0 = (uint64_t)b*blockSize; const uint64_t i1 = std::min(i0 + blockSize, _n); uint64_t smallCur = blockOffsets[b*2]; uint64_t bigCur = blockOffsets[b*2+1]; - g.copy(onesAcc[b], g.zero()); + g.copy(ones[b], g.zero()); for (uint64_t i = i0; i < i1; i++) { switch (classes[i]) { case CLS_ONE: - g.add(onesAcc[b], onesAcc[b], _bases[i]); + g.add(ones[b], ones[b], _bases[i]); break; case CLS_SMALL: - std::memcpy(&smallScalars[smallCur], _scalars + i*_scalarSize, sizeof(uint64_t)); - smallBases[smallCur] = _bases[i]; + std::memcpy(&small->ownScalars64[smallCur], _scalars + i*_scalarSize, sizeof(uint64_t)); + small->ownBases[smallCur] = _bases[i]; smallCur++; break; case CLS_BIG: - std::memcpy(&bigScalars[bigCur*_scalarSize], _scalars + i*_scalarSize, _scalarSize); - bigBases[bigCur] = _bases[i]; + std::memcpy(&big->ownScalars[bigCur*_scalarSize], _scalars + i*_scalarSize, _scalarSize); + big->ownBases[bigCur] = _bases[i]; bigCur++; break; default: @@ -150,132 +258,157 @@ void MSM::run(typename Curve::Point &r, } }); - typename Curve::Point acc; - - g.copy(acc, onesAcc[0]); - for (uint64_t b = 1; b < nBlocks; b++) { - g.add(acc, acc, onesAcc[b]); - } - - if (nSmall > 0) { - typename Curve::Point rSmall; + if (small) preparePartition(*small, nThreads); + if (big) preparePartition(*big, nThreads); +} - runPartition(rSmall, smallBases.get(), (uint8_t *)smallScalars.get(), - sizeof(uint64_t), maxSmallBits + 2, nSmall); - g.add(acc, acc, rSmall); - } - if (nBig > 0) { - typename Curve::Point rBig; +template +uint64_t MSM::maxBuckets() const +{ + uint64_t m = 0; - runPartition(rBig, bigBases.get(), bigScalars.get(), - _scalarSize, maxBigBits + 2, nBig); - g.add(acc, acc, rBig); + for (const Partition &p : partitions) { + if (p.nBuckets > m) m = p.nBuckets; } - - g.copy(r, acc); + return m; } template -void MSM::runPartition(typename Curve::Point &r, - typename Curve::PointAffine *_bases, - uint8_t* _scalars, - uint64_t _scalarSize, - uint64_t _nBits, - uint64_t _n) +void MSM::collectTasks(std::vector &tasks, + typename Curve::Point *bucketArena, + uint64_t bucketsPerThread) { - ThreadPool &threadPool = ThreadPool::defaultPool(); + for (Partition &part : partitions) { + Partition *p = ∂ - const uint64_t nThreads = threadPool.getThreadCount(); - const uint64_t nPoints = _n; - - scalars = _scalars; - scalarSize = _scalarSize; + for (uint64_t s = 0; s < p->nSlices; s++) { + const uint64_t i0 = p->n * s / p->nSlices; + const uint64_t i1 = p->n * (s+1) / p->nSlices; -#ifdef MSM_BITS_PER_CHUNK - bitsPerChunk = MSM_BITS_PER_CHUNK; -#else - bitsPerChunk = calcBitsPerChunk(nPoints, _nBits, nThreads); -#endif + for (uint64_t j = 0; j < p->nChunks; j++) { + tasks.push_back([this, p, s, j, i0, i1, bucketArena, bucketsPerThread] (uint64_t threadId) { + typename Curve::Point *buckets = &bucketArena[threadId*bucketsPerThread]; + const int32_t *digits = &p->digits[j*p->n]; + typename Curve::PointAffine *bases = p->bases; + const uint64_t nBuckets = p->nBuckets; - if (nPoints == 0) { - g.copy(r, g.zero()); - return; - } - if (nPoints == 1) { - g.mulByScalar(r, _bases[0], scalars, scalarSize); - return; - } + for (uint64_t i = 0; i < nBuckets; i++) { + g.copy(buckets[i], g.zero()); + } - const uint64_t nChunks = calcChunkCount(_nBits, bitsPerChunk); - const uint64_t nBuckets = calcBucketCount(bitsPerChunk); - const uint64_t matrixSize = nThreads * nBuckets; - const uint64_t nSlices = nChunks*nPoints; + for (uint64_t i = i0; i < i1; i++) { + const int32_t bucketIndex = digits[i]; - std::unique_ptr bucketMatrix(new typename Curve::Point[matrixSize]); - std::unique_ptr chunks(new typename Curve::Point[nChunks]); - std::unique_ptr slicedScalars(new int32_t[nSlices]); + if (bucketIndex > 0) { + g.add(buckets[bucketIndex-1], buckets[bucketIndex-1], bases[i]); - threadPool.parallelFor(0, nPoints, [&] (int begin, int end, int numThread) { + } else if (bucketIndex < 0) { + g.sub(buckets[-bucketIndex-1], buckets[-bucketIndex-1], bases[i]); + } + } - for (int i = begin; i < end; i++) { - int carry = 0; + typename Curve::Point t, tmp; - for (int j = 0; j < nChunks; j++) { - int bucketIndex = getBucketIndex(i, j) + carry; + g.copy(t, buckets[nBuckets - 1]); + g.copy(tmp, t); - if (bucketIndex >= nBuckets) { - bucketIndex -= nBuckets*2; - carry = 1; - } else { - carry = 0; - } + for (int64_t i = nBuckets - 2; i >= 0 ; i--) { + g.add(tmp, tmp, buckets[i]); + g.add(t, t, tmp); + } - slicedScalars[j*nPoints + i] = bucketIndex; + p->partials[s*p->nChunks + j] = t; + }); } } - }); + } +} - threadPool.parallelFor(0, nChunks, [&] (int begin, int end, int numThread) { +template +void MSM::reducePartition(Partition &p, typename Curve::Point &r) +{ + typename Curve::Point chunkSum; - for (int j = begin; j < end; j++) { + for (int64_t j = p.nChunks - 1; j >= 0; j--) { + g.copy(chunkSum, p.partials[j]); + for (uint64_t s = 1; s < p.nSlices; s++) { + g.add(chunkSum, chunkSum, p.partials[s*p.nChunks + j]); + } - typename Curve::Point *buckets = &bucketMatrix[numThread*nBuckets]; + if (j == (int64_t)p.nChunks - 1) { + g.copy(r, chunkSum); + } else { + g.add(r, r, chunkSum); + } - for (int i = 0; i < nBuckets; i++) { - g.copy(buckets[i], g.zero()); + if (j > 0) { + for (uint64_t b = 0; b < p.bitsPerChunk; b++) { + g.dbl(r, r); } + } + } +} - for (int i = 0; i < nPoints; i++) { - const int bucketIndex = slicedScalars[j*nPoints + i]; +template +void MSM::finish(typename Curve::Point &r) +{ + if (trivial) { + g.copy(r, trivialResult); + prepared = false; + return; + } - if (bucketIndex > 0) { - g.add(buckets[bucketIndex-1], buckets[bucketIndex-1], _bases[i]); + typename Curve::Point acc, part; - } else if (bucketIndex < 0) { - g.sub(buckets[-bucketIndex-1], buckets[-bucketIndex-1], _bases[i]); - } - } + g.copy(acc, g.zero()); - typename Curve::Point t, tmp; + for (Partition &p : partitions) { + reducePartition(p, part); + g.add(acc, acc, part); + } - g.copy(t, buckets[nBuckets - 1]); - g.copy(tmp, t); + for (uint64_t b = 0; b < nOnesBlocks; b++) { + g.add(acc, acc, onesAcc[b]); + } - for (int i = nBuckets - 2; i >= 0 ; i--) { - g.add(tmp, tmp, buckets[i]); - g.add(t, t, tmp); - } + g.copy(r, acc); - chunks[j] = t; - } - }); + partitions.clear(); + onesAcc.reset(); + nOnesBlocks = 0; + prepared = false; +} - g.copy(r, chunks[nChunks - 1]); +template +void MSM::run(typename Curve::Point &r, + typename Curve::PointAffine *_bases, + uint8_t* _scalars, + uint64_t _scalarSize, + uint64_t _n, + uint64_t _nThreads) +{ + ThreadPool &threadPool = ThreadPool::defaultPool(); + + prepare(_bases, _scalars, _scalarSize, _n); + + if (!trivial) { + const uint64_t nThreads = threadPool.getThreadCount(); + const uint64_t bucketsPerThread = maxBuckets(); + + std::unique_ptr arena( + new typename Curve::Point[nThreads * bucketsPerThread]); - for (int j = nChunks - 2; j >= 0; j--) { - for (int i = 0; i < bitsPerChunk; i++) { - g.dbl(r, r); + std::vector tasks; + collectTasks(tasks, arena.get(), bucketsPerThread); + + if (!tasks.empty()) { + threadPool.parallelFor(0, tasks.size(), [&] (int begin, int end, int numThread) { + for (int t = begin; t < end; t++) { + tasks[t]((uint64_t)numThread); + } + }); } - g.add(r, r, chunks[j]); } + + finish(r); } diff --git a/c/msm.hpp b/c/msm.hpp index f5ffd19..145de74 100644 --- a/c/msm.hpp +++ b/c/msm.hpp @@ -2,9 +2,29 @@ #define MSM_HPP #include - +#include +#include +#include + +// Pippenger bucket-method MSM with scalar-size partitioning and a +// task-based execution model. +// +// Usage (single MSM): run() — prepare, execute and reduce in one call. +// +// Usage (batched, e.g. the prover's A/B1/B2/C phase): +// msm.prepare(bases, scalars, size, n, parallelismShare); +// msm.collectTasks(tasks, arena); // arena: nThreads*maxBuckets() Points +// ... run all MSMs' tasks in one parallel region ... +// msm.finish(r); +// Tasks from several MSM instances (even over different curves) can be +// mixed in one region; each task gets the executing thread id and uses +// that thread's row of its own curve's bucket arena. template class MSM { +public: + typedef std::function Task; + +private: const uint64_t MIN_CHUNK_SIZE_BITS = 3; const uint64_t MAX_CHUNK_SIZE_BITS = 16; @@ -14,40 +34,49 @@ class MSM { // still: they need no scalar multiplication at all. static const uint64_t SMALL_SCALAR_BITS = 64; + // Don't point-split below this many points per slice: the per-slice + // running-sum cost (2^c bucket additions) would dominate. + static const uint64_t MIN_POINTS_PER_SLICE = 4096; + + // One scalar-size class of the input, ready for bucket accumulation. + struct Partition { + typename Curve::PointAffine *bases; // points (caller's or gathered) + uint8_t *scalars; // scalars (caller's or gathered) + uint64_t scalarSize; + uint64_t n; + uint64_t nBits; // significant bits + carry headroom + uint64_t bitsPerChunk; + uint64_t nChunks; + uint64_t nBuckets; + uint64_t nSlices; // point-split factor + std::unique_ptr digits; // chunk-major [nChunks][n] + std::unique_ptr partials; // [nSlices][nChunks] + + // backing storage when the class was gathered + std::unique_ptr ownBases; + std::unique_ptr ownScalars64; + std::unique_ptr ownScalars; + }; + Curve &g; + + // significantBits()/getBucketIndex() context for the classification and + // recode passes; set before each pass. uint8_t *scalars; uint64_t scalarSize; uint64_t bitsPerChunk; -private: - // Estimated wall-clock cost in point additions. Chunks run in parallel, - // so the elapsed time is the per-chunk cost times the number of waves of - // nThreads chunks; the second term charges 1/8 of the total work so that - // among near-equal wall costs the one burning fewer total additions wins. - uint64_t calcCost(uint64_t nPoints, uint64_t nBits, uint64_t bitsPerChunk, - uint64_t nThreads) const { - const uint64_t chunkCost = nPoints + ((uint64_t)1 << bitsPerChunk) + bitsPerChunk + 1; - const uint64_t nChunks = calcChunkCount(nBits, bitsPerChunk); - const uint64_t waves = (nChunks + nThreads - 1) / nThreads; - - return waves*chunkCost + nChunks*chunkCost/(8*nThreads); - } - - uint64_t calcBitsPerChunk(uint64_t n, uint64_t nBits, uint64_t nThreads) const { - uint64_t bitsPerChunk = MIN_CHUNK_SIZE_BITS; - uint64_t minCost = calcCost(n, nBits, bitsPerChunk, nThreads); + std::vector partitions; + std::unique_ptr onesAcc; // per-block partial sums of 1-scalar points + uint64_t nOnesBlocks; + bool prepared; - for (uint64_t k = MIN_CHUNK_SIZE_BITS + 1; k <= MAX_CHUNK_SIZE_BITS; k++) { - const uint64_t curCost = calcCost(n, nBits, k, nThreads); - - if (curCost < minCost) { - minCost = curCost; - bitsPerChunk = k; - } - } - return bitsPerChunk; - } + // Set when prepare() resolved the whole MSM without bucket work + // (n==0, n==1, or every scalar in {0,1}). + bool trivial; + typename Curve::Point trivialResult; +private: uint64_t calcChunkCount(uint64_t nBits, uint64_t bitsPerChunk) const { return ((nBits - 1) / bitsPerChunk) + 1; } @@ -56,6 +85,42 @@ class MSM { return ((uint64_t)1 << (bitsPerChunk-1)); } + // Estimated wall-clock cost in point additions of one partition executed + // as nSlices*nChunks tasks on nThreads threads. Tasks run in waves; the + // second term charges 1/8 of the total work so that among near-equal + // wall costs the one burning fewer total additions wins. + uint64_t calcCost(uint64_t n, uint64_t nBits, uint64_t bitsPerChunk, + uint64_t nSlices, uint64_t nThreads) const { + const uint64_t sliceCost = n/nSlices + ((uint64_t)1 << bitsPerChunk) + bitsPerChunk + 1; + const uint64_t nTasks = nSlices * calcChunkCount(nBits, bitsPerChunk); + const uint64_t waves = (nTasks + nThreads - 1) / nThreads; + + return waves*sliceCost + nTasks*sliceCost/(8*nThreads); + } + + // Pick window size and point-split factor minimizing estimated wall cost. + void calcChunkConfig(uint64_t n, uint64_t nBits, uint64_t nThreads, + uint64_t &bestC, uint64_t &bestSlices) const { + const uint64_t maxSlices = std::max(1, std::min( + 2*nThreads, n / MIN_POINTS_PER_SLICE)); + + bestC = MIN_CHUNK_SIZE_BITS; + bestSlices = 1; + uint64_t minCost = calcCost(n, nBits, bestC, 1, nThreads); + + for (uint64_t k = MIN_CHUNK_SIZE_BITS; k <= MAX_CHUNK_SIZE_BITS; k++) { + for (uint64_t s = 1; s <= maxSlices; s *= 2) { + const uint64_t curCost = calcCost(n, nBits, k, s, nThreads); + + if (curCost < minCost) { + minCost = curCost; + bestC = k; + bestSlices = s; + } + } + } + } + uint64_t getBucketIndex(uint64_t scalarIdx, uint64_t chunkIdx) const { uint64_t bitStart = chunkIdx*bitsPerChunk; @@ -88,29 +153,46 @@ class MSM { return 0; } -public: - MSM(Curve &_g): g(_g) {} + // Recode a partition's scalars into signed digits and size its partials. + void preparePartition(Partition &p, uint64_t nThreads); - // Partitions the scalars by significant bits (0, 1, up to 64 bits, wider) - // and runs the bucket method separately per partition, so the mostly-0/1 - // scalars of a circom witness don't pay full-width window costs. + // Reduce one partition's task partials into a single point. + void reducePartition(Partition &p, typename Curve::Point &r); + +public: + MSM(Curve &_g): g(_g), prepared(false), trivial(false) {} + + // Classify scalars by size, gather the classes and recode digits. + // parallelismShare: number of threads this MSM should assume it has for + // itself when sizing windows/slices — pass the pool size when the MSM + // runs alone, or roughly poolSize/nMSMs when batched with others. + void prepare(typename Curve::PointAffine *bases, + uint8_t *scalars, + uint64_t scalarSize, + uint64_t n, + uint64_t parallelismShare = 0); + + // Largest bucket row any of this MSM's tasks needs; the caller provides + // an arena of nThreads*maxBuckets() Points to collectTasks(). + uint64_t maxBuckets() const; + + // Append one task per (partition, slice, chunk). Tasks only touch their + // own partials and bucketArena[threadId*bucketsPerThread..]. When several + // MSMs share an arena, bucketsPerThread is the max of their maxBuckets(). + void collectTasks(std::vector &tasks, + typename Curve::Point *bucketArena, + uint64_t bucketsPerThread); + + // Reduce all partials into the final result. Call after every task ran. + void finish(typename Curve::Point &r); + + // Single-MSM convenience: prepare + run own tasks + finish. void run(typename Curve::Point &r, typename Curve::PointAffine *_bases, uint8_t* _scalars, uint64_t _scalarSize, uint64_t _n, uint64_t _nThreads=0); - - // One bucket-method pass over all n points using only the lowest nBits - // bits of each scalar. nBits must exceed the largest scalar's significant - // bit count by at least 2, so a signed-digit carry can never propagate - // out of the top chunk. - void runPartition(typename Curve::Point &r, - typename Curve::PointAffine *bases, - uint8_t* scalars, - uint64_t scalarSize, - uint64_t nBits, - uint64_t n); }; #include "msm.cpp" From 4dbb440728921daa95ad9381eccded6dc08debc9 Mon Sep 17 00:00:00 2001 From: Oleksandr Brezhniev Date: Fri, 10 Jul 2026 20:13:44 +0100 Subject: [PATCH 3/6] perf: batch-affine bucket accumulation Bucket additions now run in affine coordinates in batches of up to 512, sharing one field inversion via Montgomery's trick (~5M+1S per addition instead of 9M+2S for an XYZZ mixed add) and halving bucket cache footprint. Additions that conflict with the pending batch and equal-x cases (doubling, cancellation) divert to an XYZZ shadow bucket array; infinity input bases are skipped explicitly. All-zero bytes encode infinity in both representations, so bucket init is a memset. Gated per partition: window >= 10 bits and fill density >= 1 point per bucket; sparser work keeps the XYZZ path. The task arena becomes byte-based (arenaBytesPerThread) to carve affine buckets, shadow buckets and batch scratch per thread. Measured on 20 cores, interleaved: uniform-scalar MSMs -22% (G1) and -24% (G2), realistic iden3 witness -17/-22%; sha256 proof -16% end to end on top of the previous waves. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01JpMJeCdh75Shs3hBRFFLkq --- c/msm.cpp | 241 +++++++++++++++++++++++++++++++++++++++++++++--------- c/msm.hpp | 46 +++++++++-- 2 files changed, 241 insertions(+), 46 deletions(-) diff --git a/c/msm.cpp b/c/msm.cpp index 8118976..2e3f113 100644 --- a/c/msm.cpp +++ b/c/msm.cpp @@ -21,6 +21,13 @@ void MSM::preparePartition(Partition &p, uint64_t nThreads) p.digits.reset(new int32_t[p.nChunks * p.n]); p.partials.reset(new typename Curve::Point[p.nSlices * p.nChunks]); + // Batch-affine pays off only when the bucket array is large (the batch + // stays conflict-free) and densely filled (its two bucket arrays get + // amortized over many additions). + p.batchAffine = (p.bitsPerChunk >= MIN_BATCH_AFFINE_CHUNK_BITS) + && (p.n / p.nSlices >= p.nBuckets); + p.batchSize = std::min(BATCH_SIZE, p.nBuckets/8); + // recode context for getBucketIndex scalars = p.scalars; scalarSize = p.scalarSize; @@ -263,20 +270,201 @@ void MSM::prepare(typename Curve::PointAffine *_bases, } template -uint64_t MSM::maxBuckets() const +uint64_t MSM::arenaBytesPerThread() const { uint64_t m = 0; for (const Partition &p : partitions) { - if (p.nBuckets > m) m = p.nBuckets; + const uint64_t bytes = partitionArenaBytes(p); + if (bytes > m) m = bytes; } - return m; + return (m + 63) & ~(uint64_t)63; +} + +template +void MSM::fillChunkXYZZ(Partition &p, uint64_t j, + uint64_t i0, uint64_t i1, + uint64_t sliceIdx, uint8_t *taskArena) +{ + typename Curve::Point *buckets = (typename Curve::Point *)taskArena; + const int32_t *digits = &p.digits[j*p.n]; + typename Curve::PointAffine *bases = p.bases; + const uint64_t nBuckets = p.nBuckets; + + for (uint64_t i = 0; i < nBuckets; i++) { + g.copy(buckets[i], g.zero()); + } + + for (uint64_t i = i0; i < i1; i++) { + const int32_t bucketIndex = digits[i]; + + if (bucketIndex > 0) { + g.add(buckets[bucketIndex-1], buckets[bucketIndex-1], bases[i]); + + } else if (bucketIndex < 0) { + g.sub(buckets[-bucketIndex-1], buckets[-bucketIndex-1], bases[i]); + } + } + + typename Curve::Point t, tmp; + + g.copy(t, buckets[nBuckets - 1]); + g.copy(tmp, t); + + for (int64_t i = nBuckets - 2; i >= 0 ; i--) { + g.add(tmp, tmp, buckets[i]); + g.add(t, t, tmp); + } + + p.partials[sliceIdx*p.nChunks + j] = t; +} + +template +void MSM::fillChunkBatchAffine(Partition &p, uint64_t j, + uint64_t i0, uint64_t i1, + uint64_t sliceIdx, uint8_t *taskArena) +{ + typedef typename Curve::PointAffine PointAffine; + typedef typename Curve::Point Point; + typedef typename BaseField::Element Element; + + const uint64_t nBuckets = p.nBuckets; + const uint64_t batchSize = p.batchSize; + const int32_t *digits = &p.digits[j*p.n]; + PointAffine *bases = p.bases; + BaseField &F = g.F; + + uint8_t *cur = taskArena; + PointAffine *buckets = (PointAffine *)cur; cur += nBuckets*sizeof(PointAffine); + Point *shadow = (Point *)cur; cur += nBuckets*sizeof(Point); + PointAffine *batchP = (PointAffine *)cur; cur += batchSize*sizeof(PointAffine); + Element *dx = (Element *)cur; cur += batchSize*sizeof(Element); + Element *prod = (Element *)cur; cur += batchSize*sizeof(Element); + uint32_t *batchB = (uint32_t *)cur; cur += batchSize*sizeof(uint32_t); + uint8_t *inBatch = cur; // nBuckets bytes + + // all-zero bytes encode infinity in both representations + std::memset(buckets, 0, nBuckets*sizeof(PointAffine)); + std::memset(shadow, 0, nBuckets*sizeof(Point)); + std::memset(inBatch, 0, nBuckets); + + uint64_t count = 0; + + // Execute the pending independent affine additions, amortizing one + // inversion over the whole batch (Montgomery's trick). + auto executeBatch = [&] () { + if (count == 0) return; + + for (uint64_t k = 0; k < count; k++) { + F.sub(dx[k], batchP[k].x, buckets[batchB[k]].x); + + if (k == 0) { + F.copy(prod[0], dx[0]); + } else { + F.mul(prod[k], prod[k-1], dx[k]); + } + } + + Element invAll, invK, lambda, t1, x3; + + F.inv(invAll, prod[count-1]); + + for (int64_t k = count - 1; k >= 0; k--) { + PointAffine &B = buckets[batchB[k]]; + + if (k > 0) { + F.mul(invK, invAll, prod[k-1]); + F.mul(invAll, invAll, dx[k]); + } else { + F.copy(invK, invAll); + } + + // chord addition: B = B + P + F.sub(t1, batchP[k].y, B.y); + F.mul(lambda, t1, invK); + + F.square(x3, lambda); + F.sub(x3, x3, B.x); + F.sub(x3, x3, batchP[k].x); + + F.sub(t1, B.x, x3); + F.mul(t1, t1, lambda); + F.sub(B.y, t1, B.y); + F.copy(B.x, x3); + + inBatch[batchB[k]] = 0; + } + count = 0; + }; + + for (uint64_t i = i0; i < i1; i++) { + const int32_t d = digits[i]; + + if (d == 0) continue; + if (g.isZero(bases[i])) continue; + + const uint32_t b = (uint32_t)(d > 0 ? d : -d) - 1; + + PointAffine P; + F.copy(P.x, bases[i].x); + if (d > 0) { + F.copy(P.y, bases[i].y); + } else { + F.neg(P.y, bases[i].y); + } + + if (inBatch[b]) { + // the bucket has a pending addition: divert to its shadow + g.add(shadow[b], shadow[b], P); + continue; + } + if (F.isZero(buckets[b].x) && F.isZero(buckets[b].y)) { + buckets[b] = P; + continue; + } + if (F.eq(buckets[b].x, P.x)) { + if (F.eq(buckets[b].y, P.y)) { + // doubling: fold 2P into the shadow bucket + Point t2; + g.dbl(t2, P); + g.add(shadow[b], shadow[b], t2); + } + // else P == -bucket: they cancel + std::memset(&buckets[b], 0, sizeof(PointAffine)); + continue; + } + + batchB[count] = b; + batchP[count] = P; + inBatch[b] = 1; + count++; + + if (count == batchSize) executeBatch(); + } + executeBatch(); + + typename Curve::Point t, tmp; + + g.copy(t, g.zero()); + g.copy(tmp, g.zero()); + + for (int64_t b = nBuckets - 1; b >= 0; b--) { + if (!(F.isZero(buckets[b].x) && F.isZero(buckets[b].y))) { + g.add(tmp, tmp, buckets[b]); + } + if (!g.isZero(shadow[b])) { + g.add(tmp, tmp, shadow[b]); + } + g.add(t, t, tmp); + } + + p.partials[sliceIdx*p.nChunks + j] = t; } template void MSM::collectTasks(std::vector &tasks, - typename Curve::Point *bucketArena, - uint64_t bucketsPerThread) + uint8_t *bucketArena, + uint64_t bytesPerThread) { for (Partition &part : partitions) { Partition *p = ∂ @@ -286,38 +474,14 @@ void MSM::collectTasks(std::vector &tasks, const uint64_t i1 = p->n * (s+1) / p->nSlices; for (uint64_t j = 0; j < p->nChunks; j++) { - tasks.push_back([this, p, s, j, i0, i1, bucketArena, bucketsPerThread] (uint64_t threadId) { - typename Curve::Point *buckets = &bucketArena[threadId*bucketsPerThread]; - const int32_t *digits = &p->digits[j*p->n]; - typename Curve::PointAffine *bases = p->bases; - const uint64_t nBuckets = p->nBuckets; - - for (uint64_t i = 0; i < nBuckets; i++) { - g.copy(buckets[i], g.zero()); - } - - for (uint64_t i = i0; i < i1; i++) { - const int32_t bucketIndex = digits[i]; - - if (bucketIndex > 0) { - g.add(buckets[bucketIndex-1], buckets[bucketIndex-1], bases[i]); - - } else if (bucketIndex < 0) { - g.sub(buckets[-bucketIndex-1], buckets[-bucketIndex-1], bases[i]); - } - } - - typename Curve::Point t, tmp; + tasks.push_back([this, p, s, j, i0, i1, bucketArena, bytesPerThread] (uint64_t threadId) { + uint8_t *taskArena = bucketArena + threadId*bytesPerThread; - g.copy(t, buckets[nBuckets - 1]); - g.copy(tmp, t); - - for (int64_t i = nBuckets - 2; i >= 0 ; i--) { - g.add(tmp, tmp, buckets[i]); - g.add(t, t, tmp); + if (p->batchAffine) { + fillChunkBatchAffine(*p, j, i0, i1, s, taskArena); + } else { + fillChunkXYZZ(*p, j, i0, i1, s, taskArena); } - - p->partials[s*p->nChunks + j] = t; }); } } @@ -393,13 +557,12 @@ void MSM::run(typename Curve::Point &r, if (!trivial) { const uint64_t nThreads = threadPool.getThreadCount(); - const uint64_t bucketsPerThread = maxBuckets(); + const uint64_t bytesPerThread = arenaBytesPerThread(); - std::unique_ptr arena( - new typename Curve::Point[nThreads * bucketsPerThread]); + std::unique_ptr arena(new uint8_t[nThreads * bytesPerThread]); std::vector tasks; - collectTasks(tasks, arena.get(), bucketsPerThread); + collectTasks(tasks, arena.get(), bytesPerThread); if (!tasks.empty()) { threadPool.parallelFor(0, tasks.size(), [&] (int begin, int end, int numThread) { diff --git a/c/msm.hpp b/c/msm.hpp index 145de74..204f188 100644 --- a/c/msm.hpp +++ b/c/msm.hpp @@ -38,6 +38,15 @@ class MSM { // running-sum cost (2^c bucket additions) would dominate. static const uint64_t MIN_POINTS_PER_SLICE = 4096; + // Batch-affine accumulation: buckets live in affine coordinates and + // additions are executed in batches sharing one field inversion + // (~5M+1S per addition instead of 9M+2S for an XYZZ mixed add). + // Additions that conflict with the pending batch, doublings and + // cancellations go to an XYZZ shadow bucket instead. Used only when + // the bucket array is large and densely filled enough. + static const uint64_t BATCH_SIZE = 512; + static const uint64_t MIN_BATCH_AFFINE_CHUNK_BITS = 10; + // One scalar-size class of the input, ready for bucket accumulation. struct Partition { typename Curve::PointAffine *bases; // points (caller's or gathered) @@ -49,6 +58,8 @@ class MSM { uint64_t nChunks; uint64_t nBuckets; uint64_t nSlices; // point-split factor + bool batchAffine; // bucket accumulation strategy + uint64_t batchSize; std::unique_ptr digits; // chunk-major [nChunks][n] std::unique_ptr partials; // [nSlices][nChunks] @@ -159,6 +170,25 @@ class MSM { // Reduce one partition's task partials into a single point. void reducePartition(Partition &p, typename Curve::Point &r); + // Bucket accumulation + running sum for points [i0,i1) of chunk j, + // writing the result into the (slice, chunk) partial. + void fillChunkXYZZ(Partition &p, uint64_t j, uint64_t i0, uint64_t i1, + uint64_t sliceIdx, uint8_t *taskArena); + void fillChunkBatchAffine(Partition &p, uint64_t j, uint64_t i0, uint64_t i1, + uint64_t sliceIdx, uint8_t *taskArena); + + uint64_t partitionArenaBytes(const Partition &p) const { + if (!p.batchAffine) { + return p.nBuckets * sizeof(typename Curve::Point); + } + return p.nBuckets * (sizeof(typename Curve::PointAffine) + + sizeof(typename Curve::Point) + 1) + + p.batchSize * (sizeof(typename Curve::PointAffine) + + 2*sizeof(typename BaseField::Element) + + sizeof(uint32_t)) + + 64; // alignment slack + } + public: MSM(Curve &_g): g(_g), prepared(false), trivial(false) {} @@ -172,16 +202,18 @@ class MSM { uint64_t n, uint64_t parallelismShare = 0); - // Largest bucket row any of this MSM's tasks needs; the caller provides - // an arena of nThreads*maxBuckets() Points to collectTasks(). - uint64_t maxBuckets() const; + // Largest per-thread scratch any of this MSM's tasks needs; the caller + // provides an arena of nThreads*arenaBytesPerThread() bytes to + // collectTasks() (8-byte aligned, e.g. from new uint8_t[]). + uint64_t arenaBytesPerThread() const; // Append one task per (partition, slice, chunk). Tasks only touch their - // own partials and bucketArena[threadId*bucketsPerThread..]. When several - // MSMs share an arena, bucketsPerThread is the max of their maxBuckets(). + // own partials and bucketArena[threadId*bytesPerThread..]. When several + // MSMs share an arena, bytesPerThread is the max of their + // arenaBytesPerThread(). void collectTasks(std::vector &tasks, - typename Curve::Point *bucketArena, - uint64_t bucketsPerThread); + uint8_t *bucketArena, + uint64_t bytesPerThread); // Reduce all partials into the final result. Call after every task ran. void finish(typename Curve::Point &r); From 0a8075aba7d1a0bb613b8a53a8b84d80dab2ddae Mon Sep 17 00:00:00 2001 From: Oleksandr Brezhniev Date: Fri, 10 Jul 2026 21:27:34 +0100 Subject: [PATCH 4/6] perf: permutation-free FFT transform pair Adds ifftDIFNatToRev (decimation in frequency with inverse twiddles, natural order in, bit-reversed out, 1/n left to the caller) and fftDITRevToNat (the existing DIT butterfly core without the input permutation), plus rootInv/nInv accessors. Chaining DIF, a pointwise pass indexed through BR, and DIT computes coset evaluations without any bit-reversal permutation of the data. fft()/ifft() are unchanged. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01JpMJeCdh75Shs3hBRFFLkq --- c/fft.cpp | 54 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ c/fft.hpp | 13 +++++++++++++ 2 files changed, 67 insertions(+) diff --git a/c/fft.cpp b/c/fft.cpp index c557ae0..6c1c34e 100644 --- a/c/fft.cpp +++ b/c/fft.cpp @@ -198,6 +198,60 @@ void FFT::fft(Element *a, u_int64_t n) { } } +template +void FFT::fftDITRevToNat(Element *a, u_int64_t n) { + u_int64_t domainPow = log2(n); + assert(((u_int64_t)1 << domainPow) == n); + + for (u_int32_t s=1; s<=domainPow; s++) { + u_int64_t m = 1 << s; + u_int64_t mdiv2 = m >> 1; + + threadPool.parallelFor(0, (n>>1), [&] (int begin, int end, int numThread) { + for (u_int64_t i=begin; i< end; i++) { + Element t; + Element u; + u_int64_t k=(i/mdiv2)*m; + u_int64_t j=i%mdiv2; + + f.mul(t, root(s, j), a[k+j+mdiv2]); + f.copy(u,a[k+j]); + f.add(a[k+j], t, u); + f.sub(a[k+j+mdiv2], u, t); + } + }); + } +} + +// Inverse of fftDITRevToNat run backwards: decimation in frequency with +// inverse twiddles. Leaves the result scaled by n; the caller folds 1/n +// into its next pointwise pass. +template +void FFT::ifftDIFNatToRev(Element *a, u_int64_t n) { + u_int64_t domainPow = log2(n); + assert(((u_int64_t)1 << domainPow) == n); + + for (u_int32_t s=domainPow; s>=1; s--) { + u_int64_t m = 1 << s; + u_int64_t mdiv2 = m >> 1; + + threadPool.parallelFor(0, (n>>1), [&] (int begin, int end, int numThread) { + for (u_int64_t i=begin; i< end; i++) { + Element t; + Element u; + u_int64_t k=(i/mdiv2)*m; + u_int64_t j=i%mdiv2; + + f.copy(u, a[k+j]); + f.copy(t, a[k+j+mdiv2]); + f.add(a[k+j], u, t); + f.sub(t, u, t); + f.mul(a[k+j+mdiv2], t, rootInv(s, j)); + } + }); + } +} + template void FFT::ifft(Element *a, u_int64_t n ) { fft(a, n); diff --git a/c/fft.hpp b/c/fft.hpp index ab3c63b..c88fc0d 100644 --- a/c/fft.hpp +++ b/c/fft.hpp @@ -24,8 +24,21 @@ class FFT { void fft(Element *a, u_int64_t n ); void ifft(Element *a, u_int64_t n ); + // Permutation-free pair: ifftDIFNatToRev takes natural order and leaves + // the (unscaled by 1/n!) inverse transform in bit-reversed order; + // fftDITRevToNat takes bit-reversed order and leaves the forward + // transform in natural order. Chaining them with a pointwise pass in + // between (indexed through BR) round-trips to natural order without any + // bit-reversal permutation of the data. + void ifftDIFNatToRev(Element *a, u_int64_t n); + void fftDITRevToNat(Element *a, u_int64_t n); + u_int32_t log2(u_int64_t n); inline Element &root(u_int32_t domainPow, u_int64_t idx) { return roots[ idx << (s-domainPow)]; } + inline Element &rootInv(u_int32_t domainPow, u_int64_t idx) { + return roots[ idx == 0 ? 0 : ((((u_int64_t)1 << domainPow) - idx) << (s-domainPow)) ]; + } + inline Element &nInv(u_int32_t domainPow) { return powTwoInv[domainPow]; } void printVector(Element *a, u_int64_t n ); From 9e42fc79c712c53213d838e72b0d40056987a67f Mon Sep 17 00:00:00 2001 From: Oleksandr Brezhniev Date: Sat, 11 Jul 2026 16:53:19 +0100 Subject: [PATCH 5/6] fix: undefined MSM constants when linked without optimization std::min bound BATCH_SIZE by const reference, ODR-using the in-class static const without an out-of-class definition; Apple clang's -O0 dylib link (iOS simulator Debug) failed with an undefined symbol. Turn the class constants into a typed enum, which is never ODR-used. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01JpMJeCdh75Shs3hBRFFLkq --- c/msm.cpp | 2 +- c/msm.hpp | 40 ++++++++++++++++++++++------------------ 2 files changed, 23 insertions(+), 19 deletions(-) diff --git a/c/msm.cpp b/c/msm.cpp index 2e3f113..af85eaf 100644 --- a/c/msm.cpp +++ b/c/msm.cpp @@ -26,7 +26,7 @@ void MSM::preparePartition(Partition &p, uint64_t nThreads) // amortized over many additions). p.batchAffine = (p.bitsPerChunk >= MIN_BATCH_AFFINE_CHUNK_BITS) && (p.n / p.nSlices >= p.nBuckets); - p.batchSize = std::min(BATCH_SIZE, p.nBuckets/8); + p.batchSize = std::min(BATCH_SIZE, p.nBuckets/8); // recode context for getBucketIndex scalars = p.scalars; diff --git a/c/msm.hpp b/c/msm.hpp index 204f188..b4db2e3 100644 --- a/c/msm.hpp +++ b/c/msm.hpp @@ -28,24 +28,28 @@ class MSM { const uint64_t MIN_CHUNK_SIZE_BITS = 3; const uint64_t MAX_CHUNK_SIZE_BITS = 16; - // Scalars of at most this many significant bits go to the small - // partition, which runs the bucket method over a single 64-bit word - // and so pays ~4 windows instead of ~16. Scalars 0 and 1 are cheaper - // still: they need no scalar multiplication at all. - static const uint64_t SMALL_SCALAR_BITS = 64; - - // Don't point-split below this many points per slice: the per-slice - // running-sum cost (2^c bucket additions) would dominate. - static const uint64_t MIN_POINTS_PER_SLICE = 4096; - - // Batch-affine accumulation: buckets live in affine coordinates and - // additions are executed in batches sharing one field inversion - // (~5M+1S per addition instead of 9M+2S for an XYZZ mixed add). - // Additions that conflict with the pending batch, doublings and - // cancellations go to an XYZZ shadow bucket instead. Used only when - // the bucket array is large and densely filled enough. - static const uint64_t BATCH_SIZE = 512; - static const uint64_t MIN_BATCH_AFFINE_CHUNK_BITS = 10; + // (enum members: in-class integral constants that are never ODR-used, + // so no out-of-class definitions are needed) + enum : uint64_t { + // Scalars of at most this many significant bits go to the small + // partition, which runs the bucket method over a single 64-bit word + // and so pays ~4 windows instead of ~16. Scalars 0 and 1 are cheaper + // still: they need no scalar multiplication at all. + SMALL_SCALAR_BITS = 64, + + // Don't point-split below this many points per slice: the per-slice + // running-sum cost (2^c bucket additions) would dominate. + MIN_POINTS_PER_SLICE = 4096, + + // Batch-affine accumulation: buckets live in affine coordinates and + // additions are executed in batches sharing one field inversion + // (~5M+1S per addition instead of 9M+2S for an XYZZ mixed add). + // Additions that conflict with the pending batch, doublings and + // cancellations go to an XYZZ shadow bucket instead. Used only when + // the bucket array is large and densely filled enough. + BATCH_SIZE = 512, + MIN_BATCH_AFFINE_CHUNK_BITS = 10 + }; // One scalar-size class of the input, ready for bucket accumulation. struct Partition { From b9a62615428e74cc151249d4b556df7126a21391 Mon Sep 17 00:00:00 2001 From: Oleksandr Brezhniev Date: Sun, 12 Jul 2026 01:50:11 +0100 Subject: [PATCH 6/6] perf: reduce MSM and FFT working memory Three memory reductions, time-neutral or slightly faster: - MSM signed digits stored as int16 instead of int32: they fit [-2^15, 2^15-1] exactly for windows up to 16 bits. Halves the largest transient allocation and the digit memory traffic. - Scalar-size partitions address the caller's bases and scalars through an ascending uint32 index list instead of gathering copies (4 bytes per point instead of a point + scalar copy per MSM). - New FFT::higherRootOfUnity(extraPow) derives a root of unity finer than the table (same nqr derivation as the constructor), so callers needing omega_2n for coset shifts no longer construct the FFT at twice the transform size, halving the roots table. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01JpMJeCdh75Shs3hBRFFLkq --- c/fft.cpp | 30 ++++++++++++++++++++++++++ c/fft.hpp | 6 ++++++ c/msm.cpp | 63 ++++++++++++++++++++++++++++++------------------------- c/msm.hpp | 13 ++++++++---- 4 files changed, 80 insertions(+), 32 deletions(-) diff --git a/c/fft.cpp b/c/fft.cpp index 6c1c34e..cdb9f0f 100644 --- a/c/fft.cpp +++ b/c/fft.cpp @@ -114,6 +114,36 @@ FFT::FFT(u_int64_t maxDomainSize, uint32_t _nThreads) mpz_clear(m_aux); } +template +void FFT::higherRootOfUnity(Element &r, u_int32_t extraPow) { + mpz_t m_q, m_aux, m_nqr; + + mpz_init(m_q); + mpz_init(m_aux); + mpz_init(m_nqr); + + f.toMpz(m_aux, f.negOne()); + mpz_add_ui(m_q, m_aux, 1); + + // (q-1) / 2^(s+extraPow); the primitive root exists iff the division + // is exact, i.e. s+extraPow is within the field's 2-adicity + if (mpz_scan1(m_aux, 0) < s + extraPow) { + mpz_clear(m_q); + mpz_clear(m_aux); + mpz_clear(m_nqr); + throw std::range_error("Root order exceeds the field's 2-adicity"); + } + mpz_fdiv_q_2exp(m_aux, m_aux, s + extraPow); + + f.toMpz(m_nqr, nqr); + mpz_powm(m_aux, m_nqr, m_aux, m_q); + f.fromMpz(r, m_aux); + + mpz_clear(m_q); + mpz_clear(m_aux); + mpz_clear(m_nqr); +} + template FFT::~FFT() { delete[] roots; diff --git a/c/fft.hpp b/c/fft.hpp index c88fc0d..21e5d49 100644 --- a/c/fft.hpp +++ b/c/fft.hpp @@ -35,6 +35,12 @@ class FFT { u_int32_t log2(u_int64_t n); inline Element &root(u_int32_t domainPow, u_int64_t idx) { return roots[ idx << (s-domainPow)]; } + + // Primitive 2^(s+extraPow)-th root of unity — an order finer than the + // table covers (e.g. the omega_2n of a coset shift, without paying for + // a table twice the transform size). Requires s+extraPow within the + // field's 2-adicity. + void higherRootOfUnity(Element &r, u_int32_t extraPow); inline Element &rootInv(u_int32_t domainPow, u_int64_t idx) { return roots[ idx == 0 ? 0 : ((((u_int64_t)1 << domainPow) - idx) << (s-domainPow)) ]; } diff --git a/c/msm.cpp b/c/msm.cpp index af85eaf..3460f1f 100644 --- a/c/msm.cpp +++ b/c/msm.cpp @@ -18,7 +18,7 @@ void MSM::preparePartition(Partition &p, uint64_t nThreads) p.nChunks = calcChunkCount(p.nBits, p.bitsPerChunk); p.nBuckets = calcBucketCount(p.bitsPerChunk); - p.digits.reset(new int32_t[p.nChunks * p.n]); + p.digits.reset(new int16_t[p.nChunks * p.n]); p.partials.reset(new typename Curve::Point[p.nSlices * p.nChunks]); // Batch-affine pays off only when the bucket array is large (the batch @@ -36,15 +36,17 @@ void MSM::preparePartition(Partition &p, uint64_t nThreads) const uint64_t nChunks = p.nChunks; const uint64_t nBuckets = p.nBuckets; const uint64_t nPoints = p.n; - int32_t *digits = p.digits.get(); + const uint32_t *indices = p.indices; + int16_t *digits = p.digits.get(); - threadPool.parallelFor(0, nPoints, [&, nChunks, nBuckets, nPoints] (int begin, int end, int numThread) { + threadPool.parallelFor(0, nPoints, [&, nChunks, nBuckets, nPoints, indices] (int begin, int end, int numThread) { for (int i = begin; i < end; i++) { int carry = 0; + const uint64_t scalarIdx = indices ? indices[i] : (uint64_t)i; for (uint64_t j = 0; j < nChunks; j++) { - int bucketIndex = getBucketIndex(i, j) + carry; + int bucketIndex = getBucketIndex(scalarIdx, j) + carry; if (bucketIndex >= (int)nBuckets) { bucketIndex -= nBuckets*2; @@ -53,7 +55,7 @@ void MSM::preparePartition(Partition &p, uint64_t nThreads) carry = 0; } - digits[j*nPoints + i] = bucketIndex; + digits[j*nPoints + i] = (int16_t)bucketIndex; } } }); @@ -95,6 +97,7 @@ void MSM::prepare(typename Curve::PointAffine *_bases, if (_scalarSize < 8) { partitions.emplace_back(); Partition &p = partitions.back(); + p.indices = NULL; p.bases = _bases; p.scalars = _scalars; p.scalarSize = _scalarSize; @@ -171,6 +174,7 @@ void MSM::prepare(typename Curve::PointAffine *_bases, if (nBig >= _n - _n/16) { partitions.emplace_back(); Partition &p = partitions.back(); + p.indices = NULL; p.bases = _bases; p.scalars = _scalars; p.scalarSize = _scalarSize; @@ -194,6 +198,7 @@ void MSM::prepare(typename Curve::PointAffine *_bases, } }); + p.indices = NULL; p.bases = _bases; p.scalars = (uint8_t *)s64; p.scalarSize = sizeof(uint64_t); @@ -209,21 +214,21 @@ void MSM::prepare(typename Curve::PointAffine *_bases, if (nSmall > 0) { partitions.emplace_back(); small = &partitions.back(); - small->ownScalars64.reset(new uint64_t[nSmall]); - small->ownBases.reset(new typename Curve::PointAffine[nSmall]); - small->bases = small->ownBases.get(); - small->scalars = (uint8_t *)small->ownScalars64.get(); - small->scalarSize = sizeof(uint64_t); + small->ownIndices.reset(new uint32_t[nSmall]); + small->indices = small->ownIndices.get(); + small->bases = _bases; + small->scalars = _scalars; + small->scalarSize = _scalarSize; small->n = nSmall; small->nBits = maxSmallBits + 2; } if (nBig > 0) { partitions.emplace_back(); big = &partitions.back(); - big->ownScalars.reset(new uint8_t[nBig*_scalarSize]); - big->ownBases.reset(new typename Curve::PointAffine[nBig]); - big->bases = big->ownBases.get(); - big->scalars = big->ownScalars.get(); + big->ownIndices.reset(new uint32_t[nBig]); + big->indices = big->ownIndices.get(); + big->bases = _bases; + big->scalars = _scalars; big->scalarSize = _scalarSize; big->n = nBig; big->nBits = maxBigBits + 2; @@ -249,14 +254,10 @@ void MSM::prepare(typename Curve::PointAffine *_bases, g.add(ones[b], ones[b], _bases[i]); break; case CLS_SMALL: - std::memcpy(&small->ownScalars64[smallCur], _scalars + i*_scalarSize, sizeof(uint64_t)); - small->ownBases[smallCur] = _bases[i]; - smallCur++; + small->ownIndices[smallCur++] = (uint32_t)i; break; case CLS_BIG: - std::memcpy(&big->ownScalars[bigCur*_scalarSize], _scalars + i*_scalarSize, _scalarSize); - big->ownBases[bigCur] = _bases[i]; - bigCur++; + big->ownIndices[bigCur++] = (uint32_t)i; break; default: break; @@ -287,8 +288,9 @@ void MSM::fillChunkXYZZ(Partition &p, uint64_t j, uint64_t sliceIdx, uint8_t *taskArena) { typename Curve::Point *buckets = (typename Curve::Point *)taskArena; - const int32_t *digits = &p.digits[j*p.n]; + const int16_t *digits = &p.digits[j*p.n]; typename Curve::PointAffine *bases = p.bases; + const uint32_t *indices = p.indices; const uint64_t nBuckets = p.nBuckets; for (uint64_t i = 0; i < nBuckets; i++) { @@ -297,12 +299,13 @@ void MSM::fillChunkXYZZ(Partition &p, uint64_t j, for (uint64_t i = i0; i < i1; i++) { const int32_t bucketIndex = digits[i]; + typename Curve::PointAffine &base = bases[indices ? indices[i] : i]; if (bucketIndex > 0) { - g.add(buckets[bucketIndex-1], buckets[bucketIndex-1], bases[i]); + g.add(buckets[bucketIndex-1], buckets[bucketIndex-1], base); } else if (bucketIndex < 0) { - g.sub(buckets[-bucketIndex-1], buckets[-bucketIndex-1], bases[i]); + g.sub(buckets[-bucketIndex-1], buckets[-bucketIndex-1], base); } } @@ -330,8 +333,9 @@ void MSM::fillChunkBatchAffine(Partition &p, uint64_t j, const uint64_t nBuckets = p.nBuckets; const uint64_t batchSize = p.batchSize; - const int32_t *digits = &p.digits[j*p.n]; + const int16_t *digits = &p.digits[j*p.n]; PointAffine *bases = p.bases; + const uint32_t *indices = p.indices; BaseField &F = g.F; uint8_t *cur = taskArena; @@ -401,16 +405,19 @@ void MSM::fillChunkBatchAffine(Partition &p, uint64_t j, const int32_t d = digits[i]; if (d == 0) continue; - if (g.isZero(bases[i])) continue; + + PointAffine &base = bases[indices ? indices[i] : i]; + + if (g.isZero(base)) continue; const uint32_t b = (uint32_t)(d > 0 ? d : -d) - 1; PointAffine P; - F.copy(P.x, bases[i].x); + F.copy(P.x, base.x); if (d > 0) { - F.copy(P.y, bases[i].y); + F.copy(P.y, base.y); } else { - F.neg(P.y, bases[i].y); + F.neg(P.y, base.y); } if (inBatch[b]) { diff --git a/c/msm.hpp b/c/msm.hpp index b4db2e3..efb2d71 100644 --- a/c/msm.hpp +++ b/c/msm.hpp @@ -64,13 +64,18 @@ class MSM { uint64_t nSlices; // point-split factor bool batchAffine; // bucket accumulation strategy uint64_t batchSize; - std::unique_ptr digits; // chunk-major [nChunks][n] + std::unique_ptr digits; // chunk-major [nChunks][n]; signed digits + // fit int16 exactly: |digit| <= 2^(c-1) <= 2^15 for c <= 16 std::unique_ptr partials; // [nSlices][nChunks] - // backing storage when the class was gathered - std::unique_ptr ownBases; + // When set, the partition addresses the caller's bases/scalars + // through this ascending index list instead of gathered copies + // (4 bytes per point instead of a point + scalar copy). + const uint32_t *indices; + + // backing storage for indexed/gathered classes + std::unique_ptr ownIndices; std::unique_ptr ownScalars64; - std::unique_ptr ownScalars; }; Curve &g;