#include <cmath>

//*************************************************************************
float polinom(int n,float x,float *k)
{
float s=1;
for(int i=n-1;i>=0;i--)
s=s*x+k[i];
return s;
}//polinom

float dihot(int degree,float edgeNegativ,float edgePositiv,float *kf)
{
for(;;)
{
float x=0.5*(edgeNegativ+edgePositiv);
if(std::abs(x-edgeNegativ) < 1e-3 || std::abs(x-edgePositiv) < 1e-3)return x;
if(polinom(degree,x,kf)<0)edgeNegativ=x;
else edgePositiv=x;
}
}//dihot

void stepUp(
int level,
float **A,
float **B,
int *currentRootsCount
)
{

float major=0;
for(int i=0;i<level;i++)
{
float s=fabs(A[level][i]);
if(s>major)major=s;
}
major+=1.0;

currentRootsCount[level]=0;

for(int i=0;i<=currentRootsCount[level-1];i++)
{
int signLeft,signRight;
float edgeLeft,edgeRight;
float edgeNegativ,edgePositiv;

if(i==0)edgeLeft=-major;
else edgeLeft=B[level-1][i-1];

float rb=polinom(level,edgeLeft,A[level]);

if(rb==0)
{
B[level][currentRootsCount[level]]=edgeLeft;
currentRootsCount[level]++;
continue;
}

if(rb>0)signLeft=1;else signLeft=-1;

if(i==currentRootsCount[level-1])edgeRight=major;
else edgeRight=B[level-1][i];

rb=polinom(level,edgeRight,A[level]);

if(rb==0)
{
B[level][currentRootsCount[level]]=edgeRight;
currentRootsCount[level]++;
continue;
}

if(rb>0)signRight=1;else signRight=-1;
if(signLeft==signRight)continue;

if(signLeft<0){edgeNegativ=edgeLeft;edgePositiv=edgeRight;}
else {edgeNegativ=edgeRight;edgePositiv=edgeLeft;}

B[level][currentRootsCount[level]]=dihot(level,edgeNegativ,edgePositiv,A[level]);
currentRootsCount[level]++;
}
return;
}//stepUp

void polynomRealRoots(
float *rootsArray,
int n,
float *kf_,
int &rootsCount
)
{

float *kf = new float [n+1];

float **A=new float *[n+1];
float **B=new float *[n+1];

int *currentRootsCount=new int[n+1];

for(int i=1;i<=n;i++)
{
A[i]=new float[i];
B[i]=new float[i];
}

for(int i=0;i<=n;i++)kf[i]=kf_[n-i];

for(int i=0;i<n;i++)A[n][i]=kf[i]/kf[n];

for(int i1=n,i=n-1;i>0;i1=i,i--)
{
for(int j1=i,j=i-1;j>=0;j1=j,j--)
{
A[i][j]=A[i1][j1]*j1/i1;
}
}

currentRootsCount[1]=1;
B[1][0]=-A[1][0];

for(int i=2;i<=n;i++)stepUp(i,A,B,currentRootsCount);

rootsCount=currentRootsCount[n];
for(int i=0;i<rootsCount;i++)rootsArray[i]=B[n][i];

for(int i=1;i<=n;i++)
{
delete[]A[i];
delete[]B[i];
}
delete[]A;
delete[]B;
delete[]currentRootsCount;
delete[]kf;
}//polynomRealRoots

//*************************************************************************