From d364b13c539f92fc4485bc7d6bd4746cbd3e6c93 Mon Sep 17 00:00:00 2001
From: Michael Boyle <michael.oliver.boyle@gmail.com>
Date: Wed, 7 Apr 2021 12:38:52 -0400
Subject: [PATCH 01/10] Reorganize some tests

---
 tests/test_algebra.py | 429 ++++++++++++++++++++++++++++++++++++++++++
 tests/test_array.py   | 427 -----------------------------------------
 2 files changed, 429 insertions(+), 427 deletions(-)

diff --git a/tests/test_algebra.py b/tests/test_algebra.py
index 0b3f59e..fa64c62 100644
--- a/tests/test_algebra.py
+++ b/tests/test_algebra.py
@@ -36,3 +36,432 @@ def test_basis_multiplication():
     assert j*j == -one
     assert k*k == -one
     assert i*j*k == -one
+
+
+def test_array_ufunc(array):
+    np.random.seed(1234)
+
+    q = array(np.random.rand(1, 3, 4))
+    with pytest.raises(NotImplementedError):
+        np.exp(q, extra_arg=True)
+
+    with pytest.raises(NotImplementedError):
+        np.negative.at(q, [0, 1])
+
+    p = array(np.random.rand(17, 3, 4))
+    q = array(np.random.rand(1, 3, 4))
+    pq1 = np.add(p, q)
+    assert isinstance(pq1, array)
+    assert pq1.shape == (17, 3, 4)
+    assert np.array_equal(np.add(p.ndarray, q.ndarray), pq1.ndarray)
+    pq2 = array(np.empty((17, 3, 4)))
+    np.add(p, q, out=pq2)
+    assert np.array_equal(pq1, pq2)
+    assert isinstance(pq2, array)
+
+    p = array(np.random.rand(17, 3, 4))
+    q = np.random.rand(1, 3)
+    pq1 = np.multiply(p, q)
+    assert isinstance(pq1, array)
+    assert pq1.shape == (17, 3, 4)
+    pq2 = array(np.empty((17, 3, 4)))
+    np.multiply(p, q, out=pq2)
+    assert np.array_equal(pq1, pq2)
+    assert isinstance(pq2, array)
+
+    p = np.random.rand(1, 3)
+    q = array(np.random.rand(17, 3, 4))
+    pq1 = np.multiply(p, q)
+    assert isinstance(pq1, array)
+    assert pq1.shape == (17, 3, 4)
+    pq2 = array(np.empty((17, 3, 4)))
+    np.multiply(p, q, out=pq2)
+    assert np.array_equal(pq1, pq2)
+    assert isinstance(pq2, array)
+
+    p = np.random.rand(1, 3)
+    q = np.random.rand(17, 3, 4)
+    s = np.random.rand(17, 3)
+    pq1 = array(q).__array_ufunc__(np.multiply, "__call__", p, q)
+    assert pq1 == NotImplemented
+    qneg = array(q).__array_ufunc__(np.negative, "__call__", q)
+    assert qneg == NotImplemented
+    qabs = array(q).__array_ufunc__(np.absolute, "__call__", q)
+    assert qabs == NotImplemented
+    qs = array(q).__array_ufunc__(np.float_power, "__call__", q, s)
+    assert qs == NotImplemented
+    pq1 = array(q).__array_ufunc__(np.equal, "__call__", p, q)
+    assert pq1 == NotImplemented
+    qfin = array(q).__array_ufunc__(np.isfinite, "__call__", q)
+    assert qfin == NotImplemented
+
+    q = array(np.random.rand(17, 3, 4))
+    qneg = np.negative(q)
+    assert isinstance(qneg, array)
+    assert qneg.shape == q.shape
+    assert np.array_equal(np.negative(q.ndarray), qneg.ndarray)
+    qneg2 = np.empty(q.shape)
+    np.negative(q, out=qneg2)
+    assert np.array_equal(qneg, qneg2)
+    assert isinstance(qneg2, np.ndarray)
+    qneg2 = array(np.empty(q.shape))
+    np.negative(q, out=qneg2.ndarray)
+    assert np.array_equal(qneg, qneg2)
+    assert isinstance(qneg2, array)
+    qneg2 = array(np.empty(q.shape))
+    np.negative(q, out=qneg2)
+    assert np.array_equal(qneg, qneg2)
+    assert isinstance(qneg2, array)
+
+    p = np.random.rand(1, 3)
+    q = array(np.random.rand(17, 3, 4))
+    qp1 = np.float_power(q, p)
+    assert isinstance(qp1, array)
+    assert qp1.shape == (17, 3, 4)
+    qp2 = array(np.empty((17, 3, 4)))
+    np.float_power(q, p, out=qp2)
+    assert np.array_equal(qp1, qp2)
+    assert isinstance(qp2, array)
+
+    q = array(np.random.rand(17, 3, 4))
+    qabs = np.absolute(q)
+    assert isinstance(qabs, np.ndarray) and not isinstance(qabs, array)
+    assert qabs.shape == (17, 3)
+    qabs2 = np.empty((17, 3))
+    np.absolute(q, out=qabs2)
+    assert np.array_equal(qabs, qabs2)
+    assert isinstance(qabs2, np.ndarray) and not isinstance(qabs, array)
+    q = array(np.random.rand(17, 3, 4, 4))
+    qabs = array(np.empty((17, 3, 4)))
+    np.absolute(q, out=qabs)
+    assert np.array_equal(qabs, np.absolute(q))
+    assert isinstance(qabs2, np.ndarray) and isinstance(qabs, array)
+
+    p = array(np.random.rand(17, 3, 4))
+    q = array(np.random.rand(1, 3, 4))
+    pq1 = np.equal(p, q)
+    assert isinstance(pq1, np.ndarray) and not isinstance(pq1, array)
+    assert pq1.shape == (17, 3)
+    assert np.array_equal(np.all(np.equal(p.ndarray, q.ndarray), axis=-1), pq1)
+    pq2 = np.empty((17, 3), dtype=bool)
+    np.equal(p, q, out=pq2)
+    assert np.array_equal(pq1, pq2)
+    assert isinstance(pq2, np.ndarray) and not isinstance(pq2, array)
+    assert pq2.shape == (17, 3)
+    p = array(np.random.rand(17, 3, 4, 4))
+    q = array(np.random.rand(17, 3, 4, 4))
+    pq = array(np.empty((17, 3, 4)))
+    np.equal(p, q, out=pq)
+    assert isinstance(pq, np.ndarray) and isinstance(pq, array)
+
+    q = array(np.random.rand(17, 3, 4))
+    qfin = np.isfinite(q)
+    assert isinstance(qfin, np.ndarray) and not isinstance(qfin, array)
+    assert qfin.shape == (17, 3)
+    assert np.array_equal(np.all(np.isfinite(q.ndarray), axis=-1), qfin)
+    qfin2 = np.empty((17, 3), dtype=bool)
+    np.isfinite(q, out=qfin2)
+    assert np.array_equal(qfin, qfin2)
+    assert isinstance(qfin2, np.ndarray) and not isinstance(qfin2, array)
+    assert qfin2.shape == (17, 3)
+    q = array(np.random.rand(17, 3, 4, 4))
+    qfin = array(np.empty((17, 3, 4)))
+    np.isfinite(q, out=qfin)
+    assert isinstance(qfin, np.ndarray) and isinstance(qfin, array)
+
+    # Test the NotImplemented ufuncs
+    implemented_ufuncs = [
+        np.add, np.subtract, np.multiply, np.divide, np.true_divide,
+        np.bitwise_or, np.bitwise_xor, np.right_shift, np.left_shift,
+        np.negative, np.positive, np.conj, np.conjugate, np.invert,
+        np.exp, np.log, np.sqrt, np.square, np.reciprocal,
+        np.float_power,
+        np.absolute,
+        np.not_equal, np.equal, np.logical_and, np.logical_or,
+        np.isfinite, np.isinf, np.isnan,
+    ]
+    for k in dir(np):
+        attr = getattr(np, k)
+        if isinstance(attr, np.ufunc) and attr not in implemented_ufuncs:
+            p = array(np.random.rand(17, 3, 4))
+            if attr.nin == 1:
+                with pytest.raises(TypeError):
+                    attr(p)
+            elif attr.nin == 2:
+                q = array(np.random.rand(17, 3, 4))
+                with pytest.raises(TypeError):
+                    attr(p, q)
+            else:  # pragma: no cover
+                raise ValueError(f"Unexpected number of input arguments for {k}")
+
+
+# Unary bool returners
+def test_quaternion_nonzero(Qs, Q_names, Q_conditions, array):
+    Qs = array(Qs.ndarray)
+    assert np.nonzero(Qs[Q_names.q_0])[0].size == 0  # Do this one explicitly, to not use circular logic
+    assert np.nonzero(Qs[Q_names.q_1])[0].size > 0  # Do this one explicitly, to not use circular logic
+    for q in Qs[Q_conditions.zero]:
+        assert np.nonzero(q)[0].size == 0
+    for q in Qs[Q_conditions.nonzero]:
+        assert np.nonzero(q)[0].size > 0
+
+
+def test_quaternion_isnan(Qs, Q_names, Q_conditions, array):
+    Qs = array(Qs.ndarray)
+    assert not np.isnan(Qs[Q_names.q_0])  # Do this one explicitly, to not use circular logic
+    assert not np.isnan(Qs[Q_names.q_1])  # Do this one explicitly, to not use circular logic
+    assert np.isnan(Qs[Q_names.q_nan1])  # Do this one explicitly, to not use circular logic
+    for q in Qs[Q_conditions.nan]:
+        assert np.isnan(q)
+    for q in Qs[Q_conditions.nonnan]:
+        assert not np.isnan(q)
+
+
+def test_quaternion_isinf(Qs, Q_names, Q_conditions, array):
+    Qs = array(Qs.ndarray)
+    assert not np.isinf(Qs[Q_names.q_0])  # Do this one explicitly, to not use circular logic
+    assert not np.isinf(Qs[Q_names.q_1])  # Do this one explicitly, to not use circular logic
+    assert np.isinf(Qs[Q_names.q_inf1])  # Do this one explicitly, to not use circular logic
+    assert np.isinf(Qs[Q_names.q_minf1])  # Do this one explicitly, to not use circular logic
+    for q in Qs[Q_conditions.inf]:
+        assert np.isinf(q)
+    for q in Qs[Q_conditions.noninf]:
+        assert not np.isinf(q)
+
+
+def test_quaternion_isfinite(Qs, Q_names, Q_conditions, array):
+    Qs = array(Qs.ndarray)
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        assert not np.isfinite(Qs[Q_names.q_inf1])  # Do this one explicitly, to not use circular logic
+        assert not np.isfinite(Qs[Q_names.q_minf1])  # Do this one explicitly, to not use circular logic
+    assert not np.isfinite(Qs[Q_names.q_nan1])  # Do this one explicitly, to not use circular logic
+    assert np.isfinite(Qs[Q_names.q_0])  # Do this one explicitly, to not use circular logic
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        for q in Qs[Q_conditions.nonfinite]:
+            assert not np.isfinite(q)
+    for q in Qs[Q_conditions.finite]:
+        assert np.isfinite(q)
+
+
+# Binary bool returners
+def test_quaternion_equal(Qs, Q_names, Q_conditions, array):
+    Qs = array(Qs.ndarray)
+    for j in Q_conditions.nonnan:
+        assert Qs[j] == Qs[j]  # self equality
+        for k in range(len(Qs)):  # non-self inequality
+            assert (j == k) or (not (Qs[j] == Qs[k]))
+    for q in Qs:
+        for p in Qs[Q_conditions.nan]:
+            assert not q == p  # nan should never equal anything
+
+
+def test_quaternion_not_equal(Qs, Q_names, Q_conditions, array):
+    Qs = array(Qs.ndarray)
+    for j in Q_conditions.nonnan:
+        assert not (Qs[j] != Qs[j])  # self non-not_equality
+        for k in Q_conditions.nonnan:  # non-self not_equality
+            assert (j == k) or (Qs[j] != Qs[k])
+    for q in Qs:
+        for p in Qs[Q_conditions.nan]:
+            assert q != p  # nan should never equal anything
+
+
+# Unary float returners
+def test_quaternion_absolute(Qs, Q_names, Q_conditions, on_windows, array):
+    Qs = array(Qs.ndarray)
+    for abs_func in [np.abs, lambda q: q.abs]:
+        for q in Qs[Q_conditions.nan]:
+            assert np.isnan(abs_func(q))
+        for q in Qs[Q_conditions.inf]:
+            if on_windows:  # pragma: no cover
+                assert np.isinf(abs_func(q)) or np.isnan(abs_func(q))
+            else:
+                assert np.isinf(abs_func(q))
+        for q, a in [(Qs[Q_names.q_0], 0.0), (Qs[Q_names.q_1], 1.0), (Qs[Q_names.x], 1.0), (Qs[Q_names.y], 1.0), (Qs[Q_names.z], 1.0),
+                     (Qs[Q_names.Q], np.sqrt(Qs[Q_names.Q].w ** 2 + Qs[Q_names.Q].x ** 2 + Qs[Q_names.Q].y ** 2 + Qs[Q_names.Q].z ** 2)),
+                     (Qs[Q_names.Qbar], np.sqrt(Qs[Q_names.Q].w ** 2 + Qs[Q_names.Q].x ** 2 + Qs[Q_names.Q].y ** 2 + Qs[Q_names.Q].z ** 2))]:
+            assert np.allclose(abs_func(q), a)
+
+
+def test_quaternion_norm(Qs, Q_names, Q_conditions, on_windows, array):
+    Qs = array(Qs.ndarray)
+    for q in Qs[Q_conditions.nan]:
+        assert np.isnan(q.norm)
+    for q in Qs[Q_conditions.inf]:
+        if on_windows:  # pragma: no cover
+            assert np.isinf(q.norm) or np.isnan(q.norm)
+        else:
+            assert np.isinf(q.norm)
+    for q, a in [(Qs[Q_names.q_0], 0.0), (Qs[Q_names.q_1], 1.0), (Qs[Q_names.x], 1.0), (Qs[Q_names.y], 1.0), (Qs[Q_names.z], 1.0),
+                 (Qs[Q_names.Q], Qs[Q_names.Q].w ** 2 + Qs[Q_names.Q].x ** 2 + Qs[Q_names.Q].y ** 2 + Qs[Q_names.Q].z ** 2),
+                 (Qs[Q_names.Qbar], Qs[Q_names.Q].w ** 2 + Qs[Q_names.Q].x ** 2 + Qs[Q_names.Q].y ** 2 + Qs[Q_names.Q].z ** 2)]:
+        assert np.allclose(q.norm, a)
+
+
+# Unary quaternion returners
+def test_quaternion_negative(Qs, Q_names, Q_conditions, array):
+    Qs = array(Qs.ndarray)
+    assert -Qs[Q_names.Q] == Qs[Q_names.Qneg]
+    for q in Qs[Q_conditions.finite]:
+        assert -q == -1.0 * q
+    for q in Qs[Q_conditions.nonnan]:
+        assert -(-q) == q
+
+
+def test_quaternion_npconjugate(Qs, Q_names, Q_conditions, array):
+    Qs = array(Qs.ndarray)
+    assert np.conjugate(Qs[Q_names.Q]) == Qs[Q_names.Qbar]
+    for q in Qs[Q_conditions.nonnan]:
+        assert np.conjugate(q) == np.conj(q)
+        assert np.conjugate(np.conjugate(q)) == q
+        c = np.conjugate(q)
+        assert c.w == q.w
+        assert c.x == -q.x
+        assert c.y == -q.y
+        assert c.z == -q.z
+
+
+def test_quaternion_sqrt(Qs, Q_names, Q_conditions, array):
+    Qs = array(Qs.ndarray)
+    sqrt_precision = 2.e-15
+    # Test sqrt of basis elements
+    sqrthalf = np.sqrt(0.5)
+    assert np.array_equal(
+        np.sqrt(array(np.eye(4))).ndarray,
+        np.array([
+            [1, 0, 0, 0],
+            [sqrthalf, sqrthalf, 0, 0],
+            [sqrthalf, 0, sqrthalf, 0],
+            [sqrthalf, 0, 0, sqrthalf],
+        ])
+    )
+    # Test all my samples
+    one, i, j, k = tuple(array(np.eye(4)))
+    for q in Qs[Q_conditions.finitenonzero]:
+        assert np.allclose(np.sqrt(q) * np.sqrt(q), q, rtol=sqrt_precision)
+        # Ensure that non-unit quaternions are handled correctly
+        for s in [1, -1, 2, -2, 3.4, -3.4]:
+            for r in [one, i, j, k]:
+                srq = s*r*q
+                assert np.allclose(np.sqrt(srq) * np.sqrt(srq), srq, rtol=sqrt_precision)
+    # Test a huge batch of random quaternions
+    np.random.seed(1234)
+    a = array(np.random.uniform(-10, 10, size=10000*4).reshape((-1, 4)))
+    assert np.allclose(a, np.square(np.sqrt(a)), rtol=10*sqrt_precision, atol=0)
+    # Test some edge cases
+    _quaternion_resolution = 10 * np.finfo(float).resolution
+    for s in [-0.1, -1, -2, -3.4]:
+        q = np.sqrt(one * s)
+        assert np.array_equal(q, np.sqrt(-s)*i)
+        q = np.sqrt(one * s + array(0, *(0.05*np.sqrt(-s*_quaternion_resolution)*np.random.rand(3))))
+        assert np.array_equal(q, np.sqrt(-s)*i)
+
+
+def test_quaternion_square(Qs, Q_names, Q_conditions, array):
+    Qs = array(Qs.ndarray)
+    square_precision = 1.e-15
+    for q in Qs[Q_conditions.finite]:
+        assert (q*q - q**2).norm < square_precision
+        a = array([q])
+        assert (a**2 - array([q**2])).norm < square_precision
+
+
+def test_quaternion_log_exp(Qs, Q_names, Q_conditions, array):
+    Qs = array(Qs.ndarray)
+    qlogexp_precision = 4.e-15
+    zero = array([0, 0, 0, 0])
+    one, i, j, k = tuple(array(np.eye(4)))
+    assert (np.log(Qs[Q_names.Q]) - Qs[Q_names.Qlog]).abs < qlogexp_precision
+    assert (np.exp(Qs[Q_names.Q]) - Qs[Q_names.Qexp]).abs < qlogexp_precision
+    assert (np.exp(np.log(Qs[Q_names.Q])) - Qs[Q_names.Q]).abs < qlogexp_precision
+    assert (np.log(np.exp(Qs[Q_names.Q])) - Qs[Q_names.Q]).abs > qlogexp_precision  # Note order of operations!
+    assert np.log(one) == zero
+    assert np.log(i) == (np.pi / 2) * i
+    assert np.log(j) == (np.pi / 2) * j
+    assert np.log(k) == (np.pi / 2) * k
+    assert np.log(-one) == (np.pi) * i
+
+
+# Binary quat returners
+def test_quaternion_conjugate(Qs, Q_names, Q_conditions, array):
+    Qs = array(Qs.ndarray)
+    for conj_func in [np.conj, np.conjugate, lambda q: q.conj(), lambda q: q.conjugate()]:
+        assert conj_func(Qs[Q_names.Q]) == Qs[Q_names.Qbar]
+        for q in Qs[Q_conditions.nonnan]:
+            assert conj_func(conj_func(q)) == q
+            c = conj_func(q)
+            assert c.w == q.w
+            assert c.x == -q.x
+            assert c.y == -q.y
+            assert c.z == -q.z
+    for q in Qs[Q_conditions.nonnan]:
+        assert q.conjugate() == q.conj()
+        assert np.conjugate(q) == np.conj(q)
+
+
+def test_xor(array):
+    basis = one, i, j, k = tuple(array(np.eye(4)))
+    zero = 0 * one
+    assert one ^ one == one
+    assert one ^ i == i
+    assert one ^ j == j
+    assert one ^ k == k
+    assert i ^ one == i
+    assert i ^ i == zero
+    assert i ^ j == zero
+    assert i ^ k == zero
+    assert j ^ one == j
+    assert j ^ i == zero
+    assert j ^ j == zero
+    assert j ^ k == zero
+    assert k ^ one == k
+    assert k ^ i == zero
+    assert k ^ j == zero
+    assert k ^ k == zero
+
+
+def test_contractions(array):
+    basis = one, i, j, k = tuple(array(np.eye(4)))
+    zero = 0 * one
+
+    assert one << one == one
+    assert one << i == i
+    assert one << j == j
+    assert one << k == k
+    assert i << one == zero
+    assert i << i == -one
+    assert i << j == zero
+    assert i << k == zero
+    assert j << one == zero
+    assert j << i == zero
+    assert j << j == -one
+    assert j << k == zero
+    assert k << one == zero
+    assert k << i == zero
+    assert k << j == zero
+    assert k << k == -one
+
+    assert one >> one == one
+    assert one >> i == zero
+    assert one >> j == zero
+    assert one >> k == zero
+    assert i >> one == i
+    assert i >> i == -one
+    assert i >> j == zero
+    assert i >> k == zero
+    assert j >> one == j
+    assert j >> i == zero
+    assert j >> j == -one
+    assert j >> k == zero
+    assert k >> one == k
+    assert k >> i == zero
+    assert k >> j == zero
+    assert k >> k == -one
+
+    for a in basis:
+        for b in basis:
+            for c in basis:
+                assert np.allclose((a ^ b) | c, a | (b << c))
+                assert np.allclose(c | (b ^ a), (c >> b) | a)
diff --git a/tests/test_array.py b/tests/test_array.py
index f7a055f..65b505a 100644
--- a/tests/test_array.py
+++ b/tests/test_array.py
@@ -42,163 +42,6 @@ def test_array_finalize(array):
         q[1:3]
 
 
