Is there any way to get the step-by-step solution in SymPy?

Question

Is there any way to get the step-by-step solution in SymPy? For example:

x**2-5 = 4
  step 1 x**2-5+5=4+5
  step 2 : x**2=9
  step 3 :x = 3 or x= -3

Answer 1

(this is more a comment as answer)

There are some Google's-soc ideas for implementing

Step-by-step Expression Visualization

GSoC 2014 Ideas: Many times, people ask how they can tell what some functions are doing. For example, they want to know step by step...... For the former, the best you can do is to follow the code; for the latter, the algorithm doesn't work at all like you would do it by hand, so there's really no way...

GSoC 2015 Ideas:

A Strategy for step-by-step

The logic behind many SymPy operations is separated into several small methods. For example objects like sin or exp have _eval_derivative methods that are called as SymPy evaluates the derivative of a complex expression like sin(exp(x)). By capturing the inputs and outputs of all of these small methods we can collect a great quantity of information about the steps that SymPy takes. We can see that exp._eval_derivative took in exp(x) and returned exp(x) and that sin._eval_derivative took in sin(exp(x)) and returned cos(exp(x))*exp(x). These input-output pairs for each method are probably sufficient to illustrate how SymPy solves problems in many domains.

This approach of capturing the inputs of many internal functions is similar to logging systems traditionally used to analyze large codebases. We should investigate how they work and if they cause any problems with normal operation.

Once this source of information is available we can then think about interesting ways to visualize and to interact with it. A good solution will not irrevocably tie the data stream to a particular visualization technique.

This approach is straightforward intellectually but may require the student to interact with a lot of the codebase. Approaches like _eval_derivative are ubiquitous throughout SymPy but often have small variations in different modules.

here a online solution SymPy Gamma

Answer 2

A temporary solution for simple cases might be based on expession trees. At the beginning you can create an expession tree of expression with no evaluation. After that you can modify results level by level:

Source code:

class TraverseSolver:
    def __init__(self, expr):
        self.expr = expr
        
    def _set_graph(self):
        self.G = nx.nx_agraph.from_agraph(pygraphviz.AGraph(sp.dotprint(self.expr))) 
        
    def _set_map(self):
        self._map = dict(zip(self.G.nodes, sp.preorder_traversal(self.expr))) 
        
    def _set_baseNode(self):
        self._baseNode = next(iter(self.G.nodes))
        
    def get_levels(self, mode='draw'):
        if mode == 'draw':
            d = nx.single_source_shortest_path_length(self.G, self._baseNode)
            u, idx = np.unique(list(d.values()), return_index=True)
            levels = [[str(m) for m in n] for n in reversed(np.split(np.array(list(d.keys())), idx[1:]))]
            return levels
        elif mode == 'traverse':
            print(self.G)
    
    def set_color(self, node, color):
        self.G.nodes[node]['color'] = color

    def display_graph(self, fig, n, nshape=(2, 3)):      
        ax = fig.add_subplot(*nshape, n)
        pos = graphviz_layout(self.G, prog='dot')
        colors = nx.get_node_attributes(self.G, 'color')    
        nx.draw(self.G, pos = pos, nodelist=[])
        # draw self.G bbox by bbox:
        for i, n in enumerate(self.G.nodes()):
            nx.draw(nx.subgraph(self.G, [n]), pos={n:pos[n]}, labels = {n:f'${sp.latex(self._map[n])}$'}, nodelist=[],
                    bbox=dict(facecolor=colors[n], edgecolor='black', boxstyle='round,pad=0.7'))
            
    def solve(self, display_graph=True, nshape=(2, 3)):
        self._set_graph() #store sp.srepr+code in each node
        self._set_map() #sp.srepr+code -> expression (without evaluation)
        self._set_baseNode() #sp.srepr+code of self.
        solutionSteps = [self._map[self._baseNode]] #first step that contains initial expression
        levels = self.get_levels(mode='draw')
        if display_graph:
            fig = plt.figure(figsize=(20,10))
        #Step forward
        for i in range(len(levels)):
            if display_graph:
                for node in self.G.nodes(): 
                    self.set_color(node, 'lightblue')
            anyChanges = False
            for activeNode in levels[i]:
                beforeEval = self._map[activeNode]
                if display_graph:
                    self.set_color(activeNode, 'yellow')
                if not beforeEval.is_Atom:
                    afterEval = beforeEval.func(*beforeEval.args, evaluate=True) #is beforeEval different with afterEval
                    if beforeEval != afterEval:
                        self._map[activeNode] = afterEval
                        if display_graph:
                            self.set_color(activeNode, 'lime')
                        anyChanges = True
            # Calculate value of baseNode() using changes, no evaluation      
            if anyChanges:
                for j in range(i+1, len(levels)):
                    for editNode in levels[j]:
                        args = [self._map[node] for node in self.G[editNode]] #each ancestor
                        if not self._map[editNode].is_Atom:
                            self._map[editNode] = self._map[editNode].func(*args, evaluate=False)
                solutionSteps.append(self._map[self._baseNode])
            if display_graph:
                self.display_graph(fig, n=len(solutionSteps), nshape=nshape)
        plt.show()
        return solutionSteps

expr = sp.sympify('-1*(2*3-5*7)', evaluate=False)
steps = TraverseSolver(expr).solve(display_graph=True, nshape=(2, 3))

print('INPUT:', sp.srepr(expr))
print('SOLUTION 1:', ' = '. join([str(step) for step in steps]))
print('SOLUTION 2:', ' = '. join([sp.StrPrinter(dict(order='none'))._print(step) for step in steps])) 

Output:

INPUT: Mul(Integer(-1), Add(Mul(Integer(-1), Mul(Integer(5), Integer(7))), Mul(Integer(2), Integer(3))))
SOLUTION 1: -(-5*7 + 2*3) = -(-1*35 + 2*3) = -(-35 + 6) = -1*(-29) = 29
SOLUTION 2: -(2*3 - 5*7) = -(2*3 - 1*35) = -(6 - 35) = -1*(-29) = 29

Requirements: networkx, matplotlib, python-graphviz

Remark: this is an initial version of my research, I'm planning to refactor this dirty code quite soon...