Repository: sympy/sympy
Pull request: 14547: avoid evaluation of (Dimension|Quantity - 1)
Issue #14547 has, at its root, a dimensional quantity being subtracted from 1 to see how it compares to 1. This raises an error for dimensional mismatch of addends. 

Addition of numbers (and expressions containing free symbols) to quantities or dimensions is now ignored but can be tested with the added `check_dimensions` function.

closes #14886 by replacement
fixes #14547
fixes #14932 

Better unification of units with the `quantity_simplify` function is now given, too.

#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. If there is no release notes entry for this PR,
write "NO ENTRY". The bot will check your release notes automatically to see
if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* physics.units
  * checking for consistency of terms is no longer done at instantiation
  * `check_dimensions` was added to units/util.py to check expressions for dimensional consistency
  * `__radd__` and `__rmul__` are defined for Dimensions
  * `quantity_simplify` now reduces all Quantities in an Expr with the same Dimension to the same internally canonical Quantity
<!-- END RELEASE NOTES -->

TODO
- [x] see if `__rsub__` and `__rdiv__` need to be defined for Dimension
`__rsub__` is not needed because negation turns `-Dim` into `Dim`
- [x] keep sign of other when subtracting?
This is being left unchanged

Problem statement:
false error message with Eq() and quantities in log()
I believe the below is a bug:
```python
from sympy import Eq, log
from sympy.physics.units import Quantity, mass, gram, kilogram, length, meter
m1 = Quantity('m_1', mass, kilogram)
m2 = Quantity('m_2', mass, kilogram)
m3 = Quantity('m_3', mass, kilogram)
print log(m1 + m2) - log(m3)
try:
    eq1 = Eq(log(m3), log(m1 + m2))
    print eq1
except Exception, error:
    print error
```
returns:
```python
-log(m_3) + log(m_1 + m_2)
summation of quantities of incompatible dimensions
```
`log(m_3)` has the same dimensions as `log(m1 + m2)`, so the error message is wrong. Interestingly, the error message does not pop up in the summation, so it might be something related to `Eq()`.


simplify(log(foot) + foo) raises an error
```
>>> a=log(x)
>>> from sympy.physics.units import foot, inch
>>> b=log(foot)
>>> a+b
log(foot) + log(x)
>>> simplify(_)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "sympy\simplify\simplify.py", line 525, in simplify

    original_expr = expr = signsimp(expr)
  File "sympy\simplify\simplify.py", line 375, in signsimp

    e = sub_post(sub_pre(expr))
  File "sympy\simplify\cse_opts.py", line 15, in sub_pre
    reps = [a for a in e.atoms(Add) if a.could_extract_minus_sign()]
  File "sympy\simplify\cse_opts.py", line 15, in <listcomp
>
    reps = [a for a in e.atoms(Add) if a.could_extract_minus_sign()]
  File "sympy\core\expr.py", line 2272, in could_extract_m
inus_sign
    (negative_self).extract_multiplicatively(-1) is not None)
  File "sympy\core\expr.py", line 2037, in extract_multipl
icatively
    quotient = self / c
  File "sympy\core\decorators.py", line 91, in __sympifyit
_wrapper
    return func(a, b)
  File "sympy\core\decorators.py", line 132, in binary_op_
wrapper
    return func(self, other)
  File "sympy\core\expr.py", line 177, in __div__
    return Mul(self, Pow(other, S.NegativeOne))
  File "sympy\core\cache.py", line 93, in wrapper
    retval = cfunc(*args, **kwargs)
  File "sympy\core\operations.py", line 48, in __new__
    obj = cls._exec_constructor_postprocessors(obj)
  File "sympy\core\basic.py", line 1724, in _exec_construc
tor_postprocessors
    obj = f(obj)
  File "sympy\physics\units\quantities.py", line 261, in _
Quantity_constructor_postprocessor_Add
    for i in expr.args
  File "sympy\physics\units\quantities.py", line 262, in <
setcomp>
    if i.free_symbols == set()  # do not raise if there are symbols
  File "sympy\physics\units\dimensions.py", line 442, in g
et_dimensional_dependencies
    return {str(i): j for i, j in dimdep.items()}
AttributeError: 'NoneType' object has no attribute 'items'
```

