bigint-mul

bigint-mul multiplies arbitrary-precision integers on the GPU by number-theoretic transform, loadable through kernels, returning an exact Python int. The reference baseline is GMP's mpz_mul via gmpy2, which it overtakes near one million bits and beats by up to 2.8x at four million.

Multiplying numbers with millions of digits is the substrate of computational number theory and record arithmetic, and it has been CPU territory: GMP switches through Karatsuba, Toom-Cook, and FFT multiplication as operands grow. This kernel moves the FFT regime onto the GPU: operands become limb sequences, the product becomes a convolution computed as a pointwise product of number-theoretic transforms, and carries settle in logarithmic depth, so a multiplication that costs GMP ten milliseconds returns in three and a half.

Two 4,194,304-bit integers multiplied live: limb waveforms, their number-theoretic spectra, the pointwise product, the inverse transform, and the carry prefix-scan settling the 2,525,223-digit result, digit-for-digit equal to GMP; 3.5 ms against GMP's 9.8 ms on the same host.

Usage

import torch
from kernels import get_kernel

bm = get_kernel("phanerozoic/bigint-mul", version=1, trust_remote_code=True)

a = 3 ** 2_000_000
b = 5 ** 2_000_000
product = bm.mul(a, b)          # exact Python int
assert product == a * b

version selects the release branch; trust_remote_code is required by kernels for publishers without the trusted-publisher mark.

API

Operands up to 2^25 base-2^16 limbs, about half a billion bits each. Signed operands and zero are handled.

Method

Operands are split into base-2^16 limbs, making the product a linear convolution of the limb sequences. The convolution is a pointwise product of forward NTTs modulo the primes 15 * 2^27 + 1 and 27 * 2^26 + 1, both congruent to 1 modulo a large power of two, giving transform lengths to 2^26. Each coefficient is bounded by n * (2^16 - 1)^2 < 2^58, below the product of the primes, so CRT determines it uniquely.

Carry normalization is a prefix scan, not an iteration. A carry entering a run of maximal limbs advances one position per pass, so iterating "keep the low half, carry the rest" costs a pass per limb. Four coarse passes reduce every coefficient to at most 2^16 + 2, leaving single-bit carries; position i generates a carry when its high bit is set and propagates one when its low half is 0xffff, and composing generate-propagate pairs is associative, so a Hillis-Steele scan resolves all carries in logarithmic depth.

Limbs cross PCIe packed two per 32-bit word, an integer's own little-endian byte layout, so neither direction converts format.

Measured

Against GMP 6.x via gmpy2 and CPython's native multiplication, same host. Square operands of the stated bit length, best of several runs, complete call including transfer and integer assembly.

Crossover against GMP is near 10^6 bits; below it, launch latency and the host round trip dominate a job GMP finishes in microseconds. Phase split at 16.8 x 10^6 bits: transform 2.6 ms, operand export 3.8 ms, device-to-host 0.3 ms, integer assembly 3.0 ms.

Verification

The suite compares against CPython's native integer multiplication at every size, digit for digit: random operands from 17 to 2^20 bits, asymmetric operand sizes, exact powers of two, all four sign combinations, and the all-maximal-limb pattern that is the worst case for carry propagation.

Requirements and limits

• NVIDIA GPU with compute capability 8.0+.
• Operands to 2^25 base-2^16 limbs (~5 x 10^8 bits) each.
• GMP remains faster below roughly one million bits; this kernel is for the regime where a single multiplication is already milliseconds of CPU time.

References

Schoenhage-Strassen FFT multiplication; number-theoretic transforms over word-size primes; GMP (mpz_mul) as the reference implementation.

License

Apache-2.0.