Решение линейных уравнений 3-й и 4-й степеней в комплексных числах на C++

!!! Внимание !!! Предложенный код не был использован в реальных приложениях - также не проводились исследования на корректность источников формул на основании которых был построен код.

    Немного о причине появления этого кода: Как-то в условиях нестабильного подключения к интернету потребовалось решать в общем виде уравнение четвертой степени (искал максимумы отклонений направлений на спутники на геостационарной орбите). Не имея у себя нужного математического справочника и ленясь рыться в старых конспектах (вообще то я им не очень доверяю), учитывая что в периоды до обрыва связи не успел найти нужный код но все-же успев найти http://www.n-t.ru/tp/ns/oam.htm - метод решения уравнений четвертой степени и http://ru.wikipedia.org/wiki/Формула_Кардано -формулы для решения уравнений третьей степени в порыве нетерпения начал их  реализовывать на C++. Не умея писать без ошибок и промучившись недельку исправляя ошибки получил правдоподобные результаты полученных функций - но и  их не смог применить - получались (в моих уравнениях) только комплексные корни - в итоге обошелся без уравнений четвертой степени и соответственно без их решения. Но код остался - не стирать же его с винчестера. Я конечно не гарантирую что код безусловно верен - но те ошибки, что сам нашел, - исправил. Код конечно не идеален - к примеру для решения уравнений третьей степени я позже наше более оптимальный код -http://algolist.manual.ru/maths/findroot/cubic.php - но возможно кому-то потребуется - хотя бы сдать лабораторную.

     Здесь далее приведены тексты файла кодов на С++ и его заголовочного файла (RadicalsX4.cpp и RadicalsX4.h соответственно).

    Итак, скопировав нижеприведенный текст в какой нибудь "блокнот" и сохранив с указанными именами, получаем файлы которые можно подключать к своим проектам. Возможно придется убрать в строке " #include <complex.h>" расширение файла ".h" - по крайней мере в справке к среде разработки (я пользовался "Borland C++Builder 4") имя библиотечного файла было  указано без расширения - но у меня оно потребовалось.

    Функции решают уравнения 2, 3, и 4-й степеней в комплексных числах. Исходные коэффициенты уравнений a, b, c и d - задаются в том виде - как они находятся в закомментированных пояснениях. Корни уравнений содержатся в массивах x которые должны быть обязательно предварительно объявлены. Необязательная переменная ToI -указывает количество нулей после запятой когда комплексная и действительная части округляются до нуля.
    Получить действительную и комплексную составляющие чисел можно при помощи функций real и imag например real(x[0]) и imag(x[0]).
     Также имеется функция нахождения всех трех кубических корней из комплексного числа - просто применяется для решения по формуле Кардано.

    Комментарий окончен - далее код:

Содержимое файла RadicalsX4.cpp
// Есьман Юрий Григорьевич 2006г
//---------------------------------------------------------------------------

#pragma hdrstop
#include <complex.h>
#include "RadicalsX4.h"