I get similar errors if I replace `log(x)` with `log(2)`, `2` and `log(inch)`.

Hints from issue discussion:
I am working on this issue. I found out that the problem begins from `fin = L, R = [i.is_finite for i in (lhs, rhs)]`. Now I am unable to figure how it propagates to  `File "/home/gagandeep/anaconda3/lib/python3.6/site-packages/sympy/physics/units/quantities.py", line 173, in _Quantity_constructor_postprocessor_Add
    raise ValueError("summation of quantities of incompatible dimensions")`.
Can you please help me in fixing this issue? Any hint will be welcomed. :-)
`_Quantity_constructor_postprocessor_Add` should raise an exception if the dimensions of summands are incompatible. I have no idea why it raises an exception in this case, as the dimensions of `m_1`, `m_2` and `m_3` are the same. If you type `log(m1 + m2) - log(m3)`, no exception is raised, so it works correctly. The weird thing is that it raises an exception if we put the same expression in `Eq()`. Can you see what is done differently in `Eq(log(m3), log(m1 + m2))`?
Hopefully I have found out the problem. Going to fix it as soon as possible. :-)
I think there is an issue with `log(m1)`. According to concepts of dimensional analysis, `log(x)`, x should be dimensionless. I think the weird behavior is due to this passing of dimensional quantity into log. I have also noticed that sometimes, `eq1 = Eq(log(m3), log(m1 + m2))` executes correctly and sometimes throws errors which begin at `dif.is_zero` and sometimes at `n.is_zero`.  
Good point about log(unit) and dimensional analysis! According to https://en.wikipedia.org/wiki/Dimensional_analysis, the basic rule ` log(a/b) = log a − log b` is only valid if both a and b are non-dimensional. Otherwise, even if a and be have the same unit, the lhs would be dimensionless, but the rhs would have units of log(unit). See also this discussion: https://physics.stackexchange.com/questions/13060/what-is-the-logarithm-of-a-kilometer-is-it-a-dimensionless-number
Does this mean that variables in log() and exponents have to be unit-less? This would also mean that derivations resulting from log-log plots do not have a physical meaning. I must have misunderstood something.
Yes they do have as far as the knowledge I have gained from [this](https://cdn-pubs.acs.org/doi/pdfplus/10.1021/ed072p954) paper. Actually what this document says is that we should divide the quantity x by its unit dimensions so that it becomes dimensionless, like log(x) is actually log(x/units(x)), units(x) denotes the units of x. Now, what I am trying to do is to convert `arg` in `log(arg)` to symbol whenever it is an instance of class of `Quantity`. I am not sure that it would fix the issue but it would definitely strengthen the conceptuality of `log`, `exp` and `trig` methods in sympy. The problem I faced is that, when `arg = a+c+d*e/f` then it becomes tricky to convert each of a,c,d,e, and f to symbols and then reusing them again in the expression of `arg`. Can you please help me by providing hints or solutions to this problem. 
In addition, I tried to make `arg` dimensionless by redefining it as `arg = Quantity(arg.name, Dimension(1))` but `arg` returned the dimension it originally had. That's why I resorted to converting it to symbols.
Waiting for your response.
Awesome paper, thanks for sharing! In my understanding, the paper does not say that we should divide x by its unit, but that fundamentally terms in logarithms are unit-less, and if they are not, such as e.g. `mu_i = mu_i0 + R*T*log(p_i)` this is because we mistakenly omitted division by a standard value, i.e. this should in fact be: `mu_i = mu_i0 + R*T*log(p_i/p0)`
Thus, all we should do is allow log(Quantity) as long as Quantity is dimensionless. I would not convert to symbols, as this could get very confusing.
Sure, I will implement the method you have agreed upon and will try to submit a PR as soon as possible. Thanks for helping me out with my queries. :-)
One of the challenges will be to make `print log(m1 + m2) - log(m3)` and `Eq(log(m3), log(m1 + m2))` return the same (informative) error if m1, m2 and m3 are not non-dimensional. I don't understand why only one returns an error. Thanks for looking into this!
I have created the PR, https://github.com/sympy/sympy/pull/14886, for solving the issue. My code converts any dimensional expression passed to log into dimensionless expression and then the further processing is done. However, there is no effect on the original quantities. I have verified on a local program and `Eq(log(m3), log(m1 + m2))` was executed successfully after applying the above mechanism. I can confirm by adding an appropriate test, though I am finding it quite tricky to find a good test for this issue. Can you please suggest such tests? 
Here is a suggested fix:

