Even the simplest better-than-schoolbook (O(n^2)) algorithms like Karatsuba are only useful in practice for large n
. But what is n
? It's not single bits, and it's not decimal digits. (Posting this tangent as requested in comments.)
Software implementations of an extended-precision multiply algorithm work in integer chunks as wide as the hardware provides. On a 64-bit CPU, that's usually 64x64 => 128-bit integer multiplication, e.g. the x86-64 mul
instruction. (@fgrieu's answer has more detail on this, including the term "limb" for such a chunk.)
That fixed-width CPU instruction runs in fixed time (regardless of the value on most CPUs; division is the only instruction that's slow enough to justify variable latency in a modern pipelined CPU, and in the most recent x86-64 CPUs even it's constant). e.g. on modern Intel and AMD CPUs, mul r64
or mulx
have a throughput of 1 per cycle and a latency of 3 to 4 cycles (for the low and high halves of the output, respectively: https://www.uops.info/html-instr/MUL_R64.html).
Hardware doesn't "know" it's doing one big multiply, it's just doing each fixed-width part separately. Hardware can easily be parallel (adding partial products) if you can throw enough transistors at the problem. HW multipliers in CPUs use the Dadda tree design. This is simpler than doing 63 additions of shifted versions of the other 64-bit input (or 0
where this input has a 0 bit) using normal adders: carry propagation can be deferred. Hardware tricks like that are AFAIK unrelated to any of the sub-N^2 algorithmic tricks.
Such a multiply instruction, and add-with-carry, are the building blocks for schoolbook multiplication's O(n^2) time complexity. e.g. 128-bit multiplication (producing a 128-bit result) takes 3 multiplies on x86-64: https://godbolt.org/z/qBAbfQ. To also produce the high half, all of those multiplies would have to be "full" 64x64=>128 instead of only 64x64 => 64 for the low x high and high x low cross products, and we'd need to do the high x high product, for a total of 4 mul
instructions.
e.g. this SO answer shows 32x32 => 64-bit multiply using 16-bit x86 so each input is 2 limbs, and the output is 2+2 = 4 limbs, requiring 2*2 = 4 multiplies of 16x16 => 32 bits each. Exactly the same pattern would apply for 64x64 => 128 on a 32-bit machine, or 128x128 => 256 on a 64-bit machine.
Since that building block is opaque to software, and/or shuffling individual bits around would be much more expensive than it's worth, n
is only 64 for 4096-bit integer multiply.
To allow better instruction-level parallelism (letting superscalar CPUs do the same work in less time) and reducing overhead of mov
instructions, Intel introduced (in Broadwell) the ADX extension that allows two parallel dependency chains of add-with-carry. This whitepaper shows the advantages it gives for small problems (like 512-bit x 512-bit multiplication (8 x 8 limbs)).
For floating-point, an FP multiplier involves an integer multiplier for the 53x53-bit => 53-bit correctly rounded mantissa (the most significant 53 bits of the full integer product) plus hardware to add the exponents, and check for / handle overflow / underflow and NaN. See Why does Intel's Haswell chip allow floating point multiplication to be twice as fast as addition? for some info about how FP ALUs are designed, and the barely-related question of why Intel made the design choices they did in Haswell and Skylake.
To get extra FP precision, one technique is so-called "double-double": wide mantissa using two double
s, but only the exponent from one of them. Using that only takes a handful of double-precision math operations, like 6 to 20 depending on which operation and whether FMA (fused multiply-add without intermediate rounding) is available. The relevant width is n=2 doubles, not n=36 decimal digits. (And IEEE FP is a binary format, not decimal, although there are decimal FP formats that exist, with some CPUs even having hardware support for them, such as PowerPC.)
Note that a SIMD multiplier just replicates that for each SIMD element. double-double can SIMD efficiently if you store separate vectors of lo / hi halves so you don't need to shuffle to line up the corresponding halves of a single number. e.g. this Q&A.
Other extended-precision number representations
You could store numbers as an array of bytes, each byte holding a single decimal digit. But that's pretty terrible. Historically, it was not uncommon to use a simplistic format like that, especially for a score counter in a game that gets printed on screen in decimal format constantly. Or BCD (2 decimal digits per 8-bit byte, each in a separate 4-bit nibble).
But this is pretty bad, especially for multiplying numbers stored in this format, because then n
becomes large and complexity scales with N^2 (for the simple schoolbook algorithm).
@davidbak commented:
w.r.t. "nobody uses decimal digits as an extended-precision format" - is that true? I know there used to be implementations of multi precision integer arithmetic that used the largest power of 10 that would fit in a word as the base - e.g., 10^9 for 32-bit machines. Made conversions to<->from a human-readable base 10 notation much easier and cost only a "reasonable" overhead for some definition of reasonable that might depend on your use case. Is that not done anymore? (Although strictly speaking those aren't decimal digits, just power-of-ten digits...)
Indeed, larger powers of 10 could be sane when you need frequent conversion to/from a decimal string, or multiply/divide by powers of 10. But then a 36-digit number is 4 chunks of 9, not 36 chunks of 1. e.g. one use-case was printing the first 1000 decimal digits of Fib(10^9)
(x86-64 asm code-golf) where it's handy to have right shift by 1 limb be division by a power of 10, and for conversion to decimal to only need to consider the current limb, converting that to 9 decimal digits without having to do extended-precision division where the remainder depends on all higher bits.
See also this code-review answer about an implementation based on single decimal digits. I included some details about what CPython does, and some other links. It's not rare for beginners to come up with that as an idea, but non-toy libraries use at least 10^9 as the base for "limbs", unless we're talking about BCD.
Or more commonly binary extended precision using all 32 bits per 32-bit integer, or sometimes only 2^30 to leave room for high-level language handling of carry in/out (like in CPython) without access to an asm carry flag.
Another advantage of leaving some spare bits per limb is to allow deferred carry normalization, making SIMD for addition efficiently possible. See @Mysticial's answer on Can long integer routines benefit from SSE?. Especially for extended-precision addition, leaving some slack in each limb is actually interesting if you design around that format with awareness of when to normalize as an extra step. (@Mysticial is the author of y-cruncher and also works on Prime95; he implemented its use of FP-FMA to take advantage of the FP mantissa multipliers for bit-exact integer work.)
That answer also points out that "really large bignum" multiplications can be done as an FFT.
Normally (with standard techniques) it's very hard to take advantage of SIMD for extended-precision; within one operation, there's a serial dependency between each element: you don't know if there's carry-in to this element until you process the previous element (for addition).
For multiplication, it's usually even worse: SIMD doesn't usually have very wide multipliers, and with the result being twice as wide as the inputs it's a problem where to put them.
The amount of work done by one building block should be measured as the "product bits" you compute per cycle, e.g. 64x64 => 128-bit full multiply does 64x64 = 4096 units of work. But a 4x 32x32=>64-bit SIMD multiply (like AVX2 vpmuludq
) does 32^2
= 1024 units of work per element, with 4 elements, for a total of 4096 units of multiply work. And it leaves more adding of partial products not done. So even in theory, ignoring other factors, AVX2 vpmuludq
on a 256-bit vector is break-even with scalar.
AVX512 has 64x64 => 64-bit multiply (but still no way to get the upper-half of the full result so it's no more helpful for BigInteger than 32x32 => 64, I think). AVX512IFMA more directly exposes what the FP mantissa multipliers can do, providing separate low and high half 52x52 => 104-bit multiply.
(Other SIMD integer multiply instructions like vpmulld
that do 32x32 => 32-bit usually decode to two separate uops for the vector-ALU ports, so they can use the same per-element multipliers as FP mantissas. But those multipliers are only 52x52 or 24x24-bit. Making them wider would cost significantly more for these wide SIMD ALUs, and only help the fairly rarely used SIMD-integer multiply instructions.)