vdr-math 0.1.0

Exact arithmetic with Value, Denominator, Remainder triples

vdr-math

Exact Arithmetic with Value, Denominator, Remainder Triples

pip install vdr-math

from vdr import VDR
x = VDR(1, 3)    # exact 1/3 — three integers, zero truncation
y = VDR(2, 7)    # exact 2/7
z = x + y        # exact 13/21

Remainder is first-class. Never error. Never residue.

vdr-math is an exact arithmetic library where every value is an ordered triple [V, D, R] — Value, Denominator, Remainder. The remainder slot carries the exact structure that conventional systems discard. No float. No truncation. No drift.

Validated across 921 tests in 38 mathematical and computational domains with zero VDR computation errors.

---

Table of Contents

1. The Problem VDR Solves
2. The Triple: V, D, R
3. Core Concepts
   • Closed and Active Objects
   • Normalization
   • Lift and Rebase
   • The divmod Rule: D Never Explodes
   • Functional Remainders
   • Two Equality Relations
4. Installation and Setup
5. Core Modules
   • vdr.core — The Triple
   • vdr.active — Active Multiplication and Division
   • vdr.fn — Functional Remainders and Discrete Calculus
   • vdr.linalg — Exact Linear Algebra
   • vdr.export — Lossy Boundary
   • vdr.basis — D-Frame Management and Q335
6. Math Domains
   • Number Theory
   • Continued Fractions
   • Combinatorics
   • Sequences
   • Polynomial Algebra
   • Symbolic Algebra
   • Probability
   • Geometry
   • Optimization
   • Differential Equations
   • Graph Theory
   • Game Theory
   • Cryptographic Primitives
   • Coding Theory
   • Algebraic Topology
   • Tropical and Lattice Algebra
   • Control Theory
   • Wavelets
   • Chaos and Sensitivity
   • Transcendental Arithmetic
7. Signal Processing
   • Convolution
   • Discrete Fourier Transform
   • Digital Filters
   • Noise Schedules
8. Physics
   • QED Electron g-2
   • Quantum Mechanics
   • Orbital Mechanics
   • Paraxial Optics
   • Structural Mechanics
   • Thermodynamics
   • Crystallography
   • Geodesy
9. Machine Learning
   • Softmax and Exponential
   • Neural Network Layers
   • Automatic Differentiation
   • Optimizers
   • Attention and Transformers
   • Sampling and Decoding
   • Training Pipeline
   • Datasets and Checkpoints
10. Diffusion Models
    • Noise Schedules
    • Forward Process
    • Reverse Process
    • Verification and Drift Testing
11. Q335 Basis — The Configurable D-Frame
12. Design Principles
13. Honest Boundaries
14. Float Comparison Table
15. Validation Summary
16. API Quick Reference
17. Example: Toy LLM

---

The Problem VDR Solves

Floating-point arithmetic silently discards information at every operation. The mantissa is finite, the rounding is invisible, and the error compounds.

Start at 1/7.
Add 1/13 one hundred times.
Subtract 1/13 one hundred times.
Float64 result: 0.14285714285714282 (error ≈ 2.78e-16)
VDR result:     [1, 7, 0] (error = 0)

For a single operation the error is negligible. For a chain of thousands — video generation frame by frame, Kalman filter cycle by cycle, financial risk position by position — the errors accumulate and become indistinguishable from the signal.

VDR eliminates per-step arithmetic error entirely. The remainder slot catches what the denominator frame can't absorb, exactly, with zero loss. Chain length becomes irrelevant to accumulated error.

---

The Triple: V, D, R

Every VDR object is an ordered triple [V, D, R]:

V — Value slot. Arbitrary-precision integer. The settled numerator
    in the current denominator frame.

D — Denominator slot. Nonzero arbitrary-precision integer. The frame
    in which V and R are interpreted.

R — Remainder slot. Exact unresolved structure. Three forms:
    Atomic:     a single integer r
    Composite:  integer base r plus finite list of child VDR objects
    Functional: a Python callable f(depth) → VDR

VDR(3)                    # [3, 1, 0]  — integer 3
VDR(1, 2)                 # [1, 2, 0]  — rational 1/2
VDR(1, 3, Remainder(1))   # [1, 3, 1]  — active: (1+1)/3 = 2/3

V and D are always integers. Recursion exists only in R. Every valid object has finite depth, finite branching, finite total node count. No limits. No approximation. No infinity.

---

Core Concepts

Closed and Active Objects

A closed object has R = 0. It behaves exactly as the rational V/D. This is the bridge to conventional mathematics.

x = VDR(3, 7)      # closed: 3/7
x.is_closed         # True
x.to_fraction()     # Fraction(3, 7)

An active object has R ≠ 0. It carries exact structure that the denominator frame couldn't absorb. The remainder is not error — it is the part of the value that lives outside the current frame.

x = VDR(1, 3, Remainder(1))  # active: (1+1)/3 = 2/3
x.is_active                    # True
x.to_fraction()                # Fraction(2, 3)

Denominator-sensitive completion: R is interpreted within the parent's D frame — divided by D, not added externally. [3, 7, 1] means (3 + 1)/7 = 4/7, not 3/7 + 1.

Normalization

Normalization produces a deterministic canonical form. It is idempotent — normalizing twice gives the same result as normalizing once.

The steps, in order:

1. Sign convention: if D < 0, negate both V and D
2. GCD reduction: divide V and D by their GCD
3. Child normalization: normalize all children bottom-up
4. Atomic consolidation: combine integer contributions at same level
5. Canonical ordering: sort children by |D|, then V, then R structure
6. Same-denominator merge: add closed children sharing D; remove zero-sum pairs
7. Closed-form preference (N7): if the entire remainder tree is value-equivalent to zero, collapse to closed form and GCD-reduce

Rule N7 solves the Newton perfect-square problem: sqrt(4) via Newton iteration produces [2k, k, 0] for large k, which N7 collapses to [2, 1, 0].

VDR(4, 8).normalize()        # VDR(1, 2) — GCD reduced
VDR(1, -3).normalize()       # VDR(-1, 3) — sign convention
VDR(2000000, 1000000).normalize()  # VDR(2, 1) — large GCD