```diff
diff --git a/sympy/physics/units/dimensions.py b/sympy/physics/units/dimensions.py
index 2aabd59..561f06d 100644
--- a/sympy/physics/units/dimensions.py
+++ b/sympy/physics/units/dimensions.py
@@ -427,8 +427,10 @@ def _get_dimensional_dependencies_for_name(self, name):

             if isinstance(result, Dimension):
                 return dimsys_default.get_dimensional_dependencies(result)
-            # TODO shall we consider a result that is not a dimension?
-            # return Dimension._get_dimensional_dependencies_for_name(result)
+            elif result.func == name.func:
+                return {}
+            else:
+                return dimsys_default._get_dimensional_dependencies_for_name(result)

     def get_dimensional_dependencies(self, name, mark_dimensionless=False):
         if isinstance(name, Dimension):
diff --git a/sympy/physics/units/tests/test_quantities.py b/sympy/physics/units/
tests/test_quantities.py
index 8649301..af5b610 100644
--- a/sympy/physics/units/tests/test_quantities.py
+++ b/sympy/physics/units/tests/test_quantities.py
@@ -459,3 +459,7 @@ def test_eval_subs():
     expr2 = force/mass
     units = {force:gravitational_constant*kilogram**2/meter**2, mass:kilogram}
     assert expr2.subs(units) == gravitational_constant*kilogram/meter**2
+
+
+def test_issue_14923():
+    assert (log(inch) - log(2)).simplify() == log(inch/2)
```

Commit message: remove Dimension postprocessing for Add

Patch to review:
diff --git a/sympy/physics/units/dimensions.py b/sympy/physics/units/dimensions.py
index 2aabd5923875..ec3b165c11a1 100644
--- a/sympy/physics/units/dimensions.py
+++ b/sympy/physics/units/dimensions.py
@@ -18,7 +18,9 @@
 
 from sympy import Integer, Matrix, S, Symbol, sympify, Basic, Tuple, Dict, default_sort_key
 from sympy.core.compatibility import reduce, string_types
+from sympy.core.basic import Basic
 from sympy.core.expr import Expr
+from sympy.core.power import Pow
 from sympy.utilities.exceptions import SymPyDeprecationWarning
 
 
@@ -150,26 +152,31 @@ def _register_as_base_dim(self):
         dimsys_default._args = (dimsys_default.args[:2] + (Dict(d),))
 
     def __add__(self, other):
-        """Define the addition for Dimension.
-
-        Addition of dimension has a sense only if the second object is
-        the same dimension (we don't add length to time).
-        """
-
-        if not isinstance(other, Dimension):
-            raise TypeError("Only dimension can be added; '%s' is not valid"
-                            % type(other))
-        elif isinstance(other, Dimension) and self != other:
-            raise ValueError("Only dimension which are equal can be added; "
-                             "'%s' and '%s' are different" % (self, other))
-
+        from sympy.physics.units.quantities import Quantity
+        other = sympify(other)
+        if isinstance(other, Basic):
+            if other.has(Quantity):
+                other = Dimension(Quantity.get_dimensional_expr(other))
+            if isinstance(other, Dimension):
+                if self == other:
+                    return self
+                raise ValueError('mismatched dimensions of addends')
+            return super(Dimension, self).__add__(other)
         return self
 
+    def __radd__(self, other):
+        return self + other
+
     def __sub__(self, other):
         # there is no notion of ordering (or magnitude) among dimension,
         # subtraction is equivalent to addition when the operation is legal
         return self + other
 
+    def __rsub__(self, other):
+        # there is no notion of ordering (or magnitude) among dimension,
+        # subtraction is equivalent to addition when the operation is legal
+        return self + other
+
     def __pow__(self, other):
         return self._eval_power(other)
 
