The substitution *must* be done *after* the `bs**e` step. That's kind of the whole point of this extra condition.
Initially, we were having problems with expressions like:
`(x**n + O(x**m))**1.2 + O(x**p)`
The problem what that to add two order terms, you must know the order of both of them. So, here, we need to find out the order of the term raised to the power of 1.2. To do that, SymPy called `calculate_series` in gruntz which after a couple of calculations ended up again at this step (where you need to get the sum of terms contaning order). And that led to recursion.
Note that there was no problem if we need to add something like:
`(x + x**3 + O(x**4))**1.2 + O(x**8)`
that is when the entity being raised to a power doesn't have only one other term except the order term. So, for the very specific case where `bs` is a two term sum raised to a power, this method has been employed. This is just a different method to get the order of `(x**n + O(x**m))**1.2`, where `bs = x**n + O(x**m)` and `e = 1.2`.
Getting back to the point, that is why we need to first do a subs, then `bs**e` and then for the dummy, we need to subs back `O(1)`. Hope this clarifies the issue.