Lift and Rebase

Lift is the remainder transport operator. When a parent denominator frame changes by factor k, the remainder must be rescaled.

lift(r, k) = k * r                              (atomic)
lift([V, D, R], k) = [k*V, D, lift(R, k)]       (child VDR: D preserved, V and R scaled)

Lift composes multiplicatively: lift(lift(R, a), b) = lift(R, a*b).

Rebase changes the top-level denominator while preserving exact value.

x = VDR(1, 2)
y = x.rebase(6)    # VDR(3, 6) — closed rebase, clean

x = VDR(1, 3)
y = x.rebase(7)    # active rebase — mismatch witness in R
y.to_fraction()     # still Fraction(1, 3)

Active rebase uses divmod: V*B = Q*D + S. The result is [Q, B, [S, D, 0] + lift(R, B)]. The mismatch witness [S, D, 0] is the exact part the target denominator couldn't absorb.

The divmod Rule: D Never Explodes

This is the operational rule governing the entire system. When any operation would produce a larger denominator, you divmod and put the overflow in R.

Example: multiplying two Q335 values.

A = [p1, 2^335, 0]
B = [p2, 2^335, 0]

# Naive multiplication would give [p1*p2, 2^670, 0] — D explosion. Never do this.
# Instead:

Q, S = divmod(p1 * p2, 2**335)

# Result: [Q, 2^335, [S, 2^335, 0]]
# D stays at 2^335. V absorbed what fit. R caught what didn't. Exactly. Zero loss.

This is what R exists for. The remainder slot is the pressure valve that keeps D stable while preserving every bit of exactness. Multiplication doesn't change D. Division doesn't change D.

Functional Remainders

A functional remainder is a Python callable stored in the R slot that returns a VDR at any requested depth. This is how VDR handles irrational-like quantities.

from vdr.fn import make_newton_fn, resolve

# sqrt(2) via Newton: x_{n+1} = (x + 2/x) / 2
sqrt2_fn = make_newton_fn("sqrt2", lambda x: (x + VDR(2)/x) / VDR(2))

# wrap in a VDR object
obj = VDR(0, 1, sqrt2_fn)

# resolve at any depth
val = resolve(obj, depth=10)   # >100 correct digits of sqrt(2)

Each depth is a complete exact rational value — not an approximation of a limit. The residual x^2 - 2 is an exact inspectable rational. You know precisely how far from sqrt(2) you are, not approximately.

Quadratic convergence: digits double per step. Depth 1 gives 1 digit. Depth 7 gives >100 digits.

The precision knob: by default resolve fully. But if you want performance, resolve to depth 5 for ~30 digits. You choose where approximation enters, you know exactly how much, and it doesn't compound through subsequent exact rational operations.

Two Equality Relations

Structural equality (≡s): slot-by-slot recursive match. VDR(1, 2) is not structurally equal to VDR(2, 4).

Value equality (≡n): match after normalization. VDR(1, 2) == VDR(2, 4) is True. Python == uses value equality.

a = VDR(1, 2)
b = VDR(2, 4)

a.structural_eq(b)   # False — different V and D
a == b                # True  — same value after normalization
hash(a) == hash(b)    # True  — works as dict keys

---

Installation and Setup

pip install vdr-math

Python 3.8+ required. No external dependencies. Optional: mpmath for arbitrary-precision decimal export.

pip install vdr-math[decimal]   # includes mpmath

from vdr import VDR, Remainder, Vec, Mat

Active multiplication and functional remainder awareness are installed automatically on import. No manual setup required.

---

Core Modules

vdr.core — The Triple

The foundation. Everything else imports from here.

from vdr.core import VDR, Remainder, ZERO, ONE, NEG_ONE

Construction:

VDR(3)                          # integer 3
VDR(1, 2)                       # rational 1/2
VDR(1, 3, 5)                    # active: R = atomic 5
VDR(1, 3, Remainder(1, [VDR(1, 7)]))  # composite R
VDR.from_fraction(Fraction(5, 6))      # from stdlib Fraction
VDR.from_int(42)                       # from int

Arithmetic operators: +, -, *, /, -x, abs(x)

Comparison operators: ==, !=, <, <=, >, >=

Coercion: VDR operators accept int and Fraction on either side.

VDR(1, 2) + 1          # VDR(3, 2)
3 * VDR(1, 7)           # VDR(3, 7)
VDR(1, 2) == Fraction(1, 2)  # True

Projection:

x = VDR(3, 7)
x.to_fraction()     # Fraction(3, 7) — exact for closed, legacy-flattened for active
float(x)            # 0.4285714285714286 — lossy, loss belongs to float

Structural metrics:

x.depth()            # nesting depth of remainder tree
x.size()             # total node count
x.den_complexity()   # (distinct_denoms, sum_denoms, node_count)

Module constants:

ZERO = VDR(0, 1)     # additive identity
ONE = VDR(1, 1)       # multiplicative identity
NEG_ONE = VDR(-1, 1)

vdr.active — Active Multiplication and Division

Extends arithmetic to active objects (R ≠ 0).

from vdr.active import active_mul, active_div

For [V1, D1, R1] * [V2, D2, R2]:

• Frame: D1 * D2
• Closed part: V1 * V2
• Remainder captures three cross-terms: V1·R2, V2·R1, R1·R2

a = VDR(1, 2, Remainder(1))   # value = 1
b = VDR(1, 3, Remainder(1))   # value = 2/3
c = a * b                      # value preserved exactly
c.to_fraction() == a.to_fraction() * b.to_fraction()  # True

Division by active objects: v1 compromise. The divisor is projected to an exact rational via scalar projection, inverted, then multiplied. The divisor's remainder structure is lost. This is acknowledged, not hidden.

Installed automatically by import vdr. Can be uninstalled for testing:

from vdr.active import uninstall, install
uninstall()   # * and / on active objects raise ArithmeticFailure
install()     # restore

vdr.fn — Functional Remainders and Discrete Calculus

from vdr.fn import (
    FnRemainder, resolve, resolve_recursive, is_functional,
    make_newton_fn, make_series_fn, make_iterative_fn, make_constant_fn,
    discrete_derivative, discrete_derivative_nth,
    discrete_integral, discrete_integral_trapz,
)