@@ -178,16 +185,22 @@ def _eval_power(self, other):
         return Dimension(self.name**other)
 
     def __mul__(self, other):
-        if not isinstance(other, Dimension):
-            return self
+        from sympy.physics.units.quantities import Quantity
+        if isinstance(other, Basic):
+            if other.has(Quantity):
+                other = Dimension(Quantity.get_dimensional_expr(other))
+            if isinstance(other, Dimension):
+                return Dimension(self.name*other.name)
+            if not other.free_symbols:
+                return self
+            return super(Dimension, self).__mul__(other)
+        return self
 
-        return Dimension(self.name*other.name)
+    def __rmul__(self, other):
+        return self*other
 
     def __div__(self, other):
-        if not isinstance(other, Dimension):
-            return self
-
-        return Dimension(self.name/other.name)
+        return self*Pow(other, -1)
 
     def __rdiv__(self, other):
         return other * pow(self, -1)
@@ -407,28 +420,31 @@ def _get_dimensional_dependencies_for_name(self, name):
         if name.is_Number:
             return {}
 
+        get_for_name = dimsys_default._get_dimensional_dependencies_for_name
+
         if name.is_Mul:
             ret = collections.defaultdict(int)
-            dicts = [dimsys_default._get_dimensional_dependencies_for_name(i) for i in name.args]
+            dicts = [get_for_name(i) for i in name.args]
             for d in dicts:
                 for k, v in d.items():
                     ret[k] += v
             return {k: v for (k, v) in ret.items() if v != 0}
 
         if name.is_Pow:
-            dim = dimsys_default._get_dimensional_dependencies_for_name(name.base)
+            dim = get_for_name(name.base)
             return {k: v*name.exp for (k, v) in dim.items()}
 
         if name.is_Function:
             args = (Dimension._from_dimensional_dependencies(
-                dimsys_default._get_dimensional_dependencies_for_name(arg)
-            ) for arg in name.args)
+                get_for_name(arg)) for arg in name.args)
             result = name.func(*args)
 
             if isinstance(result, Dimension):
                 return dimsys_default.get_dimensional_dependencies(result)
-            # TODO shall we consider a result that is not a dimension?
-            # return Dimension._get_dimensional_dependencies_for_name(result)
+            elif result.func == name.func:
+                return {}
+            else:
+                return get_for_name(result)
 
     def get_dimensional_dependencies(self, name, mark_dimensionless=False):
         if isinstance(name, Dimension):
diff --git a/sympy/physics/units/quantities.py b/sympy/physics/units/quantities.py
index 7c6bb7c92f0b..8677f7ce57bb 100644
--- a/sympy/physics/units/quantities.py
+++ b/sympy/physics/units/quantities.py
@@ -247,28 +247,3 @@ def convert_to(self, other):
     def free_symbols(self):
         """Return free symbols from quantity."""
         return self.scale_factor.free_symbols
-
-
-def _Quantity_constructor_postprocessor_Add(expr):
-    # Construction postprocessor for the addition,
-    # checks for dimension mismatches of the addends, thus preventing
-    # expressions like `meter + second` to be created.
-
-    deset = {
-        tuple(sorted(dimsys_default.get_dimensional_dependencies(
-            Dimension(Quantity.get_dimensional_expr(i) if not i.is_number else 1
-        )).items()))
-        for i in expr.args
-        if i.free_symbols == set()  # do not raise if there are symbols
-                    # (free symbols could contain the units corrections)
-    }
-    # If `deset` has more than one element, then some dimensions do not
-    # match in the sum:
-    if len(deset) > 1:
-        raise ValueError("summation of quantities of incompatible dimensions")
-    return expr
-
-
-Basic._constructor_postprocessor_mapping[Quantity] = {
-    "Add" : [_Quantity_constructor_postprocessor_Add],
-}
diff --git a/sympy/physics/units/tests/test_dimensions.py b/sympy/physics/units/tests/test_dimensions.py
index 986ff26dcda5..27b846402b6a 100644
--- a/sympy/physics/units/tests/test_dimensions.py
+++ b/sympy/physics/units/tests/test_dimensions.py
@@ -5,6 +5,7 @@
 
 from sympy import S, Symbol, sqrt
 from sympy.physics.units.dimensions import Dimension, length, time, dimsys_default