-def test_array_ufunc(array):
-    np.random.seed(1234)
-
-    q = array(np.random.rand(1, 3, 4))
-    with pytest.raises(NotImplementedError):
-        np.exp(q, extra_arg=True)
-
-    with pytest.raises(NotImplementedError):
-        np.negative.at(q, [0, 1])
-
-    p = array(np.random.rand(17, 3, 4))
-    q = array(np.random.rand(1, 3, 4))
-    pq1 = np.add(p, q)
-    assert isinstance(pq1, array)
-    assert pq1.shape == (17, 3, 4)
-    assert np.array_equal(np.add(p.ndarray, q.ndarray), pq1.ndarray)
-    pq2 = array(np.empty((17, 3, 4)))
-    np.add(p, q, out=pq2)
-    assert np.array_equal(pq1, pq2)
-    assert isinstance(pq2, array)
-
-    p = array(np.random.rand(17, 3, 4))
-    q = np.random.rand(1, 3)
-    pq1 = np.multiply(p, q)
-    assert isinstance(pq1, array)
-    assert pq1.shape == (17, 3, 4)
-    pq2 = array(np.empty((17, 3, 4)))
-    np.multiply(p, q, out=pq2)
-    assert np.array_equal(pq1, pq2)
-    assert isinstance(pq2, array)
-
-    p = np.random.rand(1, 3)
-    q = array(np.random.rand(17, 3, 4))
-    pq1 = np.multiply(p, q)
-    assert isinstance(pq1, array)
-    assert pq1.shape == (17, 3, 4)
-    pq2 = array(np.empty((17, 3, 4)))
-    np.multiply(p, q, out=pq2)
-    assert np.array_equal(pq1, pq2)
-    assert isinstance(pq2, array)
-
-    p = np.random.rand(1, 3)
-    q = np.random.rand(17, 3, 4)
-    s = np.random.rand(17, 3)
-    pq1 = array(q).__array_ufunc__(np.multiply, "__call__", p, q)
-    assert pq1 == NotImplemented
-    qneg = array(q).__array_ufunc__(np.negative, "__call__", q)
-    assert qneg == NotImplemented
-    qabs = array(q).__array_ufunc__(np.absolute, "__call__", q)
-    assert qabs == NotImplemented
-    qs = array(q).__array_ufunc__(np.float_power, "__call__", q, s)
-    assert qs == NotImplemented
-    pq1 = array(q).__array_ufunc__(np.equal, "__call__", p, q)
-    assert pq1 == NotImplemented
-    qfin = array(q).__array_ufunc__(np.isfinite, "__call__", q)
-    assert qfin == NotImplemented
-
-    q = array(np.random.rand(17, 3, 4))
-    qneg = np.negative(q)
-    assert isinstance(qneg, array)
-    assert qneg.shape == q.shape
-    assert np.array_equal(np.negative(q.ndarray), qneg.ndarray)
-    qneg2 = np.empty(q.shape)
-    np.negative(q, out=qneg2)
-    assert np.array_equal(qneg, qneg2)
-    assert isinstance(qneg2, np.ndarray)
-    qneg2 = array(np.empty(q.shape))
-    np.negative(q, out=qneg2.ndarray)
-    assert np.array_equal(qneg, qneg2)
-    assert isinstance(qneg2, array)
-    qneg2 = array(np.empty(q.shape))
-    np.negative(q, out=qneg2)
-    assert np.array_equal(qneg, qneg2)
-    assert isinstance(qneg2, array)
-
-    p = np.random.rand(1, 3)
-    q = array(np.random.rand(17, 3, 4))
-    qp1 = np.float_power(q, p)
-    assert isinstance(qp1, array)
-    assert qp1.shape == (17, 3, 4)
-    qp2 = array(np.empty((17, 3, 4)))
-    np.float_power(q, p, out=qp2)
-    assert np.array_equal(qp1, qp2)
-    assert isinstance(qp2, array)
-
-    q = array(np.random.rand(17, 3, 4))
-    qabs = np.absolute(q)
-    assert isinstance(qabs, np.ndarray) and not isinstance(qabs, array)
-    assert qabs.shape == (17, 3)
-    qabs2 = np.empty((17, 3))
-    np.absolute(q, out=qabs2)
-    assert np.array_equal(qabs, qabs2)
-    assert isinstance(qabs2, np.ndarray) and not isinstance(qabs, array)
-    q = array(np.random.rand(17, 3, 4, 4))
-    qabs = array(np.empty((17, 3, 4)))
-    np.absolute(q, out=qabs)
-    assert np.array_equal(qabs, np.absolute(q))
-    assert isinstance(qabs2, np.ndarray) and isinstance(qabs, array)
-
-    p = array(np.random.rand(17, 3, 4))
-    q = array(np.random.rand(1, 3, 4))
-    pq1 = np.equal(p, q)
-    assert isinstance(pq1, np.ndarray) and not isinstance(pq1, array)
-    assert pq1.shape == (17, 3)
-    assert np.array_equal(np.all(np.equal(p.ndarray, q.ndarray), axis=-1), pq1)
-    pq2 = np.empty((17, 3), dtype=bool)
-    np.equal(p, q, out=pq2)
-    assert np.array_equal(pq1, pq2)
-    assert isinstance(pq2, np.ndarray) and not isinstance(pq2, array)
-    assert pq2.shape == (17, 3)
-    p = array(np.random.rand(17, 3, 4, 4))
-    q = array(np.random.rand(17, 3, 4, 4))
-    pq = array(np.empty((17, 3, 4)))
-    np.equal(p, q, out=pq)
-    assert isinstance(pq, np.ndarray) and isinstance(pq, array)
-
-    q = array(np.random.rand(17, 3, 4))
-    qfin = np.isfinite(q)
-    assert isinstance(qfin, np.ndarray) and not isinstance(qfin, array)
-    assert qfin.shape == (17, 3)
-    assert np.array_equal(np.all(np.isfinite(q.ndarray), axis=-1), qfin)
-    qfin2 = np.empty((17, 3), dtype=bool)
-    np.isfinite(q, out=qfin2)
-    assert np.array_equal(qfin, qfin2)
-    assert isinstance(qfin2, np.ndarray) and not isinstance(qfin2, array)
-    assert qfin2.shape == (17, 3)
-    q = array(np.random.rand(17, 3, 4, 4))
-    qfin = array(np.empty((17, 3, 4)))
-    np.isfinite(q, out=qfin)
-    assert isinstance(qfin, np.ndarray) and isinstance(qfin, array)
-
-    # Test the NotImplemented ufuncs
-    implemented_ufuncs = [
-        np.add, np.subtract, np.multiply, np.divide, np.true_divide,
-        np.bitwise_or, np.bitwise_xor, np.right_shift, np.left_shift,
-        np.negative, np.positive, np.conj, np.conjugate, np.invert,
-        np.exp, np.log, np.sqrt, np.square, np.reciprocal,
-        np.float_power,
-        np.absolute,
-        np.not_equal, np.equal, np.logical_and, np.logical_or,
-        np.isfinite, np.isinf, np.isnan,
-    ]
-    for k in dir(np):
-        attr = getattr(np, k)
-        if isinstance(attr, np.ufunc) and attr not in implemented_ufuncs:
-            p = array(np.random.rand(17, 3, 4))
-            if attr.nin == 1:
-                with pytest.raises(TypeError):
-                    attr(p)
-            elif attr.nin == 2:
-                q = array(np.random.rand(17, 3, 4))
-                with pytest.raises(TypeError):
-                    attr(p, q)
-            else:  # pragma: no cover
-                raise ValueError(f"Unexpected number of input arguments for {k}")
-
-
 def test_repr(array):
     q = array(np.random.rand(17, 3, 4))
     assert repr(q) == 'quaternionic.' + repr(q.ndarray)
