Cross product of vectors gives zero erroneously with expand(), but correct results without

[rho-novatron]

Given that r is a sympy.vector.CoordSys3D instance, expanding r.i+r.i gives a Mul instance in sympy >= 1.8, but a VectorMul in earlier versions. This leads to very surprising behavior, such as this:

>>> sympy.vector.cross(r.i+r.i, r.j+r.k)
(-2)*r.j + 2*r.k
>>> sympy.vector.cross((r.i+r.i).expand(), r.j+r.k)
0
>>>

Sympy 1.7.1 and earlier give the same result for both the original expression and the expanded expression, but later versions (including current git) give the above.

[rho-novatron]

This is the test I ran:

import sympy.vector as sv

r = sv.CoordSys3D("r")
print(f"type(r.i+r.i): {type(r.i+r.i)}")
print(f"type((r.i+r.i).expand()): {type((r.i+r.i).expand())}")
print(f"sv.cross((r.i+r.i), r.j+r.k)) = {sv.cross((r.i+r.i), r.j+r.k)}")
print(f"sv.cross((r.i+r.i).expand(), r.j+r.k)) = {sv.cross((r.i+r.i).expand(), r.j+r.k)}")

On sympy 1.7.1, it gives this:

(sympy-bug) rho@louie:~/git/sympy-bug-report$ pip install sympy==1.7.1
[...]
Successfully installed sympy-1.7.1
(sympy-bug) rho@louie:~/git/sympy-bug-report$ python test.py
type(r.i+r.i): <class 'sympy.vector.vector.VectorMul'>
type((r.i+r.i).expand()): <class 'sympy.vector.vector.VectorMul'>
sv.cross((r.i+r.i), r.j+r.k)) = (-2)*r.j + 2*r.k
sv.cross((r.i+r.i).expand(), r.j+r.k)) = (-2)*r.j + 2*r.k

But on 1.8 and newer, it shows

(sympy-bug) rho@louie:~/git/sympy-bug-report$ pip install sympy==1.8
[...]
Successfully installed sympy-1.8
(sympy-bug) rho@louie:~/git/sympy-bug-report$ python test.py
type(r.i+r.i): <class 'sympy.vector.vector.VectorMul'>
type((r.i+r.i).expand()): <class 'sympy.core.mul.Mul'>
sv.cross((r.i+r.i), r.j+r.k)) = (-2)*r.j + 2*r.k
sv.cross((r.i+r.i).expand(), r.j+r.k)) = 0

[rho-novatron]

With git bisect, I identified the breaking commit as 29a1a1d .

[rho-novatron]

Ok, I got a bit further. It seems the root cause might be that sympy.vector.Cross does not have a components member, and is thus ignored in BasisDependentAdd.__new__.

[rho-novatron]

Ok, so these are my thoughts now:

Since the breaking commit above, expand sometimes creates Mul instances instead of VectorMul instances. I don't understand exactly why VectorMul is a special case.

That would probably be fine, if vector.cross handled Mul instances. Currently, it has a special case for VectorMul instances, but fails for regular Mul.

Finally, BasisDependentAdd.__new__ requires that the arguments can provide components, and ignores anything passed that doesn't. Since vector.Cross does not provide components, it is ignored and this is why the cross product finally comes out as zero.

The linked PR adds another special case to vector.cross. I have no idea if this is the desired solution.

[oscarbenjamin]

The basic problem is really just this:

In [7]: from sympy.vector import CoordSys3D

In [8]: r = CoordSys3D('r')

In [9]: type(Mul(2, r.i))
Out[9]: sympy.core.mul.Mul

In [10]: type(2*r.i)
Out[10]: sympy.vector.vector.VectorMul

It isn't enough in SymPy to override __mul__ because you also need Mul to work (in fact it's a bad idea to override __mul__ in Basic subclasses because it only pretends to work).

The Matrix types have a way of making this work:

In [13]: M = MatrixSymbol('M', 2, 2)

In [14]: type(Mul(2, M))
Out[14]: sympy.matrices.expressions.matmul.MatMul

In [15]: type(2 * M)
Out[15]: sympy.matrices.expressions.matmul.MatMul

I think this is handled here:

sympy/sympy/core/operations.py

Line 101 in 64e187e

obj = cls._exec_constructor_postprocessors(obj)

Which leads through to here:

sympy/sympy/matrices/expressions/matexpr.py

Lines 474 to 516 in 64e187e

	def get_postprocessor(cls):
	def _postprocessor(expr):
	# To avoid circular imports, we can't have MatMul/MatAdd on the top level
	mat_class = {Mul: MatMul, Add: MatAdd}[cls]
	nonmatrices = []
	matrices = []
	for term in expr.args:
	if isinstance(term, MatrixExpr):
	matrices.append(term)
	else:
	nonmatrices.append(term)
	if not matrices:
	return cls._from_args(nonmatrices)
	if nonmatrices:
	if cls == Mul:
	for i in range(len(matrices)):
	if not matrices[i].is_MatrixExpr:
	# If one of the matrices explicit, absorb the scalar into it
	# (doit will combine all explicit matrices into one, so it
	# doesn't matter which)
	matrices[i] = matrices[i].__mul__(cls._from_args(nonmatrices))
	nonmatrices = []
	break
	else:
	# Maintain the ability to create Add(scalar, matrix) without
	# raising an exception. That way different algorithms can
	# replace matrix expressions with non-commutative symbols to
	# manipulate them like non-commutative scalars.
	return cls._from_args(nonmatrices + [mat_class(*matrices).doit(deep=False)])
	if mat_class == MatAdd:
	return mat_class(*matrices).doit(deep=False)
	return mat_class(cls._from_args(nonmatrices), *matrices).doit(deep=False)
	return _postprocessor
	Basic._constructor_postprocessor_mapping[MatrixExpr] = {
	"Mul": [get_postprocessor(Mul)],
	"Add": [get_postprocessor(Add)],
	}

I'm not a fan of that mechanism but BasisDependent should probably use the same mechanism as MatrixExpr for this.