Device Functions#

The following are __device__ member functions of a descriptor’s bigint type (see Types). Its double-width counterpart, bigint_wide, carries almost none of these directly—it exposes only construction, access, div_rem, and reduce_montgomery—and otherwise carries the lo/hi result of mul_wide, operator*, and square, which Barrett reduction, Montgomery reduction, and its own div_rem then consume.

For a warp-executed instance (TPI > 1), every cooperating lane must call these functions together, in lockstep, with identical arguments, including the same address for any pointer argument. Some operations are available only for thread execution; see Warp Execution Coverage.

Unless an entry says otherwise, both operands of a binary operation must be the same type—same bit width, threads per instance, and error policy—with no implicit conversion between widths. div_rem is the sole exception; see Operand Widths.

Construction and Access#

The semantics below are written in terms of limb[i], the i-th 32-bit limb, so that the value is limb[0] + limb[1] * 2^32 + ... + limb[num_limbs-1] * 2^(32 * (num_limbs-1)). Each assignment holds for every i from 0 to num_limbs - 1.

__device__ bigint::bigint()#

Default-constructs a zero-initialized value. Useful for declaring an output such as a quotient, remainder, or struct member before it is assigned by another operation (for example bigint qv, rv; before calling div_rem).

limb[i] = 0

__device__ explicit bigint::bigint(uint32_t value)#

Constructs from a 32-bit scalar.

limb[0] = value
limb[i] = 0                    for i > 0

__device__ explicit bigint::bigint(uint64_t value)#

Constructs from a 64-bit scalar. Deleted when num_limbs == 1 (a single 32-bit limb cannot hold a 64-bit value).

limb[0] = value & 0xffffffff
limb[1] = value >> 32
limb[i] = 0                    for i > 1

A value can also be loaded from memory on construction and written back with store. Under warp execution every lane passes the same pointer, and each lane transfers only the limbs it owns.

__device__ explicit bigint::bigint(const uint32_t *ptr)#

Loads num_limbs limbs starting at ptr.

limb[i] = ptr[i]

__device__ bigint::bigint(const uint32_t *ptr, unsigned int index)#

Loads instance index of a batched array at ptr, which stores num_limbs limbs per instance. index addresses a big-integer-sized element, not a limb.

limb[i] = ptr[index * num_limbs + i]

__device__ void bigint::store(uint32_t *ptr) const#

Stores num_limbs limbs starting at ptr.

ptr[i] = limb[i]

__device__ void bigint::store(uint32_t *ptr, unsigned int index) const#

Stores to instance index of a batched array at ptr, addressed as in the indexed constructor.

ptr[index * num_limbs + i] = limb[i]

__device__ uint32_t bigint::to_uint32() const#

Returns the low 32 bits of the value. Under warp execution the limb is broadcast, so every lane of the group receives the same result.

result = limb[0]

__device__ uint64_t bigint::to_uint64() const#

Returns the low 64 bits of the value, broadcast to every lane in the same way.

result = limb[0] + limb[1] * 2^32      if num_limbs > 1
result = limb[0]                       if num_limbs == 1

__device__ uint32_t &bigint::operator[](unsigned int index)#

__device__ const uint32_t &bigint::operator[]( unsigned int index, ) const#

Direct access to limb index, least significant first. The non-const form returns a mutable reference, so a limb can be written in place. Here index counts limbs, unlike the index parameter of the batched constructor and store, which counts whole big integers.

result = limb[index]

Available only for Thread() execution (tpi == 1). Under Warp() no single lane holds the whole instance, so this overload does not exist and a call fails to compile.

bigint_wide provides the same construction and access surface over 2 * num_limbs limbs, where limb[i] spans both halves: lo holds limb[0] through limb[num_limbs - 1] and hi holds the rest. It has no scalar constructors and no to_uint32 or to_uint64.

__device__ bigint_wide::bigint_wide()#

__device__ explicit bigint_wide::bigint_wide(const uint32_t *ptr)#

__device__ explicit bigint_wide::bigint_wide( const uint32_t *ptr, unsigned int index, )#

__device__ void bigint_wide::store(uint32_t *ptr) const#

__device__ void bigint_wide::store( uint32_t *ptr, unsigned int index, ) const#: Zero-initialization, pointer and batched loads, and the matching stores, each behaving as the bigint form of the same name with 2 * num_limbs limbs per instance.

