rns_engine

Fast, exact integer arithmetic on a four-rail Residue Number System.

No floating point. No approximation. Arithmetic is exact modulo the engine's fixed dynamic range.

Recent Colab benchmark: https://colab.research.google.com/drive/18RIQSLf4Vf5xUnwEmDwYNPJDLSM_Stqa?usp=sharing

What it does

rns_engine decomposes each integer across four coprime rails:

• 127
• 8191
• 65536
• 524287

Arithmetic runs independently on the rails, then Garner-style CRT reconstructs the canonical result.

M = 127 × 8191 × 65536 × 524287
  = 35,742,890,181,197,824

Unsigned results live in:

[0, M)

The centered signed interpretation lives in:

[-M/2, M/2 - 1]
=
[-17,871,445,090,598,912, 17,871,445,090,598,911]

Install

pip install rns_engine

Supported Python versions: 3.10–3.13.

AVX2 note

HAS_AVX2 is a build-time property of the compiled extension, not runtime CPU detection. Supported x86-64 wheels can use AVX2 acceleration. macOS wheels intentionally avoid an external OpenMP dependency.

Unsigned quick start

import numpy as np
import rns_engine as rns

x = np.array([123456789, 999999999], dtype=np.uint64)
y = np.array([987654321, 111111111], dtype=np.uint64)

ex = rns.encode(x)
ey = rns.encode(y)

out_add = rns.decode(*rns.add(*ex, *ey))
out_mul = rns.decode(*rns.mul(*ex, *ey))

# Stay in residue space across several exact operations and decode once.
s1 = rns.add(*ex, *ey)
s2 = rns.mul(*s1, *ey)
s3 = rns.sub(*s2, *ex)
out = rns.decode(*s3)

Session API

import numpy as np
import rns_engine as rns

session = rns.Session(cache_capacity=32)
x = np.array([1, 2, 3, 4], dtype=np.uint64)

encoded = session.encode(x)
result = session.mul(session.add(encoded, encoded), encoded)
out = session.decode(result)

one_step = session.one_shot_affine(
    x,
    multiplier=1_000_003,
    addend=7,
)

hot_loop = session.hot_loop_affine(
    x,
    multiplier=1_000_003,
    addend=7,
    iterations=1000,
)

Centered signed arithmetic

import numpy as np
import rns_engine as rns

values = np.array([-5, 0, 7], dtype=np.int64)
rails = rns.encode_signed(values)
round_trip = rns.decode_signed(*rails)

assert np.array_equal(round_trip, values)

encode_signed(...) accepts only integers in [-M/2, M/2 - 1]. It rejects out-of-range inputs instead of silently wrapping them modulo M.

decode_signed(...) returns the canonical centered representative. That representation alone does not prove that an unknown mathematical result did not wrap. Unique signed reconstruction requires independent range evidence:

|result| < M/2

certificate = rns.certify_signed_bound(max_abs_bound=1_000_000)
certificate.require_unique()

print(certificate.headroom)
print(certificate.minimum_required_modulus)

For a signed dot product—or one output entry of a GEMM—the engine can construct the conservative exact bound directly:

certificate = rns.certify_signed_dot_bound(
    inner_dimension=5120,
    left_abs_bound=127,
    right_abs_bound=127,
)

assert certificate.max_abs_bound == 82_580_480
assert certificate.unique
certificate.require_unique()

The law is:

output bound
=
inner_dimension × left bound × right bound + addend bound

SignedSession exposes the existing rail arithmetic with signed external values:

session = rns.SignedSession()
a = session.encode_signed([-5, 7])
b = session.encode_signed([3, -10])
out = session.decode_signed(session.add(a, b))

Weighted INT32 partial accumulation

Version 0.7 introduced the fused direct bridge for digit-plane, chunked-GEMM, and Tensor Core pipelines.

Give the engine signed INT32 partial outputs with shape:

(terms, *output_shape)

and one exact integer positional weight per term:

import numpy as np
import rns_engine as rns

partials = np.array(
    [
        [[1, -2], [3, -4]],
        [[5, 6], [-7, -8]],
        [[-9, 10], [11, -12]],
    ],
    dtype=np.int32,
)

receipt = rns.accumulate_weighted_int32(
    partials,
    weights=[1, -3, 8],
    require_unique=True,
)

exact = receipt.decode_signed()
print(receipt.certificate.max_abs_bound)
print(receipt.certificate.headroom)

For partials P[t] and full integer weights w[t], the receipt uses the exact conservative law:

max_abs_bound
=
sum(abs(w[t]) × max_abs(P[t]))

Important details:

• The original arbitrary-precision Python weights are retained in the proof receipt.
• Only the execution weights are reduced modulo M, which is lawful for rail arithmetic.
• Negative weights and weights wider than uint64 are supported.
• INT32_MIN is measured with its correct magnitude, 2,147,483,648.
• decode_signed() requires a unique range certificate by default.
• decode_modular() always returns the canonical residue in [0, M).
• Empty term sets produce an exact zero result.

