As promised in this post! 
It will find all the real roots as long as the coefficients are reasonable (a.k.a. not infinity) and it will sometimes even finds roots when they’re unreasonable. I’ve found that it works pretty well even with polynomials of degree 100. It does not work well with polynomials of degree 200.
Example of an unreasonable polynomial that worked anyway:
0.66853180993545+1.7558181104719 x-0.42926850647086 x^2+0.69164854420848 x^3+1.3273258012241 x^4+0.43469577626489 x^5+1.4850991614606e-006 x^6=0
x=-292701.82949164,-0.34498050186085
The issue here is that the leading coefficient is super small compared to the coefficient of the next degree down, so there is a root way out around -292702. Oh well, the finder found it anyway.
I just expect it to work under normal conditions, like finding the closest point on an ellipsoid to some other point (my next post), or getting eigenvalues and eigenvectors of a matrix.
Woohoo free code
--By xXxMoNkEyMaNxXx
local select=select
local abs=math.abs
local sqrt=math.sqrt
local cos,sin=math.cos,math.sin
local atan2=math.atan2
local sort=table.sort
local insert=table.insert
local remove=table.remove
local unpack=unpack or table.unpack
local sign=function(x)
if x==0 then
return 0
elseif x>0 then
return 1
else
return -1
end
end
local function cpow(re,im,p) --Raise a complex number to an arbitrary real power
local magnitude=(re*re+im*im)^(p/2)
local theta=atan2(im,re)*p
return magnitude*cos(theta),magnitude*sin(theta)
end
local function csqrt(re,im) -- Square root of a complex number with no trig! :>
if im and im~=0 then
local m=sqrt(re*re+im*im)
return sqrt((m+re)/2),re>=0 and (sqrt(m+im)-sqrt(m-im))/2 or im>=0 and (sqrt(m+im)+sqrt(m-im))/2 or (sqrt(m+im)+sqrt(m-im))/-2 --shut up this is the fastest way
else
return re>0 and sqrt(re) or 0,re<0 and sqrt(-re) or 0
end
end
local function cbrt(x) --wow. just wow.
return x<0 and -(-x)^(1/3) or x^(1/3)
end
local acc=1e-10 --floating point error tolerance
local r3=0.75^0.5 -- common constant
local function Zeroes(a0,a1,a2,a3,a4)
local z={}
if not a4 or abs(a4)<acc then
if not a3 or abs(a3)<acc then
if not a2 or abs(a2)<acc then
z[1]=-a0/a1
else
local p1=a1/(-2*a2)
local p2=a1*a1-4*a0*a2
if p2>acc then
local p3=sqrt(p2)/(2*a2)
if p3>0 then
z[1]=p1-p3
z[2]=p1+p3
else
z[1]=p1+p3
z[2]=p1-p3
end
elseif p2>=0 then
z[1]=p1
end
end
else -- No guarantees for sorting from now on
local p1=a2*a2-3*a3*a1
local p2=a2*a2*a2-4.5*a3*a2*a1+13.5*a3*a3*a0
local p3=p2*p2-p1*p1*p1
if abs(p3)<acc then
if abs(p1)<acc then
z[1]=a2/(-3*a3)
else
z[1]=(9*a3*a0-a2*a1)/(2*p1)
z[2]=(4*a2*a1-9*a3*a0-a2*a2*a2/a3)/p1
end
else
local p5r,p5i
local p6r,p6i
if p3<0 then
local p4i=sqrt(-p3)
p5r,p5i=cpow(p2,p4i,1/3)
p6r,p6i=cpow(p2,-p4i,1/3)
else
local p4=sqrt(p3)
p5r,p5i=cbrt(p2+p4),0
p6r,p6i=cbrt(p2-p4),0
end
local z1i,z2i,z3i=(p5i+p6i)/(-3*a3),((p5r-p6r)*r3-(p5i+p6i)/2)/(-3*a3),((p6r-p5r)*r3-(p5i+p6i)/2)/(-3*a3)
if abs(z1i)<acc then
z[#z+1]=(a2+p5r+p6r)/(-3*a3)
end
if abs(z2i)<acc then -- wow I actually have to do this?????? ok FPU
z[#z+1]=(a2-(p5r+p6r)/2+(p6i-p5i)*r3)/(-3*a3)
end
if abs(z3i)<acc then
z[#z+1]=(a2-(p5r+p6r)/2+(p5i-p6i)*r3)/(-3*a3)
end
end
end
else
local p1=a2*a2-3*a3*a1+12*a4*a0
local p2=a2*a2*a2-4.5*a3*a2*a1+13.5*(a3*a3*a0+a4*a1*a1)-36*a4*a2*a0
local p3=2*a2/a4-3*a3*a3/(4*a4*a4)
local p4=(a3*a3*a3-4*a4*a3*a2+8*a4*a4*a1)/(4*a4*a4*a4)
local p5=p2*p2-p1*p1*p1
--print(p1,p2,p3,p4,p5)
local p6
if abs(p1)<acc then
p6=cbrt(2*p2)
elseif p5<0 then
p6=2*(p2*p2-p5)^(1/6)*cos(atan2(sqrt(-p5),p2)/3) -- r u mad
else
local p5r=sqrt(p5)
p6=cbrt(p2+p5r)+cbrt(p2-p5r)
end
local p7=(p6/a4-p3)/3
if p7<-acc then --This part does not seem to ever run, probably needs more imagination
local p7i=sqrt(-p7)
local p8r,p8i=csqrt(-p7-p3,p4*p7i)
--local p9r,p9i=csqrt(-p7-p3,p4*p7i)
if abs(-p8i-p7i)/2<acc then
z[#z+1]=a3/(-4*a4)+p8r/2
end
if abs(p8i-p7i)/2<acc then
z[#z+1]=a3/(-4*a4)-p8r/2
end
if abs(p7i+p8i)/2<acc then
z[#z+1]=a3/(-4*a4)+p8r/2
end
if abs(p7i-p8i)/2<acc then
z[#z+1]=a3/(-4*a4)-p8r/2
end
elseif p7>acc then
local p7r=sqrt(p7)
local p8=p4/p7r-p7-p3
local p9=-p4/p7r-p7-p3
if p8>acc then
local p8r=sqrt(p8)/2
z[#z+1]=a3/(-4*a4)-p7r/2+p8r
z[#z+1]=a3/(-4*a4)-p7r/2-p8r
elseif p8>=0 then
z[#z+1]=a3/(-4*a4)-p7r/2
end
if p9>acc then
local p9r=sqrt(p9)/2
z[#z+1]=a3/(-4*a4)+p7r/2+p9r
z[#z+1]=a3/(-4*a4)+p7r/2-p9r
elseif p9>=0 then
z[#z+1]=a3/(-4*a4)+p7r/2
end
else
local p1=a2/(-2*a4)
local p2=a2*a2-4*a0*a4
if p2>acc then
local p3=sqrt(p2)/(2*a4)
local p4=p1+p3
local p5=p1-p3
if p4>acc then
local p4r=sqrt(p4)
z[#z+1]=p4r
z[#z+1]=-p4r
elseif p4>=0 then
z[#z+1]=0
end
if p5>acc then
local p5r=sqrt(p5)
z[#z+1]=p5r
z[#z+1]=-p5r
elseif p5>=0 then
z[#z+1]=0
end
elseif abs(p2)<acc then
if p1>acc then
local p1r=sqrt(p1)
z[#z+1]=p1r
z[#z+1]=-p1r
elseif p1>=0 then
z[#z+1]=0
end
end
end
end
--sort(z)
return z
end
local function Evaluate(coefficients,x)
local Result=coefficients[0]
local Multiplier=1
for Coefficient=1,#coefficients do
Multiplier=Multiplier*x
Result=Result+coefficients[Coefficient]*Multiplier
end
return Result
end
local function EvaluateWithDerivative(coefficients,x)
local Result,Derivative=coefficients[0],0
local Multiplier=1
for Coefficient=1,#coefficients do
local c=coefficients[Coefficient]
Derivative=Derivative+Coefficient*c*Multiplier
Multiplier=Multiplier*x
Result=Result+c*Multiplier
end
return Result,Derivative
end
--Searches for 0 starting from <value> in the direction of <interval>. Pass unBounded as true if the zero may be out of the interval. This is pretty slow compared to the other functions...
local function BinarySearch(coefficients,value,interval,unBounded)
local LastResult=Evaluate(coefficients,value)
local InitialSign=sign(LastResult)
local Guess=0.5
local Multiplier=0.5
for i=1,75 do
local Result=Evaluate(coefficients,value+Guess*interval)
local ResultSign=sign(Result)
if ResultSign==0 or Result==LastResult then
return value+Guess*interval,Result
else
if unBounded then
Guess=Guess+InitialSign*ResultSign*Multiplier
if InitialSign==ResultSign then
Multiplier=Multiplier*2
else
unBounded=false
end
else
Multiplier=Multiplier/2
Guess=Guess+InitialSign*ResultSign*Multiplier
end
LastResult=Result
end
end
return value+Guess*interval,LastResult
end
--Uses Newton's method
local function IntervalSearch(coefficients,lowerBound,upperBound)
local LowerResult,UpperResult=Evaluate(coefficients,lowerBound),Evaluate(coefficients,upperBound)
local AverageSlope=(UpperResult-LowerResult)/(upperBound-lowerBound)
local LastGuess=lowerBound-LowerResult/AverageSlope
local LastResult,LastDerivative=EvaluateWithDerivative(coefficients,LastGuess)
for i=1,25 do
local Guess=LastGuess-LastResult/LastDerivative
if Guess<=lowerBound or upperBound<=Guess then
Guess=LastGuess-LastResult/AverageSlope -- Should be guaranteed to move the guess towards the right answer...
end
local Result,Derivative=EvaluateWithDerivative(coefficients,Guess)
if Result==0 or Result==LastResult then
return Guess,Result
else
LastGuess,LastResult,LastDerivative=Guess,Result,Derivative
end
end
return BinarySearch(coefficients,lowerBound,upperBound-lowerBound)
end
local function OneSidedSearch(coefficients,bound,side)
local LastGuess=bound+side*abs(Evaluate(coefficients,bound)/coefficients[#coefficients])^(1/#coefficients)
local LastResult,LastDerivative=EvaluateWithDerivative(coefficients,LastGuess)
for i=1,25 do
local Guess=LastGuess-LastResult/LastDerivative
if (Guess-bound)/side<0 then
Guess=LastGuess-LastResult/side -- Not sure if it will be in bounds for certain...
if (Guess-bound)/side<0 then
return BinarySearch(coefficients,bound,side,true)
end
end
local Result,Derivative=EvaluateWithDerivative(coefficients,Guess)
if Result==0 or Result==LastResult then
return Guess,Result
else
LastGuess,LastResult,LastDerivative=Guess,Result,Derivative
end
end
return BinarySearch(coefficients,bound,side,true)
end
local function FindRoots(...) --Found roots are guaranteed to be in ascending order.
local Offset=1
while select(Offset,...)==0 do
Offset=Offset+1
end
local Coefficients={[0]=select(Offset,...),select(Offset+1,...)}
local Degree=#Coefficients
while Coefficients[Degree]==0 do
remove(Coefficients)
Degree=Degree-1
end
local Roots
if Degree>4 then
Roots={}
local DifferentiatedCoefficients={}
for Coefficient=1,Degree do
DifferentiatedCoefficients[Coefficient-1]=Coefficient*Coefficients[Coefficient]
end
local Extremities=FindRoots(unpack(DifferentiatedCoefficients,0))
if #Extremities>0 then
local LeadingCoefficient=Coefficients[Degree]
local LastExtremity=Extremities[1]
local LastExtremityResult=Evaluate(Coefficients,LastExtremity)
if LastExtremityResult==0 then
Roots[#Roots+1]=LastExtremity
elseif LastExtremityResult*LeadingCoefficient*(1-Degree%2*2)<0 then
Roots[#Roots+1]=OneSidedSearch(Coefficients,LastExtremity,-1)
end
for ExtremityIndex=2,#Extremities do
local Extremity=Extremities[ExtremityIndex]
local ExtremityResult=Evaluate(Coefficients,Extremity)
if ExtremityResult==0 then
Roots[#Roots+1]=Extremity
elseif ExtremityResult*LastExtremityResult<0 then
Roots[#Roots+1]=IntervalSearch(Coefficients,LastExtremity,Extremity)
end
LastExtremity=Extremity
LastExtremityResult=ExtremityResult
end
if LastExtremityResult*LeadingCoefficient<0 then
Roots[#Roots+1]=OneSidedSearch(Coefficients,LastExtremity,1)
end
else
local Intercept=Coefficients[0]
if Intercept==0 then
Roots[#Roots+1]=0
else
Roots[#Roots+1]=OneSidedSearch(Coefficients,0,-sign(Intercept)*sign(Coefficients[Degree]))
end
end
elseif Degree>2 then
Roots=Zeroes(unpack(Coefficients,0))
sort(Roots)
elseif Degree>0 then
Roots=Zeroes(unpack(Coefficients,0))
end
if Offset>1 then --Polynomial was degenerate with a root at 0
local BooleanValueThatTellsYouWetherOrNotTheZeroAtZeroHasBeenInsertedIntoTheRootsTable=true
for i=1,#Roots do
if Roots[i]>0 then
insert(Roots,i,0)
BooleanValueThatTellsYouWetherOrNotTheZeroAtZeroHasBeenInsertedIntoTheRootsTable=false
break
end
end
if BooleanValueThatTellsYouWetherOrNotTheZeroAtZeroHasBeenInsertedIntoTheRootsTable then
Roots[#Roots+1]=0
end
end
return Roots
end
return FindRootsMy algorithm is to go from extremity to extremity (by finding the roots of the derivative) and see if two extremities are on opposite sides of 0. If so, there must be a root in between, and the curve is not flat anywhere between the extremities. First I try newton’s method, but if that fails, I use a binary search. There are special cases for multiplicities and roots that don’t lie between two extremities.
I made it into a module because modules are nice. I hope someone actually uses this because it was fun to make ._.