Где почитать про метод Лобачевского?

Question

[kykyryky]

Метод решения алгебраических уравнений.
Нашел метод у Демидовича, но, во-первых там не совсем то - метод Лобачевского-Греффе. Да и слишком сложно написано, не осилил. Никаких нормальных примеров, ничего.
Так же он есть в книге Березина "Методы вычислений", но это копипаст Демидовича.

В общем, нужны еще варианты. Спасибо.

Answer 1

[Армянское Радио]

Я тоже когда-то брал материал из того же Демидовича.
Идея метода состоит в том, чтобы повторными преобразованиями получать многочлен, корнями которого будут квадраты корней многочлена, с которого процесс начали.

При этом коэффициенты многочлена начинают принимать огромные значения, для которых может не хватить точности 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.