The implementation uses one fused native _weighted call to read signed INT32 partials, collect exact per-term magnitude bounds, apply positional weights modulo M, and accumulate all four rails. The pre-fusion staged encode → scale → add body remains available internally as an exact reference for A/B tests. Local AVX2 diagnostics measured roughly 2.1×–2.9× speedup at one million outputs and larger gains on tiny many-term cases; those figures are hardware-specific CPU evidence, not CUDA, Tensor Core, or universal performance claims.

A precomputed receipt can also be created without materializing partial arrays:

certificate = rns.certify_weighted_sum_bound(
    weights=[1, 128, 128**2],
    term_abs_bounds=[100, 200, 300],
)
certificate.require_unique()

Exact GEMM capacity planning

Version 0.8 separates two different safety questions that must not be conflated:

1. Does each grouped raw digit-plane coefficient fit the native signed accumulator?
2. Does the final radix-weighted result fit uniquely inside the global RNS modulus product?

import rns_engine as rns

local = rns.plan_grouped_coefficient_capacity(
    inner_dimension=5120,
    left_digit_abs_bounds=[127] * 8,
    right_digit_abs_bounds=[127] * 8,
    accumulator_bits=32,
)

global_plan = rns.plan_digit_plane_gemm_capacity(
    inner_dimension=5120,
    radix=128,
    left_digit_abs_bounds=[127] * 8,
    right_digit_abs_bounds=[127] * 8,
)

assert local.max_abs_bound == 660_643_840
assert local.safe
assert local.minimum_signed_accumulator_bits == 31

assert global_plan.capacity.required_modulus_bits == 126
assert global_plan.capacity.additional_bits_required == 71
assert not global_plan.current_unique

That workload is locally safe in INT32 but globally modular-only in the current four-rail range. The planner can also validate explicit proposed CRT moduli and rejects candidates that are not pairwise coprime with the existing product.

See CAPACITY_PLANNING.md for the exact laws, examples, and limits of what the receipts prove.

Pre-CUDA exact shared-scale GEMM

Version 0.9 adds the complete CPU-side oracle and integration contract for the future CUDA/T4 body. Before using GPU credits, run:

rns-pre-cuda --report readiness.json --fixture t4_fixture.json

A green receipt proves:

• both seven-rail profiles are pairwise coprime and close the 126-bit target;
• radix-128 signed digit planes round-trip exactly;
• grouped plane GEMMs fit INT32 for the frozen witness;
• weighted seven-rail reconstruction equals direct Python-big-integer GEMM;
• shared rational scales multiply exactly with no denominator GEMM;
• the portable CUDA fixture and full backend verifier agree.

The future CUDA backend must implement both methods in the contract:

class ExactPipelineBackend:
    def grouped_partials(self, left_planes, right_planes):
        ...  # INT8 Tensor Core / cuBLASLt grouped coefficients

    def weighted_rails(self, grouped_partials, weights, moduli):
        ...  # GPU seven-rail modular weighting and accumulation

The verifier refuses a backend that produces correct GEMM partials but silently falls back to CPU RNS. See PRE_CUDA_READINESS.md, CUDA_PARALLEL_LANES.md, and colab/T4_EXACT_GEMM_BRINGUP.ipynb.

import numpy as np
import rns_engine as rns

left = rns.SharedScaleMatrix(
    np.array([[2**40 + 3, -17], [29, 2**35 - 5]], dtype=object),
    scale=6,
)
right = rns.SharedScaleMatrix(
    np.array([[31, -(2**30) + 9], [2**33 + 1, 43]], dtype=object),
    scale=35,
)

receipt = rns.exact_shared_scale_gemm(
    left,
    right,
    config=rns.WideRNSConfig.balanced_seven_rail(),
    left_plane_count=8,
    right_plane_count=8,
)

assert receipt.exact_match
assert receipt.output_scale == 210
result = receipt.as_matrix()

Core API

Rail operations

• encode(x) / decode(r0, r1, r2, r3)
• encode_signed(x) / decode_signed(r0, r1, r2, r3)
• add(...), sub(...), mul(...), div_(...), fma(...)
• op(..., code) where 0=add, 1=mul, 2=sub, 3=div

Signed range receipts

• certify_signed_bound(max_abs_bound)
• certify_signed_dot_bound(inner_dimension, left_abs_bound, right_abs_bound, addend_abs_bound=0)
• SignedRangeCertificate.unique
• SignedRangeCertificate.headroom
• SignedRangeCertificate.minimum_required_modulus
• SignedRangeCertificate.require_unique()

Weighted partial reconstruction

• accumulate_weighted_int32(partials, weights, require_unique=False)
• certify_weighted_sum_bound(weights, term_abs_bounds, addend_abs_bound=0)
• WeightedInt32Result.decode_signed(require_unique=True)
• WeightedInt32Result.decode_modular()

Capacity planning

• plan_signed_capacity(max_abs_bound, additional_moduli=())
• plan_weighted_sum_capacity(weights, term_abs_bounds, addend_abs_bound=0, additional_moduli=())
• plan_digit_plane_gemm_capacity(inner_dimension, radix, left_digit_abs_bounds, right_digit_abs_bounds, ...)
• plan_grouped_coefficient_capacity(inner_dimension, left_digit_abs_bounds, right_digit_abs_bounds, accumulator_bits=32)
• SignedCapacityPlan
• DigitPlaneGemmCapacityPlan
• GroupedCoefficientCapacityPlan

