Skip to content

Instantly share code, notes, and snippets.

@micycle1
Created February 7, 2021 22:40
Show Gist options
  • Save micycle1/6bdae2f84e1d2f9a6e590151a4c135cd to your computer and use it in GitHub Desktop.
Save micycle1/6bdae2f84e1d2f9a6e590151a4c135cd to your computer and use it in GitHub Desktop.
A fast solution to cubic equations with three real roots
http://momentsingraphics.de/CubicRoots.html#_Blinn06a
Returns the three real roots of the cubic polynomial (of form 1 + 2x + 3x^2 + 4x^3)
float3 SolveCubic(float4 Coefficient){
// Normalize the polynomial
Coefficient.xyz/=Coefficient.w;
// Divide middle coefficients by three
Coefficient.yz/=3.0f;
// Compute the Hessian and the discrimant
float3 Delta=float3(
mad(-Coefficient.z,Coefficient.z,Coefficient.y),
mad(-Coefficient.y,Coefficient.z,Coefficient.x),
dot(float2(Coefficient.z,-Coefficient.y),Coefficient.xy)
);
float Discriminant=dot(float2(4.0f*Delta.x,-Delta.y),Delta.zy);
// Compute coefficients of the depressed cubic
// (third is zero, fourth is one)
float2 Depressed=float2(
mad(-2.0f*Coefficient.z,Delta.x,Delta.y),
Delta.x
);
// Take the cubic root of a normalized complex number
float Theta=atan2(sqrt(Discriminant),-Depressed.x)/3.0f;
float2 CubicRoot;
sincos(Theta,CubicRoot.y,CubicRoot.x);
// Compute the three roots, scale appropriately and
// revert the depression transform
float3 Root=float3(
CubicRoot.x,
dot(float2(-0.5f,-0.5f*sqrt(3.0f)),CubicRoot),
dot(float2(-0.5f, 0.5f*sqrt(3.0f)),CubicRoot)
);
Root=mad(2.0f*sqrt(-Depressed.y),Root,-Coefficient.z);
return Root;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment