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

//*************************************************************************
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 <math.h>
#include <polynomComplexRoot.h>
#include <polynomRealRoots.h>

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<degree;n=np1,np1++)
for(int k=0;k<=np1;k++)
{//реализация рекуррентных формул
//ck[n+1,k]=ck[n,k-1]+sk[n,k-2]-sk[n,k];
//sk[n+1,k]=sk[n,k-1]+ck[n,k];
if(k>=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(mc<minmod)
{//точка up лучше точки dn
minmod=mc;
mindn=false;
}//точка up лучше точки dn
}//контрольный расчёт в точках up и dn

//формирование результата
if(mindn)
{//выбор точки dn в качестве результата
rootRe=rdn*cos(fidn);
rootIm=rdn*sin(fidn);
}//выбор точки dn в качестве результата
else
{//выбор точки up в качестве результата
rootRe=rup*cos(fiup);
rootIm=rup*sin(fiup);
}//выбор точки up в качестве результата

return errorCode;
}//run

//основная рабочая процедура
//при заданном r находит минимальное значение вещественной части
//исследуемого полинома от аргумента r*exp(i*fi)
//при том значении fi, которое обеспечивает равенство нулю его мнимой части
//это значение fi будет присвоено второму аргументу процедуры по окончании вызова
double polinomComplexRoot::minAxeXcrossVal(double r,double& fi)
{
//предварительно формируем вспомогательный массив rpw[k]=r^k
double rb=1;
for(int i=0;i<=degree;i++){rpw[i]=rb;rb*=r;}

//значение исследуемого полинома от вещественного аргумента r
rb=kf[0];for(int i=1;i<=degree;i++)rb+=rpw[i]*kf[i];
double rez=rb;fi=0;

//значение исследуемого полинома от вещественного аргумента -r
rb=kf[0];
for(int i=1,od=1;i<=degree;i++,od=1-od)if(od)rb-=rpw[i]*kf[i];else rb+=rpw[i]*kf[i];

//мы ищем минимальное значение абсциссы эпициклоиды
if(rb<rez){rez=rb;fi=cPI;}

//мнимая часть исследуемого полинома от комплексного аргумента,
//представленного параметрами r и fi, выражается так
//im=r*kf[1]*sin(fi)+r^2*kf[2]*sin(2*fi)+...+r^n*kf[n]*sin(n*fi)
//она будет представлена с участием некоторого полинома от cos(fi)
//коэффициенты этого полинома обозначены идентификатором ki
//im=sin(fi)*(ki[0]+ki[1]*cos(fi)+ki[2]*cos^2(fi)+...+ki[n]*cos^n(fi))
//коэффициенты ki выражаются через коэффициенты исследуемого полинома kf
//и коэффициенты ck и sk, вычисленные в конструкторе класса по формулам
//ki[0]=r*kf[1]*sk[1][0]+r^2*kf[2]*sk[2][0]+...+r^n*kf[n]*sk[n][0]
//ki[1]=r*kf[1]*sk[1][1]+r^2*kf[2]*sk[2][1]+...+r^n*kf[n]*sk[n][1]
//...
//ki[n]=r*kf[1]*sk[1][n]+r^2*kf[2]*sk[2][n]+...+r^n*kf[n]*sk[n][n]
for(int i=0;i<=degree;i++)
{//вычисление коэффициентов ki
rb=0;
for(int j=i+1;j<=degree;j++)rb+=(rpw[j]*kf[j]*sk[j][i]);
ki[i]=rb;
}//вычисление коэффициентов ki

//cosDegree это степень косинусного полинома, представленного коэффициентами ki
//страхуемся от возможности равенства нулю старших коэффициентов
int cosDegree=0,cosRootsCount=0;
for(int i=degree-1;i>0;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;i<cosRootsCount;i++)if(fabs(cosRootsArray[i])<1)
{//расчёт fi и коррекция rez
double x=acos(cosRootsArray[i]),
im=0,re=kf[0];

//расчёт вещественной (re) и мнимой (im) части исследуемого полинома
//при очередном значении найденного корня
for(int j=1;j<=degree;j++)
{
re+=(kf[j]*rpw[j]*cos(j*x));
im+=(kf[j]*rpw[j]*sin(j*x));
}

//существенно ненулевое значение im интерпретируем как внутреннюю ошибку
if(fabs(im)>1e-6)errorCode+=6;

//мы ищем минимальное значение абсциссы эпициклоиды
if(re<rez){rez=re;fi=x;}
}//расчёт fi и коррекция rez

//интерпретируем невероятный случай fi=0 или fi=cPI
//как внутреннюю ошибку
if(fi==0||fi==cPI)errorCode+=7;
return rez;
}//minAxeXcrossVal
//*************************************************************************

