Where to read about the Lobachevsky method?

Question

Mathematics

Method for solving algebraic equations.
I found a method from Demidovich, but, firstly, it’s not quite right there - the Lobachevsky-Greffe method. Yes, and it is too difficult to write, I did not master it. No normal examples, nothing.
It is also in Berezin's book "Methods of Computation", but this is Demidovich's copy-paste.
In general, more options are needed. Thank you.

Answer 1

I also once took material from the same Demidovich.
The idea of ​​the method is to obtain a polynomial by repeated transformations, the roots of which will be the squares of the roots of the polynomial from which the process began.
At the same time, the coefficients of the polynomial begin to take on huge values, for which the double precision may not be enough. In this situation, I used the GMP library .
My implementation is on github . Works only with real coefficients and roots!
The polynomial is given by a class derived from PolyCoeffProxy:

class PolyCoeffProxy
{
public:
    PolyCoeffProxy();
    virtual mpf_class getCoeff(int idx)=0;
    virtual int power()=0;
    void print();
};

Two methods need to be implemented - getCoeff(int idx) should return the coefficients (numbering from 0 to powers of x)
power() should return the power:

mpf_class SamplePoly::getCoeff(int idx)
{
    static const mpf_class coeff[]={12,-8,1};
    return(coeff[idx]);
}

int SamplePoly::power()
{
    return(2);
}

Then the solver is used like this:

int main()
{
    SamplePoly sample;
    Point* ret=LobachevskySolver::solve(&sample,1e-10,100);
    ret->print();
    delete ret;
    return 0;
}

Here 1e-10 is the deviation of the roots from the previous iteration when squaring (less is more precise), 100 is the maximum number of iterations.

Solver itself:

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);
}

The program outputs 6 and 2 if given the Maxima polynomial:
it also outputs 6 and 2.