Newton factory — quadratic convergence:

# sqrt(n) for any n
sqrt2_fn = make_newton_fn("sqrt2", lambda x: (x + VDR(2)/x) / VDR(2))
result = sqrt2_fn.expand(10)  # >100 correct digits

Series factory — partial sums:

# Leibniz series for pi/4
def leibniz_term(n):
    sign = 1 if n % 2 == 0 else -1
    return VDR(sign, 2*n + 1)

pi4_fn = make_series_fn("leibniz_pi4", leibniz_term)
result = pi4_fn.expand(1000)  # ~3 digits at 1000 terms (slow convergence)

Iterative factory — general iteration:

fn = make_iterative_fn("contract", lambda x: x / VDR(2) + VDR(1), VDR(100))
fn.expand(20)  # converges to 2

Discrete calculus — exact, no limits:

# Discrete derivative: Dh(f)(x) = (f(x+h) - f(x)) / h
f = lambda x: x * x
df = discrete_derivative(f, VDR(1, 1000))
df(VDR(3))     # VDR(6001, 1000) — exact

# Discrete integral: left Riemann sum
result = discrete_integral(lambda x: x*x, VDR(0), VDR(1), 10)
# VDR(57, 200) — exact

# Trapezoidal rule — more accurate, still exact
result = discrete_integral_trapz(lambda x: x*x, VDR(0), VDR(1), 100)

# Finite difference tables
d3f = discrete_derivative_nth(lambda x: x*x*x, VDR(1), order=3)
d3f(VDR(0))    # VDR(6) — exactly 3! for cubic
d4f = discrete_derivative_nth(lambda x: x*x*x, VDR(1), order=4)
d4f(VDR(0))    # VDR(0) — exactly zero, no float noise floor

These are explicitly not approximations of continuous calculus. They are a separate exact system where each step size h gives a complete exact answer.

vdr.linalg — Exact Linear Algebra

from vdr.linalg import Vec, Mat

Vec — exact vector:

v = Vec.from_ints([1, 2, 3])
w = Vec.from_fracs([(1, 2), (1, 3), (1, 7)])
v + w, v - w, v * VDR(2), v.dot(w), v.norm_sq()

Mat — exact matrix:

m = Mat.from_ints([[1, 2], [3, 4]])
m.det()          # VDR(-2)
m.inv()          # exact inverse
m.solve(b)       # exact solution to Ax = b
m.rank()         # exact integer rank
m.T              # transpose
m.trace()        # sum of diagonal
m.rref()         # reduced row echelon form
m.matmul(other)  # matrix multiplication
m.matvec(v)      # matrix-vector product

Automatic dispatch: for n ≤ 4, uses cofactor/adjugate/Cramer (simple). For n ≥ 5, uses Gaussian elimination O(n³) (practical). Both available explicitly:

m.det_cofactor()   # O(n!) — always available
m.det_gauss()      # O(n³) — always available
m.inv_adjugate()
m.inv_gauss()
m.solve_cramer(b)
m.solve_gauss(b)

Headline result — Hilbert matrices:

# Hilbert matrix: H[i,j] = 1/(i+j+1). Notoriously ill-conditioned.
H = Mat([[VDR(1, i+j+1) for j in range(5)] for i in range(5)])
H_inv = H.inv()
H.matmul(H_inv) == Mat.identity(5)  # True — exactly zero off-diagonal
# Float gives residual ~1e-9 for H5, meaningless for H10. VDR: exact at any size.

Serialization:

from vdr.linalg import parse_vdr, vdr_to_dict, vdr_from_dict, vdr_to_latex

parse_vdr("[1, 2, 0]")           # VDR(1, 2)
parse_vdr("[1, 3, [1, 6, 0]]")   # nested active

d = vdr_to_dict(VDR(1, 2))       # {"v": 1, "d": 2, "r": {"base": 0, "children": []}}
vdr_from_dict(d)                  # VDR(1, 2) — exact roundtrip

vdr_to_latex(VDR(1, 2))          # "\\frac{1}{2}"

vdr.export — Lossy Boundary

This is where exact VDR precision ends and target-format precision begins. Any loss belongs to the target format, not to VDR.

from vdr.export import to_fraction, to_float, to_decimal

x = VDR(1, 7)
to_fraction(x)          # Fraction(1, 7) — exact
to_float(x)             # 0.14285714285714285 — lossy (64-bit IEEE 754)
to_decimal(x, digits=50) # "0.14285714285714285714285714285714285714285714285714" — via long division or mpmath

vdr.basis — D-Frame Management and Q335

The configurable denominator frame. Like mpmath lets you set mp.dps, vdr-math lets you set a default D.

from vdr.basis import set_default, get_default, to_qbasis, qb_mul, qb_div

set_default(bits=335)    # Q335 — the default, proven across ~1000 tests
set_default(bits=668)    # 200-digit precision
set_default(bits=3322)   # 1000-digit precision

# project values onto the basis
a = to_qbasis(VDR(22, 7), bits=335)   # 22/7 on the 2^335 grid
b = to_qbasis(VDR(1, 3), bits=335)    # 1/3 on the grid

# multiply — D stays fixed, overflow in R
c = qb_mul(a, b, bits=335)
c.d == 2**335   # True — D never changed

Q335 is the default because it's proven, not because it's special. Any D is valid. D = 7 for working in sevenths. D = 2^16 for binary fixed-point. D = 1 for plain integers.

---

Math Domains

Every domain is a built-in library module, ready to use on import. These are not examples — they are active library components.

Number Theory

from vdr.math.number_theory import (
    harmonic, euler_totient, farey, egyptian_fractions,
    stern_brocot, vdr_gcd, vdr_lcm, vdr_pow, convergents, e_cf,
)

harmonic(10)              # VDR(7381, 2520) — exact
harmonic(100)             # exact, ~85-digit denominator
euler_totient(100)        # 40
farey(5)                  # 11 fractions, mediant property |ad-bc|=1 for all pairs