Wide reference and pre-CUDA APIs

• WideRNSConfig.current()
• WideRNSConfig.smallest_product_seven_rail()
• WideRNSConfig.balanced_seven_rail()
• accumulate_weighted_int32_wide(...)
• SharedScaleMatrix
• exact_integer_matmul(...)
• decompose_signed_radix(...) / reconstruct_signed_radix(...)
• grouped_plane_gemm(...) / reconstruct_grouped_partials(...)
• exact_shared_scale_gemm(...)
• CudaGemmFixture
• CpuExactPipelineBackend
• build_cuda_gemm_fixture(...)
• verify_backend(...)
• run_pre_cuda_readiness()
• build_default_pre_cuda_fixture()

Parallel rail-lane planning

• build_parallel_rail_prior(capacity)
• RailLearningLane
• ParallelRailPrior
• search_mersenne_rails_for_capacity(...)

Scalar-broadcast encoded operations

• mul_u64(...)
• fma_u64(...)
• affine_repeat_u64(...)

Fused raw uint64 operations

• add_u64_io(...)
• sub_u64_io(...)
• mul_u64_io(...)
• fma_u64_io(...)
• affine_repeat_u64_io(...)

OpenMP forms use the _omp suffix. Auto-dispatch forms use the _auto suffix.

High-level objects

• Session
• SignedSession
• SessionCache
• EncodedArray
• WeightedInt32Result

Division constraint

Division requires the divisor to be invertible on every rail:

• b % 127 != 0
• b % 8191 != 0
• b is odd for modulus 65536
• b % 524287 != 0

Data model

• Unsigned encode(...) treats inputs as uint64 and reduces values modulo M.
• encode_signed(...) accepts only integers in the centered interval and refuses silent wrapping.
• decode_signed(...) returns the canonical centered int64 representative.
• accumulate_weighted_int32(...) requires an array whose dtype is exactly int32.
• Encoded rails use:
  • r0: uint16
  • r1: uint16
  • r2: uint16
  • r3: uint32
• EncodedArray stores four read-only rails.

Existing verified benchmark

Verified on Google Colab Linux x86-64 with AVX2=True, two OpenMP processors, and two OpenMP threads. Workloads used the installed wheel.

Median throughput over five runs on 1,000,000 uint64 values:

Fused single-step affine

Workload: fma_u64_io(x, 1_000_003, 7)

• fma_u64_io: 47.8 million values/sec
• fma_u64_io_omp with one thread: 80.7 million values/sec
• fma_u64_io_omp with two threads: 84.6 million values/sec

Repeated affine loop

Workload: affine_repeat_u64_io(x, 1_000_003, 7, iterations=1000)

• scalar fused path: 61.19 billion ops/sec
• OpenMP with one thread: 82.80 billion ops/sec
• OpenMP with two threads: 94.86 billion ops/sec

These measurements cover the existing affine kernels, not weighted INT32 reconstruction or the pure-Python capacity planners.

Build from source

git clone https://github.com/playfularchitect/rns_engine.git
cd rns_engine
pip install -e .
python -m pytest tests/ -v

Requirements:

• Python 3.10–3.13
• C++17 compiler
• Access to the build dependencies declared in pyproject.toml

The release workflow compiles and tests the package on Linux, macOS, and Windows, then installs and tests every produced wheel before publication.

Introspection

import rns_engine as rns

rns.info()

rns.M
rns.HALF_M
rns.SIGNED_MIN
rns.SIGNED_MAX
rns.M0
rns.M1
rns.M2
rns.M3
rns.HAS_AVX2

Current release

v0.9.0

• configurable correctness-first CRT oracle for arbitrary pairwise-coprime rails
• two lawful seven-rail profiles closing the 126-bit eight-plane target
• exact signed radix-128 decomposition and reconstruction
• exact grouped digit-plane GEMM coefficients with strict INT32 overflow refusal
• exact shared-scale rational GEMM with integer numerator and scale receipts
• portable JSON CUDA fixture containing inputs, planes, partials, rails, and answer
• full backend contract requiring both grouped partials and GPU rail accumulation
• deterministic pre-CUDA readiness command and frozen witness
• parallel rail-lane shared-prior plan with hard correctness firewall
• Colab/T4 bring-up notebook with one isolated CUDA integration point
• exact Mersenne rail-set search and local/global capacity planners
• fused native four-rail weighted signed INT32 accumulation
• canonical centered signed encoding and no-wrap certificates

The seven-rail body in 0.9 is a CPU correctness oracle, not a native CUDA or production-speed wide RNS kernel.

License

AGPL-3.0-only

If this software is used in a network service, the AGPL requires modified source to be made available to users. Commercial licensing for proprietary or closed-source use is available. Inquiries: ewesley541@gmail.com

Copyright 2026 Evan Wesley