+from sympy.physics.units import foot
 from sympy.utilities.pytest import raises
 
 
@@ -65,17 +66,33 @@ def test_Dimension_properties():
 
 
 def test_Dimension_add_sub():
-    assert length + length == length
+    assert length + length == length + foot == foot + length == length
     assert length - length == length
     assert -length == length
 
-    raises(TypeError, lambda: length + 1)
-    raises(TypeError, lambda: length - 1)
     raises(ValueError, lambda: length + time)
     raises(ValueError, lambda: length - time)
 
+    # issue 14547 - only raise error for dimensional args; allow
+    # others to pass
+    x = Symbol('x')
+    e = length + x
+    assert e == x + length and e.is_Add and set(e.args) == {length, x}
+    e = length + 1
+    assert e == 1 + length == 1 - length and e.is_Add and set(e.args) == {length, 1}
+
 
 def test_Dimension_mul_div_exp():
+    assert 2*length == length*2 == length/2 == length
+    assert 2/length == 1/length
+    x = Symbol('x')
+    m = x*length
+    assert m == length*x and m.is_Mul and set(m.args) == {x, length}
+    d = x/length
+    assert d == x*length**-1 and d.is_Mul and set(d.args) == {x, 1/length}
+    d = length/x
+    assert d == length*x**-1 and d.is_Mul and set(d.args) == {1/x, length}
+
     velo = length / time
 
     assert (length * length) == length ** 2
diff --git a/sympy/physics/units/tests/test_quantities.py b/sympy/physics/units/tests/test_quantities.py
index 8649301a1a46..4b94674f7e91 100644
--- a/sympy/physics/units/tests/test_quantities.py
+++ b/sympy/physics/units/tests/test_quantities.py
@@ -283,8 +283,7 @@ def test_issue_quart():
 
 
 def test_issue_5565():
-    raises(ValueError, lambda: m < s)
-    assert (m < km).is_Relational
+    assert (m < s).is_Relational
 
 
 def test_find_unit():
@@ -317,20 +316,6 @@ def test_Quantity_derivative():
     assert diff(meter**2, meter) == 2*meter
 
 
-def test_sum_of_incompatible_quantities():
-    raises(ValueError, lambda: meter + 1)
-    raises(ValueError, lambda: meter + second)
-    raises(ValueError, lambda: 2 * meter + second)
-    raises(ValueError, lambda: 2 * meter + 3 * second)
-    raises(ValueError, lambda: 1 / second + 1 / meter)
-    raises(ValueError, lambda: 2 * meter*(mile + centimeter) + km)
-
-    expr = 2 * (mile + centimeter)/second + km/hour
-    assert expr in Basic._constructor_postprocessor_mapping
-    for i in expr.args:
-        assert i in Basic._constructor_postprocessor_mapping
-
-
 def test_quantity_postprocessing():
     q1 = Quantity('q1')
     q2 = Quantity('q2')
@@ -459,3 +444,32 @@ def test_eval_subs():
     expr2 = force/mass
     units = {force:gravitational_constant*kilogram**2/meter**2, mass:kilogram}
     assert expr2.subs(units) == gravitational_constant*kilogram/meter**2
+
+
+def test_issue_14932():
+    assert (log(inch) - log(2)).simplify() == log(inch/2)
+    assert (log(inch) - log(foot)).simplify() == -log(12)
+    p = symbols('p', positive=True)
+    assert (log(inch) - log(p)).simplify() == log(inch/p)
+
+
+def test_issue_14547():
+    # the root issue is that an argument with dimensions should
+    # not raise an error when the the `arg - 1` calculation is
+    # performed in the assumptions system
+    from sympy.physics.units import foot, inch
+    from sympy import Eq
+    assert log(foot).is_zero is None
+    assert log(foot).is_positive is None
+    assert log(foot).is_nonnegative is None
+    assert log(foot).is_negative is None
+    assert log(foot).is_algebraic is None
+    assert log(foot).is_rational is None
+    # doesn't raise error
+    assert Eq(log(foot), log(inch)) is not None  # might be False or unevaluated
+
+    x = Symbol('x')
+    e = foot + x
+    assert e.is_Add and set(e.args) == {foot, x}
+    e = foot + 1
+    assert e.is_Add and set(e.args) == {foot, 1}
diff --git a/sympy/physics/units/tests/test_util.py b/sympy/physics/units/tests/test_util.py
index 3d961f1b0f98..079404aed207 100644
--- a/sympy/physics/units/tests/test_util.py
+++ b/sympy/physics/units/tests/test_util.py
@@ -6,14 +6,17 @@
 
 from sympy.utilities.exceptions import SymPyDeprecationWarning
 
