bignum

Batched arbitrary-precision integer arithmetic on the GPU, loadable through kernels: addition, subtraction, multiplication, signed comparison, and modular reduction over a batch of wide integers, exact at every width. The reference baseline is Python's native integer arithmetic, compared digit for digit, with GMP as the throughput yardstick.

Exact algorithms over the integers, Hermite and Smith normal forms, lattice basis reduction, integer relation detection, carry coefficients that grow past 64 bits and are not decomposable by CRT the way determinants are. They need wide arithmetic on the device itself, and GPUs have no such type. This kernel supplies it as a batch primitive: thousands of 8192-bit values computed at once, so the wide-integer inner loops of exact linear algebra run where the rest of the computation lives.

4,096 lanes of 2048-bit a*b mod M in flight: 2.5 million exact modular multiplies per second, 2.6x one GMP core, every lane digit-for-digit equal to GMP.

Usage

import torch
from kernels import get_kernel

bn = get_kernel("phanerozoic/bignum", version=1, trust_remote_code=True)

a = bn.to_device([3**400, -(5**300), 2**1000])
b = bn.to_device([7**200, 11**150, 1], limbs=a.size(0))

bn.from_device(bn.add(a, b))
bn.from_device(bn.mul(a, b))
bn.compare(a, b)                       # -1, 0 or 1 per lane
bn.from_device(bn.mod(a, (1 << 127) - 1))

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

API

Method

A batch of N integers, each L limbs of 64 bits, is an [L, N] int64 tensor: limb i of value t sits at i * N + t, so consecutive threads touch consecutive addresses on every limb access and the carry chains that are sequential within one value amortize across the batch. Values are little-endian two's complement at a fixed width, up to 128 limbs (8192 bits).

Addition, subtraction, and comparison run directly on the two's-complement representation, one thread per value. Multiplication converts to sign and magnitude, accumulates schoolbook partial products through __umul64hi, and reapplies the sign. Modular reduction is Barrett: q = floor(a * mu / 2^(64 * La)) with mu precomputed on the host, then a - q * M corrected by at most a few subtractions. Sizing mu to the operand width rather than the modulus width keeps the quotient estimate valid for any dividend below 2^(64 * La), so wide values reduce against a small modulus without a narrowing step. Working buffers live in a host-allocated global scratch region indexed by thread.

Measured

Batched a*b mod M at 2048 bits, 4,096 lanes, against GMP via gmpy2 on one CPU core of the same host, results compared digit for digit:

The batch is the regime: per-op latency favors GMP, aggregate throughput and residency next to tensors favor the device.

Verification

Every operation is compared against Python's native integer arithmetic, digit for digit: round-trip packing at the signed word boundaries, addition and subtraction at 8 to 700 bits, multiplication across all four sign combinations, comparison at the boundaries, and reduction against moduli of 61 to 256 bits with dividends far exceeding the modulus squared.

Requirements and limits

• NVIDIA GPU with compute capability 8.0+.
• Fixed width per batch, up to 128 limbs (8192 bits); mod/mulmod take one shared modulus per call.
• Single-value latency is not the target; batch width is what fills the device.

References

Barrett, "Implementing the Rivest Shamir and Adleman public key encryption algorithm on a standard digital signal processor" (1986); GMP as the CPU reference implementation.

License

Apache-2.0.