diff --git a/kingdon/numerical.py b/kingdon/numerical.py
new file mode 100644
index 0000000..a699ea4
--- /dev/null
+++ b/kingdon/numerical.py
@@ -0,0 +1,39 @@
+def exp(x, n=10):
+    """
+    Compute the exponential function using Taylor series expansion.
+
+    Uses incremental term computation for improved numerical stability and performance.
+
+    :param x: The input value (scalar or multivector).
+    :param n: Number of terms in the Taylor series (default: 10).
+    :return: Approximation of exp(x).
+    """
+    result = 1.0
+    term = 1.0
+    for k in range(1, n):
+        term *= x / k
+        result += term
+    return result
+
+
+def is_close(a, b, tol=1e-6):
+    """
+    Check if two multivectors are close within tolerance.
+    
+    Handles comparison of multivectors with potentially different sparse representations.
+    
+    :param a: First multivector or scalar value.
+    :param b: Second multivector or scalar value.
+    :param tol: Tolerance for comparison (default: 1e-6).
+    :return: True if a and b are close within tolerance, False otherwise.
+    """
+    if hasattr(a, 'keys') and hasattr(b, 'keys'):
+        all_keys = set(a.keys()) | set(b.keys())
+        for key in all_keys:
+            val_a = a.e if key == 0 else getattr(a, a.algebra.bin2canon[key], 0) if key in a.keys() else 0
+            val_b = b.e if key == 0 else getattr(b, b.algebra.bin2canon[key], 0) if key in b.keys() else 0
+            if abs(val_a - val_b) >= tol:
+                return False
+        return True
+    else:
+        return abs(a - b) < tol
diff --git a/tests/test_kingdon.py b/tests/test_kingdon.py
index 781bffc..1befc6d 100644
--- a/tests/test_kingdon.py
+++ b/tests/test_kingdon.py
@@ -12,6 +12,7 @@
 from sympy import Symbol, expand, sympify, cos, sin, sinc
 from kingdon import Algebra, MultiVector, symbols
 from kingdon.operator_dict import UnaryOperatorDict
+from kingdon.numerical import exp as numerical_exp, is_close
 
 import timeit
 
@@ -881,6 +882,19 @@ def test_simple_exp():
     l = (-(B | B).e) ** 0.5
     assert R == cos(l) + sinc(l) * B
 
+def test_numerical_exp_simple_bivector(sta):
+    """Test numerical.exp() on a simple bivector matches MultiVector.exp()"""
+    B = sta.bivector([1.0, 0, 0, 0, 0, 0])
+    mv_exp = B.exp()
+    num_exp = numerical_exp(B, n=20)
+    assert is_close(mv_exp, num_exp, tol=1e-10)
+
+def test_numerical_exp_non_simple_bivector(sta):
+    """Test numerical.exp() on a non-simple bivector doesn't raise an error"""
+    B = sta.bivector([1.0, 2.0, 0, 0, 0, 0])
+    result = numerical_exp(B, n=20)
+    assert result is not None
+
 def test_dual_numbers():
     alg = Algebra(r=1)
     x = alg.multivector(name='x')