Way to find positive roots by treating complex numbers as negatives? [SymEngine + Roots]

djsegal · June 15, 2018, 7:25am

Let’s say you had some hideous nonlinear equation like so:

cur_equation = -1.0 + 32.5*(0.237334965250543*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) + 8.02921230974917e-05*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258) + 2.38498439733623e-10*((I_P^0.96)^9.6774193548387*(0.907736861673663*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)*(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)) + (0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 0.82398621004115*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)^2)*(1.2 + 994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 9.67240598811777e-05*(I_P^0.96)^3.2258064516129))^0.5)/(I_P*(0.771372182862621*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) - 0.000133194034060396*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258)))

However the value of I_P stands for a real, positive quantity – i.e. a current.

Because of the funkiness in the equation, when you try to get roots from:

Roots.find_zeros(cur_equation, 0.1, 100.0)

You get complex solutions and the root solver breaks (see the following error):

MethodError: no method matching isless(::Complex{Float64}, ::Int64)
Closest candidates are:
  isless(::Missings.Missing, ::Any) at /Users/dan/.julia/v0.6/Missings/src/Missings.jl:74
  isless(::AbstractFloat, ::Real) at operators.jl:98
  isless(::ForwardDiff.Dual{Tx,V,N} where N where V<:Real, ::Integer) where Tx at /Users/dan/.julia/v0.6/ForwardDiff/src/dual.jl:104

The way I’ve worked around this is doing some lambdify-foo a la:

cur_lambda = lambdify(cur_equation)

cur_func = function (work_I_P)
  cur_value = cur_lambda(complex(float(work_I_P)))
  iszero(imag(cur_value)) || return float(-100.0)
  return Real(cur_value)
end

cur_I_P_list = Roots.find_zeros(cur_equation, 0.1, 100.0)

And this has rewarded me with:

cur_I_P_list == [22.2564, 27.8827, 29.8268]

The problem is that this is really slow.

Has anyone had experience with a similar problem?

How could this be sped up?

// I don’t have a second per root solve

djsegal · June 15, 2018, 7:43am

For reference, here is an image of the function with the complex values replaced with the value -5

The complex values happen at low values of current. (and are representative of impossible to achieve plasmas)

Tamas_Papp · June 15, 2018, 8:28am

I would consider introducing some

\zeta = (I_P^{0.96})^{3.2258064516129}

then the equation seems to be in integer powers (particularly -1, 1, 2, there may be others I missed) of \zeta.

EDIT: fixed mistake pointed out by @mzaffalon.

mzaffalon · June 15, 2018, 11:22am

Not all I_P have the numerical factor 0.96.

mzaffalon · June 15, 2018, 11:24am

Except for the numerical factor, @Tamas_Papp’ suggestion is a sound one in that a change of variables may help. Indeed there are recurrent factors (I_P^{3.09677419354838})^{(0.3226I3.09677419354838P+6363.67829013947)}.

And do you really need to keep all the terms?

Tamas_Papp · June 15, 2018, 11:25am

I now see that it was actually an exponent, corrected it, thanks for pointing it out.

dpsanders · June 15, 2018, 12:52pm

Here is a session using IntervalRootFinding.jl on the original function, although I agree that it would be better to simplify the expression.

f(I_P) = -1.0 + 32.5*(0.237334965250543*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) + 8.02921230974917e-05*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258) + 2.38498439733623e-10*((I_P^0.96)^9.6774193548387*(0.907736861673663*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)*(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)) + (0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 0.82398621004115*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)^2)*(1.2 + 994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 9.67240598811777e-05*(I_P^0.96)^3.2258064516129))^0.5)/(I_P*(0.771372182862621*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) - 0.000133194034060396*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258)))


using IntervalArithmetic, IntervalRootFinding

rts = roots(f, -100..101, Bisection, 1e-1)

