Conversation
Multiplies two unsigned integers of kDigits digits into 2*kDigits digits,
built on MulAdd52Lo/Hi. Digits are 52 bits ("nails"), the largest width for
which the two ops are exact (for inputs < 2^52 they return the full 104-bit
product), so no separate fallback is needed: the ops exist on every target
(native HWY_NATIVE_MULADD52, generic otherwise).
Layout is SoA: digit i of each operand is passed in its own vector, so the
lanes of those vectors hold independent multiplications - the crypto use
case (several RSA/ECC instances in parallel).
Adds the CMake wiring (HWY_CONTRIB_SOURCES + test list) and a test that
checks kDigits = 1..4 against a scalar reference on all targets.
Part of google#3352. Depends on the MulAdd52 ops from google#3365.
- WideMulBits<kBits>: bit-width form of the SoA API (kDigits = ceil(kBits/52)). - WideMulLimbs<kBits>: scalar convenience that multiplies two kBits-bit integers given as little-endian 64-bit limbs, writing the 2*kBits/64-limb product (e.g. u128/u192/u256). Converts to/from the 52-bit digits and reuses the same schoolbook. - Test the limb path (kBits = 128/192/256) against an independent ripple-carry reference, plus maximum inputs.
MulRec now uses Karatsuba's 3-way split for even digit counts of at least kWideMulKaratsubaMinDigits = 8 digits: z0 = a_lo*b_lo z2 = a_hi*b_hi z1 = (a_lo+a_hi)*(b_lo+b_hi) - z0 - z2 // = a_lo*b_hi + a_hi*b_lo i.e. 3 recursive multiplications instead of 4, giving ~O(kDigits^1.585) instead of O(kDigits^2). Below the threshold, and for odd digit counts (whose split would be unbalanced), the schoolbook kernel is kept. The z1 subtraction reuses the z0/z2 already accumulated in `out`, so only s, t, st and one temporary buffer per split are needed. Adds the helpers AddDigitsCarry, AddInto and SubInto (the latter needs no separate sign handling because the subtracted result is non-negative). Test kDigits = 8 and 16 as well, still checked against the scalar reference.
|
Nice, this sounds sensible but I'm not experienced in this field. Can we first write a short 'market research' of existing approaches and both their APIs and their implementations, to help make sure our interface meets user needs? On naming: do we intend to provide more such operations, or will it stay mul-only? If the latter, we should probably rename the directory. Also, the test has arrays of vectors which does not compile on RVV/SVE. Let's replace with dynamically allocated arrays of T. |
|
Here is the market research you asked for (§1-3 are the survey, §4-5 the
Naming: I would keep The RVV/SVE test fix (arrays of vectors replaced by dynamically allocated arrays Market research: fixed-width multi-precision multiplyWhat existing approaches do, what their APIs and implementations look like, and 1. The questionThe interface has to serve a specific user: RSA/ECC-style code that multiplies
2. Prior art2.1 High-level, allocating bignums
2.2 Low-level limb arrays, caller-managed memoryThis is the classic shape, and the closest relative of what we are adding.
2.3 Fixed-size, formally verified kernelsThe modern high-assurance school converges on a very specific API shape: one
2.4 Const-generic fixed-width types
The takeaway: the width belongs in the type (or a compile-time constant), and 2.5 Field-specialized, hand-written with nails
2.6 SIMD / multibuffer libraries
Our 3. Comparison
4. What this means for our interfacePoints of agreement with the field (keep them):
Gaps worth closing before users arrive:
5. On namingThe question was whether we intend to provide more such operations or stay
6. Sources
|
|
One more finding while preparing the test fix, because I do not think the test
// multiprec-inl.h, SchoolbookMul
Vec<D> low[2 * kDigits - 1];
Vec<D> high[2 * kDigits - 1];and the public entry point takes them by array: // multiprec-inl.h
template <size_t kDigits, class D>
HWY_INLINE void SchoolbookMul(D d, const Vec<D>* HWY_RESTRICT a, ...)So replacing the test's arrays with arrays of Two ways out, and I would like your preference:
My preference is (1) for this PR (so the target-restriction is explicit and |
jan-wassenberg
left a comment
jan-wassenberg
left a comment
There was a problem hiding this comment.
Thanks for the survey. I have heard of Montgomery reduction being central in this space, and agree the multiprec name makes sense.
| // low[k] is the sum of the low 52 bits of a[i]*b[j] over all i+j == k; | ||
| // high[k] is the sum of bits 52..103 of the same products. Both fit in 64 | ||
| // bits because there are at most kDigits terms, each < 2^52. | ||
| Vec<D> low[2 * kDigits - 1]; |
There was a problem hiding this comment.
In addition to the SVE issue (no arrays of vectors), this may spill even for non-scalable targets for large kDigits. Comba's method loops over k and sums i+j = k and would only use 3 vectors.
| } | ||
|
|
||
| // s = a_lo + a_hi and t = b_lo + b_hi, each kHalf + 1 digits. | ||
| Vec<D> s[kHalf + 1]; |
There was a problem hiding this comment.
Also array of T with MaxLanes * kNumVectors.
|
|
||
| // Fixed-width unsigned integer multiplication for arbitrary precision. | ||
| // | ||
| // Digits are 52 bits ("nails"), the largest width for which the IFMA |
There was a problem hiding this comment.
Nails are the unused bits in the word/limb. Digits seems like confusing terminology, I'd suggest just using limbs and saying whether they are 52 or 64-bit.
| // `out` at digits 0 and 2*kHalf and nothing else has been added yet, so | ||
| // subtract them from there. | ||
| Vec<D> st[kSumDigits]; | ||
| MulRec<kHalf + 1>(d, s, t, st); |
There was a problem hiding this comment.
Do we care about power of two kDigits? If so, that's often going to hit the schoolbook codepath because this is always odd.
I think we can either handle odd cases inside Karatsuba, or keep the recursion at kHalf by handling carries separately.
|
|
||
| // Multiplies two kDigits-digit unsigned integers, each digit < 2^52, into | ||
| // 2*kDigits digits. `d` must be a u64 vector descriptor; `out` must have room | ||
| // for 2*kDigits digits. Input digits at or above 2^52 are ignored (as in |
There was a problem hiding this comment.
Actually these inputs are no longer ignored, they must be < 2^52.
| } | ||
| } | ||
|
|
||
| // dst[0..kTotal) -= src[kOffset..kOffset + kCount), zero-extending src, with |
There was a problem hiding this comment.
Behavior of kOffset differs between Add and Sub - out vs src. Should we be consistent in one or the other?
| const Vec<D> sub = Add(sv, borrow); | ||
| const Mask<D> underflow = Lt(dst[i], sub); | ||
| dst[i] = And(Sub(dst[i], sub), mask); | ||
| borrow = IfThenElse(underflow, one, zero); |
There was a problem hiding this comment.
Can't we compute borrow as the upper bit shifted right 63?
| static_assert(kDigits <= 2048, "accumulators may overflow"); | ||
| static constexpr size_t kNumResultDigits = 2 * kDigits; | ||
|
|
||
| template <class D> |
There was a problem hiding this comment.
Probably want HWY_IF_U64_D(D) here.
|
Round 2 in 732c73c, one commit.
Arrays of T, in the API and the test. Karatsuba removed. Two reasons, both from your comments: its temporaries were digits -> limbs, with the 12 unused bits of each 64-bit word called nails, as
Build. BUILD and meson.build did not register the new contrib at all, which Extra coverage while in there: the vector test runs an all-ones round as well Local: clean under |
732c73c to
98504f6
Compare
|
Correction to my previous comment: I did not end up removing Karatsuba. It is back
The test now exercises both paths for every size where For the record, the review points that were about code that no longer exists |
98504f6 to
6fe22ac
Compare
|
Montgomery is in (same single commit, 6fe22ac). Design notes, since you know this space better than I do:
Local: 15/15 (TestAllWideMul, TestAllWideMulLimbs, TestAllMontgomery) on SSE2, SSSE3, SSE4, AVX2 and EMU128, clean under -Wall -Wextra -Wconversion -Wsign-conversion. Still open, and I would rather ask than guess: conversion helpers. Multiplying in Montgomery form needs R^2 mod N to get in and out, and computing that needs modular arithmetic we do not have yet. Do you want that in this PR (a generic ModMul/reduction would be the building block), or is supplying R^2 mod N the caller a reasonable contract for now? I have documented it as the latter. |
…ication
Addresses the second review round, without dropping any of the previous work.
- SchoolbookMul accumulated all 2*kNumLimbs-1 columns at once, which needs an
array of vectors (not allowed on SVE/RVV) and can spill for large kNumLimbs.
It is now Comba's method: for each output limb, sum the products a[i]*b[k-i],
keeping only three vectors live.
- Karatsuba is kept, and now recurses only at kHalf instead of at kHalf + 1, by
tracking the carries of a_lo + a_hi and b_lo + b_hi as separate 0/1 vectors
and folding them back in with branch-free masked adds. All of its temporaries
(z0, z2, s, t, st) are arrays of T in caller-provided scratch, so it no longer
needs arrays of vectors and no longer falls back to schoolbook at every other
level. The scratch size is exposed as WideMul<>::kScratchLimbs; the no-scratch
Mul remains Comba.
- AddInto, AddIntoIf and SubInto all take kOffset indexing the destination, and
SubInto derives the borrow from bit 63 of the wrapped difference instead of
comparing (both operands are below 2^52, so the bit is exact).
- Added Montgomery modular multiplication. MontgomeryN0 computes -N[0]^{-1} mod
2^52 per lane by Newton iteration, and Montgomery<>::Mul is a CIOS (Koc) kernel
that alternates one row of the product with one reduction step, so the working
array is only kNumLimbs + 2 limbs and no separate 2*kNumLimbs product is
needed. It is branch-free: the final conditional subtraction is a per-lane
select, run twice so that moduli above R/2 are covered too, and the borrow
comes from bit 63 of the wrapped difference. Moduli and n0 may differ per lane,
which is the use case (several RSA/ECC instances in parallel). Scratch is
kScratchLimbs * Lanes(d) limbs.
- The API takes arrays of T (TFromD<D>) instead of arrays of Vec<D>, and the
tests do the same, so everything works where vector types are sizeless.
- Renamed digits to limbs (52-bit) and noted that the 12 unused bits of each
64-bit word are what GMP calls nails. Corrected the docs: inputs must be below
2^52 (they are no longer ignored above that), out and scratch must not overlap
the inputs, and the top result limb may be partial.
- Constrained the descriptor with HWY_IF_U64_D, and tightened the accumulator
bound to 1024 limbs with the arithmetic spelled out (2048 could overflow with
the carry from the previous limb).
- BUILD and meson.build now register the contrib (header library and test);
CMakeLists.txt already listed both files.
- Tests: T-based API with no arrays of vectors, an all-ones round in addition to
random inputs, the Karatsuba path checked against the same reference wherever
kScratchLimbs is nonzero, static_asserts for the bit-width alias, and for
Montgomery an independent reference (modular multiplication by binary
double-and-add, with no Montgomery and no R) over random operands, 0, 1 and
n-1, a modulus above R/2, and a modulus whose low limb is 1 so n0 is -1.
- Added MontgomeryR2, which computes R^2 mod n - the constant needed to enter
Montgomery form - as 104*kNumLimbs modular doublings from 1, branch-free in the
same style as the reduction: the carry out of the top limb and the borrow of the
conditional subtraction become one per-lane select. With it callers need no
modular arithmetic of their own:
ToMontgomery(x) = Montgomery<>::Mul(x, r2, ...)
FromMontgomery(y) = Montgomery<>::Mul(y, one, ...) // one = {1, 0, ...}
Scratch is 2*kNumLimbs limbs per lane (MontgomeryR2ScratchLimbs).
- Tests for that path: R^2 against the reference, and the round trip
x -> Mul(x, R2) -> Mul(., one) == x, for every size where it applies. The
kernels were checked separately by a standalone sweep over 1, 2, 4 and 8-limb
moduli, which also showed that the failures I hit were in the test's operand
reduction rather than in the code under test: a single conditional subtraction
is not enough for small moduli, and a reference that assumes a < n is circular.
The operands are now reduced properly.
Local: compiles with -Wall -Wextra -Wconversion -Wsign-conversion, and 20/20
tests pass (TestAllWideMul, TestAllWideMulLimbs, TestAllMontgomery,
TestAllMontgomeryForms) on SSE2,
SSSE3, SSE4, AVX2 and EMU128. The Montgomery test caught two real bugs in the
first draft of the kernel: a dropped carry out of the cancelled low limb, and
unmasked subtraction limbs leaking bits above 2^52.
6fe22ac to
2f44c3f
Compare
|
R^2 is in (same single commit, 2f44c3f), so callers no longer have to supply it:
Worth recording, because it cost me a detour: the first version of these tests failed, and the failures were in the test rather than in the kernel. It was reducing operands with a single conditional subtraction, which is only valid when n exceeds half the limb range, and my first replacement reference was circular - it assumed a < n while being used to establish a < n. I ended up checking the kernels separately with a standalone sweep over 1, 2, 4 and 8-limb moduli, 300 random cases each, which passed and pointed back at the test. Operand reduction is proper now, and the suite is 20/20 (it was 15/15). |
2 participants
Adds
hwy/contrib/multiprec: a fixed-width unsigned multiply layer built on the MulAdd52 ops from #3365.Design
MulAdd52Lo/Hiare exact: for inputs < 2^52 the two together return the full 104-bit product. They exist on every target (nativeHWY_NATIVE_MULADD52, generic fallback otherwise), so there is no separate code path.iof each operand is passed in its own vector, and the lanes of those vectors hold independent multiplications — the crypto use case (several RSA/ECC instances in parallel).API
WideMul<kDigits>::Mul(d, a, b, out)— SoA vector form;outhas2*kDigitsdigits.WideMulBits<kBits>— bit-width form (kDigits = ceil(kBits/52)).WideMulLimbs<kBits>(a, b, out)— scalar convenience over little-endian 64-bit limbs (e.g. u128/u192/u256).Tests
WideMul<kDigits>forkDigits = 1..4, checked against a scalar reference on all targets.WideMulLimbs<128/192/256>checked against an independent ripple-carry reference, plus maximum inputs.Dependency
Depends on #3365 (
MulAdd52); this branch is stacked onfeat/muladd52, so the diff includes it until #3365 merges. Part of #3352.