local function IsZero(a)
	return a<0.0001
end

function SolveQuadric(c0, c1, c2)
	    local s0 = 0
        local s1 = 0

        local p, q, D;

        p = c1 / (2 * c0);
        q = c2 / c0;

        D = p * p - q;

        if (IsZero(D)) then
	        s0 = -p;
	        return 1, s0, s1;
        elseif (D < 0) then
	        return 0, s0, s1
        else
	        local sqrt_D = math.sqrt(D);

	        s0 =   sqrt_D - p;
	        s1 = -sqrt_D - p;
	        return 2, s0, s1
        end
end

function SolveCubic(c0, c1, c2, c3)
	    local s0 = 0
        local s1 = 0
        local s2 = 0

        local     num;
        local  sub;
        local  A, B, C;
        local  sq_A, p, q;
        local  cb_p, D;

        --/* normal form: x^3 + Ax^2 + Bx + C = 0 */
        A = c1 / c0;
        B = c2 / c0;
        C = c3 / c0;

       -- /*  substitute x = y - A/3 to eliminate quadric term:  x^3 +px + q = 0 */
        sq_A = A * A;
        p = 1.0/3 * (- 1.0/3 * sq_A + B);
        q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);

        --/* use Cardano's formula */
        cb_p = p * p * p;
        D = q * q + cb_p;

        if (IsZero(D)) then
	        if (IsZero(q)) then
	            s0 = 0;
	            num = 1;
	        else
	            local u = math.pow(-q, 1.0/3.0);
	            s0 = 2 * u;
	            s1 = - u;
	            num = 2;
	        end
        elseif (D < 0) then
	        local phi = 1.0/3 * math.acos(-q / math.sqrt(-cb_p))
	        local t = 2 * math.sqrt(-p);

	        s0 =   t * math.cos(phi);
	        s1 = - t * math.cos.Cos(phi + math.pi / 3)
	        s2 = - t * math.cos.Cos(phi - math.pi / 3)
	        num = 3;
        else --/* one real solution */ {
	        local sqrt_D = math.sqrt(D);
	        local u = math.pow(sqrt_D - q, 1.0/3.0)
	        local v = - math.pow(sqrt_D + q, 1.0/3.0)

	        s0 = u + v;
	        num = 1;
        end

        sub = 1.0/3 * A;

        if (num > 0) then    s0 = s0-sub end
		if (num > 1) then    s1 = s1-sub end
		if (num > 2) then    s2 = s2-sub end

        return num, s0, s1, s2
end

function SolveQuartic(c0,c1,c2,c3,c4)
	 	local s0 = 0
        local s1 = 0
        local s2 = 0
        local s3 = 0

        local coeffs = {}
        local  z, u, v, sub
        local  A, B, C, D
        local  sq_A, p, q, r
        local     num;

        local A = c1 / c0;
        local B = c2 / c0;
        local C = c3 / c0;
        local D = c4 / c0;

        local sq_A = A * A;
        local p = - 3.0/8 * sq_A + B;
        local q = 1.0/8 * sq_A * A - 1.0/2 * A * B + C;
        local r = - 3.0/256*sq_A*sq_A + 1.0/16*sq_A*B - 1.0/4*A*C + D;

        if (r < 0.001) then
	        coeffs[3] = q;
	        coeffs[2] = p;
	        coeffs[1] = 0;
	        coeffs[0] = 1;

	        num, s0, s1, s3 = SolveCubic(coeffs[0], coeffs[1], coeffs[2], coeffs[3])

        else
	        coeffs[ 3 ] = 1.0/2 * r * p - 1.0/8 * q * q;
	        coeffs[ 2 ] = - r;
	        coeffs[ 1 ] = - 1.0/2 * p;
	        coeffs[ 0 ] = 1;

            num, s0, s1, s3 = SolveCubic(coeffs[0], coeffs[1], coeffs[2], coeffs[3])

	        z = s0;

	        u = z * z - r;
	        v = 2 * z - p;

	        if (IsZero(u)) then
	            u = 0
	        elseif (u > 0) then
	            u = math.sqrt(u)
	        else
	            return 0
			end

	        if (IsZero(v)) then
	            v = 0;
	        elseif (q > 0) then
	            v = math.sqrt(v)
	        else
	            return 0;
			end

	        coeffs[ 2 ] = z - u
			if q < 0 then
				coeffs[ 1 ] = -v
			else
				coeffs[ 1 ] = v
			end
	        coeffs[ 0 ] = 1;

	        num, s0, s1 = SolveQuadric(coeffs[0], coeffs[1], coeffs[2]);

	        coeffs[ 2 ]= z + u
	        if q < 0 then
				coeffs[ 1 ] = v
			else
				coeffs[ 1 ] = -v
			end
	        coeffs[ 0 ] = 1;


            if (num == 0) then
				local temp, s0, s1 = SolveQuadric(coeffs[0], coeffs[1], coeffs[2])
				num = num + temp
			end
            if (num == 1) then
				local temp, s1, s2 = SolveQuadric(coeffs[0], coeffs[1], coeffs[2])
				num = num + temp
			end
            if (num == 2) then
				local temp, s2, s3 = SolveQuadric(coeffs[0], coeffs[1], coeffs[2])
				num = num + temp
			end
        end

        sub = 1.0/4 * A;

        if (num > 0) then   s0 = s0 - sub end
        if (num > 1) then   s1 = s1 - sub end
        if (num > 2) then   s2 = s2 - sub end
        if (num > 3) then   s3 = s3 - sub end

        return num, s0, s1, s2, s3
end

local function Solve(x, z, target_velocity, proj_speed)

		local G = 196.2

        local A = x
        local B = 0
        local C = z
        local M = 0
        local N = 0
        local O = 0
        local P = target_velocity.x
        local Q = target_velocity.y
        local R = target_velocity.z
        local S = proj_speed

        local H = M - A;
        local J = O - C;
        local K = N - B;
        local L = -.5 * G;

        local c0 = L*L;
        local c1 = 2*Q*L;
        local c2 = Q*Q + 2*K*L - S*S + P*P + R*R;
        local c3 = 2*K*Q + 2*H*P + 2*J*R;
        local c4 = K*K + H*H + J*J;

		local numTimes, t1, t2, t3, t4 = SolveQuartic(c0, c1, c2, c3, c4)

		print(numTimes)

		local solutions = {}
		local times = {t2, t3, t4}
		times[0] = t1

		local numSolutions = 0

		for i=0, 3 do
			local t = times[i]
			if t == nil then break end
			--if (t>0) == false then break end
			solutions[numSolutions] = Vector3.new(((H+P*t)/t), ((K+Q*t-L*t*t)/ t), ((J+R*t)/t))
			numSolutions = numSolutions + 1
		end

		local s0, s1

		if numSolutions > 0 then
			s0 = solutions[0]
		end
		if numSolutions > 0 then
			s1 = solutions[1]
		end

		return numSolutions, s0, s1, solutions

end

return Solve