julia> rts
19-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([29.8768, 29.9757], :unknown)
 Root([29.7779, 29.8769], :unknown)
 Root([29.6805, 29.778], :unknown)
 Root([27.9235, 28.0195], :unknown)
 Root([27.8277, 27.9236], :unknown)
 Root([27.7334, 27.8278], :unknown)
 Root([22.633, 22.732], :unknown)
 Root([22.4337, 22.5327], :unknown)
 Root([22.3348, 22.4338], :unknown)
 Root([22.2374, 22.3349], :unknown)
 Root([22.1868, 22.2375], :unknown)
 Root([22.137, 22.1869], :unknown)
 Root([22.0381, 22.1371], :unknown)
 Root([21.9392, 22.0382], :unknown)
 Root([21.8419, 21.9393], :unknown)
 Root([21.743, 21.842], :unknown)
 Root([19.2784, 19.3759], :unknown)
 Root([19.1811, 19.2785], :unknown)
 Root([19.0852, 19.1812], :unknown)


julia> rts = roots(f, rts, Bisection, 1e-2)
20-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([29.833, 29.8393], :unknown)
 Root([29.8269, 29.8331], :unknown)
 Root([29.8207, 29.827], :unknown)
 Root([29.8145, 29.8208], :unknown)
 Root([27.8871, 27.8932], :unknown)
 Root([27.8811, 27.8872], :unknown)
 Root([27.8752, 27.8812], :unknown)
 Root([27.8692, 27.8753], :unknown)
 Root([22.2796, 22.2858], :unknown)
 Root([22.2735, 22.2797], :unknown)
 Root([22.2674, 22.2736], :unknown)
 Root([22.2614, 22.2675], :unknown)
 Root([22.2553, 22.2615], :unknown)
 Root([22.2493, 22.2554], :unknown)
 Root([22.2433, 22.2494], :unknown)
 Root([22.2374, 22.2434], :unknown)
 Root([22.231, 22.2375], :unknown)
 Root([22.2246, 22.2311], :unknown)
 Root([19.1811, 19.1871], :unknown)
 Root([19.1749, 19.1812], :unknown)

julia> roots(f, rts, Newton)
4-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([29.8267, 29.8269], :unique)
 Root([27.8826, 27.8828], :unique)
 Root([22.2562, 22.2566], :unique)
 Root([19.1825, 19.1833], :unknown)

dpsanders · June 15, 2018, 12:55pm

julia> roots(f, rts, Newton, 1e-10)
3-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([29.8267, 29.8268], :unique)
 Root([27.8826, 27.8827], :unique)
 Root([22.2564, 22.2565], :unique)

So I agree with the Roots package that these are the only 3 real roots of the equation in the interval -100…100.

j_verzani · June 15, 2018, 2:03pm

If you replace log with safer_log(x) = x <= 0 ? NaN : log(x) and make the whole thing a function of I_P and not a symbolic expression it runs through find_zeros just fine and is reasonably fast. However, converting with lambdify makes this trick harder. If you had to, you could get the expression to lambdify with convert(Expr, cur_equation), splice in such a substitute, then call find_zeros. Using MacroTools, something like this seems to work, which might help in general, but isn’t very friendly or speedy:

u = convert(Expr, cur_equation)
ex = MacroTools.postwalk(x -> @capture(x, log(xs__)) ? :(safer_log($(xs...))) : x, u) # slowish
fn = eval(Expr(:function, Expr(:call, gensym(), :(I_P)), ex))
find_zeros(fn, 1.0, 100.0)

djsegal · June 15, 2018, 11:19pm

After trying to do,

using IntervalArithmetic, IntervalRootFinding

function test_sanders()
    f(I_P) = -1.0 + 32.5*(0.237334965250543*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) + 8.02921230974917e-05*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258) + 2.38498439733623e-10*((I_P^0.96)^9.6774193548387*(0.907736861673663*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)*(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)) + (0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 0.82398621004115*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)^2)*(1.2 + 994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 9.67240598811777e-05*(I_P^0.96)^3.2258064516129))^0.5)/(I_P*(0.771372182862621*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) - 0.000133194034060396*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258)))

    @time rts = roots(f, -100..101, Bisection, 1e-1)
    @time rts = roots(f, rts, Bisection, 1e-2)
    @time roots(f, rts, Newton)
end

test_sanders()

I get,

  6.549414 seconds (21.20 M allocations: 797.128 MiB, 6.75% gc time)
  5.895622 seconds (20.01 M allocations: 752.655 MiB, 6.61% gc time)