ef = egyptian_fractions(3, 7)
sum(f.to_fraction() for f in ef) == Fraction(3, 7)  # True

convergents([3, 7, 15, 1])  # [3/1, 22/7, 333/106, 355/113]
e_cf(10)                     # first 10 CF coefficients of e

Continued Fractions

from vdr.math.continued_fractions import to_cf, from_cf, sqrt_cf_period, best_rational

to_cf(355, 113)           # [3, 7, 16]
from_cf([3, 7, 16])       # VDR(355, 113) — roundtrip exact

a0, period = sqrt_cf_period(2)   # (1, [2])
a0, period = sqrt_cf_period(7)   # (2, [1, 1, 1, 4])

best_rational(VDR(314159, 100000), 1000)  # best approximation with denom ≤ 1000

Combinatorics

from vdr.math.combinatorics import binom, bell, derangement, multinomial, catalan, stirling2

binom(20, 10)              # VDR(184756)
bell(5)                    # VDR(52)
derangement(7)             # VDR(1854)
multinomial(10, [3, 3, 4]) # VDR(4200)
catalan(5)                 # VDR(42)
stirling2(5, 3)            # VDR(25)

# Pascal's rule verified for all n, k:
binom(10, 3) == binom(9, 2) + binom(9, 3)   # True

Sequences

from vdr.math.sequences import fibonacci, lucas, bernoulli, tribonacci, rational_recurrence

fibonacci(30)       # VDR(832040) — fast doubling algorithm
lucas(10)           # VDR(123)
bernoulli(12)       # VDR(-691, 2730) — exact, cached

# Cassini identity: F(n-1)*F(n+1) - F(n)^2 = (-1)^n
# Verified exactly for all n — no float noise floor

# General rational recurrence
seq = rational_recurrence(
    [VDR(3, 2), VDR(-1, 2)],    # coefficients
    [VDR(1), VDR(2)],            # initial values
    14                            # terms
)

Polynomial Algebra

from vdr.math.polynomial import (
    poly_eval, poly_add, poly_mul, poly_divmod, poly_gcd,
    lagrange_interpolate, definite_integral, rational_roots,
)

# Polynomials as coefficient lists [a0, a1, a2, ...] (constant first)
p = [VDR(1), VDR(1), VDR(1)]   # 1 + x + x^2
poly_eval(p, VDR(2))            # VDR(7)

# Lagrange interpolation
points = [(VDR(0), VDR(1)), (VDR(1), VDR(3)), (VDR(2), VDR(7))]
p = lagrange_interpolate(points)  # recovers [1, 1, 1] exactly

# Polynomial GCD — exact zero-testing
gcd = poly_gcd([VDR(-1), VDR(0), VDR(1)], [VDR(1), VDR(2), VDR(1)])
# (x+1) — exact. Decimal cannot do this at any precision.

# Exact definite integral via antiderivative
definite_integral([VDR(0), VDR(0), VDR(1)], VDR(0), VDR(1))  # VDR(1, 3)

Symbolic Algebra

from vdr.math.symbolic import partial_fractions_simple, ratfun_add, power_sum

# Partial fraction decomposition
pf = partial_fractions_simple([VDR(1)], [VDR(1), VDR(2)])
# 1/((x-1)(x-2)) = -1/(x-1) + 1/(x-2)

# Power sums
power_sum(3, 100)   # 1^3 + 2^3 + ... + 100^3 = 25502500
# Nicomachus' theorem: S1(n)^2 == S3(n) — verified exactly

Probability

from vdr.math.probability import (
    binom_pmf_full, bayes_update, bayes_sequential,
    markov_steady_state, expected_value, variance,
)

# Binomial PMF sums to exactly 1
pmf = binom_pmf_full(10, VDR(1, 3))
total = VDR(0)
for p in pmf:
    total = total + p
assert total == VDR(1)   # exactly 1, not 0.9999999999999998

# Bayesian updating — exact posteriors
post = bayes_update(VDR(1, 2), VDR(3))   # VDR(3, 4)

# Markov steady state — sums to exactly 1
P = Mat.from_fracs([[(1, 2), (1, 2)], [(1, 4), (3, 4)]])
ss = markov_steady_state(P)   # [VDR(1, 3), VDR(2, 3)]

Geometry

from vdr.math.geometry import (
    line_intersect, polygon_area, barycentric, point_in_triangle,
    circumcenter, dist_sq,
)

# Line intersection — exact
pt = line_intersect((VDR(0), VDR(0)), (VDR(2), VDR(2)),
                    (VDR(0), VDR(2)), (VDR(2), VDR(0)))
# (VDR(1), VDR(1)) — no epsilon needed

# Point-in-triangle — exact boolean
point_in_triangle((VDR(1, 3), VDR(0)), 
                  (VDR(0), VDR(0)), (VDR(1), VDR(0)), (VDR(0), VDR(1)))
# True — on edge, no "almost on edge" ambiguity

# Barycentric coordinates sum to exactly 1
l1, l2, l3 = barycentric(p, a, b, c)
assert l1 + l2 + l3 == VDR(1)

# Circumcenter — equidistant from all three vertices exactly
cx, cy = circumcenter(a, b, c)
assert dist_sq((cx, cy), a) == dist_sq((cx, cy), b) == dist_sq((cx, cy), c)

Optimization

from vdr.math.optimization import newton_optimize, bisection, simplex_2d

# Newton for optimization — converges in 1 step for quadratic
x = newton_optimize(
    lambda x: VDR(2)*x - VDR(4),   # f'
    lambda x: VDR(2),               # f''
    VDR(0), 1
)   # VDR(2)

# Bisection — exact rational midpoints
x = bisection(lambda x: x*x - VDR(2), VDR(1), VDR(2), 30)
# |x^2 - 2| < 10^-8

# Simplex — exact rational vertex enumeration
result = simplex_2d([VDR(1), VDR(1)], [[VDR(1), VDR(1)]], [VDR(4)])

Differential Equations

from vdr.math.differential_eq import euler_solve, rk4_solve, lotka_volterra_solve

# Euler method — exact at each step
traj = euler_solve(lambda x, y: y, VDR(1), VDR(0), VDR(1, 10), 10)
# y(1) = (11/10)^10 exact

