@pytest.mark.parametrize('dims', [2, 3, 4, 5])
def test_geometric_algebra(dims):
    pytest.importorskip('numpy')
    import numpy as np
    from pymbolic.geometric_algebra import MultiVector as MV
    vec1 = MV(np.random.randn(dims))
    vec2 = MV(np.random.randn(dims))
    vec3 = MV(np.random.randn(dims))
    vec4 = MV(np.random.randn(dims))
    vec5 = MV(np.random.randn(dims))
    assert ((vec1 ^ vec2) + (vec1 | vec2)).close_to((vec1 * vec2))
    assert ((vec1 ^ vec2) ^ vec3).close_to((((- vec2) ^ vec1) ^ vec3))
    vecs = [vec1, vec2, vec3, vec4, vec5]
    if (len(vecs) > dims):
        from operator import xor as outer
        assert reduce(outer, vecs).close_to(0)
    assert (vec1.inv() * vec1).close_to(1)
    assert (vec1 * vec1.inv()).close_to(1)
    assert ((1 / vec1) * vec1).close_to(1)
    assert (vec1 / vec1).close_to(1)
    for (a, b, c) in [(vec1, vec2, vec3), ((vec1 * vec2), vec3, vec4), (vec1, (vec2 * vec3), vec4), (vec1, vec2, (vec3 * vec4)), (vec1, vec2, ((vec3 * vec4) * vec5)), (vec1, (vec2 * vec1), ((vec3 * vec4) * vec5))]:
        assert ((a * b) * c).close_to((a * (b * c)))
        assert ((a ^ b) ^ c).close_to((a ^ (b ^ c)))
        assert (c * b).project(0).close_to(b.scalar_product(c))
        assert (c.rev() * b).project(0).close_to(b.rev().scalar_product(c))
        assert (b.rev() * b).project(0).close_to(b.norm_squared())
        assert (b.norm_squared() >= 0)
        assert (c.norm_squared() >= 0)
        assert (b.scalar_product(c) <= ((abs(b) * abs(c)) + 1e-13))
        assert (abs((b.scalar_product((a ^ c)) - (b >> a).scalar_product(c))) < 1e-13)
        assert ((a ^ b) << c).close_to((a << (b << c)))
        assert (c << c.I.rev()).close_to((c | c.I.rev()))
        for div in (list(b.gen_blades()) + [vec1, vec1.I]):
            assert (div.inv() * div).close_to(1)
            assert (div * div.inv()).close_to(1)
            assert ((1 / div) * div).close_to(1)
            assert (div / div).close_to(1)
            assert ((c / div) * div).close_to(c)
            assert ((c * div) / div).close_to(c)
        assert (c.rev().rev() == c)
        assert (b ^ c).rev().close_to((c.rev() ^ b.rev()))
        assert c.dual().close_to((c | c.I.rev()))
        assert c.dual().close_to((c * c.I.rev()))
        assert (c.invol().invol() == c)
        assert (b ^ c).invol().close_to((b.invol() ^ c.invol()))
        assert ((a.x(b.x(c)) + b.x(c.x(a))) + c.x(a.x(b))).close_to(0)
        assert a.x((b * c)).close_to(((a.x(b) * c) + (b * a.x(c))))
@pytest.mark.parametrize('dims', [2, 3, 4, 5])
def test_geometric_algebra(dims):
    pytest.importorskip('numpy')
    import numpy as np
    from pymbolic.geometric_algebra import MultiVector as MV
    vec1 = MV(np.random.randn(dims))
    vec2 = MV(np.random.randn(dims))
    vec3 = MV(np.random.randn(dims))
    vec4 = MV(np.random.randn(dims))
    vec5 = MV(np.random.randn(dims))
    assert ((vec1 ^ vec2) + (vec1 | vec2)).close_to((vec1 * vec2))
    assert ((vec1 ^ vec2) ^ vec3).close_to((((- vec2) ^ vec1) ^ vec3))
    vecs = [vec1, vec2, vec3, vec4, vec5]
    if (len(vecs) > dims):
        from operator import xor as outer
        assert reduce(outer, vecs).close_to(0)
    assert (vec1.inv() * vec1).close_to(1)
    assert (vec1 * vec1.inv()).close_to(1)
    assert ((1 / vec1) * vec1).close_to(1)
    assert (vec1 / vec1).close_to(1)
    for (a, b, c) in [(vec1, vec2, vec3), ((vec1 * vec2), vec3, vec4), (vec1, (vec2 * vec3), vec4), (vec1, vec2, (vec3 * vec4)), (vec1, vec2, ((vec3 * vec4) * vec5)), (vec1, (vec2 * vec1), ((vec3 * vec4) * vec5))]:
        assert ((a * b) * c).close_to((a * (b * c)))
        assert ((a ^ b) ^ c).close_to((a ^ (b ^ c)))
        assert (c * b).project(0).close_to(b.scalar_product(c))
        assert (c.rev() * b).project(0).close_to(b.rev().scalar_product(c))
        assert (b.rev() * b).project(0).close_to(b.norm_squared())
        assert (b.norm_squared() >= 0)
        assert (c.norm_squared() >= 0)
        assert (b.scalar_product(c) <= ((abs(b) * abs(c)) + 1e-13))
        assert (abs((b.scalar_product((a ^ c)) - (b >> a).scalar_product(c))) < 1e-13)
        assert ((a ^ b) << c).close_to((a << (b << c)))
        assert (c << c.I.rev()).close_to((c | c.I.rev()))
        for div in (list(b.gen_blades()) + [vec1, vec1.I]):
            assert (div.inv() * div).close_to(1)
            assert (div * div.inv()).close_to(1)
            assert ((1 / div) * div).close_to(1)
            assert (div / div).close_to(1)
            assert ((c / div) * div).close_to(c)
            assert ((c * div) / div).close_to(c)
        assert (c.rev().rev() == c)
        assert (b ^ c).rev().close_to((c.rev() ^ b.rev()))
        assert c.dual().close_to((c | c.I.rev()))
        assert c.dual().close_to((c * c.I.rev()))
        assert (c.invol().invol() == c)
        assert (b ^ c).invol().close_to((b.invol() ^ c.invol()))
        assert ((a.x(b.x(c)) + b.x(c.x(a))) + c.x(a.x(b))).close_to(0)
        assert a.x((b * c)).close_to(((a.x(b) * c) + (b * a.x(c))))