//*************************************************************************
//демонстрация расчёта комплексного корня полинома
#include <stdio.h>
#include <math.h>
#include <polynomComplexRoot.h>

int main()
{
//задание степени и коэффициентов исходного полинома
int degree=4;
//double kf[5]={24,24,12,4,1};
//double kf[5]={2,-2,3,-2,1};
double kf[5]={4,0,0,0,1};

//запуск процесса расчёта комплексного корня
double re,im;
int ret=polinomComplexRoot(degree,kf).run(re,im);

//распечатка результата расчёта корня
if(ret==0||ret>4)
{//успешное завершение процедуры расчёта корня
//распечатка результата расчёта корня
if(ret==0)printf("успешное решениеn");
else 
{
printf("ПРЕДУПРЕЖДЕНИЕnбыли странные ситуацииn");
printf("возможно из-за недостаточной точности аппаратного представления вещественных чиселn");
}
printf("кореньВещт=%fnкореньМним=%fn",re,im);
//проверка достоверности найденного корня
//используется следующая формула для полинома от суммы аргументов:
//p0(x+y)=p0(x)+p1(x)*y+p2(x)*y^2/2!+...+pn*y^n/n!
//представляющая собой конечную сумму ряда Тейлора по приращению y
//для исходного полинома p0 в точке x
//здесь p1, p2, ... , pn - производные 1-го, 2-го, ... , n-го порядка от исходного полинома
//в тексте программы далее использован обратный порядок индексации исходного и производных полиномов
//там индекс полинома равен его степени
double** kfx=new double*[degree+1];
for(int i=0;i<=degree;i++)kfx[i]=new double[degree+1];
double* kfy=new double[degree+1];

for(int i=degree;i>=0;i--)
{//исходные присвоения
for(int j=0;j<=degree;j++)kfx[i][j]=0;
kfy[i]=0;
kfx[degree][i]=kf[i];
}//исходные присвоения

//расчёт коэффициентов производных полиномов
for(int i=degree;i>0;i--)
for(int j=i;j>0;j--)kfx[i-1][j-1]=kfx[i][j]*j;

double fact=1;
for(int i=0,dmi=degree;i<=degree;i++,dmi--)
{//расчёт коэффициентов полинома по y
//вычисляется по схеме Горнера значение производного полинома от аргумента re
kfy[i]=kfx[dmi][dmi];
for(int j=dmi-1;j>=0;j--)kfy[i]=re*kfy[i]+kfx[dmi][j];
kfy[i]/=fact;
fact*=(i+1);
}//расчёт коэффициентов полинома по y

//массив степеней мнимой единицы
const int ipw[4]={1,1,-1,-1};

double ypw=1;

//вещественная и мнимая части значения исходного полинома
//при значении комплексного аргумента re+i*im
double fre=0,fim=0;
for(int i=0,od=0;i<=degree;i++,od=1-od)
{//расчёт значения исходного полинома
if(od==0)fre+=(ypw*kfy[i]*ipw[i%4]);
else fim+=(ypw*kfy[i]*ipw[i%4]);
ypw*=im;
}//расчёт значения исходного полинома

//печать погрешности расчёта корня исходного полинома
//это модуль полинома в найденном корне,
//нормированный на значение свободного члена
printf("погрешность=%7.1en",sqrt(fre*fre+fim*fim)/kf[0]);

//заключительное освобождение занимаемой памяти
for(int i=0;i<=degree;i++)delete[] kfx[i];
delete[] kfx;
delete[] kfy;
}//успешное решение
if(ret>0)switch(ret)
{//печать краткой диагностики произошедшей ошибки
case 1:
printf("#ошибка 1:нулевое значение свободного члена исходного полиномаn");
break;
case 2:
printf("#ошибка 2:нулевое значение старшего коэффициента полиномаn");
break;
case 3:
printf("#ошибка 3:степень исходного полинома слишком малаn");
break;
case 4:
printf("#ошибка 4:исходный полином имеет вещественные корниn");
break;
default:printf("#код ошибки=%dn",ret);
}//печать краткой диагностики произошедшей ошибки
return 0;
}//main
//*************************************************************************