# RK4 — ~140x more accurate than Euler, still exact per step
traj = rk4_solve(lambda x, y: y, VDR(1), VDR(0), VDR(1, 10), 10)

# Lotka-Volterra — 200 steps, all exact
traj = lotka_volterra_solve(
    (VDR(10), VDR(5)), VDR(1, 100),
    VDR(1, 10), VDR(1, 100), VDR(1, 10), VDR(1, 100), 200
)

VDR eliminates arithmetic error but not method error. Euler at h=1/10 gives (11/10)^10 regardless of exact or approximate arithmetic. The advantage is zero drift and exact reproducibility.

Graph Theory

from vdr.math.graph import dijkstra, bellman_ford, floyd_warshall, pagerank, prim_mst

# Dijkstra with exact rational weights
dist = dijkstra({0: [(1, VDR(1, 3))], 1: [(2, VDR(1, 2))], 2: []}, 0)

# PageRank — sums to exactly 1 via Cramer's rule
pr = pagerank(transition_matrix)

# Floyd-Warshall — all-pairs shortest paths, exact
result = floyd_warshall(3, distance_matrix)

Game Theory

from vdr.math.game_theory import minimax_2x2, shapley_values, cournot_duopoly

# Minimax — exact mixed strategies
p, q, val = minimax_2x2(Mat.from_ints([[3, -1], [-2, 4]]))
# p* = 3/5, q* = 1/2, value = 1

# Shapley values — sum to exactly v(N)
phi = shapley_values(v_func, 3)

# Cournot duopoly — exact equilibrium
q1, q2, pi1, pi2 = cournot_duopoly(VDR(100), VDR(1), VDR(10), VDR(20))

Cryptographic Primitives

from vdr.math.cryptographic import rsa_keygen, rsa_encrypt, rsa_decrypt, baby_giant

# RSA — exact modular integer arithmetic (float categorically excluded)
n, e, d = rsa_keygen(61, 53, 17)   # (3233, 17, 2753)
c = rsa_encrypt(42, e, n)
m = rsa_decrypt(c, d, n)            # 42 exactly

# Discrete logarithm
baby_giant(2, 8, 13)   # 3 (because 2^3 = 8 mod 13)

Coding Theory

from vdr.math.coding_theory import (
    gf_inv, hamming74_encode, hamming74_correct, min_distance,
)

# Galois field arithmetic
gf_inv(3, 7)   # 5 (because 3*5 = 15 = 1 mod 7)

# Hamming(7,4) — all errors corrected exactly
codeword = hamming74_encode([1, 0, 1, 1])
received = list(codeword)
received[2] ^= 1                       # inject error
corrected = hamming74_correct(received) # corrected == codeword

Algebraic Topology

from vdr.math.topology import betti_numbers, euler_characteristic, simplicial_complex_boundaries

simplices = {
    0: [(0,), (1,), (2,)],
    1: [(0,1), (0,2), (1,2)],
    2: [(0,1,2)],
}
boundaries = simplicial_complex_boundaries(simplices)
betti = betti_numbers(boundaries)   # [1, 0] — filled triangle
euler_characteristic(betti)          # 1

# d∘d = 0 verified as exact zero matrix

Tropical and Lattice Algebra

from vdr.math.tropical import (
    trop_add, trop_mul, trop_mat_mul,
    gram_schmidt_exact, lll_reduce, lovasz_condition,
)

# Tropical: min-plus algebra
trop_add(VDR(3), VDR(5))   # VDR(3) — min
trop_mul(VDR(3), VDR(5))   # VDR(8) — ordinary add

# Exact Gram-Schmidt — cross-dot products are exactly 0
ortho, mu = gram_schmidt_exact([v1, v2, v3])
assert ortho[0].dot(ortho[1]) == VDR(0)

# LLL — exact Lovász condition, no float rounding in threshold
reduced = lll_reduce(basis)
# mu = 0.500000000000001 vs 0.499999999999999: float makes wrong decision. VDR: exact.

Control Theory

from vdr.math.control import (
    state_evolve, is_controllable, cayley_hamilton_verify,
    transfer_function_eval, mat_pow,
)

# State evolution — zero drift after 100+ steps
traj = state_evolve(A, B, x0, inputs)

# Cayley-Hamilton — p(A) = exact zero matrix
result = cayley_hamilton_verify(A)
assert result == Mat.zero(n, n)   # not ≈ 10^-15, exactly zero

# Transfer function at complex frequency
H = transfer_function_eval([VDR(1)], [VDR(2), VDR(3), VDR(1)], (VDR(0), VDR(1)))
# (VDR(1, 10), VDR(-3, 10)) — exact complex value

Wavelets

from vdr.math.wavelets import haar_forward, haar_inverse, haar_multilevel, parseval_verify

# Perfect reconstruction — exact
signal = [VDR(1), VDR(3), VDR(5), VDR(7)]
avgs, dets = haar_forward(signal)
recovered = haar_inverse(avgs, dets)
assert recovered == signal   # exact

# Multi-level: 64 samples, 6 levels, perfect reconstruction
signal64 = [VDR(i) for i in range(64)]
decomp = haar_multilevel(signal64, 6)
recovered = haar_reconstruct_multilevel(decomp)
assert recovered == signal64

# Parseval energy identity — exact
assert parseval_verify(signal, decomp)

Chaos and Sensitivity

from vdr.math.chaos import tent_map, iterate_map, detect_period, logistic_map

# Tent map on 1/7: period 3, exact forever
orbit = iterate_map(tent_map, VDR(1, 7), 20)
period = detect_period(orbit)   # 3
# Float diverges at ~25 steps. VDR: period 3 forever.

# Logistic map: exact but exponential denominator growth at r=4
x = logistic_map(VDR(1, 3), VDR(4))   # VDR(8, 9) exact
# NOTE: In flat Fraction, denominator digits grow as ~2^n (step 10: ~258,
# step 20: ~260,000). VDR with Q335 frame avoids this — D stays 2^335,
# tree depth grows by 1 per step. Step 20 in Q335: ~4,080 tree digits,
# not ~260,000. Cost is tree depth, not denominator explosion.

