Алгоритм вычисления комплексного корня полинома произвольной степени16.06.2016 21:48

//*************************************************************************
class polinomComplexRoot
{
public:
//тело класса определяется степенью исследуемого полинома и набором его коэффициентов
polinomComplexRoot(int _degree,double* _kf);
~polinomComplexRoot();

//основная процедура внешнего интерфейса класса
//возвращает код ошибки
//при успешном завершении код ошибки равен нулю,
//а значения параметров rootRe и rootIm станут равны соответственно
//вещественной и мнимой части комплексного корня исследуемого полинома
int run(double& rootRe,double& rootIm);

private:
int degree,errorCode;
double *kf;
double *ki;
double *cosRootsArray;
double *rpw;
int **ck;
int **sk;

//основная рабочая процедура
//при заданном r находит минимальное значение вещественной части
//исследуемого полинома от аргумента r*exp(i*fi)
//при том значении fi, которое обеспечивает равенство нулю
//его мнимой части. Это значение fi 
//будет присвоено второму аргументу процедуры по окончании вызова
double minAxeXcrossVal(double r,double& fi);
};//class polinomComplexRoot
//*************************************************************************

//*************************************************************************
#include 
#include 
#include 

const double cPI=3.14159265358979323846;

polinomComplexRoot::polinomComplexRoot(int _degree,double* _kf)
{
degree=_degree;errorCode=0;

//предварительная инициализация
kf=0;ck=0;sk=0;rpw=0;ki=0;cosRootsArray=0;

//отбраковка неприемлемых значений коэффициентов исходного полинома
if(_kf[0]==0){errorCode=1;return;}
if(_kf[degree]==0){errorCode=2;return;}
if(degree<2){errorCode=3;return;}

kf=new double[1+degree];
ki=new double[1+degree];
cosRootsArray=new double[1+degree];
rpw=new double[degree+1];
ck=new int*[1+degree];
sk=new int*[1+degree];
for(int i=0;i<=degree;i++)
{
ck[i]=new int[1+degree];
sk[i]=new int[1+degree];
}
//отбраковка исходного полинома, имеющего вещественные корни
polynomRealRoots(degree,_kf,kf,errorCode);
if(errorCode>0){errorCode=4;return;}

for(int i=0;i<=degree;i++)kf[i]=_kf[i]/_kf[0];


for(int i=0;i<=degree;i++)
for(int j=0;j<=degree;j++)
{ck[i][j]=0;sk[i][j]=0;}

//коэффициенты косинусных полиномов для представления
//кратных синусов и косинусов по формулам
//cos(n*x)=ck[n,0]+ck[n,1]*cos(x)+ck[n,2]*cos^2(x)+...+ck[n,n]*cos^n(x)
//sin(n*x)=sin(x)*(sk[n,0]+sk[n,1]*cos(x)+sk[n,2]*cos^2(x)+...+sk[n,n]*cos^n(x))
ck[1][1]=1;
sk[1][0]=1;

//расчёт ck и sk по рекуррентным формулам
//здесь np1=n+1
for(int n=1,np1=2;n=1)ck[np1][k]+=ck[n][k-1];
if(k>=2)ck[np1][k]+=sk[n][k-2];
if(k>=1)sk[np1][k]+=sk[n][k-1];
if(k<=n){ck[np1][k]-=sk[n][k];sk[np1][k]+=ck[n][k];}
}//реализация рекуррентных формул
return;
}//constructor

polinomComplexRoot::~polinomComplexRoot()
{
if(kf)delete[] kf;
if(ki)delete[] ki;
if(cosRootsArray)delete[] cosRootsArray;
if(rpw)delete[] rpw;
if(ck)
{
for(int i=0;i<=degree;i++)delete[] ck[i];
delete[] ck;
}

if(sk)
{
for(int i=0;i<=degree;i++)delete[] sk[i];
delete[] sk;
}

}//destructor

int polinomComplexRoot::run(double& rootRe,double& rootIm)
{
rootRe=0;rootIm=0;
if(errorCode>0)return errorCode;

//верхняя rup и нижняя rdn границы поиска для r 
double fidn=0,rdn=0,fiup=0,rup=1;

//удваиваем rup пока оно недостаточно велико
for(;minAxeXcrossVal(rup,fiup)>0;)rup+=rup;

for(;;)
{//цикл деления пополам интервала (rdn, rup)
double fit,rt=0.5*(rdn+rup);
if(rt==rdn||rt==rup)break;
if(minAxeXcrossVal(rt,fit)>0){rdn=rt;fidn=fit;}else {rup=rt;fiup=fit;}
}//цикл деления пополам интервала (rdn, rup)

//формальный выбор для выдачи результата одного из двух
//практически одинаковых решений (rdn, fidn) или (rup, fiup)
//будет вычислен модуль исследуемого полинома в каждой из этих двух точек
//и в качестве результата решения будет выдана та из этих точек,
//которой соответствует вычислительно меньшее значение этого модуля
double minmod; //минимальное значение модуля полинома
bool mindn=true; //минимум достигнут в точке с суффиксом dn

for(int c=0;c<2;c++)
{//контрольный расчёт в точках up и dn
double rc,fic;
if(c==0){rc=rdn;fic=fidn;}
else {rc=rup;fic=fiup;}

//rec - вещественная часть исследуемого полинома в пробной точке
//imc - мнимая часть исследуемого полинома в пробной точке
double rec=kf[degree]*cos(degree*fic),
imc=kf[degree]*sin(degree*fic);
for(int i=degree-1;i>0;i--)
{//gorner
rec=rc*rec+kf[i]*cos(i*fic);
imc=rc*imc+kf[i]*sin(i*fic);
}//gorner
rec=rc*rec+kf[0];
imc=rc*imc;

//mc - квадрат модуля исследуемого полинома в пробной точке
double mc=rec*rec+imc*imc;
if(c==0)minmod=mc;
else if(mc0;i--)if(fabs(ki[i])>0){cosDegree=i;break;}

//интерпретируем ситуацию вырождения ki-полинома как внутреннюю ошибку
if(cosDegree<1){errorCode=5;return rez;}

//находим все вещественные корни ki-полинома
polynomRealRoots(cosDegree,ki,cosRootsArray,cosRootsCount);

//обследование найденных корней ki-полинома
for(int i=0;i1e-6)errorCode+=6;

//мы ищем минимальное значение абсциссы эпициклоиды
if(re