@@ -209,273 +52,3 @@ def test_str(array):
     assert str(q) == str(q.ndarray)
 
 
-# Unary bool returners
-def test_quaternion_nonzero(Qs, Q_names, Q_conditions, array):
-    Qs = array(Qs.ndarray)
-    assert np.nonzero(Qs[Q_names.q_0])[0].size == 0  # Do this one explicitly, to not use circular logic
-    assert np.nonzero(Qs[Q_names.q_1])[0].size > 0  # Do this one explicitly, to not use circular logic
-    for q in Qs[Q_conditions.zero]:
-        assert np.nonzero(q)[0].size == 0
-    for q in Qs[Q_conditions.nonzero]:
-        assert np.nonzero(q)[0].size > 0
-
-
-def test_quaternion_isnan(Qs, Q_names, Q_conditions, array):
-    Qs = array(Qs.ndarray)
-    assert not np.isnan(Qs[Q_names.q_0])  # Do this one explicitly, to not use circular logic
-    assert not np.isnan(Qs[Q_names.q_1])  # Do this one explicitly, to not use circular logic
-    assert np.isnan(Qs[Q_names.q_nan1])  # Do this one explicitly, to not use circular logic
-    for q in Qs[Q_conditions.nan]:
-        assert np.isnan(q)
-    for q in Qs[Q_conditions.nonnan]:
-        assert not np.isnan(q)
-
-
-def test_quaternion_isinf(Qs, Q_names, Q_conditions, array):
-    Qs = array(Qs.ndarray)
-    assert not np.isinf(Qs[Q_names.q_0])  # Do this one explicitly, to not use circular logic
-    assert not np.isinf(Qs[Q_names.q_1])  # Do this one explicitly, to not use circular logic
-    assert np.isinf(Qs[Q_names.q_inf1])  # Do this one explicitly, to not use circular logic
-    assert np.isinf(Qs[Q_names.q_minf1])  # Do this one explicitly, to not use circular logic
-    for q in Qs[Q_conditions.inf]:
-        assert np.isinf(q)
-    for q in Qs[Q_conditions.noninf]:
-        assert not np.isinf(q)
-
-
-def test_quaternion_isfinite(Qs, Q_names, Q_conditions, array):
-    Qs = array(Qs.ndarray)
-    with warnings.catch_warnings():
-        warnings.simplefilter("ignore")
-        assert not np.isfinite(Qs[Q_names.q_inf1])  # Do this one explicitly, to not use circular logic
-        assert not np.isfinite(Qs[Q_names.q_minf1])  # Do this one explicitly, to not use circular logic
-    assert not np.isfinite(Qs[Q_names.q_nan1])  # Do this one explicitly, to not use circular logic
-    assert np.isfinite(Qs[Q_names.q_0])  # Do this one explicitly, to not use circular logic
-    with warnings.catch_warnings():
-        warnings.simplefilter("ignore")
-        for q in Qs[Q_conditions.nonfinite]:
-            assert not np.isfinite(q)
-    for q in Qs[Q_conditions.finite]:
-        assert np.isfinite(q)
-
-
-# Binary bool returners
-def test_quaternion_equal(Qs, Q_names, Q_conditions, array):
-    Qs = array(Qs.ndarray)
-    for j in Q_conditions.nonnan:
-        assert Qs[j] == Qs[j]  # self equality
-        for k in range(len(Qs)):  # non-self inequality
-            assert (j == k) or (not (Qs[j] == Qs[k]))
-    for q in Qs:
-        for p in Qs[Q_conditions.nan]:
-            assert not q == p  # nan should never equal anything
-
-
-def test_quaternion_not_equal(Qs, Q_names, Q_conditions, array):
-    Qs = array(Qs.ndarray)
-    for j in Q_conditions.nonnan:
-        assert not (Qs[j] != Qs[j])  # self non-not_equality
-        for k in Q_conditions.nonnan:  # non-self not_equality
-            assert (j == k) or (Qs[j] != Qs[k])
-    for q in Qs:
-        for p in Qs[Q_conditions.nan]:
-            assert q != p  # nan should never equal anything
-
-
-# Unary float returners
-def test_quaternion_absolute(Qs, Q_names, Q_conditions, on_windows, array):
-    Qs = array(Qs.ndarray)
-    for abs_func in [np.abs, lambda q: q.abs]:
-        for q in Qs[Q_conditions.nan]:
-            assert np.isnan(abs_func(q))
-        for q in Qs[Q_conditions.inf]:
-            if on_windows:  # pragma: no cover
-                assert np.isinf(abs_func(q)) or np.isnan(abs_func(q))
-            else:
-                assert np.isinf(abs_func(q))
-        for q, a in [(Qs[Q_names.q_0], 0.0), (Qs[Q_names.q_1], 1.0), (Qs[Q_names.x], 1.0), (Qs[Q_names.y], 1.0), (Qs[Q_names.z], 1.0),
-                     (Qs[Q_names.Q], np.sqrt(Qs[Q_names.Q].w ** 2 + Qs[Q_names.Q].x ** 2 + Qs[Q_names.Q].y ** 2 + Qs[Q_names.Q].z ** 2)),
-                     (Qs[Q_names.Qbar], np.sqrt(Qs[Q_names.Q].w ** 2 + Qs[Q_names.Q].x ** 2 + Qs[Q_names.Q].y ** 2 + Qs[Q_names.Q].z ** 2))]:
-            assert np.allclose(abs_func(q), a)
-
-
-def test_quaternion_norm(Qs, Q_names, Q_conditions, on_windows, array):
-    Qs = array(Qs.ndarray)
-    for q in Qs[Q_conditions.nan]:
-        assert np.isnan(q.norm)
-    for q in Qs[Q_conditions.inf]:
-        if on_windows:  # pragma: no cover
-            assert np.isinf(q.norm) or np.isnan(q.norm)
-        else:
-            assert np.isinf(q.norm)
-    for q, a in [(Qs[Q_names.q_0], 0.0), (Qs[Q_names.q_1], 1.0), (Qs[Q_names.x], 1.0), (Qs[Q_names.y], 1.0), (Qs[Q_names.z], 1.0),
-                 (Qs[Q_names.Q], Qs[Q_names.Q].w ** 2 + Qs[Q_names.Q].x ** 2 + Qs[Q_names.Q].y ** 2 + Qs[Q_names.Q].z ** 2),
-                 (Qs[Q_names.Qbar], Qs[Q_names.Q].w ** 2 + Qs[Q_names.Q].x ** 2 + Qs[Q_names.Q].y ** 2 + Qs[Q_names.Q].z ** 2)]:
-        assert np.allclose(q.norm, a)
-
-
-# Unary quaternion returners
-def test_quaternion_negative(Qs, Q_names, Q_conditions, array):
-    Qs = array(Qs.ndarray)
-    assert -Qs[Q_names.Q] == Qs[Q_names.Qneg]
-    for q in Qs[Q_conditions.finite]:
-        assert -q == -1.0 * q
-    for q in Qs[Q_conditions.nonnan]:
-        assert -(-q) == q
-
-
-def test_quaternion_npconjugate(Qs, Q_names, Q_conditions, array):
-    Qs = array(Qs.ndarray)
-    assert np.conjugate(Qs[Q_names.Q]) == Qs[Q_names.Qbar]
-    for q in Qs[Q_conditions.nonnan]:
-        assert np.conjugate(q) == np.conj(q)
-        assert np.conjugate(np.conjugate(q)) == q
-        c = np.conjugate(q)
-        assert c.w == q.w
-        assert c.x == -q.x
-        assert c.y == -q.y
-        assert c.z == -q.z
-
-
-def test_quaternion_sqrt(Qs, Q_names, Q_conditions, array):
-    Qs = array(Qs.ndarray)
-    sqrt_precision = 2.e-15
-    # Test sqrt of basis elements
-    sqrthalf = np.sqrt(0.5)
-    assert np.array_equal(
-        np.sqrt(array(np.eye(4))).ndarray,
-        np.array([
-            [1, 0, 0, 0],
-            [sqrthalf, sqrthalf, 0, 0],
-            [sqrthalf, 0, sqrthalf, 0],
-            [sqrthalf, 0, 0, sqrthalf],
-        ])
-    )
-    # Test all my samples
-    one, i, j, k = tuple(array(np.eye(4)))
-    for q in Qs[Q_conditions.finitenonzero]:
-        assert np.allclose(np.sqrt(q) * np.sqrt(q), q, rtol=sqrt_precision)
-        # Ensure that non-unit quaternions are handled correctly
-        for s in [1, -1, 2, -2, 3.4, -3.4]:
-            for r in [one, i, j, k]:
-                srq = s*r*q
-                assert np.allclose(np.sqrt(srq) * np.sqrt(srq), srq, rtol=sqrt_precision)
-    # Test a huge batch of random quaternions
-    np.random.seed(1234)
-    a = array(np.random.uniform(-10, 10, size=10000*4).reshape((-1, 4)))
-    assert np.allclose(a, np.square(np.sqrt(a)), rtol=10*sqrt_precision, atol=0)
-    # Test some edge cases
-    _quaternion_resolution = 10 * np.finfo(float).resolution
-    for s in [-0.1, -1, -2, -3.4]:
-        q = np.sqrt(one * s)
-        assert np.array_equal(q, np.sqrt(-s)*i)
-        q = np.sqrt(one * s + array(0, *(0.05*np.sqrt(-s*_quaternion_resolution)*np.random.rand(3))))
-        assert np.array_equal(q, np.sqrt(-s)*i)
-
-
-def test_quaternion_square(Qs, Q_names, Q_conditions, array):
-    Qs = array(Qs.ndarray)
-    square_precision = 1.e-15
-    for q in Qs[Q_conditions.finite]:
-        assert (q*q - q**2).norm < square_precision
-        a = array([q])
-        assert (a**2 - array([q**2])).norm < square_precision
-
-
-def test_quaternion_log_exp(Qs, Q_names, Q_conditions, array):
-    Qs = array(Qs.ndarray)
-    qlogexp_precision = 4.e-15
-    zero = array([0, 0, 0, 0])
-    one, i, j, k = tuple(array(np.eye(4)))
-    assert (np.log(Qs[Q_names.Q]) - Qs[Q_names.Qlog]).abs < qlogexp_precision
-    assert (np.exp(Qs[Q_names.Q]) - Qs[Q_names.Qexp]).abs < qlogexp_precision
-    assert (np.exp(np.log(Qs[Q_names.Q])) - Qs[Q_names.Q]).abs < qlogexp_precision
-    assert (np.log(np.exp(Qs[Q_names.Q])) - Qs[Q_names.Q]).abs > qlogexp_precision  # Note order of operations!
-    assert np.log(one) == zero
-    assert np.log(i) == (np.pi / 2) * i
-    assert np.log(j) == (np.pi / 2) * j
-    assert np.log(k) == (np.pi / 2) * k
-    assert np.log(-one) == (np.pi) * i
-
-
-# Binary quat returners
-def test_quaternion_conjugate(Qs, Q_names, Q_conditions, array):
-    Qs = array(Qs.ndarray)
-    for conj_func in [np.conj, np.conjugate, lambda q: q.conj(), lambda q: q.conjugate()]:
-        assert conj_func(Qs[Q_names.Q]) == Qs[Q_names.Qbar]
-        for q in Qs[Q_conditions.nonnan]:
-            assert conj_func(conj_func(q)) == q
-            c = conj_func(q)
-            assert c.w == q.w
-            assert c.x == -q.x
-            assert c.y == -q.y
-            assert c.z == -q.z
-    for q in Qs[Q_conditions.nonnan]:
-        assert q.conjugate() == q.conj()
-        assert np.conjugate(q) == np.conj(q)
-
-
-def test_xor(array):
-    basis = one, i, j, k = tuple(array(np.eye(4)))
-    zero = 0 * one
-    assert one ^ one == one
-    assert one ^ i == i
-    assert one ^ j == j
-    assert one ^ k == k
-    assert i ^ one == i
-    assert i ^ i == zero
-    assert i ^ j == zero
-    assert i ^ k == zero
-    assert j ^ one == j
-    assert j ^ i == zero
-    assert j ^ j == zero
-    assert j ^ k == zero
-    assert k ^ one == k
-    assert k ^ i == zero
-    assert k ^ j == zero
-    assert k ^ k == zero
-
-
-def test_contractions(array):
-    basis = one, i, j, k = tuple(array(np.eye(4)))
-    zero = 0 * one
-
-    assert one << one == one
-    assert one << i == i
-    assert one << j == j
-    assert one << k == k
-    assert i << one == zero
-    assert i << i == -one
-    assert i << j == zero
-    assert i << k == zero
-    assert j << one == zero
-    assert j << i == zero
-    assert j << j == -one
-    assert j << k == zero
-    assert k << one == zero
-    assert k << i == zero
-    assert k << j == zero
-    assert k << k == -one
-
-    assert one >> one == one
-    assert one >> i == zero
-    assert one >> j == zero
-    assert one >> k == zero
-    assert i >> one == i
-    assert i >> i == -one
-    assert i >> j == zero
-    assert i >> k == zero
-    assert j >> one == j
-    assert j >> i == zero
-    assert j >> j == -one
-    assert j >> k == zero
-    assert k >> one == k
-    assert k >> i == zero
-    assert k >> j == zero
-    assert k >> k == -one
-
-    for a in basis:
-        for b in basis:
-            for c in basis:
-                assert np.allclose((a ^ b) | c, a | (b << c))
-                assert np.allclose(c | (b ^ a), (c >> b) | a)