DomainError:
log will only return a complex result if called with a complex argument. Try log(complex(x)).

Stacktrace:
 [1] nan_dom_err at ./math.jl:300 [inlined]
 [2] log at ./math.jl:419 [inlined]
 [3] (::#f#1)(::Float64) at ./In[1]:4
 [4] guarded_mid(::#f#1, ::IntervalArithmetic.Interval{Float64}) at /Users/dan/.julia/v0.6/IntervalRootFinding/src/IntervalRootFinding.jl:48
 [5] N at /Users/dan/.julia/v0.6/IntervalRootFinding/src/newton.jl:37 [inlined]

Is this the fastest I can make it to catch all three roots?
Why does your code sample not trigger the DomainError for log calls?

djsegal · June 15, 2018, 11:30pm

Why would you do this if it was slower?

Just rehashing this quick, this is simple setup:

using SymEngine, Roots, MacroTools

I_P = SymEngine.symbols("I_P")
cur_equation = -1.0 + 32.5*(0.237334965250543*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) + 8.02921230974917e-05*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258) + 2.38498439733623e-10*((I_P^0.96)^9.6774193548387*(0.907736861673663*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)*(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)) + (0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 0.82398621004115*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)^2)*(1.2 + 994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 9.67240598811777e-05*(I_P^0.96)^3.2258064516129))^0.5)/(I_P*(0.771372182862621*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) - 0.000133194034060396*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258)))

This is your suggestion:

safer_log(x) = x <= 0 ? NaN : log(x)
u = convert(Expr, cur_equation)

@time ex = MacroTools.postwalk(x -> @capture(x, log(xs__)) ? :(safer_log($(xs...))) : x, u) # slowish
@time fn = eval(Expr(:function, Expr(:call, gensym(), :(I_P)), ex)) 
@time Roots.find_zeros(fn, 0.01, 100.0)

-------
  2.028164 seconds (102.16 k allocations: 4.982 MiB)
  0.010968 seconds (1.15 k allocations: 86.956 KiB)
  0.153697 seconds (172.94 k allocations: 7.850 MiB)

And this is the existing code from the initial question,

@time cur_lambda = lambdify(cur_equation)
@time cur_func = function (work_I_P)
    cur_value = cur_lambda(complex(float(work_I_P)))
    iszero(imag(cur_value)) || return float(-100.0)
    return Real(cur_value)
end
@time Roots.find_zeros(cur_func, 0.01, 100.0)

-------
  0.020340 seconds (22.95 k allocations: 831.798 KiB)
  0.000018 seconds (19 allocations: 1.483 KiB)
  1.060586 seconds (1.24 M allocations: 40.453 MiB, 1.03% gc time)

djsegal · June 15, 2018, 11:39pm

The reason for this factor is because we substitute expressions for radius and magnetic field into the current equation (through the G_blank composite functions).

14%20PM

Where R and B equations come in pairs for various constraints:

20%20PM

Therefore simplification would have to be done either manually for every constraint or digitally – which can take a few second.

Also, if you made your substitution, there is still a factor of I_P floating around that will have a fractional exponent.

// i.e. the one that is currently sitting to a power of one:

... /(I_P*(0.771372182862621* ...

So I guess part of the question: is having one fractional exponent any better than having a bunch?

djsegal · June 17, 2018, 6:19am

Follow-up, is there a way to handle exponentiation with this trick to?

DomainError:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

And the next one after that would be negative exponents // to complete the set

(basically can we throw away type stability for symbolic root solving functions)

edit: is lambdify what’s slow? (b/c of the invokelatest call)

j_verzani · June 17, 2018, 1:11pm

There is overhead with lambdify, about twice the time for smaller cases, but in my example the slowness comes from walking the expression tree. You can speed it up by converting to a string, replacing the calls as desired, then converting back to an expr. As for complex numbers, many of the algorithms in find_zero can use complex numbers, but none of those that use bracketing (as find_zeros does).

djsegal · June 19, 2018, 12:42pm

For the interested, a robust find_zeros for complex numbers was proposed in:

https://github.com/JuliaMath/Roots.jl/pull/108

// it’s able to find all the roots I need in around 1 second !!!