Consider adding `Symbol.__str__()` overload

Currently we have `Symbol.__repr__()` implemented
```python
In [1]: from dwave.optimization import Model

In [2]: model = Model()

In [3]: i = model.integer()

In [5]: repr(i)
Out[5]: ''
```

However, we've had a few requests for something closer to `dimod`'s CQM printing
```python
In [6]: import dimod

In [7]: cqm = dimod.ConstrainedQuadraticModel()

In [8]: i = cqm.add_variable("INTEGER")

In [9]: j = cqm.add_variable("INTEGER")

In [17]: cqm.set_objective([(i, j, 1)])

In [18]: print(cqm)
Constrained quadratic model: 2 variables, 0 constraints, 3 biases

Objective
Integer(1)*Integer(0)

Constraints

Bounds
0.0 <= Integer(0) <= 9007199254740991.0
0.0 <= Integer(1) <= 9007199254740991.0
```

**Approaches Considered**

I think don't want to print everything for every node, so maybe something like 1-nesting
```python
Max(Sum(...), Constant(...))
```
or
```python
Max(Sum, Constant)
```
or something like that.

We could also take inspiration from [sympy](https://www.sympy.org) which does print the whole thing
```python
In [1]: from sympy import *

In [2]: x, y, z = symbols("x y z")

In [10]: out = x + y

In [11]: for _ in range(100):
...: out = x*out + y
...:

In [12]: out
Out[12]: x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x*(x + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y) + y)
```
Though how to handle say indexing gets a bit funny.