Полное решение кубического уравнения (формула Кардано)

Необходимо решить кубическое уравнение с действительными коэффициентами:   ax³ + bx² + cx + d = 0 .

По основной теореме алгебры оно имеет три корня (4 различных типа).

Формула была украдена у Тартальи и опубликована Джероламо Кардано (тем самым изобретателем карданного вала) в книге «Великое искусство» в 1545 году.  Описание алгоритма доступно в Интернете, однако программная реализация имеет некоторые нюансы, связанные с необходимостью анализа исходных данных. Вот что получилось:

// Решение кубического уравнения по формуле Кардано
private static void Kardano(double a, double b, double c, double d, ref int tip, ref double p1, ref double p2, ref double p3)
{
   double eps = 1E-14;
   double p = (3 * a * c - b * b) / (3 * a * a);
   double q = (2 * b * b * b - 9 * a * b * c + 27 * a * a * d) / (27 * a * a * a);
   double det = q * q / 4 + p * p * p / 27;
   if (Math.Abs(det) < eps)
      det = 0;
   if (det > 0)
   {
      tip = 1; // один вещественный, два комплексных корня
      double u = -q / 2 + Math.Sqrt(det);
      u = Math.Exp(Math.Log(u) / 3);
      double yy = u - p / (3 * u);
      p1 = yy - b / (3 * a); // первый корень
      p2 = -(u - p / (3 * u)) / 2 - b / (3 * a);
      p3 = Math.Sqrt(3) / 2 * (u + p / (3 * u));
   }
   else
   {
      if (det < 0)
      {
         tip = 2; // три вещественных корня
         double fi;
         if (Math.Abs(q) < eps) // q=0
            fi = Math.PI / 2;
         else
         {
            if (q < 0) // q<0
               fi = Math.Atan(Math.Sqrt(-det) / (-q / 2));
            else // q<0
               fi = Math.Atan(Math.Sqrt(-det) / (-q / 2)) + Math.PI;
         }
         double r = 2 * Math.Sqrt(-p / 3);
         p1 = r * Math.Cos(fi / 3) - b / (3 * a);
         p2 = r * Math.Cos((fi + 2 * Math.PI) / 3) - b / (3 * a);
         p3 = r * Math.Cos((fi + 4 * Math.PI) / 3) - b / (3 * a);
      }
      else // det=0
      {
         if (Math.Abs(q) < eps)
         {
            tip = 4; // 3-х кратный 
            p1 = -b / (3 * a); // 3-х кратный 
            p2 = -b / (3 * a);
            p3 = -b / (3 * a);
         }
         else
         {
            tip = 3; // один и два кратных
            double u = Math.Exp(Math.Log(Math.Abs(q)/2)/ 3);
            if (q < 0)
            u = -u;
            p1 = -2 * u - b / (3 * a);
            p2 = u - b / (3 * a);
            p3 = u - b / (3 * a);
         }
      }
   }
} // end Kardano

Входными параметрами метода являются коэффициенты a, b, c, d; по ссылке возвращается тип корней (tip=1,2,3,4) и их значения (p1,p2,p3):

Для типа 1 (tip=1) имеется один действительный и два комплексных корня: x1=p1; x2=p2+i*p3; x3=p2-i*p3, где i — мнимая единица.
Тип 2 — три различных действительных корня, тип 3 —  один отличающийся и два кратных действительных корня, тип 4 — три кратных действительных корня, для всех типов (tip=2,3,4) x1=p1, x2=p2, x3=p3.
Если кубическое уравнение является характеристическим уравнением исходного линейного дифференциального уравнения 3 степени, то для его решения важно знать именно тип решения (tip).

Для тестирования метода используйте следующую программу:

using System;
namespace KArdano
{
class Program
{
   static void Main(string[] args)
   {
      Console.WriteLine("Введите параметры кубического уравнения a*x^3+b*x^2+c*x+d=0 :");
      Console.Write("a=");
      double a =Convert.ToDouble(Console.ReadLine());
      Console.Write("b=");
      double b = Convert.ToDouble(Console.ReadLine());
      Console.Write("c=");
      double c = Convert.ToDouble(Console.ReadLine());
      Console.Write("d=");
      double d = Convert.ToDouble(Console.ReadLine()); 
      int tip = 0;
      double p1 = 0, p2 = 0, p3 = 0;
      Kardano(a, b, c, d, ref tip, ref p1, ref p2, ref p3);
      if (tip==1)
         Console.WriteLine("тип=1. Один вещественный и два комплексно сопряженных корня: x1={0} Re[x2,x3]={1} Im[x2,x3]={2}", p1, p2, p3);
      else
         Console.WriteLine("Вещественные корни: тип={0} p1={1} p2={2} p3={3}", tip, p1, p2, p3);
      Console.ReadKey();
   }
   // Сюда вставьте текст метода private static void Kardano( ) - см. выше
}

Минимальный набор тестов:
a, b, c, d:
187.5, 50, 10, 1 -> тип 1
1, 6, 3, -10 -> тип 2
1, 12, 36, 32 -> тип 3
3, -9, 9, -3 -> тип 4