Periodic rational orbits under chaotic maps are computationally free — denominators stay bounded.

Transcendental Arithmetic

from vdr.math.transcendental import (
    PI, E, LN2, SQRT2, PHI, ZETA3, ZETA5, CATALAN,  # Q335 constants
    PI_FN, SQRT2_FN, ZETA3_FN,                        # functional remainders
    sqrt_newton, exp_series, sin_series, cos_series, ln_series,
    borwein_zeta, elliptic_k, hypergeometric_2f1, pi_machin,
)

22 named constants at 100-digit precision:

PI       # VDR(2198864..., 2^335) — ready to use
E        # VDR(1902580..., 2^335)
SQRT2    # VDR(9898366..., 2^335)
ZETA3    # VDR(8413439..., 2^335)
# ... and 18 more

# Addition is one integer add over shared denominator
pi_plus_e = PI + E   # D stays 2^335, one integer add

Functional remainder versions for arbitrary precision:

from vdr.fn import resolve
val = resolve(VDR(0, 1, SQRT2_FN), depth=20)  # ~200 digits

Series functions:

exp_series(VDR(1), depth=20)          # e to ~18 digits
sin_series(VDR(1, 4), depth=16)       # sin(1/4) exact rational
pi_machin(terms=160)                   # pi to ~67 digits

Borwein acceleration:

borwein_zeta(5, n=210)   # zeta(5) to 100+ digits
borwein_zeta(7, n=210)   # zeta(7) to 100+ digits
# Geometric convergence 3^(-n). Universal for any s.

Elliptic integrals:

elliptic_k(VDR(1, 2), terms=500)   # K(k) at k^2 = 1/2
elliptic_e(VDR(1, 4), terms=500)   # E(k) at k^2 = 1/4

---

Signal Processing

Convolution

from vdr.signal.convolution import convolve, correlate, deconvolve

y = convolve([VDR(1), VDR(2), VDR(3)], [VDR(1), VDR(1)])
# [VDR(1), VDR(3), VDR(5), VDR(3)] — exact

# Deconvolution: exact inverse
x = deconvolve(convolve(a, b), b)
assert x == a   # exact recovery

Discrete Fourier Transform

from vdr.signal.dft import exact_dft, exact_idft, parseval_verify

x = [VDR(1), VDR(2), VDR(3), VDR(4)]
X = exact_dft(x)          # list of (real, imag) VDR pairs
recovered = exact_idft(X)
assert recovered == x      # IDFT(DFT(x)) == x exactly

parseval_verify(x, X)      # True — energy preserved exactly

Digital Filters

from vdr.signal.filters import iir_filter, moving_average

# IIR: y[n] = a*y[n-1] + x[n] — exact at every step
y = iir_filter([VDR(1)] + [VDR(0)]*19, VDR(1, 2))
# y[n] = (1/2)^n — exact forever
# Year-long operation: same precision as second 1

Noise Schedules

from vdr.signal.schedule import linear_schedule, compute_alpha_bars

betas = linear_schedule(5, VDR(1, 100), VDR(1, 20))
# All exact rationals

---

Physics

QED Electron g-2

from vdr.physics.qed import a2_coefficient, anomalous_moment

a2 = a2_coefficient()
# 197/144 + pi^2/12 + 3*zeta(3)/4 - (pi^2/2)*ln(2)
# All Q335 basis constants. Matches -0.328478965579... to 100 digits.

ae = anomalous_moment(n_loops=3)
# Perturbation series with exact coefficients

Quantum Mechanics

from vdr.physics.quantum import (
    SIGMA_X, SIGMA_Z, spin_rotation, measurement_probabilities,
    verify_unitarity, pauli_multiply,
)

# Pauli algebra — structural identity, not approximate
pauli_multiply("x", "y")   # ((0, 1), "z") — sigma_x * sigma_y = i * sigma_z
pauli_multiply("x", "x")   # ((1, 0), "I") — sigma_x^2 = I

# Spin rotation — 4x pi/2 about z returns to I exactly
U = spin_rotation("z", VDR(1, 2))
verify_unitarity(U)   # True, exact

# Measurement probabilities sum to exactly 1
probs = measurement_probabilities(Vec([VDR(3, 5), VDR(4, 5)]))
# [VDR(9, 25), VDR(16, 25)] — sum = VDR(1) exactly

Orbital Mechanics

from vdr.physics.orbital import kepler_newton, orbit_closure_verify

E = kepler_newton(VDR(1), VDR(1, 2), depth=20)
# Eccentric anomaly, >100 correct digits

err = orbit_closure_verify(positions)
# Float: ~1e-12 position error. VDR: exactly 0.

Paraxial Optics

from vdr.physics.optics import free_space, thin_lens, system_matrix, verify_symplecticity, matrix_power

M = system_matrix([free_space(VDR(1)), thin_lens(VDR(2)), free_space(VDR(1))])
verify_symplecticity(M)   # True — det(M) == 1 exactly

M1000 = matrix_power(M, 1000)
verify_symplecticity(M1000)   # Still True. Float: ~1e-12 after 1000.

Structural Mechanics

from vdr.physics.structural import solve_structure, verify_equilibrium

u = solve_structure(K, F)
verify_equilibrium(K, u, F)   # True — K@u == F exactly, not "residual < tolerance"

Thermodynamics

from vdr.physics.thermo import partition_function, ising_1d_transfer

Z = partition_function([VDR(0), VDR(1), VDR(2)], VDR(1))
Z_ising = ising_1d_transfer(VDR(1), VDR(0), VDR(1), 10)
# Transfer matrix power, exact. No Monte Carlo, no sign problems.

Crystallography

from vdr.physics.crystallography import point_group_matrix, verify_group_closure

matrices = [point_group_matrix(op) for op in ["E", "C4z", "C2z", "C4z_inv"]]
verify_group_closure(matrices)   # True — exact structural equality, not tolerance

Geodesy

from vdr.physics.geodesy import helmert_forward, helmert_roundtrip_verify

helmert_roundtrip_verify(coords, params)   # True — exact recovery
# Float: ~1 nm error. VDR: zero.