//---------------------------------------------------------------------------
#pragma package(smart_init)
//---------------------------------------------------------------------------
inline __fastcall int RadicalsX3(double a,double b,double c,double d,radicalsX3 &x,int ToI)
{//решение уравнения третьей степени a*x3+b*x2+c*x+d=0
radicalsX2 y;
radicalsX3 A,B;
int r;//счетчик действительных корней
if (a==0)//уравнение  не является уравнением третьей степени
   {r=RadicalsX2((complex<double>)b,(complex<double>)c,(complex<double>)d,y);
   x[0]=y[0];x[1]=y[1];if (r!=2)return r;//одно или нет решений
   else {if (imag(y[0])!=0) return 0;//нет действительных корней
        else return 2;}}//два действительных корня
d=d/a;c=c/a;b=b/a;//приведение к виду x3+b*x2+c*x+d=0
double p=c-b*b/3,q=d+(b*b/13.5-c/3)*b;//приведение к виду x3+p*x+q=0
complex<double> DL=-27*q*q-4*p*p*p;//дискриминант уравнения
DL=sqrt(-DL/(complex<double>)108);
double Roundoff=exp(log(0.1)*ToI);//минимальная учитываемае величина комплексной составляющей числа
try{
complex<double> ab=-(complex<double>)q/(complex<double>)2;
sqrt3(ab+DL, &A);
sqrt3(ab-DL, &B);
ab=A[0]*B[0];
if ((real(ab)<(-p/3+Roundoff))&&(real(ab)>(-p/3-Roundoff))&&(imag(ab)<Roundoff)
       &&(imag(ab)>-Roundoff)) x[0]=A[0]+B[0];
else {ab=A[0]*B[1];
 if ((real(ab)<(-p/3+Roundoff))&&(real(ab)>(-p/3-Roundoff))&&(imag(ab)<Roundoff)
        &&(imag(ab)>-Roundoff)) x[0]=A[0]+B[1];
 else {ab=A[0]*B[2];
  if ((real(ab)<(-p/3+Roundoff))&&(real(ab)>(-p/3-Roundoff))&&(imag(ab)<Roundoff)
         &&(imag(ab)>-Roundoff)) x[0]=A[0]+B[2];
  else {ab=A[1]*B[1];
   if ((real(ab)<(-p/3+Roundoff))&&(real(ab)>(-p/3-Roundoff))&&(imag(ab)<Roundoff)
          &&(imag(ab)>-Roundoff)) x[0]=A[1]+B[1];
   else {ab=A[1]*B[2];
    if ((real(ab)<(-p/3+Roundoff))&&(real(ab)>(-p/3-Roundoff))&&(imag(ab)<Roundoff)
          &&(imag(ab)>-Roundoff)) x[0]=A[1]+B[2];
    else {ab=A[2]*B[2];
    if ((real(ab)<(-p/3+Roundoff))&&(real(ab)>(-p/3-Roundoff))&&(imag(ab)<Roundoff)
          &&(imag(ab)>-Roundoff)) x[0]=A[2]+B[2];
    }}}}}
}catch ( ... ){return 0;}//уравнение не разрешилось
x[0]=x[0]-b/3;//обратная подстановка
r=RadicalsX2((complex<double>)1,(complex<double>)(b+x[0]),(complex<double>)((b+x[0])*x[0]+c),y);
x[1]=y[0]; x[2]=y[1];
if (r>2) return 999999;//решений бесконечное множество
if (r==0) return 0;//уравнение не разрешилось
r=0;
for (int i=0;i<3;i++) if (((imag(x[i]))<Roundoff)&&((imag(x[i]))>-Roundoff))
        {x[i]=real(x[i]);//округляем мелкие комплексные части
        r++;}//считаем действительные корни
return r;
}
//---------------------------------------------------------------------------
inline __fastcall int RadicalsX2(complex<double> a,complex<double> b,complex<double> c,radicalsX2 &x)
{//решение уравнения второй степени a*x2+b*x+c=0
if (a==(complex<double>)0)//линейное уравнение b*x+c=0
        {if (b==(complex<double>)0)//нет переменных
                {if (c==(complex<double>)0) return 999999;//бесконечное множиство решений
                else return 0;}//нет решений
        x[0]=x[1]=-c/b;return 1;}//едиственное решение
x[1]=sqrt(b*b-((complex<double>)4*a*c));
x[0]=(-b+x[1])/((complex<double>)2*a);
x[1]=(-b-x[1])/((complex<double>)2*a);
return 2;//уравнение действительно квадратное
}
//---------------------------------------------------------------------------
inline __fastcall void sqrt3(complex<double> a,radicalsX3 &x)
{//три корня третьей степени из комплексного числа а
double X=real(a),y=imag(a),R=sqrt(y*y+X*X);
if (R==0) {x[0]=x[1]=x[2]=0; return;}
double Q=acos(X/R);
complex<double> i(0.0, 1.0);
R=exp(log(R)/3);
if (y<0) Q=6.28318530717958647692528676655901-Q;
x[0]=R*cos(Q/3)+R*sin(Q/3)*i;
x[1]=R*cos((Q+3.1415926535897932384626433832795)/3)+R*sin((Q+3.1415926535897932384626433832795)/3)*i;
x[2]=R*cos((Q+6.28318530717958647692528676655901)/3)+R*sin((Q+6.28318530717958647692528676655901)/3)*i;
}
//---------------------------------------------------------------------------
inline __fastcall int RadicalsX4(double a,double k,double t,double p,double q,radicalsX4 &x,int ToI)
{//решение уравнения четвертой степени a*x4+k*x3+t*x2+p*x+q=0
int r;//счетчик действительных корней
radicalsX3 A;//корни уравнения третьей степени
radicalsX2 x01,x23;//корни урравнений второй степени
double Roundoff=exp(log(0.1)*ToI);//минимальная учитываемае величина комплексной составляющей числа
if (a==0)//уравнение  не является уравнением четвертой степени, решаем урвнение третьей степени
   {r=RadicalsX3(k,t,p,q,A,ToI);x[0]=A[0];x[1]=A[1];x[2]=A[2];return r;}
q=q/a;p=p/a;t=t/a;k=k/a;//приведение к виду x4+k*x3+t*x2+p*x+q=0
k=k/4; q=q+((t-k*k*3)*k-p)*k; p=p+k*(k*k*8-t*2); t=t-k*k*6;//приведение к виду x4+t*x2+p*x+q=0
r=RadicalsX3(1,t*2,t*t-q*4,-p*p,A,ToI);//решение приведенного уравнения третьей степени
if ((real(A[0])<0)&&(real(A[0])>-Roundoff))A[0]=0;//приближение малых отрицательных чисел к нулю (избегаем комплексных решений)
else A[0]=0; // первый член всегда положителен -комплексная часть результат погрешностей
complex<double> A12=sqrt(A[0]),AT12=sqrt((A[0]+t)*(A[0]+t)-q*4);
RadicalsX2(1,-A12,(A[0]+t+AT12)*0.5,x01);
RadicalsX2(1,A12,(A[0]+t-AT12)*0.5,x23);
x[0]=x01[0];x[1]=x01[1];x[2]=x23[0];x[3]=x23[1];
x[0]=x[0]-k;x[1]=x[1]-k;x[2]=x[2]-k;x[3]=x[3]-k;
//используем только первый из трех корней A - две другие четверки корней примерно равны - с точностью вычислений
r=0;
for (int i=0;i<4;i++)
     {if (((imag(x[i]))<Roundoff)&&((imag(x[i]))>-Roundoff))
        {x[i]=real(x[i]);//округляем мелкие комплексные части
        A12=x[r];x[r]=x[i];x[i]=A12;//перемещаем целые части вперед
        r++;}//считаем действительные корни
     if (((real(x[i]))<Roundoff)&&((real(x[i]))>-Roundoff))
        x[i]=x[i]-real(x[i]);//округляем мелкие целые части
     }    
return r;
}
//---------------------------------------------------------------------------

