# Lambert Solver Of Sympy (GSoC 2019 - Week 7)

15 Jul 2019This was the seventh week meeting with the GSoC mentors which was scheduled on Monday 15th July, 2019 between 6:00 - 7:00 PM (IST). Me and Yathartha were the attendees of the meeting. In this blog I will be describing the lambert equation solver for Sympy and what problems it faced before and how PR #16890 will solve the problems. It is preassumed that you know what lambert type equations are, so I will not be explaining that.

## Explaining the function `_solve_lambert`

(main function to solve lambert equations)

```
Input - f, symbol, gens
OutPut - Solution of f = 0 if its lambert type expression else NotImplementedError
```

This function separates out cases as below based on the main function present in the main equation.

```
For the first ones:
1a1) B**B = R != 0 (when 0, there is only a solution if the base is 0,
but if it is, the exp is 0 and 0**0=1
comes back as B*log(B) = log(R)
1a2) B*(a + b*log(B))**p = R or with monomial expanded or with whole
thing expanded comes back unchanged
log(B) + p*log(a + b*log(B)) = log(R)
lhs is Mul:
expand log of both sides to give:
log(B) + log(log(B)) = log(log(R))
1b) d*log(a*B + b) + c*B = R
lhs is Add:
isolate c*B and expand log of both sides:
log(c) + log(B) = log(R - d*log(a*B + b))
```

If the equation are of type 1a1, 1a2 and 1b then the mainlog of the equation is taken into concern as the deciding factor lies in the main logarithmic term of equation.

```
For the next two,
collect on main exp
2a) (b*B + c)*exp(d*B + g) = R
lhs is mul:
log to give
log(b*B + c) + d*B = log(R) - g
2b) -b*B + g*exp(d*B + h) = R
lhs is add:
add b*B
log and rearrange
log(R + b*B) - d*B = log(g) + h
```

If the equation are of type 2a and 2b then the mainexp of the equation is taken into concern as the deciding factor lies in the main exponential term of equation.

```
3) d*p**(a*B + b) + c*B = R
collect on main pow
log(R - c*B) - a*B*log(p) = log(d) + b*log(p)
```

If the equation are of type 3 then the mainpow of the equation is taken into concern as the deciding factor lies in the main power term of equation.

Eventually from all of the three cases the equation is meant to be converted to this form:-

```
f(x, a..f) = a*log(b*X + c) + d*X - f = 0 which has the
solution, X = -c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(f/a)).
```

And the solution calculation process is done by `_lambert`

function.

Everything seems flawless?? You might be thinking no modification is required. Lets see what loopholes are there in it.

## What does PR #16890 do?

There are basically two flaws present with the this approach.

- Not considering all branches of equation while taking log both sides.
- Calculation of roots should consider all roots in case having rational power.

### 1. Not considering all branches of equation while taking log both sides.

Let us consider this equation to be solved by `_solve_lambert`

function.

```
-1/x**2 + exp(x/2)/2 = 0
```

So what the old `_solve_lambert`

do is to convert this equation to following.

```
2*log(x) + x/2 = 0
```

and calculates its roots from `_lambert`

.
But it missed this branch of equation while taking log on main equation.

```
2*log(-x) + x/2 = 0
```

Yeah you can reproduce the original equation from this equation.So basically the problem was that it missed the branches of equation while taking log. And when does the main equation have more than one branch?? The terms having even powers of variable x leads to two different branches of equation.

So how it is solved?
What I has done is that before actually gets into solving I preprocess the main equation
and if it has more than one branches of equation while converting taking log then I consider
all the equations generated from them.(with the help of `_solve_even_degree_expr`

)

How I preprocess the equation? So what I do is I replace all the even powers of x present with even powers of t(dummy variable).

```
Code for targeted replacement
lhs = lhs.replace(
lambda i: # find symbol**even
i.is_Pow and i.base == symbol and i.exp.is_even,
lambda i: # replace t**even
t**i.exp)
Example:-
Main equation -> -1/x**2 + exp(x/2)/2 = 0
After replacement -> -1/t**2 + exp(x/2)/2 = 0
```

Now I take logarithms on both sides and simplify it.

```
After simplifying -> 2*log(t) + x/2 = 0
```

Now I call function `_solve_even_degree_expr`

to replace the t with +/-x to generate two equations.

```
Replacing t with +/-x
1. 2*log(x) + x/2 = 0
2. 2*log(-x) + x/2 = 0
```

And consider the solutions of both of the equations to return all lambert real solutions
of `-1/x**2 + exp(x/2)/2 = 0`

.

Hope you could understand the logic behind this work.

### 2. Calculation of roots should consider all roots in case having rational power.

This flaw is in the calculation of roots in function `_lambert`

.
Earlier the function_lambert has the working like :-

- Find all the values of a, b, c, d, e in the required loagrithmic equation
- Then it defines a solution of the form
`-c/b + (a/d)*l where l = LambertW(d/(a*b)*exp(c*d/a/b)*exp(-f/a), k)`

and then it included that solution. I agree everything seems flawless here. but try to see the step where we are defining l.

Let us suppose a hypothetical algorithm just like algorithm used in `_lambert`

in which equation to be solved is

```
x**3 - 1 = 0
```

and in which we define solution of the form

```
x = exp(I*2*pi/n) where n is the power of x in equation
```

so the algorithm will give solution

```
x = exp(I*2*pi/3) # but expected was [1, exp(I*2*pi/3), exp(-I*2*pi/3)]
```

which can be found by finding all solutions of

```
x**n - exp(2*I*pi) = 0
```

by a different correct algorithm. Thats y it was wrong.
The above algorithm would have given correct values for `x - 1 = 0`

.

And the question in your mind may arise that why only exp() because the
possiblity of having more than one roots is in exp(), because if the algorithm
would have been like `x = a`

, where a is some real constant then there is not
any possiblity of further roots rather than solution like `x = a**(1/n)`

.
And its been done in code like this:

```
code
num, den = ((c*d-b*f)/a/b).as_numer_denom()
p, den = den.as_coeff_Mul()
e = exp(num/den)
t = Dummy('t')
args = [d/(a*b)*t for t in roots(t**p - e, t).keys()]
```

Thank you! Thats all it was!!.

Follow @jmig5776