diff --git a/numga/multivector/extension/decompose.py b/numga/multivector/extension/decompose.py
index 0ea6ae3..3f963dc 100644
--- a/numga/multivector/extension/decompose.py
+++ b/numga/multivector/extension/decompose.py
@@ -71,10 +71,30 @@ def motor_rotor(m) -> "Translator":
 	""""""
 	return m.nondegenerate()
 
-mv.motor_split = SubspaceDispatch("""Split motor wrt a given origin""")
-@mv.motor_split.register(lambda m, o: m.inside.motor() and o.inside.vector())
-def motor_split(m: Motor, o: Vector) -> Tuple[Motor, Motor]:
-	o = o.dual()
+mv.motor_split = SubspaceDispatch("""
+	Split motor wrt a given origin in two components;
+	one that leaves the origin invariant and one
+	complimentary to that
+	""")
+@mv.motor_split.register(lambda m, o: m.inside.motor() and o.inside.antivector() and o.dual().is_degenerate)
+def motor_split_euclidian(m: Motor, o: AntiVector) -> Tuple[Motor, Motor]:
+	"""If the dual of o is degenerate, it represents the euclidean origin
+	In this special case there is a more optimized solution
+	m = t * r
+	"""
+	return motor_translator(m), motor_rotor(m)
+@mv.motor_split.register(lambda m, o: m.inside.motor() and o.inside.antivector() and o.is_blade)
+def motor_split_blade(m: Motor, o: AntiVector) -> Tuple[Motor, Motor]:
+	"""If the dual of o is a blade, we can also specialize
+	"""
+	# drop translational axes
+	R = m.subspace.reject(o.subspace.dual())
+	T = m.context.subspace.vector().product(o.subspace.dual())
+	t, r = motor_split(m, o)
+	return t.select_subspace(T), r.select_subspace(R)
+@mv.motor_split.register(lambda m, o: m.inside.motor() and o.inside.antivector())
+def motor_split(m: Motor, o: AntiVector) -> Tuple[Motor, Motor]:
 	# construct shortest trajectory from o to m >> o
 	t = ((m >> o) / o).motor_square_root()
-	return t, ~t * m
\ No newline at end of file
+	# the residual is a rotation that leaves o invariant
+	return t, ~t * m
diff --git a/numga/multivector/extension/norms.py b/numga/multivector/extension/norms.py
index be0b5f1..61dc558 100644
--- a/numga/multivector/extension/norms.py
+++ b/numga/multivector/extension/norms.py
@@ -68,4 +68,5 @@ def reverse_study_normalized(x):
 	"""Binds to 4d and 5d even grade, as well as many other objects"""
 	s = x.symmetric_reverse_product()
 	s = s.square_root().inverse()
+	# NOTE: order of s and x here is important!
 	return s * x