__device__ uint32_t &bigint_wide::operator[](unsigned int index)#

__device__ const uint32_t &bigint_wide::operator[]( unsigned int index, ) const#

Direct access to limb index of the combined value, so one index reaches into lo or hi without naming either. Restricted to tpi == 1 like the bigint form.

result = lo[index]                  if index < num_limbs
result = hi[index - num_limbs]      if index >= num_limbs

Basic Arithmetic#

Unsigned arithmetic on the full bit_width value. The three groups below differ mainly in what happens when a result does not fit: addition, subtraction, and scalar multiplication wrap modulo 2^bit_width; mul_wide, operator*, and square widen instead, returning the exact product as a bigint_wide; division cannot overflow but can fail, so it reports a bigint_error status.

Addition and Subtraction#

Fixed-width addition and subtraction, both bigint-bigint and bigint-scalar. Results wrap modulo 2^bit_width, so the carry or borrow out of the most significant limb is discarded and a smaller minus a larger value gives the wrapped result rather than a negative one. add_scalar(limb) and this + limb are two spellings of one function, as are sub_scalar(limb) and this - limb: each operator calls the named function. Both behave as the bigint forms with other equal to limb.

__device__ bigint bigint::add_scalar(uint32_t limb) const#

__device__ bigint bigint::operator+(uint32_t limb) const#

__device__ bigint bigint::operator+(const bigint &other) const#

result = (this + other) mod 2^bit_width

__device__ bigint &bigint::operator+=(uint32_t limb)#

__device__ bigint &bigint::operator+=(const bigint &other)#

this = (this + other) mod 2^bit_width

__device__ bigint bigint::sub_scalar(uint32_t limb) const#

__device__ bigint bigint::operator-(uint32_t limb) const#

__device__ bigint bigint::operator-(const bigint &other) const#

result = (this - other) mod 2^bit_width

__device__ bigint &bigint::operator-=(uint32_t limb)#

__device__ bigint &bigint::operator-=(const bigint &other)#

this = (this - other) mod 2^bit_width

Multiplication#

The product of two bit_width-bit values needs 2 * bit_width bits. mul_wide, operator*, and square return all of it as a bigint_wide, whose lo and hi halves are exactly what mul_low and mul_high return individually. Multiplication by a uint32_t scalar is fixed-width instead, wrapping modulo 2^bit_width like addition.

__device__ bigint bigint::mul_low(const bigint &other) const#

The low bit_width bits of the product. Costs about half of mul_wide.

result = (this * other) mod 2^bit_width

__device__ bigint bigint::mul_high(const bigint &other) const#

The high bit_width bits of the product. This is not faster than a full mul_wide; if the low half is also required, use mul_wide instead.

result = (this * other) / 2^bit_width

__device__ bigint_wide bigint::mul_wide(const bigint &other) const#

__device__ bigint_wide bigint::operator*(const bigint &other) const#

The full 2 * bit_width-bit product, the only bigint-bigint multiply that keeps every bit instead of wrapping.

a * b and a.mul_wide(b) are the same function. operator* does nothing but call mul_wide, so the two are interchangeable and cost the same.

The wide result is the usual input to a reduction step: pass it to reduce_montgomery, reduce_barrett, or the double-width div_rem to bring it back to a single-width value modulo your modulus.

result.lo = (this * other) mod 2^bit_width
result.hi = (this * other) / 2^bit_width

__device__ bigint_wide bigint::square() const#

The full 2 * bit_width-bit square. square is faster than mul_wide(this) only for TPI == 1; for warp-executed instances (TPI > 1) the two cost the same.

result.lo = (this * this) mod 2^bit_width
result.hi = (this * this) / 2^bit_width

__device__ bigint bigint::mul_scalar(uint32_t limb) const#

__device__ bigint bigint::operator*(uint32_t limb) const#

Fixed-width multiplication by a 32-bit scalar. mul_scalar(limb) and this * limb are two spellings of one function: operator* calls mul_scalar.

result = (this * limb) mod 2^bit_width

__device__ bigint &bigint::operator*=(uint32_t limb)#

this = (this * limb) mod 2^bit_width

Division#