-from sympy import Add, Mul, Pow, Tuple, pi, sin, sqrt, sstr, sympify
+from sympy import (Add, Mul, Pow, Tuple, pi, sin, sqrt, sstr, sympify,
+    symbols)
 from sympy.physics.units import (
     G, centimeter, coulomb, day, degree, gram, hbar, hour, inch, joule, kelvin,
     kilogram, kilometer, length, meter, mile, minute, newton, planck,
     planck_length, planck_mass, planck_temperature, planck_time, radians,
-    second, speed_of_light, steradian, time)
+    second, speed_of_light, steradian, time, km)
 from sympy.physics.units.dimensions import dimsys_default
-from sympy.physics.units.util import convert_to, dim_simplify
+from sympy.physics.units.util import convert_to, dim_simplify, check_dimensions
+from sympy.utilities.pytest import raises
+
 
 
 def NS(e, n=15, **options):
@@ -142,7 +145,23 @@ def test_quantity_simplify():
     x, y = symbols('x y')
 
     assert quantity_simplify(x*(8*kilo*newton*meter + y)) == x*(8000*meter*newton + y)
-    assert quantity_simplify(foot*inch*(foot + inch)) == foot**2*(foot + inch)/12
-    assert quantity_simplify(foot*inch*(foot*foot + inch*(foot + inch))) == foot**2*(foot**2 + inch*(foot + inch))/12
-    assert quantity_simplify(2**(foot/inch*kilo/1000)*inch) == 4096*inch
+    assert quantity_simplify(foot*inch*(foot + inch)) == foot**2*(foot + foot/12)/12
+    assert quantity_simplify(foot*inch*(foot*foot + inch*(foot + inch))) == foot**2*(foot**2 + foot/12*(foot + foot/12))/12
+    assert quantity_simplify(2**(foot/inch*kilo/1000)*inch) == 4096*foot/12
     assert quantity_simplify(foot**2*inch + inch**2*foot) == 13*foot**3/144
+
+
+def test_check_dimensions():
+    x = symbols('x')
+    assert check_dimensions(inch + x) == inch + x
+    assert check_dimensions(length + x) == length + x
+    # after subs we get 2*length; check will clear the constant
+    assert check_dimensions((length + x).subs(x, length)) == length
+    raises(ValueError, lambda: check_dimensions(inch + 1))
+    raises(ValueError, lambda: check_dimensions(length + 1))
+    raises(ValueError, lambda: check_dimensions(length + time))
+    raises(ValueError, lambda: check_dimensions(meter + second))
+    raises(ValueError, lambda: check_dimensions(2 * meter + second))
+    raises(ValueError, lambda: check_dimensions(2 * meter + 3 * second))
+    raises(ValueError, lambda: check_dimensions(1 / second + 1 / meter))
+    raises(ValueError, lambda: check_dimensions(2 * meter*(mile + centimeter) + km))
diff --git a/sympy/physics/units/util.py b/sympy/physics/units/util.py
index 9b520c34ee29..9175c4158b5e 100644
--- a/sympy/physics/units/util.py
+++ b/sympy/physics/units/util.py
@@ -132,7 +132,7 @@ def get_total_scale_factor(expr):
 
 def quantity_simplify(expr):
     """Return an equivalent expression in which prefixes are replaced
-    with numerical values and products of related quantities are
+    with numerical values and all units of a given dimension are the
     unified in a canonical manner.
 
     Examples
@@ -143,34 +143,70 @@ def quantity_simplify(expr):
     >>> from sympy.physics.units import foot, inch
     >>> quantity_simplify(kilo*foot*inch)
     250*foot**2/3
+    >>> quantity_simplify(foot - 6*inch)
+    foot/2
     """
 
