diff --git a/tests/core/test_algorithm.py b/tests/core/test_algorithm.py
index 03441b3c..22cc3499 100644
--- a/tests/core/test_algorithm.py
+++ b/tests/core/test_algorithm.py
@@ -17,6 +17,8 @@ def test_nnls(rng):
 
     from glass.core.algorithm import nnls as nnls_glass
 
+    # cross-check output with scipy's nnls
+
     a = rng.standard_normal((100, 20))
     b = rng.standard_normal((100,))
 
@@ -24,3 +26,12 @@ def test_nnls(rng):
     x_scipy, _ = nnls_scipy(a, b)
 
     np.testing.assert_allclose(x_glass, x_scipy)
+
+    # check matrix and vector's shape
+
+    with pytest.raises(ValueError, match="input `a` is not a matrix"):
+        nnls_glass(b, a)
+    with pytest.raises(ValueError, match="input `b` is not a vector"):
+        nnls_glass(a, a)
+    with pytest.raises(ValueError, match="the shapes of `a` and `b` do not match"):
+        nnls_glass(a.T, b)
diff --git a/tests/core/test_array.py b/tests/core/test_array.py
index 8268d803..de2aee65 100644
--- a/tests/core/test_array.py
+++ b/tests/core/test_array.py
@@ -1,5 +1,46 @@
 import numpy as np
 import numpy.testing as npt
+import pytest
+
+# check if scipy is available for testing
+try:
+    import scipy
+except ImportError:
+    HAVE_SCIPY = False
+else:
+    del scipy
+    HAVE_SCIPY = True
+
+
+def broadcast_first():
+    from glass.core.array import broadcast_first
+
+    a = np.ones((2, 3, 4))
+    b = np.ones((2, 1))
+
+    # arrays with shape ((3, 4, 2)) and ((1, 2)) are passed
+    # to np.broadcast_arrays; hence it works
+    a_a, b_a = broadcast_first(a, b)
+    assert a_a.shape == (2, 3, 4)
+    assert b_a.shape == (2, 3, 4)
+
+    # plain np.broadcast_arrays will not work
+    with pytest.raises(ValueError, match="shape mismatch"):
+        np.broadcast_arrays(a, b)
+
+    # arrays with shape ((5, 6, 4)) and ((6, 5)) are passed
+    # to np.broadcast_arrays; hence it will not work
+    a = np.ones((4, 5, 6))
+    b = np.ones((5, 6))
+
+    with pytest.raises(ValueError, match="shape mismatch"):
+        broadcast_first(a, b)
+
+    # plain np.broadcast_arrays will work
+    a_a, b_a = broadcast_first(a, b)
+
+    assert a_a.shape == (4, 5, 6)
+    assert b_a.shape == (4, 5, 6)
 
 
 def test_broadcast_leading_axes():
@@ -114,3 +155,54 @@ def test_trapz_product():
     s = trapz_product((x1, f1), (x2, f2))
 
     assert np.allclose(s, 1.0)
+
+
+@pytest.mark.skipif(not HAVE_SCIPY, reason="test requires SciPy")
+def test_cumtrapz():
+    from scipy.integrate import cumulative_trapezoid
+
+    from glass.core.array import cumtrapz
+
+    # 1D f and x
+
+    f = np.array([1, 2, 3, 4])
+    x = np.array([0, 1, 2, 3])
+
+    # default dtype (int - not supported by scipy)
+
+    glass_ct = cumtrapz(f, x)
+    npt.assert_allclose(glass_ct, np.array([0, 1, 4, 7]))
+
+    # explicit dtype (float)
+
+    glass_ct = cumtrapz(f, x, dtype=float)
+    scipy_ct = cumulative_trapezoid(f, x, initial=0)
+    npt.assert_allclose(glass_ct, scipy_ct)
+
+    # explicit return array
+
+    result = cumtrapz(f, x, dtype=float, out=np.zeros((4,)))
+    scipy_ct = cumulative_trapezoid(f, x, initial=0)
+    npt.assert_allclose(result, scipy_ct)
+
+    # 2D f and 1D x
+
+    f = np.array([[1, 4, 9, 16], [2, 3, 5, 7]])
+    x = np.array([0, 1, 2.5, 4])
+
+    # default dtype (int - not supported by scipy)
+
+    glass_ct = cumtrapz(f, x)
+    npt.assert_allclose(glass_ct, np.array([[0, 2, 12, 31], [0, 2, 8, 17]]))
+
+    # explicit dtype (float)
+
+    glass_ct = cumtrapz(f, x, dtype=float)
+    scipy_ct = cumulative_trapezoid(f, x, initial=0)
+    npt.assert_allclose(glass_ct, scipy_ct)
+
+    # explicit return array
+
+    glass_ct = cumtrapz(f, x, dtype=float, out=np.zeros((2, 4)))
+    scipy_ct = cumulative_trapezoid(f, x, initial=0)
+    npt.assert_allclose(glass_ct, scipy_ct)