Division of unsigned values, producing a quotient and a remainder together. On success the two outputs satisfy the division identity below, which is what pins them down; unless bigint_error::success is returned, treat q and r as unspecified.

this = q * divisor + r        with 0 <= r < divisor

__device__ bigint_error bigint::div_rem( const bigint &divisor, bigint &q, bigint &r, ) const#

Returns bigint_error::success, bigint_error::divide_by_zero, or bigint_error::quotient_overflow.

q = this / divisor
r = this mod divisor

template<unsigned int DIVISOR_LIMBS, unsigned int QUOTIENT_LIMBS> __device__ bigint_error bigint::div_rem( const bigint_impl<DIVISOR_LIMBS, TPI, ErrorPolicy> &divisor, bigint_impl<QUOTIENT_LIMBS, TPI, ErrorPolicy> &q, bigint_impl<DIVISOR_LIMBS, TPI, ErrorPolicy> &r, ) const#

Division by a divisor narrower than the dividend: the remainder width follows the divisor, and the quotient width is chosen through the type of q. Returns the same status codes as the same-width overload, with bigint_error::quotient_overflow reported when the true quotient needs more than QUOTIENT_LIMBS limbs.

Precondition

QUOTIENT_LIMBS must be num_limbs or half of it, enforced by a static_assert. DIVISOR_LIMBS must be half of num_limbs; that one is enforced only by which specializations the library provides, so an unsupported divisor width compiles and then fails to link. See Operand Widths.

q = this / divisor           q is QUOTIENT_LIMBS limbs wide
r = this mod divisor         r is DIVISOR_LIMBS limbs wide

__device__ bigint bigint::operator/(const bigint &divisor) const#

The quotient alone. Failures go through the configured on-error policy rather than coming back as a status (see Operators), and the returned value is zero.

result = this / divisor

template<unsigned int DIVISOR_LIMBS, unsigned int QUOTIENT_LIMBS> __device__ bigint_error bigint_wide::div_rem( const bigint_impl<DIVISOR_LIMBS, TPI, ErrorPolicy> &divisor, bigint_impl<QUOTIENT_LIMBS, TPI, ErrorPolicy> &quotient, bigint_impl<DIVISOR_LIMBS, TPI, ErrorPolicy> &remainder, ) const#

Divides a double-width bigint_wide value by divisor, so the dividend is the full 2 * bit_width-bit value rather than a single-width one. This is the path that reduces a mul_wide or square result without Barrett or Montgomery setup.

Precondition

DIVISOR_LIMBS and QUOTIENT_LIMBS must both be num_limbs, that is, half the double-width dividend. The static_assert also admits a full 2 * num_limbs quotient, but that specialization is provided only for num_limbs == 8. See Operand Widths.

quotient  = this / divisor
remainder = this mod divisor

Comparison#

Three-way comparison: returns a value less than, equal to, or greater than zero when this is respectively less than, equal to, or greater than other. Values are compared as unsigned magnitudes. Under warp execution every lane of the group receives the same result, so it is safe to branch the whole group on it.

__device__ int bigint::compare(const bigint &other) const#

__device__ int bigint::compare(uint32_t other) const#

__device__ int bigint::compare(uint64_t other) const#

Semantics

result <  0   if   this <  other
result == 0   if   this == other
result >  0   if   this >  other

Relational and equality operators built on compare(). Scalar operands must be uint32_t or uint64_t, so write an unsigned literal (base == 0u) rather than a plain 0.

__device__ bool bigint::operator==(const bigint &other) const#

result = (this->compare(other) == 0)

__device__ bool bigint::operator!=(const bigint &other) const#

result = (this->compare(other) != 0)

__device__ bool bigint::operator<(const bigint &other) const#

result = (this->compare(other) < 0)

__device__ bool bigint::operator<=(const bigint &other) const#

result = (this->compare(other) <= 0)

__device__ bool bigint::operator>=(const bigint &other) const#

result = (this->compare(other) >= 0)

__device__ bool bigint::operator>(const bigint &other) const#

result = (this->compare(other) > 0)

Each of the six also takes a uint32_t or uint64_t right-hand operand as a member overload, with the same semantics against a scalar. The mirrored spellings with the scalar on the left, such as 2u < x, are provided by free functions that swap the operands and call the member back, so the sense of the comparison is preserved.

Bitwise and Shift#