---

Machine Learning

Softmax and Exponential

from vdr.ml.softmax import softmax
from vdr.ml.exp import exp_series

probs = softmax(Vec.from_ints([1, 2, 3]))
# Sums to exactly 1. Not 0.9999999999999998. Exactly 1.

e_val = exp_series(VDR(1), depth=20)   # e to ~18 digits

Neural Network Layers

from vdr.ml.nn import Linear, ReLU, Sequential

model = Sequential([
    Linear.from_ints([[1, 2], [3, 4]], [0, 0]),
    ReLU(),
    Linear.from_ints([[1, 0], [0, 1]], [1, 1]),
])
output = model.forward(Vec.from_ints([1, 1]))

# Backward — exact gradients via chain rule
grad = mse_grad(output, target)
model.backward(grad)

Automatic Differentiation

from vdr.ml.autodiff import Node, relu, mse_loss

a = Node(VDR(3))
b = Node(VDR(4))
c = a * b + a   # builds computation graph
c.backward()
print(a.grad)    # VDR(5) = d/da(a*b + a) = b + 1. Exact.

Optimizers

from vdr.ml.optim import SGD, Momentum

opt = SGD(model.parameters(), lr=VDR(1, 100))
opt.step()   # w = w - lr * grad. Exact. No float accumulation.

opt = Momentum(model.parameters(), lr=VDR(1, 100), beta=VDR(9, 10))
opt.step()   # v = beta*v + grad; w = w - lr*v. All exact.

Attention and Transformers

from vdr.ml.attention import self_attention
from vdr.ml.transformer import TransformerLM

Q = [Vec.from_ints([1, 0]), Vec.from_ints([0, 1])]
K = [Vec.from_ints([1, 1]), Vec.from_ints([1, -1])]
V = [Vec.from_ints([10, 0]), Vec.from_ints([0, 10])]
out = self_attention(Q, K, V, causal=True)
# Every attention weight row sums to exactly 1

Sampling and Decoding

from vdr.ml.sampling import categorical_sample, top_k_probs, nucleus_probs
from vdr.ml.rng import VDRRandom

rng = VDRRandom(seed=42)   # deterministic, platform-independent
idx = categorical_sample(probs, rng)   # exact rational CDF comparison

filtered = top_k_probs(probs, k=5)    # renormalized to sum exactly 1
filtered = nucleus_probs(probs, VDR(9, 10))  # top-p, exact threshold

Training Pipeline

from vdr.ml.trainer import train_step, train_epoch, evaluate_classification

loss = train_step(model, x, y, optimizer)   # exact loss as VDR
losses = train_epoch(model, dataset, optimizer)
accuracy = evaluate_classification(model, test_set)   # exact rational

Datasets and Checkpoints

from vdr.ml.datasets import tiny_text_dataset, one_hot
from vdr.ml.checkpoint import save_model, load_model_parameters

windows, vocab, inv_vocab = tiny_text_dataset("the cat sat on the mat", seq_len=2)
vec = one_hot(2, 10)   # exact: VDR(1) at index 2, VDR(0) elsewhere

state = save_model(model)   # JSON-serializable, exact VDR triples
load_model_parameters(model, state)   # bit-identical on any platform

---

Diffusion Models

Diffusion Noise Schedules

from vdr.diffusion.schedule import DiffusionSchedule, linear_schedule, cosine_schedule

schedule = linear_schedule(T=5, beta_start=VDR(1, 100), beta_end=VDR(1, 20))
# All betas, alphas, alpha_bars exact rationals
# alpha_bar_T = 26821179/31250000 — verified bit-identical against Python Fraction

schedule.sqrt_alpha_bars[t]           # sqrt via Newton, >100 digits
schedule.sqrt_one_minus_alpha_bars[t] # same
schedule.posterior_variance(t)         # exact closed positive rational

Forward Process

from vdr.diffusion.forward import forward_sample, forward_trajectory

# x_t = sqrt(alpha_bar_t) * x0 + sqrt(1 - alpha_bar_t) * epsilon
xt = forward_sample(x0, t, schedule, epsilon)

trajectory = forward_trajectory(x0, schedule, epsilons)
assert trajectory[0] == x0   # identity, exact

Reverse Process

from vdr.diffusion.reverse import compute_x0_prediction, reverse_step_ddim

# With oracle noise predictor: prediction error = exactly 0
x0_pred = compute_x0_prediction(xt, t, schedule, eps_pred)

# DDIM deterministic step — roundtrip error = exactly 0
x_prev = reverse_step_ddim(xt, t, t_prev, schedule, eps_pred)

Verification and Drift Testing

from vdr.diffusion.sampling import (
    verify_forward_reverse_roundtrip, verify_ddim_roundtrip,
    verify_multi_step_drift, make_oracle_predictor,
)

# DDIM roundtrip: error = exactly 0 (strongest result)
err = verify_ddim_roundtrip(x0, schedule, epsilon)

# Multi-cycle drift: error at cycle N = error at cycle 1 (central result)
drift = verify_multi_step_drift(x0, schedule, epsilon, num_cycles=10)
# drift[0] ≈ drift[9] — does not grow

# Oracle predictor separates arithmetic error from model error
oracle = make_oracle_predictor(x0, schedule)

Drift comparison:

Cycles	Float64 error	VDR error
1	~1e-15	< 1e-50
1,000	~1e-12	< 1e-50
36,000 (30s video)	~1e-11	< 1e-50
8,640,000 (2hr film)	~1.9e-8	< 1e-50

---

Q335 Basis — The Configurable D-Frame

Q335 means D = 2^335. It's the default because it's proven across ~1000 tests, but it's a configuration choice, not a definition.

Why 2^335? The continued fraction convergents of e include three with power-of-two denominators: 2/2^0, 11/2^2, and 87/2^5. Extending to 2^335 gives 100-digit precision for all 22 transcendental constants in the basis. n=335 is minimal: at n=334, Catalan's G fails at digit 101.

Precision floor: 2^(-336) ≈ 10^(-101.2). The Planck length is ~10^(-35). The first divergent digit is 66 orders of magnitude below the Planck length.

