Arbitrary Degree Polynomial Root Finder (real numbers only)

--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 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 dumpPoly(coefficients)
    local s=""
    for i=0,#coefficients do
        local v=coefficients[i]
        if v~=0 then
            s=s..((v<0 or #s==0) and "" or "+")..v..(i==0 and "" or i==1 and " x" or " x^"..(i>9 and "{"..i.."}" or i))
        end
    end
    return (s:gsub("e(%-?%d%d%d)"," 10^{%1}"))
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

local flag=false
--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
    print("What the hell?",unBounded,value,unBounded and interval or value+interval,value+Guess*interval,LastResult,dumpPoly(coefficients))
    flag=true
    return value+Guess*interval,LastResult
end

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
    --print("Only the worst polynomials have to use this code.",#coefficients)
    return BinarySearch(coefficients,lowerBound,upperBound-lowerBound)
end

--[[ It was such a good idea though :c
--side = 1:Integrate[D[a x^n, x] Exp[x0 - x], {x, x0, Infinity}]/Integrate[Exp[x0 - x], {x, x0, Infinity}]
--side =-1:Integrate[D[a x^n, x] Exp[x - x0], {x,-Infinity, x0}]/Integrate[Exp[x - x0], {x,-Infinity, x0}]
local function GetWeightedSlopeIntegral(coefficients,x0,side)
    local WeightedSlope=coefficients[1]
    local GammaCoefficients={1}
    for d=2,#coefficients do
        local GammaResult=0
        local Multiplier=1
        local d1=d-1
        for i=1,d1 do
            local gci=GammaCoefficients[i]*d1*side
            GammaCoefficients[i]=gci
            GammaResult=GammaResult+gci*Multiplier
            Multiplier=Multiplier*x0
        end
        GammaCoefficients[d]=1
        WeightedSlope=WeightedSlope+d*coefficients[d]*(GammaResult+Multiplier)
    end
    return WeightedSlope
end

local function OneSidedSearch(coefficients,bound,side)
    local BoundResult=Evaluate(coefficients,bound)
    local WeightedSlope=GetWeightedSlopeIntegral(coefficients,bound,side)
    local LastGuess=bound-BoundResult/WeightedSlope
    local LastResult,LastDerivative=EvaluateWithDerivative(coefficients,LastGuess)
    for i=1,30 do
        local Guess=LastGuess-LastResult/LastDerivative
        if (Guess-bound)/side<0 then
            Guess=LastGuess-LastResult/WeightedSlope -- Not sure if it will be in bounds for certain...
            if (Guess-bound)/side<0 then
                print("Your polynomial is truly ugly.")
                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 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
                --print("Your polynomial is truly ugly.",#coefficients)
                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
    --print("Your polynomial sucks.",#coefficients)
    return BinarySearch(coefficients,bound,side,true)
end

local function FindRoots(...) --Found roots are guaranteed to be in ascending order.
    local Coefficients={[0]=(...),select(2,...)}
    local Degree=#Coefficients
    while Coefficients[Degree]==0 do
        remove(Coefficients)
        Degree=Degree-1
    end
    if Degree>4 then
        local DifferentiatedCoefficients={}
        for Coefficient=1,Degree do
            DifferentiatedCoefficients[Coefficient-1]=Coefficient*Coefficients[Coefficient]
        end
        local Roots={}
        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--if Degree%2==1 then -- Degree cannot be even if there are no extremities - this test is useless
            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
        return Roots
    elseif Degree>2 then
        local Roots=Zeroes(...)
        sort(Roots)
        return Roots
    elseif Degree>0 then
        return Zeroes(...)
    end
end

--[[
do
    local c={}
    local f=1
    for i=1,100,2 do
        c[i]=f
        c[i-1]=0
        f=-f/((i+1)*(i+2))
    end
    print(dumpPoly(c))
    print(unpack(FindRoots(unpack(c,0))))
end

print(unpack(FindRoots(0,0,3,-2)))
--]]

--[[ sad time testing -- took a long time to find little mistakes like forgot to multiply by 2
print(unpack(Zeroes(6,5,4,3,-2)))
print(unpack(Zeroes(0.5,5,-6,-1,2)))
print(unpack(Zeroes(2,0,-3,0,3)))
print(unpack(Zeroes(4.1,-1,-4,-4.6,5)))
print(unpack(Zeroes(-1,0,0,0,1)))
print(unpack(Zeroes(1.9208,-6.86,8.82,-4.9,1)))
print(unpack(Zeroes(0.5,5,-2.8,-3.5,2)))
print(unpack(Zeroes(6,5,4,3)))
print(unpack(Zeroes(5,4,3)))
print(unpack(Zeroes(4,3)))
--]]

local tau=2*math.pi
local log=math.log

math.randomseed(os.time())
math.random()

local rand=math.random

local s
local function r()
    if s then
        local r=s
        s=nil
        return r
    else
        local theta=tau*rand()
        local r2
        repeat r2=rand() until r2~=0
        local magnitude=log(r2)
        local r
        r,s=magnitude*cos(theta),magnitude*sin(theta)
        return r
    end
end

--[[
do
    local c={}
    local s=""
    for i=0,100 do
        local v=r()
        c[i]=v
        s=s..((v<0 or i==0) and "" or "+")..v..(i==0 and "" or i==1 and " x" or " x^"..(i>9 and "{"..i.."}" or i))
    end
    print(s)
    print(unpack(FindRoots(unpack(c,0))))
end
--]]

--[[
for d=1,100 do
    local c={}
    local s=""
    for i=0,d do
        local v=r()
        c[i]=v
        s=s..((v<0 or i==0) and "" or "+")..v..(i==0 and "" or i==1 and " x" or " x^"..(i>9 and "{"..i.."}" or i))
    end
    print(s)
    print(unpack(FindRoots(unpack(c,0))))
end
--]]

---[[
local n=100000
local n2=1000000
local rtime
do
    local t0=os.clock()
    for i=1,n2 do
        r()
    end
    local t1=os.clock()
    rtime=(t1-t0)/n2
    print(rtime)
end
for d=1,10 do
    local t0=os.clock()
    local q=1+math.floor(n*rand())
    for _=1,n do
        local c={}
        for i=0,d do
            c[i]=r()
        end
        local Roots=FindRoots(unpack(c,0))
        if flag or _==q then
            print(dumpPoly(c))
            print(unpack(Roots))
            flag=false
        end
    end
    local t1=os.clock()
    print(d,1/((t1-t0)/n-(d+1)*rtime))
end
--]]

--[[ Extra cases that were handled by a more general block of code
elseif p3>0 then
    local p4=sqrt(p3)/2
    local p5=p2+p4
    local p6=p2-p4
    local c1,c2=p5>=0,p6>=0
    if c1 and c2 then
        --uhhh they'll always be imaginary... ok
        local p7,p8=p5^(1/3),p6^(1/3)
        z[1]=(a2+p7+p8)/(-3*a3)
        print(r3*(p7-p8)/(-3*a3),r3*(p8-p7)/(-3*a3))--imaginary parts should not be zero
    elseif c2 then
        local p7r,p7i=cpow(p5,0,1/3)
        local p8=p6^(1/3)
        local z1r,z1i=(a2+p7r+p8)/(-3*a3),p7i/(-3*a3) -- not zero o.o
        local z2r,z2i=(a2-(p7r+p8)/2-p7i*r3)/(-3*a3),((p7r-p8)*r3-p7i/2)/(-3*a3)
        local z3r,z3i=(a2-(p7r+p8)/2+p7i*r3)/(-3*a3),((p8-p7r)*r3-p7i/2)/(-3*a3)
        print(z1i,z2i,z3i)
    elseif c1 then
        local p7=p5^(1/3)
        local p8r,p8i=cpow(p5,0,1/3)
        local z1r,z1i=(a2+p7+p8r)/(-3*a3),p8i/(-3*a3)
    else
        --
    end
    local p4=sqrt(p3)/2
    local p5r,p5i=cpow(p2+p4,0,1/3)
    local p6r,p6i=cpow(p2-p4,0,1/3)
    z[1]=(a2+p5r+p6r)/(-3*a3)
    z[2]=(a2-(p5r+p6r)/2+(p6i-p5i)*r3)/(-3*a3)
    z[3]=(a2-(p5r+p6r)/2+(p5i-p6i)*r3)/(-3*a3)
    print((p5i+p6i)/(-3*a3),((p5r-p6r)*r3-(p5i+p6i)/2)/(-3*a3),((p6r-p5r)*r3-(p5i+p6i)/2)/(-3*a3))
--]]