Bitwise AND, OR, XOR, and complement over the full bit_width value, applied limb by limb. The compound forms assign the result into this.

__device__ bigint bigint::operator&(const bigint &other) const#

result.limb[i] = this->limb[i] & other.limb[i]

__device__ bigint bigint::operator|(const bigint &other) const#

result.limb[i] = this->limb[i] | other.limb[i]

__device__ bigint bigint::operator^(const bigint &other) const#

result.limb[i] = this->limb[i] ^ other.limb[i]

__device__ bigint bigint::operator~() const#

result.limb[i] = ~this->limb[i]

__device__ bigint &bigint::operator&=(const bigint &other)#

this->limb[i] = this->limb[i] & other.limb[i]

__device__ bigint &bigint::operator|=(const bigint &other)#

this->limb[i] = this->limb[i] | other.limb[i]

__device__ bigint &bigint::operator^=(const bigint &other)#

this->limb[i] = this->limb[i] ^ other.limb[i]

Logical left and right shift by shift bits. Bits shifted out are discarded and no bits are shifted back in at the opposite end. A shift of bit_width or more is well defined and yields zero, rather than being undefined as it is for C++ built-in integers.

__device__ bigint bigint::operator<<(unsigned int shift) const#

result = (this * 2^shift) mod 2^bit_width

__device__ bigint bigint::operator>>(unsigned int shift) const#

result = this / 2^shift

__device__ bigint &bigint::operator<<=(unsigned int shift)#

this = (this * 2^shift) mod 2^bit_width

__device__ bigint &bigint::operator>>=(unsigned int shift)#

this = this / 2^shift

__device__ int bigint::clz() const#

Count of leading zero bits over the full bit_width value.

result = bit_width - floor(log2(this)) - 1      if this != 0
result = bit_width                              if this == 0

Modular Arithmetic#

Three routes to arithmetic modulo a modulus, differing in where the cost of reduction falls rather than in the results they produce. Direct needs no setup at all: one call in, a reduced result out. Barrett precomputes a reciprocal once per modulus, making repeated reductions against that modulus cheaper. Montgomery changes how values are represented so that multiplication needs no division at all, which pays off once several multiplications share a modulus. For choosing between them, see Choosing Modular Reduction.

Direct#

Single-call modular operations, each returning a result already reduced to the range [0, modulus).

The two groups differ in what they expect of their inputs. add_mod and sub_mod require both operands to be reduced already, which is what lets them bring the result back into range with at most one correction by the modulus rather than a division. mod, mul_mod, pow_mod, and inv_mod reduce as part of the operation and are correspondingly more expensive.

__device__ bigint bigint::mod(const bigint &modulus) const#

__device__ bigint bigint::operator%(const bigint &modulus) const#

The remainder of the same division performed by div_rem, with the quotient discarded. a.mod(b) and a % b are the same operation. A zero modulus yields zero and is reported through the configured on-error policy.

result = this mod modulus

__device__ bigint bigint::add_mod( const bigint &other, const bigint &modulus, ) const#

Precondition

this < modulus and other < modulus. This is not checked: with unreduced operands the result may fall outside [0, modulus) or be incorrect.

result = (this + other) mod modulus

__device__ bigint bigint::sub_mod( const bigint &other, const bigint &modulus, ) const#

Precondition

this < modulus and other < modulus. This is not checked: with unreduced operands the result may fall outside [0, modulus) or be incorrect.

result = (this - other) mod modulus

__device__ bigint bigint::mul_mod( const bigint &other, const bigint &modulus, ) const#

Computed as a mul_wide followed by the double-width div_rem, keeping the remainder and discarding the quotient — the same reduction mod performs, but on the 2 * bit_width-bit product rather than a single-width value. This is not Montgomery multiplication; when multiplying repeatedly against one modulus, prefer the Montgomery path.

Unlike add_mod and sub_mod, the operands need not be reduced beforehand. A zero modulus yields zero, reported through the configured on-error policy.

result = (this * other) mod modulus

__device__ bigint bigint::pow_mod( uint64_t exponent, const bigint &modulus, ) const#

__device__ bigint bigint::pow_mod( const bigint &exponent, const bigint &modulus, ) const#

Modular exponentiation by the right-to-left binary method. The base is reduced first, so it need not be reduced by the caller. A zero exponent yields 1 mod modulus, and a zero modulus yields zero, reported through the configured on-error policy.