Содержимое файла RadicalsX4.h
// Есьман Юрий Григорьевич 2006г
//---------------------------------------------------------------------------

#ifndef RadicalsX4H
#define RadicalsX4H

//---------------------------------------------------------------------------
typedef complex<double> radicalsX4[4];//корни уравнения 4-й степени
typedef complex<double> radicalsX3[3];//корни уравнения 3-й степени
typedef complex<double> radicalsX2[2];//корни уравнения 2-й степени
    //решение уравнения третьей степени a*x3+b*x2+c*x+d=0
inline __fastcall int RadicalsX3(double a,double b,double c,double d,radicalsX3 &x,int ToI=10);
    //решение уравнения второй степени a*x2+b*x+c=0
inline __fastcall int RadicalsX2(complex<double> a,complex<double> b,complex<double> c,radicalsX2 &x);
    //три корня третьей степени из комплексного числа а
inline __fastcall void sqrt3(complex<double> a,radicalsX3 &x);
    //решение уравнения четвертой степени a*x4+k*x3+t*x2+p*x+q=0
inline __fastcall int RadicalsX4(double a,double k,double t,double p,double q,radicalsX4 &x,int ToI=10);
//---------------------------------------------------------------------------
#endif

28 января 2006г.
Есьман Юрий Григорьевич

"Скорочей" - читай быстрей - а программа в помощь

Абгрейдь мозги - читай быстрее