diff --git a/numga/multivector/extension/optimized.py b/numga/multivector/extension/optimized.py
index 1ebb611..f2e719f 100644
--- a/numga/multivector/extension/optimized.py
+++ b/numga/multivector/extension/optimized.py
@@ -17,6 +17,8 @@
 [2] May the Forque Be with You; Dynamics in PGA
 
 """
+import numpy as np
+
 # FIXME: move this to numpy backend module?
 from numga.multivector.multivector import AbstractMultiVector as mv
 @mv.normalized.register(
@@ -54,3 +56,75 @@ def normalize_3dpga_motor_numpy(m: "Motor") -> "Motor":
 
 # reset cache
 mv.normalized.cache = {}
+@mv.exp.register(
+	lambda s:
+	s.named_str == 'xw,yw,zw' and
+	s.algebra.description.description_str == 'x+y+z+w0',
+	position=0      # try and match this before all else
+)
+def exponentiate_degenerate_bivector_numpy(b: "BiVector") -> "Motor":
+	"""Optimized 3d pga exponentiation, for the given memory layout
+	"""
+	M = np.empty(b.shape + (4,))
+	M[..., 0] = 1
+	M[..., 1:] = b.values
+	return b.context.multivector(values=M, subspace=b.context.subspace.translator())
+
+@mv.exp.register(
+	lambda s:
+	s.named_str == 'xy,xz,yz' and
+	s.algebra.description.description_str == 'x+y+z+w0',
+	position=0      # try and match this before all else
+)
+def exponentiate_nondegenerate_bivector_numpy(b: "BiVector") -> "Motor":
+	"""Optimized 3d pga exponentiation, for the given memory layout
+	"""
+	xy, xz, yz = np.moveaxis(b.values, -1, 0)
+	l = xy*xy + xz*xz + yz*yz
+	a = np.sqrt(l)
+	c = np.cos(a)
+	s = np.where(a > 1e-20, np.sin(a) / a, 1)
+
+	M = np.empty(b.shape + (4,))
+	M[..., 0] = c
+	M[..., 1] = s*xy
+	M[..., 2] = s*xz
+	M[..., 3] = s*yz
+	return b.context.multivector(values=M, subspace=b.context.subspace.rotor())
+
+@mv.exp.register(
+	lambda s:
+	s.named_str == 'xy,xz,yz,xw,yw,zw' and
+	s.algebra.description.description_str == 'x+y+z+w0',
+	position=0      # try and match this before all else
+)
+def exponentiate_bivector_numpy(b: "BiVector") -> "Motor":
+	"""Optimized 3d pga exponentiation, for the given memory layout
+
+	Notes
+	-----
+	Terms involving xz are negated compared to the reference [1],
+	since our normalized blade order is to sort the axes alphabetically,
+	and our subspace object does not currently store signs of basis blades
+	"""
+	xy, xz, yz, xw, yw, zw = np.moveaxis(b.values, -1, 0)
+	l = xy*xy + xz*xz + yz*yz
+	a = np.sqrt(l)
+	m = xy*zw - xz*yw + yz*xw
+	c = np.cos(a)
+	s = np.where(a > 1e-20, np.sin(a) / a, 1)
+	t = m * (c - s) / l
+
+	M = np.empty(b.shape + (8,))
+	M[..., 0] = c
+	M[..., 1] = s*xy
+	M[..., 2] = s*xz
+	M[..., 3] = s*yz
+	M[..., 4] = s*xw + t*yz
+	M[..., 5] = s*yw - t*xz
+	M[..., 6] = s*zw + t*xy
+	M[..., 7] = m*s
+	return b.context.multivector(values=M, subspace=b.context.subspace.motor())
+
+# reset cache
+mv.exp.cache = {}
diff --git a/numga/multivector/extension/test/test_decompose.py b/numga/multivector/extension/test/test_decompose.py
index e29d0d1..ddc1b9a 100644
--- a/numga/multivector/extension/test/test_decompose.py
+++ b/numga/multivector/extension/test/test_decompose.py
@@ -14,6 +14,7 @@
 	(5, 0, 0), (4, 0, 1), (4, 1, 0), (3, 1, 1), #(2, 2, 1),
 ])
 def test_invariant_decomposition(descr):
+	np.random.seed(1)
 	print()
 	print(descr)
 	algebra = Algebra.from_pqr(*descr)
@@ -21,50 +22,66 @@ def test_invariant_decomposition(descr):
 	for i in range(10):
 		b = random_motor(ga).select[2]
 
-		print('decompose')
 		l, r = b.decompose_invariant()
-		# if ga.algebra.pseudo_scalar_squared == 0:
-		# 	assert r.subspace.is_degenerate
-		q = l.inner(r)
-		npt.assert_allclose(q.values, 0, atol=1e-9)
-		q = l.squared() # test they square to a scalar
-		npt.assert_allclose(q.values[1:], 0, atol=1e-9)
-		q = r.squared()
-		npt.assert_allclose(q.values[1:], 0, atol=1e-9)
+		assert_close(l.inner(r), 0)
+		assert_close(l.squared().select.nonscalar(), 0)
+		assert_close(r.squared().select.nonscalar(), 0)
 
-		npt.assert_allclose(b.values, (l+r).values, atol=1e-9)
+		assert_close(b, l + r)
 
 		L = l.exp()
 		R = r.exp()
 		# test that component bivectors commute under exponentiation
-		npt.assert_allclose((L * R).values, (R * L).values)
+		assert_close(L * R, R * L)
 
 
 @pytest.mark.parametrize('descr', [
 	(1, 0, 1), (2, 0, 1), (3, 0, 1), (4, 0, 1),
 ])
 def test_motor_decompose_euclidean(descr):
-
 	ga = NumpyContext(descr)
 	m = random_motor(ga, (10,))
 	t = m.motor_translator()
-	print(t)
+	# print(t)
 	r = m.motor_rotor()
 
-	assert np.allclose((t * r - m).values, 0, atol=1e-9)
+	assert_close(t * r, m)
 	# alternate method of constructing translator
 	tt = (m * ~r).select.translator()
-	assert np.allclose((tt - t).values, 0, atol=1e-9)
-
-
-def test_motor_decompose_elliptical():
-	ga = NumpyContext((4,0,0))
-	m = random_motor(ga, (10,))
-
-	t, r = m.motor_split(ga.multivector.a)
+	assert_close(tt, t)
 
-	assert np.allclose((t * r - m).values, 0, atol=1e-9)
 
-	assert np.allclose((t.symmetric_reverse_product() - 1).values, 0, atol=1e-9)
-	assert np.allclose((r.symmetric_reverse_product() - 1).values, 0, atol=1e-9)
+@pytest.mark.parametrize('descr', [
+	(2, 0, 0), (3, 0, 0), (4, 0, 0), (5, 0, 0),
+	(1, 0, 1), (2, 0, 1), (3, 0, 1), (4, 0, 1),
+	(1, 1, 0), (2, 1, 0), (3, 1, 0), (4, 1, 0),
+	# (3, 0, 0)
+])
+def test_motor_split(descr):
+	"""test splitting of motors relative to a specific point"""
+	np.random.seed(10)
+	ga = NumpyContext(descr)
+	N = 10
+	m = random_motor(ga, (N,))
+
+	# test for both canonical origin, or some arbitrary point
+	o = ga.multivector.basis()[-1].dual()
+	os = [o, random_motor(ga, (N,)) >> o]
+
+	for o in os:
+		t, r = m.motor_split(o)
+		# print(m)
+		# print(t.subspace)
+		# print(t.subspace.is_subalgebra)
+		# print(t.subspace)
+		# print(r.subspace.is_subalgebra)
+		# check that it is a valid decomposition of the original motor
+		assert_close(t * r, m)
+		# check that r leaves o invariant
+		assert_close(r >> o, o, atol=1e-6)
+		# in d<6, both of these should be simple motors
+		assert_close(t * ~t, 1, atol=1e-6)
+		assert_close(r * ~r, 1, atol=1e-6)
+		# and their respective bivectors are orthogonal
+		assert_close(t.motor_log().inner(r.motor_log()), 0)
 
diff --git a/numga/multivector/extension/test/test_logexp.py b/numga/multivector/extension/test/test_logexp.py
index 73f8563..eff950d 100644
--- a/numga/multivector/extension/test/test_logexp.py
+++ b/numga/multivector/extension/test/test_logexp.py
@@ -1,4 +1,4 @@
-
+import numpy as np
 import pytest
 
 from numga.backend.numpy.context import NumpyContext
@@ -10,12 +10,13 @@
 
 @pytest.mark.parametrize('descr', [
 	(2, 0, 0), (1, 0, 1), (1, 1, 0), #(0, 1, 1),
-	(3, 0, 0), (2, 0, 1), (2, 1, 0), (1, 1, 1), (1, 0, 2), (1, 2, 0),
+	(3, 0, 0), (2, 0, 1), (2, 1, 0), (1, 1, 1), (1, 0, 2), #(1, 2, 0),
 	(4, 0, 0), (3, 0, 1), (3, 1, 0), (2, 1, 1), (2, 0, 2), #(2, 2, 0),
 	(5, 0, 0), (4, 0, 1), (4, 1, 0), (3, 1, 1),
 ])
 def test_bisect(descr):
 	"""Test roundtrip qualities of bisection based log and exp"""
+	np.random.seed(0)
 	print()
 	print(descr)
 	algebra = Algebra.from_pqr(*descr)
@@ -25,12 +26,12 @@ def test_bisect(descr):
 
 	# check exact inverse props
 	r = m.motor_log().exp()
-	npt.assert_allclose(m.values, r.values, atol=1e-9)
+	assert_close(m, r)
 
 	# check for low iteration count
 	from numga.multivector.extension.logexp import exp_bisect, motor_log_bisect
 	r = exp_bisect(motor_log_bisect(m, 2), 2)
-	npt.assert_allclose(m.values, r.values, atol=1e-9)
+	assert_close(m, r)
 
 
 def test_interpolate():
@@ -89,6 +90,4 @@ def test_exp_high_dim():
 		r = h.squared() * i
 	# else:
 	# 	r = h.squared().solve(q)
-	res = r.symmetric_reverse_product() - 1
-	print(res)
-	npt.assert_allclose(res.values, 0, atol=1e-9)
\ No newline at end of file
+	assert_close(r.symmetric_reverse_product(), 1)
diff --git a/numga/multivector/extension/test/test_norms.py b/numga/multivector/extension/test/test_norms.py
index 233835c..18ddc70 100644
--- a/numga/multivector/extension/test/test_norms.py
+++ b/numga/multivector/extension/test/test_norms.py
@@ -27,27 +27,25 @@ def test_study(descr):
 	s = m.symmetric_reverse_product()
 	assert s.subspace.is_study
 
+	mn = m.normalized()
+	assert_close(mn * ~mn, 1)
+
 	print('inverse')
 	Is = s.inverse()
-	r = Is * s - 1
-	npt.assert_allclose(r.values, 0, atol=1e-9)
+	assert_close(Is * s, 1)
 
 	print('square root')
 	rs = s.square_root()
-	r = rs * rs - s
-	npt.assert_allclose(r.values, 0, atol=1e-9)
+	assert_close(rs * rs, s)
 	zero = s * 0
 	rs = zero.square_root()
-	npt.assert_allclose(rs.values, 0, atol=1e-9)
+	assert_close(rs, 0)
 
 	print('inverse square root')
 	Irs = s.inverse().square_root()
-
 	Irs2 = s.square_root().inverse()
-	npt.assert_allclose(Irs2.values, Irs.values, atol=1e-6)
-
-	r = s * Irs * Irs - 1
-	npt.assert_allclose(r.values, 0, atol=1e-9)
+	assert_close(Irs2, Irs, atol=1e-6)
+	assert_close(s * Irs * Irs, 1)
 
 
 def test_motor_normalize_6d(descr=[6,0,0]):
diff --git a/numga/multivector/extension/test/test_optimized.py b/numga/multivector/extension/test/test_optimized.py
index a63c719..f4f0ffc 100644
--- a/numga/multivector/extension/test/test_optimized.py
+++ b/numga/multivector/extension/test/test_optimized.py
@@ -3,11 +3,12 @@
 import numpy as np
 import numpy.testing as npt
 from numga.backend.numpy.context import NumpyContext as Context
+from numga.multivector.test.util import *
 
 
-def test_normalize_extended():
+def test_normalize():
 	ga = Context('x+y+z+w0')
-	m = ga.multivector.motor(np.random.normal(size=(20, 8)))
+	m = random_non_motor(ga, (10,))
 	print()
 	print(m)
 	print(m.norm())
@@ -23,4 +24,29 @@ def test_normalize_extended():
 	print(n2)
 	print(n2.norm())
 
-	npt.assert_allclose(n1.values, n2.values, atol=1e-12)
+	assert_close(n1, n2)
+
+
+def test_exp():
+	np.random.seed(0)
+	from time import time
+	from numga.multivector.extension import optimized
+	from numga.multivector.extension.logexp import exp_bisect
+	from numga.backend.numpy.operator import NumpyEinsumOperator as Operator
+
+	ga = Context('x+y+z+w0', otype=Operator)
+	B = ga.subspace.bivector()
+
+	subspaces = [B, B.degenerate(), B.nondegenerate()]
+	for s in subspaces:
+		print(s)
+		b = random_subspace(ga, s, (10, 10, ))
+		b = random_subspace(ga, s, ( ))
+		t = time()
+		m1 = exp_bisect(b)
+		print(time() - t)
+		t = time()
+		m2 = b.exp()
+		print(time() - t)
+
+		assert_close(m1, m2, atol=1e-8)
diff --git a/numga/multivector/extension/test/test_roots.py b/numga/multivector/extension/test/test_roots.py
index 09a4e71..27e15bf 100644
--- a/numga/multivector/extension/test/test_roots.py
+++ b/numga/multivector/extension/test/test_roots.py
@@ -1,3 +1,4 @@
+import numpy as np
 import pytest
 
 from numga.backend.numpy.context import NumpyContext
@@ -34,10 +35,7 @@ def test_study_sqrt(descr):
 		if x.subspace.is_study:
 			print(x.subspace)
 			sqrt = x.square_root()
-			# print(sqrt)
-			r = sqrt.squared() - x
-			npt.assert_allclose(r.values, 0, atol=1e-9)
-			print('PASS')
+			assert_close(sqrt.squared(), x)
 
 
 def test_object_square_root():
@@ -46,15 +44,25 @@ def test_object_square_root():
 	algebra = Algebra.from_pqr(2, 0, 0)
 	context = NumpyContext(algebra)
 	v = random_subspace(context, algebra.subspace.bivector(), (10,))
-	sqrt = v.square_root()
-	r = sqrt.squared() - v
-	npt.assert_allclose(r.values, 0, atol=1e-9)
-
+	assert_close(v.square_root().squared(), v)
 
 	# we can take real square roots of 1-vectors in a negative sig
 	algebra = Algebra.from_pqr(0, 2, 0)
 	context = NumpyContext(algebra)
 	v = random_subspace(context, algebra.subspace.vector(), (10,))
-	sqrt = v.square_root()
-	r = sqrt.squared() - v
-	npt.assert_allclose(r.values, 0, atol=1e-9)
+	assert_close(v.square_root().squared(), v)
+
+
+def test_motor_roots():
+	np.random.seed(0)
+	algebra = Algebra.from_pqr(3, 0, 1)
+	context = NumpyContext(algebra)
+	m = random_motor(context, (10,))
+	assert_close(m.motor_square_root().squared(), m)
+
+	# construct some translators, which can have a specialized sqrt
+	m = m.select.bivector().degenerate() + 1
+	assert_close(m.motor_square_root().squared(), m)
+
+	assert_close(m.motor_geometric_mean(m), m)
+
diff --git a/numga/multivector/test/test_multivector.py b/numga/multivector/test/test_multivector.py
index 6264599..f302ab0 100644
--- a/numga/multivector/test/test_multivector.py
+++ b/numga/multivector/test/test_multivector.py
@@ -3,7 +3,7 @@
 import numpy.testing as npt
 from numga.algebra.algebra import Algebra
 from numga.backend.numpy.context import NumpyContext
-from numga.multivector.test.util import random_non_motor, random_subspace
+from numga.multivector.test.util import *
 
 
 def test_basic():
@@ -32,5 +32,5 @@ def test_5d_rotor_quality(descr):
 	r = random_non_motor(ga, (N,))
 	r = r.normalized()  # enforce R~R=1
 	q = r.full_sandwich(v)  # use full sandwich which will yield grade-5 elements
-	npt.assert_allclose(q.select[5].values, 0, atol=1e-10)
+	assert_close(q.select[5], 0)
 
diff --git a/numga/multivector/test/test_normalize.py b/numga/multivector/test/test_normalize.py
index a237028..06199a5 100644
--- a/numga/multivector/test/test_normalize.py
+++ b/numga/multivector/test/test_normalize.py
@@ -93,8 +93,7 @@ def test_multivector_normalize(descr):
 		if not s.symmetric_reverse().equals.empty():
 			x = random_subspace(context, s, (1,))
 			y = x.normalized()
-			r = y.symmetric_reverse_product() - 1
-			npt.assert_allclose(r.values, 0, atol=1e-6)
+			assert_close(y.symmetric_reverse_product(), 1, atol=1e-6)
 
 
 @pytest.mark.parametrize('descr', [
@@ -115,6 +114,8 @@ def test_motor_normalize_numerical_equivalence(descr):
 	m = v.normalized()
 	n = normalize_motor(v, inner=6, outer=1)
 
+	assert_close(m, n, atol=1e-9)
+
 	v0, v1 = motor_properties(m)
 	npt.assert_allclose(v0, 0, atol=1e-9)
 	npt.assert_allclose(v1, 0, atol=1e-9)
@@ -122,5 +123,3 @@ def test_motor_normalize_numerical_equivalence(descr):
 	v0, v1 = motor_properties(n)
 	npt.assert_allclose(v0, 0, atol=1e-9)
 	npt.assert_allclose(v1, 0, atol=1e-9)
-
-	npt.assert_allclose(m.values, n.values, atol=1e-9)
diff --git a/numga/multivector/test/test_numerical_normalize.py b/numga/multivector/test/test_numerical_normalize.py
index da6ebac..972d599 100644
--- a/numga/multivector/test/test_numerical_normalize.py
+++ b/numga/multivector/test/test_numerical_normalize.py
@@ -30,8 +30,7 @@ def test_numerical_inverse_sandwich(descr):
 	# should admit an inverse sandwich solve
 	s = e.symmetric_reverse_product()
 	isr = inverse_sandwich(s, n_iter=10)
-	r = (isr >> s) - 1
-	npt.assert_allclose(r.values, 0, atol=1e-9)
+	assert_close(isr >> s, 1)
 
 
 @pytest.mark.parametrize('descr', [
diff --git a/numga/multivector/test/util.py b/numga/multivector/test/util.py
index 4b3f785..7226929 100644
--- a/numga/multivector/test/util.py
+++ b/numga/multivector/test/util.py
@@ -84,7 +84,7 @@ def motor_properties(m):
 	v0 = np.linalg.norm(v0.values, axis=-1)
 	v1 = np.linalg.norm(v1.values, axis=-1)
 	return v0, v1
-	return np.allclose(v0, 0, atol=1e-6), np.allclose(v1, 0, atol=1e-6)
+	# return np.allclose(v0, 0, atol=1e-6), np.allclose(v1, 0, atol=1e-6)
 
 
 def all_grade_combinations(algebra):
@@ -96,5 +96,9 @@ def all_grade_combinations(algebra):
 
 
 def check_inverse(x, i, atol=1e-9):
-	assert np.allclose((x * i - 1).values, 0, atol=atol)
-	assert np.allclose((i * x - 1).values, 0, atol=atol)
+	assert_close(x*i, 1, atol)
+	assert_close(i*x, 1, atol)
+
+
+def assert_close(a, b, atol=1e-9):
+	assert np.allclose((a-b).values, 0, atol=atol)
\ No newline at end of file
diff --git a/numga/multivector/types.py b/numga/multivector/types.py
index 00d84da..d75e1ee 100644
--- a/numga/multivector/types.py
+++ b/numga/multivector/types.py
@@ -12,6 +12,10 @@ class Vector:
 	pass
 
 
+class AntiVector:
+	pass
+
+
 class BiVector:
 	pass
 
diff --git a/numga/subspace/subspace.py b/numga/subspace/subspace.py
index bdb5a9c..6944f5e 100644
--- a/numga/subspace/subspace.py
+++ b/numga/subspace/subspace.py
@@ -45,6 +45,9 @@ def grade(self) -> int:
 
 	def __len__(self):
 		return len(self.blades)
+	@property
+	def is_blade(self):
+		return len(self) == 1
 
 	def __repr__(self):
 		# return self.pretty_str
@@ -178,6 +181,8 @@ def wedge(self, other) -> "SubSpace":
 	@cache
 	def inner(self, other) -> "SubSpace":
 		return self.algebra.operator.inner(self, other).subspace
+	# def reject(self, other):
+	# 	return self.difference()
 
 	@cache
 	def squared(self) -> "SubSpace":