Is there something for implication or equivalence of equations? I was trying the example from Smaug123's sibling comment (https://news.ycombinator.com/item?id=15728598), but

  >>> import sympy as sp
  >>> x = sp.Symbol('x')
  >>> sp.simplify(sp.Implies(sp.Eq(x**2 + 2*x + 1, 0), sp.Eq(x, -1)))

  Eq(x, -1) | Ne(x**2 + 2*x + 1, 0)

  >>> sp.solve(sp.simplify(sp.Implies(sp.Eq(x**2 + 2*x + 1, 0), sp.Eq(x, -1))))

  Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/user/.local/lib/python3.5/site-packages/sympy/solvers/solvers.py", line 1065, in solve
    solution = _solve(f[0], *symbols, **flags)
  File "/home/user/.local/lib/python3.5/site-packages/sympy/solvers/solvers.py", line 1401, in _solve
    f_num, sol = solve_linear(f, symbols=symbols)
  File "/home/user/.local/lib/python3.5/site-packages/sympy/solvers/solvers.py", line 1971, in solve_linear
    eq = lhs - rhs
  TypeError: unsupported operand type(s) for -: 'Or' and 'int'

  >>> sp.solveset(sp.simplify(sp.Implies(sp.Eq(x**2 + 2*x + 1, 0), sp.Eq(x, -1))))

  Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/user/.local/lib/python3.5/site-packages/sympy/solvers/solveset.py", line 880, in solveset
    raise ValueError("%s is not a valid SymPy expression" % (f))
  ValueError: Eq(x, -1) | Ne(x**2 + 2*x + 1, 0) is not a valid SymPy expression

none of the obvious ways appear to work. Does Sympy not support this kind of equational reasoning?