diff --git a/galgebra/mv.py b/galgebra/mv.py
index e8f7ea1f..cae12bde 100644
--- a/galgebra/mv.py
+++ b/galgebra/mv.py
@@ -935,25 +935,37 @@ def is_base(self) -> bool:
 
     def is_versor(self) -> bool:
         """
-        Test for versor (geometric product of vectors)
+        Test if `self` is a versor.
 
-        This follows Leo Dorst's test for a versor.
-        Leo Dorst, 'Geometric Algebra for Computer Science,' p.533
-        Sets self.versor_flg and returns value
-        """
+        A versor is defined as a geometric product of finitely-many nonnull vectors.
+
+        According to Greg Grunberg's analysis, the following are necessary and sufficient
+        conditions for a multivector V to be a versor:
+
+        1. V*V.rev() must be a nonzero scalar
+        2. V.g_invol()*x*V.rev() must be a vector whenever x is a vector
+
+        Sets self.versor_flg and returns value.
 
+        See https://github.com/pygae/galgebra/issues/533 for more discussions.
+        """
         if self.versor_flg is not None:
             return self.versor_flg
+
         self.characterise_Mv()
-        self.versor_flg = False
         self_rev = self.rev()
-        # see if self*self.rev() is a scalar
-        test = self*self_rev
-        if not test.is_scalar():
+
+        # Test condition 1: V*V.rev() must be a nonzero scalar
+        VVrev = self * self_rev
+        if not (VVrev.is_scalar() and not VVrev.is_zero()):
+            self.versor_flg = False
             return self.versor_flg
-        # see if self*x*self.rev() returns a vector for x an arbitrary vector
-        test = self * self.Ga._XOX * self.rev()
-        self.versor_flg = test.is_vector()
+
+        # Test condition 2: V.g_invol()*x*V.rev() must be a vector
+        # where x is a generic vector
+        # and self.Ga._XOX is a generic vector in self's geometric algebra
+        VinvolXVrev = self.g_invol() * self.Ga._XOX * self_rev
+        self.versor_flg = VinvolXVrev.is_vector()
         return self.versor_flg
 
     r'''
diff --git a/justfile b/justfile
new file mode 100644
index 00000000..bf81c64e
--- /dev/null
+++ b/justfile
@@ -0,0 +1,33 @@
+# This is a Justfile. It is a file that contains commands that can be run with the `just` command.
+# To install `just`, see https://github.com/casey/just
+
+default:
+  just --list
+
+prep-uv:
+  #!/usr/bin/env bash
+  # Already installed, then exit
+  which uv && exit 0
+  # On Mac, install with Homebrew
+  which brew && brew install uv
+  # On *nix
+  which uv || (curl -LsSf https://astral.sh/uv/install.sh | sh)
+  # On Windows, see https://docs.astral.sh/uv/getting-started/installation/
+
+prep-dt:
+  #!/usr/bin/env bash
+  which uv || (echo "Please install 'uv' in your preferred way, e.g. just prep-uv" && exit 1)
+  uv venv --python=python3.11 --seed
+  source .venv/bin/activate
+  pip install -e .
+  pip install -r test_requirements.txt
+
+# Run tests during development
+# depends on prep-dt
+dt TESTS="*":
+  #!/usr/bin/env bash
+  source .venv/bin/activate
+  pytest test/test_{{TESTS}}.py
+
+lint:
+  flake8 -v
diff --git a/test/test_versors.py b/test/test_versors.py
new file mode 100644
index 00000000..43ba95d7
--- /dev/null
+++ b/test/test_versors.py
@@ -0,0 +1,68 @@
+import unittest
+from galgebra.ga import Ga
+from sympy import S
+
+
+class TestVersors(unittest.TestCase):
+    def setUp(self):
+        # Set up different geometric algebras for testing
+
+        # G(3,0) - Euclidean 3-space
+        self.g3d = Ga("e1 e2 e3", g=[1, 1, 1])
+
+        # G(1,3) - Spacetime algebra
+        self.sta = Ga("e0 e1 e2 e3", g=[1, -1, -1, -1])
+
+        # Initialize basis vectors
+        e1, e2, e3 = self.g3d.mv()
+        e0, e1, e2, e3 = self.sta.mv()
+
+    def test_basis_vectors_are_versors(self):
+        """Individual basis vectors should be versors"""
+        e1, e2, e3 = self.g3d.mv()
+        self.assertTrue(e1.is_versor())
+        self.assertTrue(e2.is_versor())
+        self.assertTrue(e3.is_versor())
+
+    def test_mixed_grades_not_versor(self):
+        """A sum of different grades cannot be a versor"""
+        e1, e2, e3 = self.g3d.mv()
+        mixed = e1 + e2 * e3  # vector + bivector
+        self.assertFalse(mixed.is_versor())
+
+    def test_dorst_counterexample(self):
+        """Test Greg1950's counterexample from the spacetime algebra"""
+        e0, e1, e2, e3 = self.sta.mv()
+        V = S.One + self.sta.I()  # 1 + I
+        self.assertFalse(V.is_versor())
+
+    def test_rotors_are_versors(self):
+        """Test that rotors (even-grade versors) are properly detected"""
+        e1, e2, e3 = self.g3d.mv()
+        rotor = S.One + e1 * e2  # 1 + e1e2
+        self.assertTrue(rotor.is_versor())
+
+    def test_null_is_not_versor(self):
+        """Test that null multivectors are not versors"""
+        e1, e2, e3 = self.g3d.mv()
+        null_mv = e1 + e1 * self.g3d.I()  # v + vI (squares to 0)
+        self.assertFalse(null_mv.is_versor())
+
+    def test_scalar_versor(self):
+        """Test that nonzero scalars are versors"""
+        self.assertTrue(self.g3d.mv(2).is_versor())
+        self.assertFalse(self.g3d.mv(0).is_versor())
+
+    def test_bivector_exponential(self):
+        """Test that exponentials of bivectors are versors"""
+        e1, e2, e3 = self.g3d.mv()
+        B = e1 * e2  # bivector
+        rotor = (B/2).exp()  # exp(B/2) is a rotor
+        self.assertTrue(rotor.is_versor())
+
+    def test_product_of_vectors(self):
+        """Test that products of vectors are versors"""
+        e1, e2, e3 = self.g3d.mv()
+        v1 = e1 + 2*e2  # vector
+        v2 = 3*e1 - e3  # vector
+        self.assertTrue((v1 * v2).is_versor())