Laguerre15m.h

Go to the documentation of this file.

 
 
int tdeg;
void bernoullid(){
    int k, k0, steps, n, kn1, be, se;
    float cla, kden,asa,qsa,i,j, s,t, pn,u1,u2,res3,d1,terr1;
 
    while(tca[k] != 0 && k < tdeg){++k;}
    k0 = k;
 
    kden = -tca[0]/(k0*tca[k0]);
 
    int *kia = new int[steps+1]; // each element equal to k0
    int *tqa = new int[steps+1]; // each element equal to tdeg*kden
    int *asa = new int[tdeg+1]; // 0
    int *qsa = new int[tdeg+1]; // 0
 
    for(int n=0; n < steps; ++n){
        kn1 = kia[n-1];
        if(kn != n){
            kia[n]=kn1;
            tqa[n]=0;
            continue;
        }
 
        k=n; be=std::max(0,n-tdeg);
 
        while(pn != 0 && k != tdeg){
            k+=1; 
            se=std::min(k,tdeg);
            for(int t=k-n+1; t <= se; ++t){
                asa[t] = tqa[k-t] != 0 ? tca[t] : 0;
                for(int s=k-t-be; n-1-be; ++s){
                    qsa[0] = (s+be == 0)*(k*kden) + (s+be != 0)*qsa[0];
                    qsa[s] = (s+be == 0)*(tqa[s+be] != 0 ? tqa[s+be] : 1);
                }
            }
 
            pn = se<=k-n ? 0 : asa[se]*qsa[k-se-be];
 
            for(int j = se; j >= k-n+2; --j){
                pn = (pn+asa[j-1])*qsa[k-j-be+1];
            }
            pn=pn+tca[k-n];
        }
 
        kia[n] =k;
        if(pn == 0) return []; // ????
        tqa[n] =-tca[0]/pn;
    }
    terr1 = err(tp, tqa[steps], true)[0];
    if(terr1 is not finite) throw error;
    if(terr < eps) return [tqa[steps],terr1];
 
    lam = (!flag_min)*lam + flag_min*terr1;
    ats = (!flag_min)*ats + flag_min*tshift;
    ass = (!flag_min)*ass + flag_min*sshift;
    tph = (!flag_min)*tph + flag_min*tp;
    stem = (!flag_min)*stem + flag_min*steps;
    sqa = (!flag_min)*sqa + flag_min*eval(tqa);
    kma = (!flag_min)*kma + flag_min*eval(kia)
}
 
void degtwoone(std::vector<std::complex<double>>& tp, int& tdeg) {
    std::complex<double> a, b, c, d, bz;
    bool flag_suc = true;
    int ladi = 0;
 
    if (tdeg == 2) {
        a = tp[2]; 
        b = -tp[1] / (static_cast<double>(2) * a); 
        c = tp[0] / a;
        d = std::sqrt(std::pow(b, 2) - c); 
        a = b + d; 
        b -= d;
 
        if (std::abs(a) > std::abs(b)) {
            finroot(bz + a); 
            finroot(bz + c / a);
        }
        if (std::abs(a) == std::abs(b)) {
            finroot(bz + a); 
            finroot(bz + b);
        }
        if (std::abs(a) < std::abs(b)) {
            finroot(bz + b); 
            finroot(bz + c / b);
        }
        tdeg -= 2;
    } else {
        finroot(-tp[0] / tp[1] + bz);
        tdeg -= 1;
    }
}
 
std::vector<double> err(double q, double x) {
    double la, ala;
    la = laguerre(q, x);
    // bool flag_lag = true;
 
    ala = degree(q) * std::fabs(la);
    if (ala < INFINITY) {
        return {ala, la};
    }
    std::cout << "No error result" << std::endl;
    return {INFINITY, false};
}
 
void minrad(){
    flag_once = false;
    double rmin = np.inf;
    for(int jl=0; j < len(rl); jl++){
        rmin = rl[jl] < rmin ? std::fabs(rl[jl]) : rmin;
    }
}
 
void findroot(std::vector<double>& p, std::vecvor<double>& rootlist, double& ro, bool& flag_suc, bool& flag_fin){
    auto nera = err(p, ro);
    auto bof = (nera[0] < inf)*bof + 0;
    if(flag_suc){
        rootlist.push_back(ro);
    }
    if(bof != 0) flag_fin = false;
    // nera[1] returned in maple
    return nera[1];
}
 
void degtwo(double bz, std::vector<double> tp, int& tdeg, std::vector<double>& rootlist, bool flag_fin){
    bool flag_suc = true;
    auto a = -0.5*tp[1] + np.sqrt(b**2 - c);
    auto b = -0.5*tp[1] - np.sqrt(b**2 - c);
    flag_f = flag_fin;
 
    finroot(tp, rootlist, bz + (abs(a) >= abs(b)*(a) + (abs(a) < abs(b))*(b)) );
    finroot(tp, rootlist, bz + (abs(a) == abs(b)*(b) + (abs(a) != abs(b))*(tp[2]/( (abs(a)>abs(b))*(a) + (abs(a)<abs(b))*(b) ))));
    tdeg -= 2;
}
 
void trial(int cl, int tdeg, std::vector<double>& tp,std::vector<double>& tp1, double bz, double bz1, double stepfactor,
                                                            double stepmin, bool& flag_suc, bool& flag_inc, bool flag_min){
    if(tp[-1]*tp[-1] < eps4 && tp[-2]*tp[-2] > eps4){
        auto qco = -tp[-1]/tp[-2];
        auto ladi = err(tp, qco)[0];
        if(ladi < eps){
            if(flag_inc = (err(tp1, qco + bz - bz1 && cl == 2)[0] >=eps && cl==2) ? true : false) return;
            flag_suc = true
            deflate(qco + bz - bz1);
            return;
        }
    }
    auto diam = 1.0 + std::max(tp[:-1]);
    auto steps = (cl == 1 || (cl==2 && flag_min))*np.ceil(stepfactor *(np.ln(diam*tdeg) + ln10n*(prec+1))) +
                 (cl >= 2 || (cl!=2 && flag_min))*(16*tdeg)
    std::vector<double> tca(tdeg, 0);
    auto _tp = list(reversed(tp));
    std::vector<double> tca = _tp;
    auto res2 = bernoullid();
    if(flag_inc) return;
    if(len(res2) == 2){
        ladi = res2[1];
        flag_suc = true;
        deflate(res2[0]+bz-bz1);
    }
    else if(len(res2) == 0){
        rad = minrad();
    }
    /*
    elif cl == 2 and flag_min:
      tqm = tqa
      kim = kia
    */
}
 
void run(double& bz, double& bzre, double& tp, int& tdeg, double& rootlist, double& stepfactor, bool flag_fin){
    flag_lag = true;
    if(tdeg < 3) degtwo(bz, tp, tdeg, rootlist, flag_fin);
    else{
        flag_suc= false;
        flag_inc = true;
        flag_min = false;
    }
 
    while(flag_inc){
        flag_inc = false;
        tp1 = tp;
        bz1 = bz;
        bz, bz1, flag_suc, flag_inc = trial(1, tp, tp1, tdeg, bz, bz1, stepfactor, flag_suc, flag_inc, flag_min);
        tp = tp1;
        bz = bz1;
    }
}