Scaling: ~3.32 bits per additional decimal digit. 200 digits → 2^668. 1000 digits → 2^3322.

from vdr.basis import set_default

set_default(bits=335)    # 100 digits (default)
set_default(bits=668)    # 200 digits
set_default(bits=3322)   # 1000 digits

---

Design Principles

1. Remainder is first-class. R carries exact unresolved structure that scalar systems discard. It is part of the value, not error.

2. D never explodes. Every operation that would increase D instead uses divmod to keep D fixed and places the exact overflow in R.

3. Preserve data, compress shape. No operation may discard remainder. Exactness is maintained through all operations.

4. Recursion only through R. V and D are integers. Nesting occurs exclusively in the remainder slot.

5. Exact as written. No valid object depends on approximation, limits, or infinite expansion.

6. The code is the specification. Every claim verified by executable tests. Implementation is normative.

7. No required external dependencies. Python 3.8+ stdlib only. mpmath optional for decimal export.

8. Honest boundaries. Limitations are documented, not hidden. Chaotic dynamics have exponential cost. Active division loses divisor structure. These are acknowledged facts, not defects.

---

Honest Boundaries

Chaotic dynamics: exact representation of chaotic orbits has real computational cost. In flat Fraction representation, logistic map at r=4 has denominator digits growing as ~2^n after n steps. VDR with Q335 fixed-frame avoids denominator explosion — D stays at 2^335 and the remainder tree grows linearly (one ~102-digit level per step). At step 30: flat Fraction needs ~10^9 digits; Q335 tree needs ~6,120 digits (~163,000× compression). The cost in VDR is tree depth, not denominator explosion. Periodic rational orbits are computationally free — denominators stay bounded.

No native irrationals or complex numbers. Functional remainders produce exact rationals approaching irrationals at any depth. Complex extension is engineering work, not mathematical obstacle.

Active division loses divisor structure. When dividing by an active object, the divisor is projected to exact rational via scalar projection, then inverted. The divisor's remainder structure is lost. This is the v1 compromise — acknowledged, not hidden.

Cofactor expansion is O(n!). Replaced by Gaussian elimination O(n³) for n ≥ 5. Practical limit ~50×50.

Computational cost. ~50-200× float per operation in Python. ~150× on Q335 GPU. Practical path: VDR for validation providing exact ground truth to verify float implementations on larger systems.

VDR does not replace real numbers. It is an exact finite alternative for domains where rational arithmetic and recursive rational refinement suffice.

---

Float Comparison Table

Computation	Float64 error	VDR error
200-op return to origin	~2.78e-16	0
Hilbert 5×5 inverse residual	~1e-9	0
Hilbert 10×10 inverse residual	~1e-2 (catastrophic)	0
IIR (1/√2)^20	~1e-16	0 (exact 1/1024)
Spin rotation ×4 = I	~1e-15 per entry	0
Binomial PMF sum n=10	~2e-16	0 (exact 1)
Orbit closure error	~1e-12	0
Cayley-Hamilton residual	~1e-15 per entry	0
det(M) after 1000 optics elements	~1e-12	0 (exact 1)
Helmert roundtrip	~1 nm	0
Tent map 1/7, 25 steps	diverged (total loss)	0 (period 3 forever)
36,000-step diffusion chain	~1e-11	< 1e-50
8.64M-step chain (2hr film)	~1.9e-8	< 1e-50

---

Validation Summary

Paper	Domains	Tests	Passed	Failed (test error)	Failed (VDR error)
VDR-1	Core	68	68	0	0
VDR-2	15 gyms	285	279	6	0
VDR-3	8 new gyms	157	152	5	0
VDR-4	LLM pipeline	198	196	2	0
VDR-13	14 physics	—	—	—	0
VDR-26	Diffusion	37	33	4	0
VDR-27	35 domains	—	—	—	0
Total	38 domains	921	903	18	0

All 18 failures traced to test-design errors (wrong expected values, too-tight thresholds, normalization presentation, precision frame mismatches). Zero VDR computation errors.

The system remains falsifiable: any test producing an incorrect exact rational from correct inputs would falsify VDR. 921 tests have not produced one.

---

API Quick Reference

Core: VDR(v, d, r), Remainder(base, children), ZERO, ONE, NEG_ONE

Arithmetic: +, -, *, /, -x, abs(x), ==, !=, <, <=, >, >=

Methods: .normalize(), .to_fraction(), .to_float(), .rebase(d), .depth(), .size(), .den_complexity(), .is_closed, .is_active, .bracket()

Constructors: VDR.from_fraction(f), VDR.from_int(n), Vec.from_ints([]), Vec.from_fracs([]), Mat.from_ints([[]]), Mat.from_fracs([[]]), Mat.identity(n), Mat.zero(r, c)

Linear algebra: Vec.dot(), Vec.norm_sq(), Mat.det(), Mat.inv(), Mat.solve(b), Mat.rank(), Mat.rref(), Mat.matmul(), Mat.matvec(), Mat.T, Mat.trace()

Functional remainder: resolve(obj, depth), make_newton_fn(name, step), make_series_fn(name, term), make_iterative_fn(name, step, start)

Discrete calculus: discrete_derivative(f, h), discrete_integral(f, a, b, n), discrete_integral_trapz(f, a, b, n)

Basis: set_default(bits), to_qbasis(x, bits), qb_mul(a, b, bits), qb_div(a, b, bits)

Export: to_fraction(x), to_float(x), to_decimal(x, digits)

Serialization: parse_vdr(text), vdr_to_dict(x), vdr_from_dict(d), vdr_to_latex(x)

---

Example: Toy LLM

Toy LLM Example — A complete single-block transformer language model running entirely in exact VDR arithmetic at D = 2^32. Trains on "the cat sat on the mat," generates text, and passes 9 verification tests including bit-identical determinism and exact softmax-sums-to-1. No floating-point anywhere — forward pass, backpropagation, attention, quadratic softmax surrogate, SGD weight updates, and categorical sampling all use basis-safe VDR operators with automatic frame management. 181 parameters, 20 epochs, under 10 seconds on a 2019 laptop. See the Toy LLM README for architecture details, results, and precision analysis.