Theoretical Foundations

RMath is built on rigorous mathematical algorithms designed for numerical stability and high-speed execution. We prioritize algorithms that eliminate floating-point error without sacrificing parallel performance.

Kahan Compensated Summation

When summing millions of floating-point numbers, the standard approach ($\sum x_i$) suffers from catastrophic cancellation because smaller numbers are swallowed by the accumulated sum's lower precision bounds.

RMath uses Kahan Summation, which maintains a running compensation variable ($c$) to accumulate these lost low-order bits and feed them back into the main sum.

Parallel Kahan

To achieve O(1) error bounds in parallel, RMath performs independent Kahan sums within chunked threads. The results are then merged using a secondary Kahan summation pass, ensuring the final result is exactly as precise as a serial execution, but orders of magnitude faster.

Welford's Algorithm Proof

Standard variance calculation $(\sum x^2 - (\sum x)^2/n)$ is prone to catastrophic cancellation. Welford's algorithm calculates a running mean and M2 sum to ensure precision even for extreme values.

Implementation Detail: RMath parallelizes Welford's by using a merge-reduction. Each thread calculates a local mean/M2, and these are then combined using the following merge rules:

$\Delta = \mu_B - \mu_A$
$\mu_{combined} = \mu_A + \Delta \frac{n_B}{n_A + n_B}$
$M2_{combined} = M2_A + M2_B + \Delta^2 \frac{n_A n_B}{n_A + n_B}$

Proven: Automatic Differentiation

Forward-mode AD in RMath is mathematically exact. We use Dual Numbers $(a + b\epsilon)$ to track derivatives through any arbitrary composition of transcendental functions.

dual_verification.py

from rmath.calculus import Dual

# f(x) = sin(x) + exp(x)
# f'(x) = cos(x) + exp(x)
x = Dual(0.0, 1.0)
y = x.sin() + x.exp()

print(f"Value f(0): {y.value}")
print(f"Derivative f'(0): {y.derivative}")

Value f(0): 1.0 Derivative f'(0): 2.0

Computational Graph Proof

The Tensor autograd engine uses a dynamic tape with reverse topological sorting. We prove its correctness by tracking node identity and verifying against analytical gradients.

autograd_proof.py

import rmath.array as ra

# f(x) = (x * 2 + 5)^2
x = ra.Tensor([3.0], requires_grad=True)
y = (x * 2 + 5)**2
y.backward()

# Analytical f'(x) = 2(x*2+5) * 2 = 4(3*2+5) = 44
print(f"Gradient: {x.grad}")

Gradient: Array([44.0000])