From c79bf7b20aed0f8fdbad062a67586b949cde2edf Mon Sep 17 00:00:00 2001
From: Michael Boyle <michael.oliver.boyle@gmail.com>
Date: Wed, 7 Apr 2021 12:39:05 -0400
Subject: [PATCH 02/10] Warn about copying

---
 README.md | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/README.md b/README.md
index 50782c0..b2c3fef 100644
--- a/README.md
+++ b/README.md
@@ -88,6 +88,14 @@ Similarly, if you multiply two quaternionic arrays, their product will be comput
 quaternion multiplication, rather than element-wise multiplication of floats as numpy usually
 performs.
 
+| :warning: WARNING                                                                                |
+|:-------------------------------------------------------------------------------------------------|
+| Because of unfortunate default values, the `np.copy` function will not preserve the quaternionic |
+| nature of an array by default; the result will just be a plain array of floats.  You can either  |
+| pass the optional argument `subok=True`, as in `q3 = np.copy(q1, subok=True)`, or you can just   |
+| use the member function: `q3 = q1.copy()`.                                                       |
+
+
 ## Algebra
 
 All the usual quaternion operations are available, including

From ee71faf7c119a45ad68984baa69a8c0b40e510cf Mon Sep 17 00:00:00 2001
From: Michael Boyle <michael.oliver.boyle@gmail.com>
Date: Wed, 7 Apr 2021 12:40:57 -0400
Subject: [PATCH 03/10] Test for failing copies

---
 tests/test_copy.py | 43 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 43 insertions(+)
 create mode 100644 tests/test_copy.py

diff --git a/tests/test_copy.py b/tests/test_copy.py
new file mode 100644
index 0000000..6201120
--- /dev/null
+++ b/tests/test_copy.py
@@ -0,0 +1,43 @@
+import warnings
+import copy
+import pickle
+import numpy as np
+import quaternionic
+import pytest
+
+
+def np_copy(q):
+    return np.copy(q)
+
+def np_array_copy(q):
+    return np.array(q, copy=True)
+
+def np_array_copy_subok(q):
+    return np.array(q, copy=True, subok=True)
+
+def ndarray_copy(q):
+    return q.copy()
+
+def pickle_roundtrip(q):
+    return pickle.loads(pickle.dumps(q))
+
+def copy_copy(q):
+    return copy.copy(q)
+
+def copy_deepcopy(q):
+    return copy.deepcopy(q)
+
+# Note that np.copy and np.array(..., copy=True) return ndarray's, and thus lose information
+copy_xfail = lambda f: pytest.param(f, marks=pytest.mark.xfail(reason="Unexpected numpy defaults"))
+local_xfail = lambda f: pytest.param(f, marks=pytest.mark.xfail(reason="Can't pickle local object"))
+
+@pytest.mark.parametrize("copier", [
+    copy_xfail(np_copy), copy_xfail(np_array_copy), np_array_copy_subok,
+    ndarray_copy, local_xfail(pickle_roundtrip), copy_copy, copy_deepcopy
+])
+def test_modes_copying_and_pickling(copier):
+    q = quaternionic.array(np.random.normal(size=(17, 4)))
+    c = copier(q)
+    assert q is not c
+    assert isinstance(c, type(q))
+    assert np.array_equal(c, q)

From bd5b71992f6f2b4a33d164bd00414126bad065ba Mon Sep 17 00:00:00 2001
From: Michael Boyle <michael.oliver.boyle@gmail.com>
Date: Wed, 7 Apr 2021 14:11:02 -0400
Subject: [PATCH 04/10] Update poetry

---
 poetry.lock | 113 ++++++++++++++++++++++++++--------------------------
 1 file changed, 57 insertions(+), 56 deletions(-)

diff --git a/poetry.lock b/poetry.lock
index e0eae75..9fae2c9 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -123,7 +123,7 @@ python-versions = ">=3.6, <3.7"
 
 [[package]]
 name = "docutils"
-version = "0.16"
+version = "0.17"
 description = "Docutils -- Python Documentation Utilities"
 category = "main"
 optional = true
@@ -300,11 +300,11 @@ python-versions = "*"
 
 [[package]]
 name = "nltk"
-version = "3.5"
+version = "3.6"
 description = "Natural Language Toolkit"
 category = "main"
 optional = true
-python-versions = "*"
+python-versions = ">=3.5.*"
 
 [package.dependencies]
 click = "*"
@@ -313,9 +313,9 @@ regex = "*"
 tqdm = "*"
 
 [package.extras]
-all = ["requests", "numpy", "python-crfsuite", "scikit-learn", "twython", "pyparsing", "scipy", "matplotlib", "gensim"]
+all = ["numpy", "matplotlib", "pyparsing", "twython", "scikit-learn", "gensim (<4.0.0)", "requests", "python-crfsuite", "scipy"]
 corenlp = ["requests"]
-machine_learning = ["gensim", "numpy", "python-crfsuite", "scikit-learn", "scipy"]
+machine_learning = ["gensim (<4.0.0)", "numpy", "python-crfsuite", "scikit-learn", "scipy"]
 plot = ["matplotlib"]
 tgrep = ["pyparsing"]
 twitter = ["twython"]
@@ -425,7 +425,7 @@ python-versions = ">=2.6, !=3.0.*, !=3.1.*, !=3.2.*"
 
 [[package]]
 name = "pytest"
-version = "6.2.2"
+version = "6.2.3"
 description = "pytest: simple powerful testing with Python"
 category = "dev"
 optional = false
@@ -478,7 +478,7 @@ python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*, !=3.5.*"
 
 [[package]]
 name = "regex"
-version = "2021.3.17"
+version = "2021.4.4"
 description = "Alternative regular expression module, to replace re."
 category = "main"
 optional = true
@@ -649,7 +649,7 @@ python-versions = ">= 3.5"
 
 [[package]]
 name = "tqdm"
-version = "4.59.0"
+version = "4.60.0"
 description = "Fast, Extensible Progress Meter"
 category = "main"
 optional = true
