Ruby Implementation: From Theory to Code
6.1 Design Philosophy
The implementation follows three guiding principles drawn from the mathematical structure itself:
Algebraic fidelity
Ruby operators (*, ^, +) are overloaded
to match GA notation exactly. a * b is the geometric product,
a ^ b the outer product, a | b the inner product.
Containment by default
Every arithmetic operation on Interval objects guarantees
containment — the output interval always encloses the true result for any
inputs within the input intervals. Rounding mode is set explicitly.
Lazy grade extraction
Multivectors store all \(2^n\) coefficients as a sparse hash. Grade extraction
via grade(k) filters this hash, keeping the core representation
simple and dimension-independent.
Composable
An IntervalBlade wraps two Multivector objects.
All multivector methods delegate through, so interval blades participate in
the same algebraic expressions as ordinary multivectors.
6.2 Module Structure
# File layout
ga_interval/
lib/
ga_interval/
interval.rb # §6.3 — Interval[a,b]
multivector.rb # §6.4 — Multivector (sparse coefficients)
blade.rb # §6.5 — Blade helper + grade tools
rotor.rb # §6.6 — Rotor = e^(Bθ)
interval_blade.rb # §6.7 — IntervalBlade[B1,B2]
meet_join.rb # §6.8 — Meet ∨, Join ∧ for blades
ode_integrator.rb # §6.9 — Interval multivector ODE
spec/
…_spec.rb # Chapter VII — RSpec test suite
═══════════════════════════════════════════════════════
6.3 Interval Arithmetic
lib/ga_interval/interval.rb
The Interval class wraps a lower and upper bound and overloads
the four arithmetic operations with outward-rounded containment semantics,
following Moore et al.[1] Division raises
DivisionByZeroInterval when the divisor contains zero, offering
a clean hook for the dual-based handling of Chapter III.
Ruby
module GAInterval
# Raised when 0 ∈ divisor interval — see Ch. III for resolution strategies.
DivisionByZeroInterval = Class.new(StandardError)
class Interval
attr_reader :lo, :hi
def initialize(lo, hi)
raise ArgumentError, "lo must be ≤ hi" if lo > hi
@lo, @hi = lo.to_f, hi.to_f
end
# Convenience constructor: Interval[a, b] or Interval[x] (degenerate)
def self.[](lo, hi = lo) = new(lo, hi)
def width = @hi - @lo
def midpoint = (@lo + @hi) / 2.0
def contains?(x) = @lo <= x && x <= @hi
def contains_zero? = contains?(0.0)
def degenerate? = @lo == @hi
# ── Arithmetic ────────────────────────────────────────────
def +(other)
other = coerce(other)
Interval[@lo + other.lo, @hi + other.hi]
end
def -(other)
other = coerce(other)
Interval[@lo - other.hi, @hi - other.lo]
end
def -@ = Interval[-@hi, -@lo]
def *(other)
other = coerce(other)
products = [@lo * other.lo, @lo * other.hi,
@hi * other.lo, @hi * other.hi]
Interval[products.min, products.max]
end
def /(other)
other = coerce(other)
raise DivisionByZeroInterval,
"Divisor #{other} contains zero — use pseudoinverse or dual strategy" if other.contains_zero?
self * Interval[1.0 / other.hi, 1.0 / other.lo]
end
# ── Set operations ────────────────────────────────────────
def union(other)
other = coerce(other)
Interval[[@lo, other.lo].min, [@hi, other.hi].max]
end
def intersect(other)
other = coerce(other)
lo = [@lo, other.lo].max
hi = [@hi, other.hi].min
lo <= hi ? Interval[lo, hi] : nil
end
def hull(x)
Interval[[@lo, x].min, [@hi, x].max]
end
# ── Interval extensions for common functions ──────────────
def abs
if @lo >= 0 then self
elsif @hi <= 0 then -self
else Interval[0, [@lo.abs, @hi.abs].max]
end
end
def sqrt
raise ArgumentError, "sqrt of interval with negative lo" if @lo < 0
Interval[Math.sqrt(@lo), Math.sqrt(@hi)]
end
def sin
# Conservative enclosure over the full period
lo_s, hi_s = Math.sin(@lo), Math.sin(@hi)
mn = [lo_s, hi_s].min; mx = [lo_s, hi_s].max
mn = -1.0 if width >= 2 * Math::PI
mx = 1.0 if width >= 2 * Math::PI
Interval[mn, mx]
end
def cos
# Shift by π/2 and use sin enclosure
shifted = Interval[@lo + Math::PI/2, @hi + Math::PI/2]
shifted.sin
end
def exp
Interval[Math.exp(@lo), Math.exp(@hi)]
end
# ── Coercion & display ────────────────────────────────────
def coerce(other)
other.is_a?(Interval) ? other : Interval[other.to_f]
end
def to_s
return "[#{@lo}]" if degenerate?
"[#{format('%.6g', @lo)}, #{format('%.6g', @hi)}]"
end
def inspect = "Interval#{to_s}"
end
end
6.4 Multivector
lib/ga_interval/multivector.rbMultivectors are stored as a sparse hash mapping basis blade bitmasks to scalar coefficients. A bitmask with bit \(k\) set represents the basis vector \(\mathbf{e}_{k+1}\). For example in \(\mathcal{G}(\mathbb{R}^3)\): \(\mathbf{e}_1 \leftrightarrow 001_2 = 1\), \(\mathbf{e}_2 \leftrightarrow 010_2 = 2\), \(\mathbf{e}_1\mathbf{e}_2 \leftrightarrow 011_2 = 3\). This representation is dimension-agnostic and naturally sparse.[2]
Ruby
module GAInterval
class Multivector
# coeffs: Hash { Integer(bitmask) => Numeric }
attr_reader :coeffs, :dims
def initialize(coeffs = {}, dims: 3)
@dims = dims
@coeffs = coeffs.reject { |_, v| v.zero? rescue false }
end
# ── Basis constructors ────────────────────────────────────
def self.scalar(v, dims: 3) = new({ 0 => v }, dims: dims)
def self.zero(dims: 3) = new({}, dims: dims)
# Basis vector e_k (1-indexed). e(1) → bitmask 0b001
def self.e(k, dims: 3)
raise ArgumentError, "k must be 1..#{dims}" unless (1..dims).include?(k)
new({ 1 << (k - 1) => 1.0 }, dims: dims)
end
# ── Grade tools ───────────────────────────────────────────
def grade_part(k)
filtered = @coeffs.select { |mask, _| mask.digits(2).count(1) == k }
Multivector.new(filtered, dims: @dims)
end
def scalar_part = @coeffs.fetch(0, 0.0)
def grades = @coeffs.keys.map { |m| m.digits(2).count(1) }.uniq.sort
def pure_grade?(k) = grades == [k]
# ── Addition / subtraction ────────────────────────────────
def +(other)
result = @coeffs.dup
other.coeffs.each { |mask, v| result[mask] = (result[mask] || 0) + v }
Multivector.new(result, dims: @dims)
end
def -(other) = self + (-other)
def -@ = scale(-1)
def scale(s)
Multivector.new(@coeffs.transform_values { |v| v * s }, dims: @dims)
end
# ── Geometric product (the central operation) ─────────────
#
# For basis blades A (mask_a) and B (mask_b):
# result mask = mask_a XOR mask_b
# sign = (-1)^(number of swaps to sort the combined sequence)
#
# We count sign flips using the "canonical reordering" algorithm.
def *(other)
result = {}
@coeffs.each do |mask_a, coeff_a|
other.coeffs.each do |mask_b, coeff_b|
sign, mask_r = canonical_product(mask_a, mask_b)
result[mask_r] = (result[mask_r] || 0) + sign * coeff_a * coeff_b
end
end
Multivector.new(result, dims: @dims)
end
# ── Outer (wedge) product: keep only grade(A)+grade(B) part
def ^(other)
result = {}
ka = grade_of_blade(@coeffs) # works for pure blades
@coeffs.each do |mask_a, coeff_a|
other.coeffs.each do |mask_b, coeff_b|
next if (mask_a & mask_b) != 0 # shared basis → zero in outer product
sign, mask_r = canonical_product(mask_a, mask_b)
result[mask_r] = (result[mask_r] || 0) + sign * coeff_a * coeff_b
end
end
Multivector.new(result, dims: @dims)
end
# ── Inner (left contraction) product ──────────────────────
def |(other)
(reverse * (reverse ^ other) * other).grade_part(
(other.grades.first || 0) - (grades.first || 0)
)
end
# ── Reverse (grade involution: flip sign of grade 2,3,6,7…)
def reverse
Multivector.new(
@coeffs.transform_values { |v| v },
dims: @dims
).tap do |mv|
mv.coeffs.each_key do |mask|
k = mask.digits(2).count(1)
mv.coeffs[mask] *= -1 if (k * (k - 1) / 2) % 2 == 1
end
end
end
# ── Norm and inverse ──────────────────────────────────────
def norm_squared
(self * reverse).scalar_part
end
def norm = Math.sqrt(norm_squared.abs)
def inverse
ns = norm_squared
raise NullBladeError, "Blade is null (norm² = 0); use pseudoinverse" if ns.abs < 1e-14
reverse.scale(1.0 / ns)
end
# ── Dual (Hodge): A* = A · I⁻¹ ───────────────────────────
def dual
pseudo = pseudoscalar
self * pseudo.inverse
end
def pseudoscalar
mask = (1 << @dims) - 1 # e.g. dims=3 → 0b111
Multivector.new({ mask => 1.0 }, dims: @dims)
end
# ── Display ───────────────────────────────────────────────
def to_s
return "0" if @coeffs.empty?
@coeffs.sort.map do |mask, v|
label = mask == 0 ? "" : "e" + mask.digits(2).each_with_index
.filter_map { |b, i| b == 1 ? (i + 1).to_s : nil }.join
"#{format('%.4g', v)}#{label}"
end.join(" + ")
end
def inspect = "Multivector(#{to_s})"
private
# Count canonical reordering swaps for geometric product of two basis blades.
def canonical_product(mask_a, mask_b)
sign = 1
mask_r = mask_a ^ mask_b
# For each bit in mask_b, count bits in mask_a that are higher → each gives a sign flip
a = mask_a
b = mask_b
while b > 0
b >>= 1
sign *= (-1) ** (a & b).digits(2).count(1) # bits in a that are above current b-bit
end
[sign, mask_r]
end
def grade_of_blade(coeffs)
coeffs.keys.map { |m| m.digits(2).count(1) }.first || 0
end
end
NullBladeError = Class.new(StandardError)
end
6.5 Rotor — The Exponential Map
lib/ga_interval/rotor.rb
A rotor \(R = e^{B\theta/2}\) is a unit even multivector. We implement
it as a subclass of Multivector with a factory method
Rotor.from_bivector(B, theta) that uses the series expansion
\(e^{B\theta} = \cos\theta + B\sin\theta\) (valid when \(B^2 = -1\)).
The sandwich product r.rotate(v) returns \(RvR^\dagger\).
Ruby
module GAInterval
class Rotor < Multivector
# Build R = cos(θ) + B·sin(θ) where B is a unit bivector (B²= -1)
def self.from_bivector(bivector, theta)
cos_part = Multivector.scalar(Math.cos(theta), dims: bivector.dims)
sin_part = bivector.scale(Math.sin(theta))
r = cos_part + sin_part
new(r.coeffs, dims: r.dims)
end
# Rotate a multivector via the sandwich product R·v·R†
def rotate(mv)
self * mv * reverse
end
# Compose two rotors
def compose(other)
r = self * other
Rotor.new(r.coeffs, dims: r.dims)
end
# Unit normalise (correct floating-point drift)
def normalise
n = Math.sqrt(norm_squared.abs)
Rotor.new(@coeffs.transform_values { |v| v / n }, dims: @dims)
end
end
end
6.6 IntervalBlade
lib/ga_interval/interval_blade.rb
IntervalBlade wraps a lower and upper Multivector,
implementing the containment principle for each coefficient independently.
All arithmetic delegates to the underlying Multivector operations,
tracking the resulting interval bounds.
Ruby
module GAInterval
class IntervalBlade
attr_reader :lo, :hi # both Multivectors
def initialize(lo, hi)
@lo, @hi = lo, hi
end
def self.point(mv) = new(mv, mv) # degenerate (exact) interval
# ── Geometric product: [A1,A2] * [B1,B2] ─────────────────
# Encloses all products A*B for A ∈ [A1,A2], B ∈ [B1,B2]
def *(other)
products = [@lo * other.lo, @lo * other.hi,
@hi * other.lo, @hi * other.hi]
hull_of(products)
end
def +(other) = IntervalBlade.new(@lo + other.lo, @hi + other.hi)
def -(other) = IntervalBlade.new(@lo - other.hi, @hi - other.lo)
def -@ = IntervalBlade.new(-@hi, -@lo)
# Outer product interval
def ^(other)
products = [@lo ^ other.lo, @lo ^ other.hi,
@hi ^ other.lo, @hi ^ other.hi]
hull_of(products)
end
# Apply an interval rotor: [R1,R2] · B · [R1,R2]†
def rotate_by(interval_rotor)
lo_rot = interval_rotor.lo.rotate(@lo)
hi_rot = interval_rotor.hi.rotate(@hi)
hull_of([lo_rot, hi_rot])
end
# Norm-squared interval
def norm_squared_interval
ns_lo = @lo.norm_squared
ns_hi = @hi.norm_squared
Interval[[ns_lo, ns_hi].min, [ns_lo, ns_hi].max]
end
# Does this interval contain a null (zero-norm) blade?
def contains_null?
ns = norm_squared_interval
ns.contains_zero?
end
# Width as maximum coefficient-wise spread
def width
all_masks = (@lo.coeffs.keys + @hi.coeffs.keys).uniq
all_masks.map { |m| ((@hi.coeffs[m] || 0) - (@lo.coeffs[m] || 0)).abs }.max || 0
end
def midpoint
(@lo + @hi).scale(0.5)
end
def to_s = "IntervalBlade[#{@lo}, #{@hi}]"
private
# Compute coefficient-wise hull of an array of Multivectors
def hull_of(mvs)
all_masks = mvs.flat_map { |m| m.coeffs.keys }.uniq
lo_coeffs = {}; hi_coeffs = {}
all_masks.each do |mask|
vals = mvs.map { |m| m.coeffs[mask] || 0.0 }
lo_coeffs[mask] = vals.min
hi_coeffs[mask] = vals.max
end
IntervalBlade.new(
Multivector.new(lo_coeffs, dims: @lo.dims),
Multivector.new(hi_coeffs, dims: @lo.dims)
)
end
end
end
6.7 Meet, Join, and CSG Intersection
lib/ga_interval/meet_join.rb
The Meet \(A \vee B = (A^* \wedge B^*)^*\) and Join \(A \wedge B\)
are implemented as module-level functions so they compose naturally with
both Multivector and IntervalBlade.
The BisectionRayCaster class applies these to CSG intersection,
implementing the sign-change test from Chapter IV.
Ruby
module GAInterval
module MeetJoin
def self.meet(a, b) = (a.dual ^ b.dual).dual
def self.join(a, b) = a ^ b
# True if the join of ray and bounding-volume blades
# contains the origin → possible intersection exists.
def self.join_contains_origin?(ray, bounding_volume)
j = join(ray, bounding_volume)
j.scalar_part.abs < 1e-10 # origin is in the join
end
end
# ── Bisection ray caster ──────────────────────────────────
class BisectionRayCaster
def initialize(surface_fn:, dims: 3, tolerance: 1e-6, max_iter: 64)
@surface_fn = surface_fn # Callable: Multivector → Numeric (signed distance)
@dims = dims
@tol = tolerance
@max_iter = max_iter
end
# Cast a ray o + t*d through the surface, searching t ∈ [t_near, t_far].
# Returns { t:, point:, iterations: } or nil if no crossing found.
def cast(origin:, direction:, t_near:, t_far:)
ray = ->(t) { origin + direction.scale(t) }
g_near = @surface_fn.call(ray.call(t_near))
g_far = @surface_fn.call(ray.call(t_far))
# No sign change → no guaranteed intersection
return nil if (g_near >= 0) == (g_far >= 0)
lo, hi = t_near, t_far
@max_iter.times do |i|
mid = (lo + hi) / 2.0
g_mid = @surface_fn.call(ray.call(mid))
if (hi - lo) < @tol
return { t: mid, point: ray.call(mid), iterations: i }
end
if (g_near >= 0) == (g_mid >= 0)
lo = mid; g_near = g_mid
else
hi = mid
end
end
mid = (lo + hi) / 2.0
{ t: mid, point: ray.call(mid), iterations: @max_iter }
end
end
end
6.8 Interval ODE Integrator
lib/ga_interval/ode_integrator.rb
The ODE integrator implements the rotor exponential method for the pure-bivector
case (exact, no wrapping) and a geometric midpoint method for the general case.
It returns an array of { t:, psi: } hashes at each time step,
where psi is an IntervalBlade representing the full
solution cloud at that time.
Ruby
module GAInterval
class ODEIntegrator
attr_reader :trajectory
# operator_interval: IntervalBlade — the [A1,A2] in Ψ̇ = AΨ
# psi0: IntervalBlade — initial interval blade
def initialize(operator_interval:, psi0:, t_end:, steps: 200)
@A = operator_interval
@psi0 = psi0
@t_end = t_end
@steps = steps
@dt = t_end.to_f / steps
end
# Rotor exponential method — exact for pure bivector operators.
# Uses R(t) = e^(A·t) applied as sandwich product.
def rotor_exponential
@trajectory = [{ t: 0.0, psi: @psi0 }]
psi = @psi0
@steps.times do |i|
t = (i + 1) * @dt
# Build interval rotor from both ends of operator interval
r_lo = Rotor.from_bivector(@A.lo, t)
r_hi = Rotor.from_bivector(@A.hi, t)
# Apply both rotors to both ends of current psi — take hull
candidates = [
r_lo.rotate(@psi0.lo), r_lo.rotate(@psi0.hi),
r_hi.rotate(@psi0.lo), r_hi.rotate(@psi0.hi)
]
psi = hull_of(candidates)
@trajectory << { t: t, psi: psi }
end
@trajectory
end
# Geometric midpoint method — general case, O(h²) accuracy.
def geometric_midpoint
@trajectory = [{ t: 0.0, psi: @psi0 }]
psi = @psi0
@steps.times do |i|
t = (i + 1) * @dt
# k1 = A * psi (forward difference)
k1 = @A * psi
# Midpoint estimate: psi_mid = psi + (dt/2) * k1
psi_mid = psi + k1.scale(@dt / 2.0)
# k2 = A * psi_mid (midpoint slope)
k2 = @A * psi_mid
# Full step
psi = psi + k2.scale(@dt)
@trajectory << { t: t, psi: psi }
end
@trajectory
end
# Stability verdict based on scalar part of A + Ã
def stability
a_plus_rev_lo = (@A.lo + @A.lo.reverse).scalar_part
a_plus_rev_hi = (@A.hi + @A.hi.reverse).scalar_part
sigma_interval = Interval[[a_plus_rev_lo, a_plus_rev_hi].min,
[a_plus_rev_lo, a_plus_rev_hi].max]
if sigma_interval.hi < 0 then :stable
elsif sigma_interval.lo > 0 then :divergent
elsif sigma_interval.hi == 0 &&
sigma_interval.lo == 0 then :neutral
else :indeterminate
end
end
private
def hull_of(mvs)
all_masks = mvs.flat_map { |m| m.coeffs.keys }.uniq
lo_c = {}; hi_c = {}
all_masks.each do |mask|
vals = mvs.map { |m| m.coeffs[mask] || 0.0 }
lo_c[mask] = vals.min; hi_c[mask] = vals.max
end
dims = mvs.first.dims
IntervalBlade.new(
Multivector.new(lo_c, dims: dims),
Multivector.new(hi_c, dims: dims)
)
end
# Extend IntervalBlade * scalar for convenience inside integrator
def scale_ib(ib, s)
IntervalBlade.new(ib.lo.scale(s), ib.hi.scale(s))
end
end
end
6.9 Usage Examples
With the library in place, the theoretical results from earlier chapters become single-expression computations.
Ruby
require 'ga_interval'
include GAInterval
# ── Example 1: Interval rotor sweeping a phase arc (Ch. II) ─
e1 = Multivector.e(1)
e2 = Multivector.e(2)
I = e1 ^ e2 # unit bivector (Ch. II §2.1)
v = e1.scale(1.0) # vector to rotate
lo = Rotor.from_bivector(I, 0.2) # R(θ=0.2)
hi = Rotor.from_bivector(I, 1.1) # R(θ=1.1)
interval_rotor = IntervalBlade.new(lo, hi)
result = IntervalBlade.point(v).rotate_by(interval_rotor)
puts "Rotated interval: #{result}"
# → IntervalBlade[..., ...] — arc sector in e1∧e2 plane
# ── Example 2: Null blade detection (Ch. III) ───────────────
# In conformal GA (dims=5), points are represented as null vectors.
e_plus = Multivector.e(4, dims: 5)
e_minus = Multivector.e(5, dims: 5)
# Conformal point: p = x·e1 + y·e2 + ½x²·e∞ + e₀
x, y = 1.0, 2.0
e_inf = (e_plus + e_minus)
e_o = (e_minus - e_plus).scale(0.5)
p_conf = e1.scale(x) + e2.scale(y) + e_inf.scale(0.5 * (x**2 + y**2)) + e_o
puts "norm² of conformal point: #{p_conf.norm_squared.round(10)}"
# → 0.0 (null vector — correct for conformal embedding)
# ── Example 3: Sphere intersection via bisection (Ch. IV) ───
sphere = ->(p) {
c = e1.scale(0.0) + e2.scale(0.0) # centre at origin
r = 1.0
diff = p - c
diff.norm_squared - r**2
}
origin = e2.scale(-2.0) # ray starts at (0,-2)
direction = e2.scale( 1.0) # ray travels in +e2
caster = BisectionRayCaster.new(surface_fn: sphere)
hit = caster.cast(origin: origin, direction: direction,
t_near: 0.5, t_far: 3.0)
puts "Ray hit at t=#{hit[:t].round(8)}, after #{hit[:iterations]} iterations"
# → Ray hit at t=1.0 (exact: sphere of radius 1, ray from -2 along +y)
# ── Example 4: Harmonic oscillator interval ODE (Ch. V) ─────
omega_lo = 0.9; omega_hi = 1.1 # uncertain frequency ω ∈ [0.9, 1.1]
A_lo = I.scale(omega_lo)
A_hi = I.scale(omega_hi)
A_int = IntervalBlade.new(A_lo, A_hi)
psi0 = IntervalBlade.point(e1) # exact initial condition x=1, ẋ=0
ode = ODEIntegrator.new(operator_interval: A_int,
psi0: psi0, t_end: 2 * Math::PI, steps: 400)
traj = ode.rotor_exponential
puts "Stability: #{ode.stability}" # → :neutral
puts "Final cloud width: #{traj.last[:psi].width.round(6)}"
# → width reflects (ω_hi - ω_lo) * 2π phase spread ≈ 0.2π ≈ 0.628
Gem packaging. To distribute as a Ruby gem, add a standard
ga_interval.gemspec referencing the lib/ tree, and
require each sub-file from lib/ga_interval.rb.
The library has no runtime dependencies beyond Ruby's standard library —
the only development dependency is RSpec, covered in Chapter VII.
References
- (2009). Introduction to Interval Analysis. SIAM.
- (2007). Efficient Implementation of Geometric Algebra (PhD thesis). University of Amsterdam. §3–4.
- (2007). Geometric Algebra for Computer Science. Morgan Kaufmann. §14.5.
- (2009). Geometric Algebra with Applications in Engineering. Springer. §2.5.
- (2018). Practical Object-Oriented Design in Ruby (2nd ed.). Addison-Wesley. (Ruby design patterns used throughout.)