-    if expr.is_Atom:
+    if expr.is_Atom or not expr.has(Prefix, Quantity):
         return expr
 
-    if isinstance(expr, Prefix):
-        return expr.scale_factor
+    # replace all prefixes with numerical values
+    p = expr.atoms(Prefix)
+    expr = expr.xreplace({p: p.scale_factor for p in p})
 
-    expr = expr.func(*map(quantity_simplify, expr.args))
+    # replace all quantities of given dimension with a canonical
+    # quantity, chosen from those in the expression
+    d = sift(expr.atoms(Quantity), lambda i: i.dimension)
+    for k in d:
+        if len(d[k]) == 1:
+            continue
+        v = list(ordered(d[k]))
+        ref = v[0]/v[0].scale_factor
+        expr = expr.xreplace({vi: ref*vi.scale_factor for vi in v[1:]})
+
+    return expr
 
-    if not expr.is_Mul or not expr.has(Quantity):
-        return expr
 
-    args_pow = [arg.as_base_exp() for arg in expr.args]
-    quantity_pow, other_pow = sift(
-        args_pow, lambda x: isinstance(x[0], Quantity), binary=True)
-    coeff = Mul.fromiter([
-        Pow(b, e, evaluate=False) for b, e in other_pow])
-    quantity_pow_by_dim = sift(quantity_pow, lambda x: x[0].dimension)
-    new_quantities = []
-    for _, bp in quantity_pow_by_dim.items():
-        if len(bp) == 1:
-            new_quantities.append(bp[0][0]**bp[0][1])
-        else:
-            # just let reference quantity be the first quantity,
-            # picked from an ordered list
-            bp = list(ordered(bp))
-            ref = bp[0][0]/bp[0][0].scale_factor
-            new_quantities.append(Mul.fromiter((
-                (ref*b.scale_factor)**p for b, p in bp)))
-    return coeff*Mul.fromiter(new_quantities)
+def check_dimensions(expr):
+    """Return expr if there are not unitless values added to
+    dimensional quantities, else raise a ValueError."""
+    from sympy.solvers.solveset import _term_factors
+    # the case of adding a number to a dimensional quantity
+    # is ignored for the sake of SymPy core routines, so this
+    # function will raise an error now if such an addend is
+    # found.
+    # Also, when doing substitutions, multiplicative constants
+    # might be introduced, so remove those now
+    adds = expr.atoms(Add)
+    DIM_OF = dimsys_default.get_dimensional_dependencies
+    for a in adds:
+        deset = set()
+        for ai in a.args:
+            if ai.is_number:
+                deset.add(())
+                continue
+            dims = []
+            skip = False
+            for i in Mul.make_args(ai):
+                if i.has(Quantity):
+                    i = Dimension(Quantity.get_dimensional_expr(i))
+                if i.has(Dimension):
+                    dims.extend(DIM_OF(i).items())
+                elif i.free_symbols:
+                    skip = True
+                    break
+            if not skip:
+                deset.add(tuple(sorted(dims)))
+                if len(deset) > 1:
+                    raise ValueError(
+                        "addends have incompatible dimensions")
+
+    # clear multiplicative constants on Dimensions which may be
+    # left after substitution
+    reps = {}
+    for m in expr.atoms(Mul):
+        if any(isinstance(i, Dimension) for i in m.args):
+            reps[m] = m.func(*[
+                i for i in m.args if not i.is_number])
+
+    return expr.xreplace(reps)


Code review comments:
diff_hunk: @@ -65,17 +66,33 @@ def test_Dimension_properties():
 
 
 def test_Dimension_add_sub():
-    assert length + length == length
+    assert length + length == length + foot == foot + length == length
     assert length - length == length
     assert -length == length
 
-    raises(TypeError, lambda: length + 1)
-    raises(TypeError, lambda: length - 1)
     raises(ValueError, lambda: length + time)
path: sympy/physics/units/tests/test_dimensions.py
text: @user1:
I thought that dimensions are not checked for consistency any more, unless specified, so `length + time` should not raise an error.

@author:
hmm...good catch. Not intentionally left behind. I will look into it.