Warning

The loop bound depends on the value of the exponent, so pow_mod is not a constant-time implementation.

result = this^exponent mod modulus

__device__ bigint_error bigint::inv_mod( const bigint &modulus, bigint &out, ) const#

The modular inverse. Supports both odd and even moduli, and both single-thread and warp-cooperative instances. Returns bigint_error::success, bigint_error::inv_mod_not_invertible, bigint_error::inv_mod_invalid_modulus, or bigint_error::inv_mod_zero_input.

Precondition

modulus > 1, and this < modulus when modulus is valid. Unlike the other preconditions in this section these are checked, and a violation is reported through the returned status rather than producing a silently wrong result.

out = this^-1 mod modulus        so that (this * out) mod modulus == 1

Barrett#

Reduction of a double-width value by a fixed modulus, split into a one-time setup_barrett per modulus and a reduce_barrett per value.

__device__ static void bigint::setup_barrett( const bigint &den, bigint &approx_out, int &den_clz_out, )#

__device__ static void bigint::reduce_barrett( bigint &rem_out, const bigint_wide &num, const bigint &den, const bigint &approx, int den_clz, )#

Barrett-reduction helpers: setup_barrett precomputes a reciprocal approximation and leading-zero count for a fixed denominator den, and reduce_barrett uses them to compute num mod den for a double-width num, avoiding recomputation of the approximation across repeated reductions by the same denominator.

Precondition

setup_barrett requires a non-zero den. reduce_barrett requires num.hi < den, and den, approx, and den_clz must all come from the same setup_barrett call.

setup_barrett:    approx_out, den_clz_out = reusable state derived from den
reduce_barrett:   rem_out = num mod den

approx and den_clz are an internal representation of the reciprocal, not a value to compute or interpret yourself; treat them as opaque state that is set up once per denominator and passed unchanged to every reduction against it.

Both functions return void, so unlike div_rem and inv_mod there is no status to inspect at the call site; errors are only reported through the configured on-error policy (see Operators):

setup_barrett reports bigint_error::barrett_divide_by_zero, leaving approx_out and den_clz_out unspecified.
reduce_barrett reports bigint_error::barrett_input_invariant_violated, leaving rem_out unspecified.

Under the default OnErrorNone a failure is therefore entirely silent, leaving a wrong value in approx_out or rem_out; use OnErrorTrap or OnErrorPrintTrap when it needs to be detectable.

Montgomery#

The Montgomery domain represents a value x as (x * R) mod modulus, where R = 2^(num_limbs * 32). Multiplication in this representation needs no division, so the cost of entering and leaving the domain is worth paying once several multiplications share a modulus:

a_mont = a.to_montgomery(modulus)
b_mont = b.to_montgomery(modulus)
c_mont = a_mont.mul_montgomery(b_mont, modulus)
c      = c_mont.from_montgomery(modulus)         c == (a * b) mod modulus

Precondition

Every function in this section requires an odd modulus; see montgomery_modulus for why. This is not checked, and an even modulus yields incorrect results.

__device__ bigint bigint::to_montgomery( const montgomery_modulus &modulus, ) const#

Converts into the Montgomery domain.

Precondition

this < modulus. This is not checked, and an unreduced input yields an incorrect result.

result = (this * R) mod modulus

__device__ bigint bigint::from_montgomery( const montgomery_modulus &modulus, ) const#

Converts out of the Montgomery domain (Montgomery reduction / REDC).

result = (this * R^-1) mod modulus

__device__ bigint bigint::mul_montgomery( const bigint &other, const montgomery_modulus &modulus, ) const#

Montgomery multiplication. The result is in the Montgomery domain.

Precondition

Both this and other must already be in the Montgomery domain, as produced by to_montgomery or a previous mul_montgomery. This is not checked, and mixing domains yields an incorrect result.

result = (this * other * R^-1) mod modulus

__device__ bigint bigint_wide::reduce_montgomery( const montgomery_modulus &modulus, ) const#

Montgomery-reduces an already double-width value (for example the output of mul_wide or square) directly against modulus, equivalent to from_montgomery applied to a double-width input without a preceding multiply.

result = (this * R^-1) mod modulus        this is 2 * bit_width bits wide