@@ -810,8 +810,8 @@ dataclasses = [
     {file = "dataclasses-0.8.tar.gz", hash = "sha256:8479067f342acf957dc82ec415d355ab5edb7e7646b90dc6e2fd1d96ad084c97"},
 ]
 docutils = [
-    {file = "docutils-0.16-py2.py3-none-any.whl", hash = "sha256:0c5b78adfbf7762415433f5515cd5c9e762339e23369dbe8000d84a4bf4ab3af"},
-    {file = "docutils-0.16.tar.gz", hash = "sha256:c2de3a60e9e7d07be26b7f2b00ca0309c207e06c100f9cc2a94931fc75a478fc"},
+    {file = "docutils-0.17-py2.py3-none-any.whl", hash = "sha256:a71042bb7207c03d5647f280427f14bfbd1a65c9eb84f4b341d85fafb6bb4bdf"},
+    {file = "docutils-0.17.tar.gz", hash = "sha256:e2ffeea817964356ba4470efba7c2f42b6b0de0b04e66378507e3e2504bbff4c"},
 ]
 future = [
     {file = "future-0.18.2.tar.gz", hash = "sha256:b1bead90b70cf6ec3f0710ae53a525360fa360d306a86583adc6bf83a4db537d"},
@@ -921,7 +921,8 @@ mypy-extensions = [
     {file = "mypy_extensions-0.4.3.tar.gz", hash = "sha256:2d82818f5bb3e369420cb3c4060a7970edba416647068eb4c5343488a6c604a8"},
 ]
 nltk = [
-    {file = "nltk-3.5.zip", hash = "sha256:845365449cd8c5f9731f7cb9f8bd6fd0767553b9d53af9eb1b3abf7700936b35"},
+    {file = "nltk-3.6-py3-none-any.whl", hash = "sha256:718de6908f538db19a77f96b9e6f5f586b0892d7de5eea32e71f2a2535ed8657"},
+    {file = "nltk-3.6.zip", hash = "sha256:5c5c0c02942cb423d38e5d66285257f53f05e420a8715bca273fb9cf28067e0e"},
 ]
 numba = [
     {file = "numba-0.53.1-cp36-cp36m-macosx_10_14_x86_64.whl", hash = "sha256:b23de6b6837c132087d06b8b92d343edb54b885873b824a037967fbd5272ebb7"},
@@ -1015,8 +1016,8 @@ pyparsing = [
     {file = "pyparsing-2.4.7.tar.gz", hash = "sha256:c203ec8783bf771a155b207279b9bccb8dea02d8f0c9e5f8ead507bc3246ecc1"},
 ]
 pytest = [
-    {file = "pytest-6.2.2-py3-none-any.whl", hash = "sha256:b574b57423e818210672e07ca1fa90aaf194a4f63f3ab909a2c67ebb22913839"},
-    {file = "pytest-6.2.2.tar.gz", hash = "sha256:9d1edf9e7d0b84d72ea3dbcdfd22b35fb543a5e8f2a60092dd578936bf63d7f9"},
+    {file = "pytest-6.2.3-py3-none-any.whl", hash = "sha256:6ad9c7bdf517a808242b998ac20063c41532a570d088d77eec1ee12b0b5574bc"},
+    {file = "pytest-6.2.3.tar.gz", hash = "sha256:671238a46e4df0f3498d1c3270e5deb9b32d25134c99b7d75370a68cfbe9b634"},
 ]
 pytest-cov = [
     {file = "pytest-cov-2.11.1.tar.gz", hash = "sha256:359952d9d39b9f822d9d29324483e7ba04a3a17dd7d05aa6beb7ea01e359e5f7"},
@@ -1050,47 +1051,47 @@ pyyaml = [
     {file = "PyYAML-5.4.1.tar.gz", hash = "sha256:607774cbba28732bfa802b54baa7484215f530991055bb562efbed5b2f20a45e"},
 ]
 regex = [
-    {file = "regex-2021.3.17-cp36-cp36m-macosx_10_9_x86_64.whl", hash = "sha256:b97ec5d299c10d96617cc851b2e0f81ba5d9d6248413cd374ef7f3a8871ee4a6"},
-    {file = "regex-2021.3.17-cp36-cp36m-manylinux1_i686.whl", hash = "sha256:cb4ee827857a5ad9b8ae34d3c8cc51151cb4a3fe082c12ec20ec73e63cc7c6f0"},
-    {file = "regex-2021.3.17-cp36-cp36m-manylinux1_x86_64.whl", hash = "sha256:633497504e2a485a70a3268d4fc403fe3063a50a50eed1039083e9471ad0101c"},
-    {file = "regex-2021.3.17-cp36-cp36m-manylinux2010_i686.whl", hash = "sha256:a59a2ee329b3de764b21495d78c92ab00b4ea79acef0f7ae8c1067f773570afa"},
-    {file = "regex-2021.3.17-cp36-cp36m-manylinux2010_x86_64.whl", hash = "sha256:f85d6f41e34f6a2d1607e312820971872944f1661a73d33e1e82d35ea3305e14"},
-    {file = "regex-2021.3.17-cp36-cp36m-manylinux2014_aarch64.whl", hash = "sha256:4651f839dbde0816798e698626af6a2469eee6d9964824bb5386091255a1694f"},
-    {file = "regex-2021.3.17-cp36-cp36m-manylinux2014_i686.whl", hash = "sha256:39c44532d0e4f1639a89e52355b949573e1e2c5116106a395642cbbae0ff9bcd"},
-    {file = "regex-2021.3.17-cp36-cp36m-manylinux2014_x86_64.whl", hash = "sha256:3d9a7e215e02bd7646a91fb8bcba30bc55fd42a719d6b35cf80e5bae31d9134e"},
-    {file = "regex-2021.3.17-cp36-cp36m-win32.whl", hash = "sha256:159fac1a4731409c830d32913f13f68346d6b8e39650ed5d704a9ce2f9ef9cb3"},
-    {file = "regex-2021.3.17-cp36-cp36m-win_amd64.whl", hash = "sha256:13f50969028e81765ed2a1c5fcfdc246c245cf8d47986d5172e82ab1a0c42ee5"},
-    {file = "regex-2021.3.17-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:b9d8d286c53fe0cbc6d20bf3d583cabcd1499d89034524e3b94c93a5ab85ca90"},
-    {file = "regex-2021.3.17-cp37-cp37m-manylinux1_i686.whl", hash = "sha256:201e2619a77b21a7780580ab7b5ce43835e242d3e20fef50f66a8df0542e437f"},
-    {file = "regex-2021.3.17-cp37-cp37m-manylinux1_x86_64.whl", hash = "sha256:d47d359545b0ccad29d572ecd52c9da945de7cd6cf9c0cfcb0269f76d3555689"},
-    {file = "regex-2021.3.17-cp37-cp37m-manylinux2010_i686.whl", hash = "sha256:ea2f41445852c660ba7c3ebf7d70b3779b20d9ca8ba54485a17740db49f46932"},
-    {file = "regex-2021.3.17-cp37-cp37m-manylinux2010_x86_64.whl", hash = "sha256:486a5f8e11e1f5bbfcad87f7c7745eb14796642323e7e1829a331f87a713daaa"},
-    {file = "regex-2021.3.17-cp37-cp37m-manylinux2014_aarch64.whl", hash = "sha256:18e25e0afe1cf0f62781a150c1454b2113785401ba285c745acf10c8ca8917df"},
-    {file = "regex-2021.3.17-cp37-cp37m-manylinux2014_i686.whl", hash = "sha256:a2ee026f4156789df8644d23ef423e6194fad0bc53575534101bb1de5d67e8ce"},
-    {file = "regex-2021.3.17-cp37-cp37m-manylinux2014_x86_64.whl", hash = "sha256:4c0788010a93ace8a174d73e7c6c9d3e6e3b7ad99a453c8ee8c975ddd9965643"},
-    {file = "regex-2021.3.17-cp37-cp37m-win32.whl", hash = "sha256:575a832e09d237ae5fedb825a7a5bc6a116090dd57d6417d4f3b75121c73e3be"},
-    {file = "regex-2021.3.17-cp37-cp37m-win_amd64.whl", hash = "sha256:8e65e3e4c6feadf6770e2ad89ad3deb524bcb03d8dc679f381d0568c024e0deb"},
-    {file = "regex-2021.3.17-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:a0df9a0ad2aad49ea3c7f65edd2ffb3d5c59589b85992a6006354f6fb109bb18"},
-    {file = "regex-2021.3.17-cp38-cp38-manylinux1_i686.whl", hash = "sha256:b98bc9db003f1079caf07b610377ed1ac2e2c11acc2bea4892e28cc5b509d8d5"},
-    {file = "regex-2021.3.17-cp38-cp38-manylinux1_x86_64.whl", hash = "sha256:808404898e9a765e4058bf3d7607d0629000e0a14a6782ccbb089296b76fa8fe"},
-    {file = "regex-2021.3.17-cp38-cp38-manylinux2010_i686.whl", hash = "sha256:5770a51180d85ea468234bc7987f5597803a4c3d7463e7323322fe4a1b181578"},
-    {file = "regex-2021.3.17-cp38-cp38-manylinux2010_x86_64.whl", hash = "sha256:976a54d44fd043d958a69b18705a910a8376196c6b6ee5f2596ffc11bff4420d"},
-    {file = "regex-2021.3.17-cp38-cp38-manylinux2014_aarch64.whl", hash = "sha256:63f3ca8451e5ff7133ffbec9eda641aeab2001be1a01878990f6c87e3c44b9d5"},
-    {file = "regex-2021.3.17-cp38-cp38-manylinux2014_i686.whl", hash = "sha256:bcd945175c29a672f13fce13a11893556cd440e37c1b643d6eeab1988c8b209c"},
-    {file = "regex-2021.3.17-cp38-cp38-manylinux2014_x86_64.whl", hash = "sha256:3d9356add82cff75413bec360c1eca3e58db4a9f5dafa1f19650958a81e3249d"},
-    {file = "regex-2021.3.17-cp38-cp38-win32.whl", hash = "sha256:f5d0c921c99297354cecc5a416ee4280bd3f20fd81b9fb671ca6be71499c3fdf"},
-    {file = "regex-2021.3.17-cp38-cp38-win_amd64.whl", hash = "sha256:14de88eda0976020528efc92d0a1f8830e2fb0de2ae6005a6fc4e062553031fa"},
-    {file = "regex-2021.3.17-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:4c2e364491406b7888c2ad4428245fc56c327e34a5dfe58fd40df272b3c3dab3"},
-    {file = "regex-2021.3.17-cp39-cp39-manylinux1_i686.whl", hash = "sha256:8bd4f91f3fb1c9b1380d6894bd5b4a519409135bec14c0c80151e58394a4e88a"},
-    {file = "regex-2021.3.17-cp39-cp39-manylinux1_x86_64.whl", hash = "sha256:882f53afe31ef0425b405a3f601c0009b44206ea7f55ee1c606aad3cc213a52c"},
-    {file = "regex-2021.3.17-cp39-cp39-manylinux2010_i686.whl", hash = "sha256:07ef35301b4484bce843831e7039a84e19d8d33b3f8b2f9aab86c376813d0139"},
-    {file = "regex-2021.3.17-cp39-cp39-manylinux2010_x86_64.whl", hash = "sha256:360a01b5fa2ad35b3113ae0c07fb544ad180603fa3b1f074f52d98c1096fa15e"},
-    {file = "regex-2021.3.17-cp39-cp39-manylinux2014_aarch64.whl", hash = "sha256:709f65bb2fa9825f09892617d01246002097f8f9b6dde8d1bb4083cf554701ba"},
-    {file = "regex-2021.3.17-cp39-cp39-manylinux2014_i686.whl", hash = "sha256:c66221e947d7207457f8b6f42b12f613b09efa9669f65a587a2a71f6a0e4d106"},
-    {file = "regex-2021.3.17-cp39-cp39-manylinux2014_x86_64.whl", hash = "sha256:c782da0e45aff131f0bed6e66fbcfa589ff2862fc719b83a88640daa01a5aff7"},
-    {file = "regex-2021.3.17-cp39-cp39-win32.whl", hash = "sha256:dc9963aacb7da5177e40874585d7407c0f93fb9d7518ec58b86e562f633f36cd"},
-    {file = "regex-2021.3.17-cp39-cp39-win_amd64.whl", hash = "sha256:a0d04128e005142260de3733591ddf476e4902c0c23c1af237d9acf3c96e1b38"},
-    {file = "regex-2021.3.17.tar.gz", hash = "sha256:4b8a1fb724904139149a43e172850f35aa6ea97fb0545244dc0b805e0154ed68"},
+    {file = "regex-2021.4.4-cp36-cp36m-macosx_10_9_x86_64.whl", hash = "sha256:619d71c59a78b84d7f18891fe914446d07edd48dc8328c8e149cbe0929b4e000"},
+    {file = "regex-2021.4.4-cp36-cp36m-manylinux1_i686.whl", hash = "sha256:47bf5bf60cf04d72bf6055ae5927a0bd9016096bf3d742fa50d9bf9f45aa0711"},
+    {file = "regex-2021.4.4-cp36-cp36m-manylinux1_x86_64.whl", hash = "sha256:281d2fd05555079448537fe108d79eb031b403dac622621c78944c235f3fcf11"},
+    {file = "regex-2021.4.4-cp36-cp36m-manylinux2010_i686.whl", hash = "sha256:bd28bc2e3a772acbb07787c6308e00d9626ff89e3bfcdebe87fa5afbfdedf968"},
+    {file = "regex-2021.4.4-cp36-cp36m-manylinux2010_x86_64.whl", hash = "sha256:7c2a1af393fcc09e898beba5dd59196edaa3116191cc7257f9224beaed3e1aa0"},
+    {file = "regex-2021.4.4-cp36-cp36m-manylinux2014_aarch64.whl", hash = "sha256:c38c71df845e2aabb7fb0b920d11a1b5ac8526005e533a8920aea97efb8ec6a4"},
+    {file = "regex-2021.4.4-cp36-cp36m-manylinux2014_i686.whl", hash = "sha256:96fcd1888ab4d03adfc9303a7b3c0bd78c5412b2bfbe76db5b56d9eae004907a"},
+    {file = "regex-2021.4.4-cp36-cp36m-manylinux2014_x86_64.whl", hash = "sha256:ade17eb5d643b7fead300a1641e9f45401c98eee23763e9ed66a43f92f20b4a7"},
+    {file = "regex-2021.4.4-cp36-cp36m-win32.whl", hash = "sha256:e8e5b509d5c2ff12f8418006d5a90e9436766133b564db0abaec92fd27fcee29"},
+    {file = "regex-2021.4.4-cp36-cp36m-win_amd64.whl", hash = "sha256:11d773d75fa650cd36f68d7ca936e3c7afaae41b863b8c387a22aaa78d3c5c79"},
+    {file = "regex-2021.4.4-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:d3029c340cfbb3ac0a71798100ccc13b97dddf373a4ae56b6a72cf70dfd53bc8"},
+    {file = "regex-2021.4.4-cp37-cp37m-manylinux1_i686.whl", hash = "sha256:18c071c3eb09c30a264879f0d310d37fe5d3a3111662438889ae2eb6fc570c31"},
+    {file = "regex-2021.4.4-cp37-cp37m-manylinux1_x86_64.whl", hash = "sha256:4c557a7b470908b1712fe27fb1ef20772b78079808c87d20a90d051660b1d69a"},
+    {file = "regex-2021.4.4-cp37-cp37m-manylinux2010_i686.whl", hash = "sha256:01afaf2ec48e196ba91b37451aa353cb7eda77efe518e481707e0515025f0cd5"},
+    {file = "regex-2021.4.4-cp37-cp37m-manylinux2010_x86_64.whl", hash = "sha256:3a9cd17e6e5c7eb328517969e0cb0c3d31fd329298dd0c04af99ebf42e904f82"},
+    {file = "regex-2021.4.4-cp37-cp37m-manylinux2014_aarch64.whl", hash = "sha256:90f11ff637fe8798933fb29f5ae1148c978cccb0452005bf4c69e13db951e765"},
+    {file = "regex-2021.4.4-cp37-cp37m-manylinux2014_i686.whl", hash = "sha256:919859aa909429fb5aa9cf8807f6045592c85ef56fdd30a9a3747e513db2536e"},
+    {file = "regex-2021.4.4-cp37-cp37m-manylinux2014_x86_64.whl", hash = "sha256:339456e7d8c06dd36a22e451d58ef72cef293112b559010db3d054d5560ef439"},
+    {file = "regex-2021.4.4-cp37-cp37m-win32.whl", hash = "sha256:67bdb9702427ceddc6ef3dc382455e90f785af4c13d495f9626861763ee13f9d"},
+    {file = "regex-2021.4.4-cp37-cp37m-win_amd64.whl", hash = "sha256:32e65442138b7b76dd8173ffa2cf67356b7bc1768851dded39a7a13bf9223da3"},
+    {file = "regex-2021.4.4-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:1e1c20e29358165242928c2de1482fb2cf4ea54a6a6dea2bd7a0e0d8ee321500"},
+    {file = "regex-2021.4.4-cp38-cp38-manylinux1_i686.whl", hash = "sha256:314d66636c494ed9c148a42731b3834496cc9a2c4251b1661e40936814542b14"},
+    {file = "regex-2021.4.4-cp38-cp38-manylinux1_x86_64.whl", hash = "sha256:6d1b01031dedf2503631d0903cb563743f397ccaf6607a5e3b19a3d76fc10480"},
+    {file = "regex-2021.4.4-cp38-cp38-manylinux2010_i686.whl", hash = "sha256:741a9647fcf2e45f3a1cf0e24f5e17febf3efe8d4ba1281dcc3aa0459ef424dc"},
+    {file = "regex-2021.4.4-cp38-cp38-manylinux2010_x86_64.whl", hash = "sha256:4c46e22a0933dd783467cf32b3516299fb98cfebd895817d685130cc50cd1093"},
+    {file = "regex-2021.4.4-cp38-cp38-manylinux2014_aarch64.whl", hash = "sha256:e512d8ef5ad7b898cdb2d8ee1cb09a8339e4f8be706d27eaa180c2f177248a10"},
+    {file = "regex-2021.4.4-cp38-cp38-manylinux2014_i686.whl", hash = "sha256:980d7be47c84979d9136328d882f67ec5e50008681d94ecc8afa8a65ed1f4a6f"},
+    {file = "regex-2021.4.4-cp38-cp38-manylinux2014_x86_64.whl", hash = "sha256:ce15b6d103daff8e9fee13cf7f0add05245a05d866e73926c358e871221eae87"},
+    {file = "regex-2021.4.4-cp38-cp38-win32.whl", hash = "sha256:a91aa8619b23b79bcbeb37abe286f2f408d2f2d6f29a17237afda55bb54e7aac"},
+    {file = "regex-2021.4.4-cp38-cp38-win_amd64.whl", hash = "sha256:c0502c0fadef0d23b128605d69b58edb2c681c25d44574fc673b0e52dce71ee2"},
+    {file = "regex-2021.4.4-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:598585c9f0af8374c28edd609eb291b5726d7cbce16be6a8b95aa074d252ee17"},
+    {file = "regex-2021.4.4-cp39-cp39-manylinux1_i686.whl", hash = "sha256:ee54ff27bf0afaf4c3b3a62bcd016c12c3fdb4ec4f413391a90bd38bc3624605"},
+    {file = "regex-2021.4.4-cp39-cp39-manylinux1_x86_64.whl", hash = "sha256:7d9884d86dd4dd489e981d94a65cd30d6f07203d90e98f6f657f05170f6324c9"},
+    {file = "regex-2021.4.4-cp39-cp39-manylinux2010_i686.whl", hash = "sha256:bf5824bfac591ddb2c1f0a5f4ab72da28994548c708d2191e3b87dd207eb3ad7"},
+    {file = "regex-2021.4.4-cp39-cp39-manylinux2010_x86_64.whl", hash = "sha256:563085e55b0d4fb8f746f6a335893bda5c2cef43b2f0258fe1020ab1dd874df8"},
+    {file = "regex-2021.4.4-cp39-cp39-manylinux2014_aarch64.whl", hash = "sha256:b9c3db21af35e3b3c05764461b262d6f05bbca08a71a7849fd79d47ba7bc33ed"},
+    {file = "regex-2021.4.4-cp39-cp39-manylinux2014_i686.whl", hash = "sha256:3916d08be28a1149fb97f7728fca1f7c15d309a9f9682d89d79db75d5e52091c"},
+    {file = "regex-2021.4.4-cp39-cp39-manylinux2014_x86_64.whl", hash = "sha256:fd45ff9293d9274c5008a2054ecef86a9bfe819a67c7be1afb65e69b405b3042"},
+    {file = "regex-2021.4.4-cp39-cp39-win32.whl", hash = "sha256:fa4537fb4a98fe8fde99626e4681cc644bdcf2a795038533f9f711513a862ae6"},
+    {file = "regex-2021.4.4-cp39-cp39-win_amd64.whl", hash = "sha256:97f29f57d5b84e73fbaf99ab3e26134e6687348e95ef6b48cfd2c06807005a07"},
+    {file = "regex-2021.4.4.tar.gz", hash = "sha256:52ba3d3f9b942c49d7e4bc105bb28551c44065f139a65062ab7912bef10c9afb"},
 ]
 requests = [
     {file = "requests-2.25.1-py2.py3-none-any.whl", hash = "sha256:c210084e36a42ae6b9219e00e48287def368a26d03a048ddad7bfee44f75871e"},
@@ -1207,8 +1208,8 @@ tornado = [
     {file = "tornado-6.1.tar.gz", hash = "sha256:33c6e81d7bd55b468d2e793517c909b139960b6c790a60b7991b9b6b76fb9791"},
 ]
 tqdm = [
-    {file = "tqdm-4.59.0-py2.py3-none-any.whl", hash = "sha256:9fdf349068d047d4cfbe24862c425883af1db29bcddf4b0eeb2524f6fbdb23c7"},
-    {file = "tqdm-4.59.0.tar.gz", hash = "sha256:d666ae29164da3e517fcf125e41d4fe96e5bb375cd87ff9763f6b38b5592fe33"},
+    {file = "tqdm-4.60.0-py2.py3-none-any.whl", hash = "sha256:daec693491c52e9498632dfbe9ccfc4882a557f5fa08982db1b4d3adbe0887c3"},
+    {file = "tqdm-4.60.0.tar.gz", hash = "sha256:ebdebdb95e3477ceea267decfc0784859aa3df3e27e22d23b83e9b272bf157ae"},
 ]
 typed-ast = [
     {file = "typed_ast-1.4.2-cp35-cp35m-manylinux1_i686.whl", hash = "sha256:7703620125e4fb79b64aa52427ec192822e9f45d37d4b6625ab37ef403e1df70"},

From 079d477a3c89ec8021a182690edd4b4eb2ed631d Mon Sep 17 00:00:00 2001
From: Michael Boyle <michael.oliver.boyle@gmail.com>
Date: Wed, 7 Apr 2021 14:11:35 -0400
Subject: [PATCH 05/10] Ignore direnv file

---
 .gitignore | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/.gitignore b/.gitignore
index 50113e5..05799c1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -115,3 +115,6 @@ notes
 
 build_artifacts
 # conda smithy ci-skeleton end
+
+# direnv file
+.envrc
\ No newline at end of file

From 763bda1db6c03e8e705f00baef9a6d2e18ba1318 Mon Sep 17 00:00:00 2001
From: Michael Boyle <michael.oliver.boyle@gmail.com>
Date: Wed, 7 Apr 2021 14:11:52 -0400
Subject: [PATCH 06/10] Favor np.random.normal

---
 README.md                   |  8 ++++----
 quaternionic/arrays.py      |  4 ++--
 quaternionic/converters.py  | 37 +++++++++++++++++++++++++++++++++++++
 tests/test_algebra.py       | 36 ++++++++++++++++++------------------
 tests/test_array.py         |  6 +++---
 tests/test_converters.py    | 34 +++++++++++++++++++++++-----------
 tests/test_interpolation.py |  4 +---
 tests/test_properties.py    | 19 +++++++++----------
 tests/test_utilities.py     | 22 +++++++++++-----------
 9 files changed, 108 insertions(+), 62 deletions(-)

diff --git a/README.md b/README.md
index b2c3fef..2d01dc5 100644
--- a/README.md
+++ b/README.md
@@ -70,7 +70,7 @@ existing arrays:
 import numpy as np
 import quaternionic
 
-a = 1.0 - np.random.rand(17, 11, 4)  # Just some random numbers; last dimension is 4
+a = np.random.normal(size=(17, 11, 4))  # Just some random numbers; last dimension is 4
 q1 = quaternionic.array(a)  # Reinterpret an existing array
 q2 = quaternionic.array([1.2, 2.3, 3.4, 4.5])  # Create a new array
 ```
@@ -91,9 +91,9 @@ performs.
 | :warning: WARNING                                                                                |
 |:-------------------------------------------------------------------------------------------------|
 | Because of unfortunate default values, the `np.copy` function will not preserve the quaternionic |
-| nature of an array by default; the result will just be a plain array of floats.  You can either  |
-| pass the optional argument `subok=True`, as in `q3 = np.copy(q1, subok=True)`, or you can just   |
-| use the member function: `q3 = q1.copy()`.                                                       |
+| nature of an array by default; the result will just be a plain array of floats.  You could pass  |
+| the optional argument `subok=True`, as in `q3 = np.copy(q1, subok=True)`, but it's easier to t   |
+| just use the member function: `q3 = q1.copy()`.                                                  |
 
 
 ## Algebra
diff --git a/quaternionic/arrays.py b/quaternionic/arrays.py
index 30c84c3..739f6d2 100644
--- a/quaternionic/arrays.py
+++ b/quaternionic/arrays.py
@@ -30,11 +30,11 @@ class QArray(QuaternionPropertiesMixin(jit), QuaternionConvertersMixin(jit), np.
         also be used as expected — as in `np.exp(q)`, `np.log(q)`, etc.
 
         Because this is a subclass of numpy's ndarray object, its constructor takes
-        anything the ndarray constructor takes, or just an ndarray to be considered
+        anything the ndarray constructor takes, or just an ndarray to be re-interpreted
         as a quaternion array:
 
             q1 = quaternionic.array([1, 2, 3, 4])  # explicit array
-            q2 = quaternionic.array(np.random.rand(10, 4))  # re-interpret ndarray
+            q2 = quaternionic.array(np.random.normal(size=(10, 4))  # re-interpret ndarray
 
         In addition to the basic numpy array features, we also have a number of
         extra properties that are particularly useful for quaternions, including
diff --git a/quaternionic/converters.py b/quaternionic/converters.py
index 2073ed1..8a3b4c6 100644
--- a/quaternionic/converters.py
+++ b/quaternionic/converters.py
@@ -841,4 +841,41 @@ def to_minimal_rotation(self, t, t_new=None, axis=0, iterations=2):
             Rγ = np.exp(z * γover2)
             return (R * Rγ).to_minimal_rotation(t_new, t_new=None, axis=axis, iterations=iterations-1)
 
+        @classmethod
+        def random(cls, shape=(4,), normalize=False):
+            """Construct random quaternions
+
+            Parameters
+            ----------
+            shape : tuple, optional
+                Shape of the array.  If this does not end with `4`, it will be appended.
+                Default is `(4,)`.
+            normalize : bool, optional
+                If True, normalize the result, so that the returned array can be
+                interpreted as rotors.  Defaults to False.
+
+            Returns
+            -------
+            q : array of quaternions
+
+            Notes
+            -----
+            This function constructs quaternions in which each component has a random value
+            drawn from a normal (Gaussian) distribution centered at 0 with scale 1.  This
+            has the nice property that the resulting distribution of quaternions is
+            isotropic — it is spherically symmetric.  If the result is normalized, these
+            are truly random rotors.
+
+            """
+            if isinstance(shape, int):
+                shape = (shape,)
+            if len(shape) == 0:
+                shape = (4,)
+            if shape[-1] != 4:
+                shape = shape + (4,)
+            q = np.random.normal(size=shape)  # Note the weird naming of this argument to `normal`
+            if normalize:
+                q /= np.linalg.norm(q, axis=-1)[..., np.newaxis]
+            return cls(q)
+
     return mixin
diff --git a/tests/test_algebra.py b/tests/test_algebra.py
index fa64c62..590551f 100644
--- a/tests/test_algebra.py
+++ b/tests/test_algebra.py
@@ -41,15 +41,15 @@ def test_basis_multiplication():
 def test_array_ufunc(array):
     np.random.seed(1234)
 
-    q = array(np.random.rand(1, 3, 4))
+    q = array(np.random.normal(size=(1, 3, 4)))
     with pytest.raises(NotImplementedError):
         np.exp(q, extra_arg=True)
 
     with pytest.raises(NotImplementedError):
         np.negative.at(q, [0, 1])
 
-    p = array(np.random.rand(17, 3, 4))
-    q = array(np.random.rand(1, 3, 4))
+    p = array(np.random.normal(size=(17, 3, 4)))
+    q = array(np.random.normal(size=(1, 3, 4)))
     pq1 = np.add(p, q)
     assert isinstance(pq1, array)
     assert pq1.shape == (17, 3, 4)
@@ -59,7 +59,7 @@ def test_array_ufunc(array):
     assert np.array_equal(pq1, pq2)
     assert isinstance(pq2, array)
 
-    p = array(np.random.rand(17, 3, 4))
+    p = array(np.random.normal(size=(17, 3, 4)))
     q = np.random.rand(1, 3)
     pq1 = np.multiply(p, q)
     assert isinstance(pq1, array)
@@ -70,7 +70,7 @@ def test_array_ufunc(array):
     assert isinstance(pq2, array)
 
     p = np.random.rand(1, 3)
-    q = array(np.random.rand(17, 3, 4))
+    q = array(np.random.normal(size=(17, 3, 4)))
     pq1 = np.multiply(p, q)
     assert isinstance(pq1, array)
     assert pq1.shape == (17, 3, 4)
@@ -80,7 +80,7 @@ def test_array_ufunc(array):
     assert isinstance(pq2, array)
 
     p = np.random.rand(1, 3)
-    q = np.random.rand(17, 3, 4)
+    q = np.random.normal(size=(17, 3, 4))
     s = np.random.rand(17, 3)
     pq1 = array(q).__array_ufunc__(np.multiply, "__call__", p, q)
     assert pq1 == NotImplemented
@@ -95,7 +95,7 @@ def test_array_ufunc(array):
     qfin = array(q).__array_ufunc__(np.isfinite, "__call__", q)
     assert qfin == NotImplemented
 
-    q = array(np.random.rand(17, 3, 4))
+    q = array(np.random.normal(size=(17, 3, 4)))
     qneg = np.negative(q)
     assert isinstance(qneg, array)
     assert qneg.shape == q.shape
@@ -114,7 +114,7 @@ def test_array_ufunc(array):
     assert isinstance(qneg2, array)
 
     p = np.random.rand(1, 3)
-    q = array(np.random.rand(17, 3, 4))
+    q = array(np.random.normal(size=(17, 3, 4)))
     qp1 = np.float_power(q, p)
     assert isinstance(qp1, array)
     assert qp1.shape == (17, 3, 4)
@@ -123,7 +123,7 @@ def test_array_ufunc(array):
     assert np.array_equal(qp1, qp2)
     assert isinstance(qp2, array)
 
-    q = array(np.random.rand(17, 3, 4))
+    q = array(np.random.normal(size=(17, 3, 4)))
     qabs = np.absolute(q)
     assert isinstance(qabs, np.ndarray) and not isinstance(qabs, array)
     assert qabs.shape == (17, 3)
@@ -131,14 +131,14 @@ def test_array_ufunc(array):
     np.absolute(q, out=qabs2)
     assert np.array_equal(qabs, qabs2)
     assert isinstance(qabs2, np.ndarray) and not isinstance(qabs, array)
-    q = array(np.random.rand(17, 3, 4, 4))
+    q = array(np.random.normal(size=(17, 3, 4, 4)))
     qabs = array(np.empty((17, 3, 4)))
     np.absolute(q, out=qabs)
     assert np.array_equal(qabs, np.absolute(q))
     assert isinstance(qabs2, np.ndarray) and isinstance(qabs, array)
 
-    p = array(np.random.rand(17, 3, 4))
-    q = array(np.random.rand(1, 3, 4))
+    p = array(np.random.normal(size=(17, 3, 4)))
+    q = array(np.random.normal(size=(1, 3, 4)))
     pq1 = np.equal(p, q)
     assert isinstance(pq1, np.ndarray) and not isinstance(pq1, array)
     assert pq1.shape == (17, 3)
@@ -148,13 +148,13 @@ def test_array_ufunc(array):
     assert np.array_equal(pq1, pq2)
     assert isinstance(pq2, np.ndarray) and not isinstance(pq2, array)
     assert pq2.shape == (17, 3)
-    p = array(np.random.rand(17, 3, 4, 4))
-    q = array(np.random.rand(17, 3, 4, 4))
+    p = array(np.random.normal(size=(17, 3, 4, 4)))
+    q = array(np.random.normal(size=(17, 3, 4, 4)))
     pq = array(np.empty((17, 3, 4)))
     np.equal(p, q, out=pq)
     assert isinstance(pq, np.ndarray) and isinstance(pq, array)
 
-    q = array(np.random.rand(17, 3, 4))
+    q = array(np.random.normal(size=(17, 3, 4)))
     qfin = np.isfinite(q)
     assert isinstance(qfin, np.ndarray) and not isinstance(qfin, array)
     assert qfin.shape == (17, 3)
@@ -164,7 +164,7 @@ def test_array_ufunc(array):
     assert np.array_equal(qfin, qfin2)
     assert isinstance(qfin2, np.ndarray) and not isinstance(qfin2, array)
     assert qfin2.shape == (17, 3)
-    q = array(np.random.rand(17, 3, 4, 4))
+    q = array(np.random.normal(size=(17, 3, 4, 4)))
     qfin = array(np.empty((17, 3, 4)))
     np.isfinite(q, out=qfin)
     assert isinstance(qfin, np.ndarray) and isinstance(qfin, array)
@@ -183,12 +183,12 @@ def test_array_ufunc(array):
     for k in dir(np):
         attr = getattr(np, k)
         if isinstance(attr, np.ufunc) and attr not in implemented_ufuncs:
-            p = array(np.random.rand(17, 3, 4))
+            p = array(np.random.normal(size=(17, 3, 4)))
             if attr.nin == 1:
                 with pytest.raises(TypeError):
                     attr(p)
             elif attr.nin == 2:
-                q = array(np.random.rand(17, 3, 4))
+                q = array(np.random.normal(size=(17, 3, 4)))
                 with pytest.raises(TypeError):
                     attr(p, q)
             else:  # pragma: no cover
diff --git a/tests/test_array.py b/tests/test_array.py
index 65b505a..f8c830e 100644
--- a/tests/test_array.py
+++ b/tests/test_array.py
@@ -27,7 +27,7 @@ def test_new(array):
 
 
 def test_getitem(array):
-    q = array(np.random.rand(17, 3, 4))
+    q = array(np.random.normal(size=(17, 3, 4)))
     p = q[1:-1]
     assert isinstance(p, array)
     assert p.shape == (q.shape[0]-2,) + q.shape[1:]
@@ -43,12 +43,12 @@ def test_array_finalize(array):
 
 
 def test_repr(array):
-    q = array(np.random.rand(17, 3, 4))
+    q = array(np.random.normal(size=(17, 3, 4)))
     assert repr(q) == 'quaternionic.' + repr(q.ndarray)
 
 
 def test_str(array):
-    q = array(np.random.rand(17, 3, 4))
+    q = array(np.random.normal(size=(17, 3, 4)))
     assert str(q) == str(q.ndarray)
 
 
diff --git a/tests/test_converters.py b/tests/test_converters.py
index 801060b..05443e0 100644
--- a/tests/test_converters.py
+++ b/tests/test_converters.py
@@ -458,14 +458,26 @@ def test_to_minimal_rotation():
     assert np.max(abs(q_minimal_rotation.to_angular_velocity(t)[:, :2])) < 1e-8
 
 
-# def test_slerp():
-#     slerp()
-
-
-# def test_slerp_pairwise():
-#     slerp_pairwise()
-
-
-# def test_squad():
-#     squad()
-
+def test_random():
+    q = quaternionic.array.random()
+    assert isinstance(q, quaternionic.array)
+    assert q.dtype == np.float64
+    assert q.shape == (4,)
+
+    q = quaternionic.array.random(17)
+    assert isinstance(q, quaternionic.array)
+    assert q.dtype == np.float64
+    assert q.shape == (17, 4)
+
+    q = quaternionic.array.random((17, 3))
+    assert isinstance(q, quaternionic.array)
+    assert q.dtype == np.float64
+    assert q.shape == (17, 3, 4)
+
+    q = quaternionic.array.random((17, 3, 4))
+    assert isinstance(q, quaternionic.array)
+    assert q.dtype == np.float64
+    assert q.shape == (17, 3, 4)
+
+    q = quaternionic.array.random((17, 3, 4), normalize=True)
+    assert np.max(np.abs(1 - q.abs)) < 4 * np.finfo(float).eps
diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py
index 9588e6f..41ce796 100644
--- a/tests/test_interpolation.py
+++ b/tests/test_interpolation.py
@@ -6,9 +6,7 @@
 def test_unflip_rotors(Rs):
     np.random.seed(12345)
     unflip_precision = 4e-16
-    f = 2 * np.random.rand(17, 1_000, 4) - 1
-    q = quaternionic.array(f)
-    q = q / abs(q)
+    q = quaternionic.array.random((17, 1_000, 4), normalize=True)
     ndim = q.ndim
     axis = -2
     inplace = False
diff --git a/tests/test_properties.py b/tests/test_properties.py
index 268ad48..8567ad4 100644
--- a/tests/test_properties.py
+++ b/tests/test_properties.py
@@ -47,8 +47,7 @@ def test_setting_components():
 
 def test_two_spinor():
     np.random.seed(1234)
-    arr = np.random.rand(17, 9, 4)
-    q = quaternionic.array(arr)
+    q = quaternionic.array.random((17, 9, 4))
     s = q.two_spinor
     a = q.w + 1j * q.z
     b = q.y + 1j * q.x
@@ -70,13 +69,13 @@ def test_rotate_vectors(Rs):
     with pytest.raises(ValueError):
         one.rotate(np.array(3.14))
     with pytest.raises(ValueError):
-        one.rotate(np.random.rand(17, 9, 4))
+        one.rotate(np.random.normal(size=(17, 9, 4)))
     with pytest.raises(ValueError):
-        one.rotate(np.random.rand(17, 9, 3), axis=1)
+        one.rotate(np.random.normal(size=(17, 9, 3)), axis=1)
 
     np.random.seed(1234)
     # Test (1)*(1)
-    vecs = np.random.rand(3)
+    vecs = np.random.normal(size=(3,))
     quats = z
     vecsprime = quats.rotate(vecs)
     assert np.allclose(vecsprime,
@@ -84,7 +83,7 @@ def test_rotate_vectors(Rs):
                        rtol=0.0, atol=0.0)
     assert quats.shape[:-1] + vecs.shape == vecsprime.shape, ("Out of shape!", quats.shape, vecs.shape, vecsprime.shape)
     # Test (1)*(5)
-    vecs = np.random.rand(5, 3)
+    vecs = np.random.normal(size=(5, 3))
     quats = z
     vecsprime = quats.rotate(vecs)
     for i, vec in enumerate(vecs):
@@ -93,7 +92,7 @@ def test_rotate_vectors(Rs):
                            rtol=0.0, atol=0.0)
     assert quats.shape[:-1] + vecs.shape == vecsprime.shape, ("Out of shape!", quats.shape, vecs.shape, vecsprime.shape)
     # Test (1)*(5) inner axis
-    vecs = np.random.rand(3, 5)
+    vecs = np.random.normal(size=(3, 5))
     quats = z
     vecsprime = quats.rotate(vecs, axis=-2)
     for i, vec in enumerate(vecs.T):
@@ -102,7 +101,7 @@ def test_rotate_vectors(Rs):
                            rtol=0.0, atol=0.0)
     assert quats.shape[:-1] + vecs.shape == vecsprime.shape, ("Out of shape!", quats.shape, vecs.shape, vecsprime.shape)
     # Test (N)*(1)
-    vecs = np.random.rand(3)
+    vecs = np.random.normal(size=(3))
     quats = Rs
     vecsprime = quats.rotate(vecs)
     assert np.allclose(vecsprime,
@@ -110,7 +109,7 @@ def test_rotate_vectors(Rs):
                        rtol=1e-15, atol=1e-15)
     assert quats.shape[:-1] + vecs.shape == vecsprime.shape, ("Out of shape!", quats.shape, vecs.shape, vecsprime.shape)
     # Test (N)*(5)
-    vecs = np.random.rand(5, 3)
+    vecs = np.random.normal(size=(5, 3))
     quats = Rs
     vecsprime = quats.rotate(vecs)
     for i, vec in enumerate(vecs):
@@ -119,7 +118,7 @@ def test_rotate_vectors(Rs):
                            rtol=1e-15, atol=1e-15)
     assert quats.shape[:-1] + vecs.shape == vecsprime.shape, ("Out of shape!", quats.shape, vecs.shape, vecsprime.shape)
     # Test (N)*(5) inner axis
-    vecs = np.random.rand(3, 5)
+    vecs = np.random.normal(size=(3, 5))
     quats = Rs
     vecsprime = quats.rotate(vecs, axis=-2)
     for i, vec in enumerate(vecs.T):
diff --git a/tests/test_utilities.py b/tests/test_utilities.py
index 1c37dd9..a988057 100644
--- a/tests/test_utilities.py
+++ b/tests/test_utilities.py
@@ -12,9 +12,9 @@ def f1(a, b, c):
         assert isinstance(c, np.ndarray) and isinstance(c, quaternionic.array)
         assert isinstance(d, np.ndarray) and not isinstance(d, quaternionic.array)
         return d
-    a = quaternionic.array(np.random.rand(17, 3, 4))
-    b = quaternionic.array(np.random.rand(13, 3, 4))
-    c = quaternionic.array(np.random.rand(11, 3, 4))
+    a = quaternionic.array.random((17, 3, 4))
+    b = quaternionic.array.random((13, 3, 4))
+    c = quaternionic.array.random((11, 3, 4))
     d1 = f1(a, b, c)
     assert isinstance(d1, np.ndarray) and not isinstance(d1, quaternionic.array)
     f2 = quaternionic.utilities.type_self_return(f1)
@@ -34,9 +34,9 @@ def f1(a, b, c):
         assert isinstance(c, np.ndarray) and not isinstance(c, quaternionic.array)
         assert isinstance(d, np.ndarray) and not isinstance(d, quaternionic.array)
         return d
-    a = quaternionic.array(np.random.rand(17, 3, 4))
-    b = quaternionic.array(np.random.rand(13, 3, 4))
-    c = quaternionic.array(np.random.rand(11, 3, 4))
+    a = quaternionic.array.random((17, 3, 4))
+    b = quaternionic.array.random((13, 3, 4))
+    c = quaternionic.array.random((11, 3, 4))
     f2 = quaternionic.utilities.ndarray_args(f1)
     d2 = f2(a, b, c)
     assert isinstance(d2, np.ndarray) and not isinstance(d2, quaternionic.array)
@@ -54,9 +54,9 @@ def f1(a, b, c):
         assert isinstance(c, np.ndarray) and not isinstance(c, quaternionic.array)
         assert isinstance(d, np.ndarray) and not isinstance(d, quaternionic.array)
         return d
-    a = quaternionic.array(np.random.rand(17, 3, 4))
-    b = quaternionic.array(np.random.rand(13, 3, 4))
-    c = quaternionic.array(np.random.rand(11, 3, 4))
+    a = quaternionic.array.random((17, 3, 4))
+    b = quaternionic.array.random((13, 3, 4))
+    c = quaternionic.array.random((11, 3, 4))
     f2 = quaternionic.utilities.ndarray_args_and_return(f1)
     d2 = f2(a, b, c)
     assert isinstance(d2, np.ndarray) and isinstance(d2, quaternionic.array)
@@ -97,8 +97,8 @@ def test_pyguvectorize():
     _quaternion_resolution = 10 * np.finfo(float).resolution
     np.random.seed(1234)
     one = quaternionic.array(1, 0, 0, 0)
-    x = quaternionic.array(np.random.rand(7, 13, 4))
-    y = quaternionic.array(np.random.rand(13, 4))
+    x = quaternionic.array.random((7, 13, 4))
+    y = quaternionic.array.random((13, 4))
     z = np.random.rand(13)
 
     arg0s = [one, -(1+2*_quaternion_resolution)*one, -one, x]

From ef5d99735e971c6bd854db3f0770ea46511d0cbd Mon Sep 17 00:00:00 2001
From: Michael Boyle <michael.oliver.boyle@gmail.com>
Date: Wed, 7 Apr 2021 14:38:09 -0400
Subject: [PATCH 07/10] Allow out arguments to overwrite inputs safely

---
 quaternionic/algebra.py | 72 +++++++++++++++++++++++++++--------------
 tests/test_algebra.py   | 49 +++++++++++++++++++++++++++-
 2 files changed, 96 insertions(+), 25 deletions(-)

diff --git a/quaternionic/algebra.py b/quaternionic/algebra.py
index a1e9624..f8ee5fe 100644
--- a/quaternionic/algebra.py
+++ b/quaternionic/algebra.py
@@ -64,20 +64,28 @@ def subtract(q1, q2, qout):
 @attach_typelist_and_signature([(float64[:], float64[:], float64[:])], '(n),(n)->(n)')
 def multiply(q1, q2, qout):
     """Multiply quaternions q1*q2"""
-    qout[0] = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3]
-    qout[1] = q1[0]*q2[1] + q1[1]*q2[0] + q1[2]*q2[3] - q1[3]*q2[2]
-    qout[2] = q1[0]*q2[2] - q1[1]*q2[3] + q1[2]*q2[0] + q1[3]*q2[1]
-    qout[3] = q1[0]*q2[3] + q1[1]*q2[2] - q1[2]*q2[1] + q1[3]*q2[0]
+    a = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3]
+    b = q1[0]*q2[1] + q1[1]*q2[0] + q1[2]*q2[3] - q1[3]*q2[2]
+    c = q1[0]*q2[2] - q1[1]*q2[3] + q1[2]*q2[0] + q1[3]*q2[1]
+    d = q1[0]*q2[3] + q1[1]*q2[2] - q1[2]*q2[1] + q1[3]*q2[0]
+    qout[0] = a
+    qout[1] = b
+    qout[2] = c
+    qout[3] = d
 
 
 @attach_typelist_and_signature([(float64[:], float64[:], float64[:])], '(n),(n)->(n)')
 def divide(q1, q2, qout):
     """Divide quaternions q1/q2 = q1 * q2.inverse"""
     q2norm = q2[0]**2 + q2[1]**2 + q2[2]**2 + q2[3]**2
-    qout[0] = (+q1[0]*q2[0] + q1[1]*q2[1] + q1[2]*q2[2] + q1[3]*q2[3]) / q2norm
-    qout[1] = (-q1[0]*q2[1] + q1[1]*q2[0] - q1[2]*q2[3] + q1[3]*q2[2]) / q2norm
-    qout[2] = (-q1[0]*q2[2] + q1[1]*q2[3] + q1[2]*q2[0] - q1[3]*q2[1]) / q2norm
-    qout[3] = (-q1[0]*q2[3] - q1[1]*q2[2] + q1[2]*q2[1] + q1[3]*q2[0]) / q2norm
+    a = (+q1[0]*q2[0] + q1[1]*q2[1] + q1[2]*q2[2] + q1[3]*q2[3]) / q2norm
+    b = (-q1[0]*q2[1] + q1[1]*q2[0] - q1[2]*q2[3] + q1[3]*q2[2]) / q2norm
+    c = (-q1[0]*q2[2] + q1[1]*q2[3] + q1[2]*q2[0] - q1[3]*q2[1]) / q2norm
+    d = (-q1[0]*q2[3] - q1[1]*q2[2] + q1[2]*q2[1] + q1[3]*q2[0]) / q2norm
+    qout[0] = a
+    qout[1] = b
+    qout[2] = c
+    qout[3] = d
 
 
 true_divide = divide
@@ -288,10 +296,14 @@ def sqrt(q, qout):
 @attach_typelist_and_signature([(float64[:], float64[:])], '(n)->(n)')
 def square(q, qout):
     """Return square of quaternion q*q"""
-    qout[0] = q[0]**2 - q[1]**2 - q[2]**2 - q[3]**2
-    qout[1] = 2*q[0]*q[1]
-    qout[2] = 2*q[0]*q[2]
-    qout[3] = 2*q[0]*q[3]
+    a = q[0]**2 - q[1]**2 - q[2]**2 - q[3]**2
+    b = 2*q[0]*q[1]
+    c = 2*q[0]*q[2]
+    d = 2*q[0]*q[3]
+    qout[0] = a
+    qout[1] = b
+    qout[2] = c
+    qout[3] = d
 
 
 @attach_typelist_and_signature([(float64[:], float64[:])], '(n)->(n)')
@@ -341,10 +353,14 @@ def bitwise_xor(q1, q2, qout):
     Note that the result may seem surprising because we sometimes think of quaternions as 
 
     """
-    qout[0] = q1[0]*q2[0]
-    qout[1] = q1[0]*q2[1] + q1[1]*q2[0]
-    qout[2] = q1[0]*q2[2] + q1[2]*q2[0]
-    qout[3] = q1[0]*q2[3] + q1[3]*q2[0]
+    a = q1[0]*q2[0]
+    b = q1[0]*q2[1] + q1[1]*q2[0]
+    c = q1[0]*q2[2] + q1[2]*q2[0]
+    d = q1[0]*q2[3] + q1[3]*q2[0]
+    qout[0] = a
+    qout[1] = b
+    qout[2] = c
+    qout[3] = d
 
 
 invert = conj  # reversion (= conjugate for quaternion algebra)
@@ -361,10 +377,14 @@ def left_shift(q1, q2, qout):
     ```
 
     """
-    qout[0] = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3]
-    qout[1] = q1[0]*q2[1]
-    qout[2] = q1[0]*q2[2]
-    qout[3] = q1[0]*q2[3]
+    a = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3]
+    b = q1[0]*q2[1]
+    c = q1[0]*q2[2]
+    d = q1[0]*q2[3]
+    qout[0] = a
+    qout[1] = b
+    qout[2] = c
+    qout[3] = d
 
 
 @attach_typelist_and_signature([(float64[:], float64[:], float64[:])], '(n),(n)->(n)')
@@ -378,10 +398,14 @@ def right_shift(q1, q2, qout):
     ```
 
     """
-    qout[0] = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3]
-    qout[1] = q1[1]*q2[0]
-    qout[2] = q1[2]*q2[0]
-    qout[3] = q1[3]*q2[0]
+    a = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3]
+    b = q1[1]*q2[0]
+    c = q1[2]*q2[0]
+    d = q1[3]*q2[0]
+    qout[0] = a
+    qout[1] = b
+    qout[2] = c
+    qout[3] = d
 
 
 @attach_typelist_and_signature([(float64[:], float64[:], boolean[:])], '(n),(n)->()')
diff --git a/tests/test_algebra.py b/tests/test_algebra.py
index 590551f..b27c526 100644
--- a/tests/test_algebra.py
+++ b/tests/test_algebra.py
@@ -48,6 +48,7 @@ def test_array_ufunc(array):
     with pytest.raises(NotImplementedError):
         np.negative.at(q, [0, 1])
 
+    # Addition
     p = array(np.random.normal(size=(17, 3, 4)))
     q = array(np.random.normal(size=(1, 3, 4)))
     pq1 = np.add(p, q)
@@ -59,6 +60,7 @@ def test_array_ufunc(array):
     assert np.array_equal(pq1, pq2)
     assert isinstance(pq2, array)
 
+    # Quaternion-scalar multiplication
     p = array(np.random.normal(size=(17, 3, 4)))
     q = np.random.rand(1, 3)
     pq1 = np.multiply(p, q)
@@ -68,7 +70,19 @@ def test_array_ufunc(array):
     np.multiply(p, q, out=pq2)
     assert np.array_equal(pq1, pq2)
     assert isinstance(pq2, array)
-
+    pq3 = p * q
+    assert np.array_equal(pq1, pq3)
+    assert isinstance(pq3, array)
+    pq4 = p.copy()
+    pq4 *= q
+    assert np.array_equal(pq1, pq4)
+    assert isinstance(pq4, array)
+    pq5 = p.copy()
+    np.multiply(pq5, q, out=pq5)
+    assert np.array_equal(pq1, pq5)
+    assert isinstance(pq5, array)
+
+    # Scalar-quaternion multiplication
     p = np.random.rand(1, 3)
     q = array(np.random.normal(size=(17, 3, 4)))
     pq1 = np.multiply(p, q)
@@ -78,6 +92,39 @@ def test_array_ufunc(array):
     np.multiply(p, q, out=pq2)
     assert np.array_equal(pq1, pq2)
     assert isinstance(pq2, array)
+    pq3 = p * q
+    assert np.array_equal(pq1, pq3)
+    assert isinstance(pq3, array)
+    pq4 = q.copy()
+    pq4 *= p
+    assert np.array_equal(pq1, pq4)
+    assert isinstance(pq4, array)
+    pq5 = q.copy()
+    np.multiply(p, pq5, out=pq5)
+    assert np.array_equal(pq1, pq5)
+    assert isinstance(pq5, array)
+
+    # Quaternion-quaternion multiplication
+    p = array(np.random.normal(size=(17, 3, 4)))
+    q = array(np.random.normal(size=(17, 3, 4)))
+    pq1 = np.multiply(p, q)
+    assert isinstance(pq1, array)
+    assert pq1.shape == (17, 3, 4)
+    pq2 = array(np.empty((17, 3, 4)))
+    np.multiply(p, q, out=pq2)
+    assert np.array_equal(pq1, pq2)
+    assert isinstance(pq2, array)
+    pq3 = p * q
+    assert np.array_equal(pq1, pq3)
+    assert isinstance(pq3, array)
+    pq4 = p.copy()
+    pq4 *= q
+    assert np.array_equal(pq1, pq4)
+    assert isinstance(pq4, array)
+    pq5 = p.copy()
+    np.multiply(pq5, q, out=pq5)
+    assert np.array_equal(pq1, pq5)
+    assert isinstance(pq5, array)
 
     p = np.random.rand(1, 3)
     q = np.random.normal(size=(17, 3, 4))

From 4f34f9b6c41b1c7ced5581448a8b2459d9518c76 Mon Sep 17 00:00:00 2001
From: Michael Boyle <michael.oliver.boyle@gmail.com>
Date: Wed, 7 Apr 2021 14:51:27 -0400
Subject: [PATCH 08/10] Fix warning syntax and note `random` method [skip ci]

---
 README.md | 14 ++++++++++----
 1 file changed, 10 insertions(+), 4 deletions(-)

diff --git a/README.md b/README.md
index 2d01dc5..19dd0ab 100644
--- a/README.md
+++ b/README.md
@@ -90,10 +90,7 @@ performs.
 
 | :warning: WARNING                                                                                |
 |:-------------------------------------------------------------------------------------------------|
-| Because of unfortunate default values, the `np.copy` function will not preserve the quaternionic |
-| nature of an array by default; the result will just be a plain array of floats.  You could pass  |
-| the optional argument `subok=True`, as in `q3 = np.copy(q1, subok=True)`, but it's easier to t   |
-| just use the member function: `q3 = q1.copy()`.                                                  |
+| Because of unfortunate default values, the `np.copy` function will not preserve the quaternionic nature of an array by default; the result will just be a plain array of floats.  You could pass the optional argument `subok=True`, as in `q3 = np.copy(q1, subok=True)`, but it's easier to just use the member function: `q3 = q1.copy()`. |
 
 
 ## Algebra
@@ -217,6 +214,15 @@ quaternions as rotations:
 np.max(quaternionic.distance.rotation.intrinsic(q1, q2))  # Typically around 1e-15
 ```
 
+Also note the classmethod
+
+   * `random`
+
+This constructs a quaternionic array in which each component is randomly selected from a normal
+(Gaussian) distribution centered at 0 with scale 1, which means that the result is isotropic
+(spherically symmetric).  It is also possible to pass the `normalize` argument to this function,
+which results in truly random unit quaternions.
+
 
 ## Distance functions
 

From d48d930962a330c420a2efb212278c91f15e1238 Mon Sep 17 00:00:00 2001
From: Michael Boyle <michael.oliver.boyle@gmail.com>
Date: Wed, 7 Apr 2021 14:57:49 -0400
Subject: [PATCH 09/10] Tweak wording

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 19dd0ab..a1fed28 100644
--- a/README.md
+++ b/README.md
@@ -90,7 +90,7 @@ performs.
 
 | :warning: WARNING                                                                                |
 |:-------------------------------------------------------------------------------------------------|
-| Because of unfortunate default values, the `np.copy` function will not preserve the quaternionic nature of an array by default; the result will just be a plain array of floats.  You could pass the optional argument `subok=True`, as in `q3 = np.copy(q1, subok=True)`, but it's easier to just use the member function: `q3 = q1.copy()`. |
+| Because of an unfortunate choice by the numpy developers, the `np.copy` function will not preserve the quaternionic nature of an array by default; the result will just be a plain array of floats.  You could pass the optional argument `subok=True`, as in `q3 = np.copy(q1, subok=True)`, but it's easier to just use the member function: `q3 = q1.copy()`. |
 
 
 ## Algebra

From 0e95370a725a8e10e40abfeac366a1641d074e7f Mon Sep 17 00:00:00 2001
From: Michael Boyle <michael.oliver.boyle@gmail.com>
Date: Wed, 7 Apr 2021 14:58:07 -0400
Subject: [PATCH 10/10] Test for empty tuple input to random

---
 tests/test_converters.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/tests/test_converters.py b/tests/test_converters.py
index 05443e0..e068fde 100644
--- a/tests/test_converters.py
+++ b/tests/test_converters.py
@@ -464,6 +464,11 @@ def test_random():
     assert q.dtype == np.float64
     assert q.shape == (4,)
 
+    q = quaternionic.array.random(tuple())
+    assert isinstance(q, quaternionic.array)
+    assert q.dtype == np.float64
+    assert q.shape == (4,)
+
     q = quaternionic.array.random(17)
     assert isinstance(q, quaternionic.array)
     assert q.dtype == np.float64