Я тоже когда-то брал материал из того же Демидовича.
Идея метода состоит в том, чтобы повторными преобразованиями получать многочлен, корнями которого будут квадраты корней многочлена, с которого процесс начали.
При этом коэффициенты многочлена начинают принимать огромные значения, для которых может не хватить точности double. В этой ситуации я использовал библиотеку
GMP.
Моя реализация на
гитхабе.
Работает только с вещественными коэффициентами и корнями!
Полином задается классом, произведенным от PolyCoeffProxy:
class PolyCoeffProxy
{
public:
PolyCoeffProxy();
virtual mpf_class getCoeff(int idx)=0;
virtual int power()=0;
void print();
};
Нужно дать реализацию двум методам - getCoeff(int idx) должен возвращать коэффициенты (нумерация с 0 по степеням x)
power() должен возвращать степень:
mpf_class SamplePoly::getCoeff(int idx)
{
static const mpf_class coeff[]={12,-8,1};
return(coeff[idx]);
}
int SamplePoly::power()
{
return(2);
}
Тогда решатель используется так:
int main()
{
SamplePoly sample;
Point* ret=LobachevskySolver::solve(&sample,1e-10,100);
ret->print();
delete ret;
return 0;
}
Здесь 1e-10 - отклонение корней от предыдущей итерации при квадрировании (меньше-точнее), 100 - максимальное количество итераций.
Сам решатель:Point* LobachevskySolver::solve(PolyCoeffProxy *prox,double meps,int nmaxiter)
{
Point* ret=new Point(prox->power());
rational* coeff[2];
rational roots[prox->power()];
for(idx i=0;i<2;i++)
{
coeff[i]=new rational[prox->power()+1];
}
for(idx i=0;i<prox->power();i++)
{
roots[i]=2;
}
for(idx i=0;i<=prox->power();i++)
{
coeff[0][i]=prox->getCoeff(i);
}
idx l=1;
idx m=0;
for(;m<nmaxiter;m++)
{
coeff[l][0]=coeff[1-l][0]*coeff[1-l][0];
coeff[l][prox->power()]=coeff[1-l][prox->power()]*coeff[1-l][prox->power()];
for(idx i=1;i<prox->power();i++)
{
const idx d=std::min(prox->power()-i,i);
coeff[l][i]=coeff[1-l][i]*coeff[1-l][i];
for(idx j=1;j<=d;j++)
{
idx k=(((1+j)%2)*2-1)<<1;
coeff[l][i]+=k*coeff[1-l][i-j]*coeff[1-l][i+j];
}
}
mpf_class eps=0;
for(idx i=0;i<prox->power();i++)
{
const mpf_class znam=coeff[l][prox->power()-i];
mpf_class tmp=coeff[l][prox->power()-1-i]/znam;
tmp=abs(tmp);
for(idx j=0;j<=m;j++)
{
tmp=sqrt(tmp);
}
const mpf_class ceps=abs(tmp-roots[i]);
eps=std::max(eps,ceps);
roots[i]=tmp;
}
// printf("%10.8lf\n",eps.get_d());
if(eps<meps)
{
break;
}
l=1-l;
}
for(int i=0;i<prox->power();i++)
{
ret->coord[i]=roots[i].get_d();
}
for(int i=0;i<2;i++)
{
delete[] coeff[i];
}
return(ret);
}
Программа выдает 6 и 2, если дать многочлен Maxima:
float(solve(x^2-8*x+12,x));
она тоже дает 6 и 2.
