################################################################################
PQDEquation := module()
option package;
export
    PolynomialSolution,
    RationalSolution,
    UniversalDenominator,
    AccurateQSummation,
    QZeilberger,
    SetOutputShort,
    SetExtended;
local
    add_to_S2,
    VarsZeroSystem,
    Extended,
    OutputShort,
    check_arguments,
    ZeroSystem,
    check_system,
    check_pair_of_systems,
    make_return,
    eval_spec,
    eval_spec2,
    Equation,
    Filter,
    ParametricGCRD2,
    FractionFree,
    Piecewise, 
    Sqrfree, 
    UpperBound, 
    UpperBoundOfCondition, 
    adjust, 
    q_DENOM;

    OutputShort := false;
    VarsZeroSystem := {};
    Extended := false;

################################################################################
SetExtended := proc(f::truefalse); 
    Extended := f;
    return NULL;
end proc:
SetOutputShort := proc(f::truefalse); 
    OutputShort := f;
    return NULL;
end proc:
################################################################################
PolynomialSolution := proc () 
local reqn, y, x, q, params, SS1, SS2, info, QEq, Res, OVZS,
      a, f, SS, wq, l, Sol, V, res, S, A, A2, cnd, eq, tmp, i, j, psi; 
global _C;

    OVZS := VarsZeroSystem;
    try reqn, y, x, q, params, SS1, SS2 := check_arguments(args[]) 
    catch: 
        error  
    end try;

    if SS1 = {1 = 0} or SS2 = {0 = 0} then 
        VarsZeroSystem := OVZS;
        return NULL 
    end if;
 
    if reqn = 0 then 
        VarsZeroSystem := OVZS;
        return make_return('ANY_POLYNOMIAL', SS1, SS2, params); 
    elif not reqn::'specfunc'('anything', 'QESol') then
        VarsZeroSystem := OVZS;
        return NULL;
    end if;

    info := eval(op(3, reqn)); 
    QEq := op([1, 1], reqn); 
    if not has(QEq, params) then 
        Res := QDifferenceEquations:-PolynomialSolution(QEq, y(x)); 
        if Res <> 0 and Res <> NULL then 
            VarsZeroSystem := OVZS;
            return make_return(collect(Res, x, radnormal), SS1, SS2, params); 
        else 
            VarsZeroSystem := OVZS;
            return NULL; 
        end if 
    end if;

    a := info['coeffs'][2 .. -1]; 
    f := -info['coeffs'][1]; 

    ## find values of parameters such that QEq is 0=0
    SS := map(ZeroSystem, {a[], f}, VarsZeroSystem, params);
    SS := check_system(map(op, SS), params);
    if SS <> {1 = 0} then
        Res := procname(reqn, 'parameters' = params, 
                                               'S1' = SS1 union SS, 'S2' = SS2);
        SS, SS2 := check_pair_of_systems(SS1, add_to_S2(SS, SS2), 
                                                                params)[1..2][];
        if SS = {1 = 0} or SS2 = {0 = 0} then
            VarsZeroSystem := OVZS;
            return Res;
        end if;
        if SS <> SS1 then 
            VarsZeroSystem := OVZS;
            return Res, procname(reqn, 'parameters' = params, 
                                                         'S1' = SS, 'S2' = SS2);
        end if;
    else
        Res := NULL;
    end if;

    wq := max(map(degree, a, q)); 
    l := max(wq, degree(f, x));

    Sol := add(_C[i]*x^i, i = 0 .. l); 
    V := Vector( [seq(_C[i], i = 0 .. l), -1]); 
    res := collect(expand(eval_spec((lhs-rhs)(QEq), y(x) = Sol)), x); 
    S := [coeffs(res, x)]; 
    S := map(primpart, S, {params[], seq(_C[i], i = 0 .. l)});
    S := eval_spec2(S, SS1, params);
    A := LinearAlgebra:-GenerateMatrix(S, [seq(_C[i], i = 0 .. l)], 
                                                                   'augmented');
    A2 := LinearAlgebra:-GaussianElimination(A);

    cnd := NULL;
    for i from LinearAlgebra:-Rank(A2) by -1 to 1 do
        eq := radnormal(A2[i].V);
        tmp := denom(eq);
        tmp := eval_spec2(tmp, SS1, params);
        if tmp = 0 then 
            tmp := ZeroSystem(denom(eq), VarsZeroSystem, params);
            SS := select(el -> lhs(el)::`*`, tmp);
            if SS <> {} then
                tmp := op(1, lhs(SS[1]));
                Res := Res, procname(reqn, 'parameters' = params, 
                                        'S1' = {tmp = 0} union SS1, 'S2' = SS2);
                Res := Res, procname(reqn, 'parameters' = params, 
                                        'S1' = SS1, 'S2' = SS2 union {tmp = 0});
                VarsZeroSystem := OVZS;
                return Res;
            end if;
            VarsZeroSystem := OVZS;
            return Res;
        end if;        
        SS := ZeroSystem(tmp, VarsZeroSystem, params);
        if SS <> {1 = 0} then
            Res := Res, procname(reqn, 'parameters' = params, 
                                               'S1' = SS1 union SS, 'S2' = SS2);
            SS, SS2 := check_pair_of_systems(SS1, add_to_S2(SS, SS2), 
                                                                params)[1..2][];
            if SS = {1 = 0} or SS2 = {0 = 0} then
                VarsZeroSystem := OVZS;
                return Res;
            end if;
            if SS <> SS1 then 
                VarsZeroSystem := OVZS;
                return Res, procname(reqn, 'parameters' = params, 
                                                         'S1' = SS, 'S2' = SS2);
            end if;
        end if;
        eq := numer(eq);
        eq := eval_spec2(eq, SS1, params);
        if eq = 0 then next end if;
        j := min(map2(op, 1, indets(eq, _C['integer'])));
        if j = infinity then
            psi := eq;
        else
            psi := coeff(eq, _C[j]);
        end if;
        SS := ZeroSystem(psi, VarsZeroSystem, params);
        if SS <> {1 = 0} then
            Res := Res, procname(reqn, 'parameters' = params, 
                                               'S1' = SS1 union SS, 'S2' = SS2);
            SS, SS2 := check_pair_of_systems(SS1, add_to_S2(SS, SS2), 
                                                                params)[1..2][];
            if SS = {1 = 0} or SS2 = {0 = 0} then
                VarsZeroSystem := OVZS;
                return Res;
            end if;
            if SS <> SS1 then 
                VarsZeroSystem := OVZS;
                return Res, procname(reqn, 'parameters' = params, 
                                                         'S1' = SS, 'S2' = SS2);
            end if;
        end if;
        if j = infinity then
            Sol := 0;
        else
            tmp := solve(eval(eq, {cnd}), _C[j]);
            tmp := eval_spec2(tmp, SS1, params);
            cnd := cnd, _C[j] = tmp;
        end if;
    end do;

    tmp := eval(Sol, {cnd});
    if tmp <> 0 then 
        Res := Res, make_return(tmp, SS1, SS2, params);
    end if;
    VarsZeroSystem := OVZS;
    return Res;
end proc:
################################################################################
RationalSolution := proc()
local reqn, info, params, SS1, SS2, y, x, q, QEq, tmp, i,j, Res, 
      W, NewQEq, z, Res1, Res2, A, t, OVZS;
global _C;

    OVZS := VarsZeroSystem;
    try reqn, y, x, q, params, SS1, SS2 := check_arguments(args[])
    catch: error 
    end try;

    info := eval(op(3, reqn));
    QEq := op([1,1], reqn);
    if not has(QEq, params) then
        Res := QDifferenceEquations:-RationalSolution(QEq, y(x));
        if Res = NULL or Res = 0 then
            VarsZeroSystem := OVZS;
            return NULL;
        else
            VarsZeroSystem := OVZS;
            return make_return(Res, SS1, SS2, params);
        end if;
    end if;

    W := UniversalDenominator(reqn, 'parameters' = params, 
                                                        'S1' = SS1, 'S2' = SS2);
    if W = NULL then 
        VarsZeroSystem := OVZS;
        return NULL;
    end if;
    if nops([W]) = 1 then
        W := {}, W;
    end if;

    Res := NULL;
    for j to nops([W])/2 do
        if W[2*j] = 'ANY_POLYNOMIAL' then
            Res := Res, W[2*j-1], 'ANY_RATIONAL';
        else
            t := W[2*j-1];
            t := `if`(t::'set', t, {t});
            NewQEq := eval_spec((lhs-rhs)(QEq), y(x) = z(x)/W[2*j]);
            tmp := PolynomialSolution(NewQEq, z(x), 'parameters' = params, 
                     'S1' = SS1 union select(type, t, `=`), 
                     'S2' = SS2 union 
                           map(el -> lhs(el) = rhs(el), select(type, t, `<>`)));
            if nops([tmp]) = 1 then
                tmp := {}, tmp;
            end if;

            Res1 := eval(seq([tmp[2*i-1],
                           map(`*`, collect(tmp[2*i]+A, x),
                                  eval_spec2(1/W[2*j], tmp[2*i-1], params)) ][],
                                                i = 1 .. nops([tmp])/2), A = 0);
            Res2 := seq([tmp[2*i-1],
                     normal(tmp[2*i]/eval_spec2(W[2*j], tmp[2*i-1], params))][], 
                                                        i = 1 .. nops([tmp])/2);

            Res := Res, seq([tmp[2*i-1],
                             `if`(length(Res1[2*i]) < length(Res2[2*i]),
                                     Res1[2*i], 
                                         Res2[2*i])][], i = 1 .. nops([tmp])/2);
        end if;
    end do;
    VarsZeroSystem := OVZS;
    return remove(`=`, [Res], {})[];
end proc:
################################################################################
UniversalDenominator := proc()
local reqn, y, x, q, params, SS1, SS2, info, QEq, psi, tmp,
      a, f, SS, wq, wx, rho, d, w, W, u, Res, i, p,
      A, B, k, c, C, j, OVZS;


    OVZS := VarsZeroSystem;
    try reqn, y, x, q, params, SS1, SS2 := check_arguments(args[])
    catch: error 
    end try;

    if SS1  = {1 = 0} or SS2 = {0 = 0} then
        VarsZeroSystem := OVZS;
        return NULL;
    end if;

    if reqn = 0 then 
        VarsZeroSystem := OVZS;
        return make_return('ANY_POLYNOMIAL', SS1, SS2, params);
    end if;

    info := eval(op(3, reqn));
    QEq := op([1,1], reqn);
    if not has(QEq, params) then
        Res := QDifferenceEquations:-UniversalDenominator(QEq, y(x));
        VarsZeroSystem := OVZS;
        return make_return(Res, SS1, SS2, params);
    end if;

    a := info['coeffs'][2 .. -1]; 
    f := -info['coeffs'][1]; 

    psi := collect(info['coeffs'][-1], VarsZeroSystem, 'distributed'); 
    SS := ZeroSystem(psi, VarsZeroSystem, params);
    if SS <> {1 = 0} then
        Res := procname(reqn, 'parameters' = params, 'S1' = SS1 union SS, 
                                                                    'S2' = SS2);
        SS, SS2 := check_pair_of_systems(SS1, add_to_S2(SS, SS2), 
                                                                params)[1..2][];
        if SS = {1 = 0} or SS2 = {0 = 0} then
            VarsZeroSystem := OVZS;
            return Res;
        end if;
        if SS <> SS1 then 
            VarsZeroSystem := OVZS;
            return Res, procname(reqn, 'parameters' = params, 
                                                         'S1' = SS, 'S2' = SS2);
        end if;
    else
        Res := NULL;
    end if;

    wq := max(map(degree, a, q)); 
    wx := max(map(degree, a, x)); 
    rho := info['order'];
    d := rho*wx^2 + 2*wx*wq;
    w := max(degree(f, x), wq);
    tmp := indets({a[-1], a[1]}, 'name') intersect params;
    if tmp = {} then
        u := q_DENOM(a[-1], a[1], rho, q, x);
        W := x^w*u;
        VarsZeroSystem := OVZS;
        return Res, make_return(W, SS1, SS2, params);
    elif nops(tmp) > 1 then
        u := mul(eval(a[-1], x = q^(-rho-i)*x), i = 0 .. d);
        W := x^w*u;
        VarsZeroSystem := OVZS;
        return Res, make_return(W, SS1, SS2, params);
    end if;

    tmp := tmp[1];
    SS := select(depends, SS1, tmp);
    if nops(SS) > 1 or nops(indets(SS) intersect params) > 1 then 
         error "unimplemented"
    end if;
    if SS = {} then p := 0 else p := lhs(SS[1]) end if;

    k := ldegree(a[-1], x);
    A := evala(Normal(eval(quo(a[-1], x^k, x), x = x*q^(-rho))));
    k := ldegree(a[1], x);
    B := evala(Normal(quo(a[1], x^k, x)));
    u := 1;
    for k from d by -1 to 0 while degree(A, x) > 0 or degree(B, x) > 0 do
        c := ParametricGCRD2(A, eval(B, x = x*q^k), x, p, tmp);
        if op(0,c) <> 'piecewise' then
            u := u*mul(eval(c, x = x*q^(-i)), i = 0 .. k);
            A := radnormal(A/c);
            c := eval(c, x = x*q^(-k));
            B := radnormal(B/c);
        else
            if nops(c)::'odd' then
                C := op(-1,c);
                u := u*mul(eval(C, x=x*q^(-i)), i = 0 .. k);
                A := radnormal(A/C);
                C := eval(C, x = x*q^(-k));
                B := radnormal(B/C);
            end if;
            for j to iquo(nops(c),2) do
                Res := Res, procname(reqn, 'parameters' = params, 
                                   'S1' = SS1 union {op(2*j-1,c)}, 'S2' = SS2 );
                SS, SS2 := check_pair_of_systems(SS1, 
                                                  add_to_S2({op(2*j-1,c)}, SS2), 
                                                                params)[1..2][];
                if SS = {1 = 0} or SS2 = {0 = 0} then
                    VarsZeroSystem := OVZS;
                    return Res;
                end if;
                if SS <> SS1 then 
                    VarsZeroSystem := OVZS;
                    return Res, procname(reqn, 'parameters' = params, 
                                                         'S1' = SS, 'S2' = SS2);
                end if;
            end do;
            if SS1  = {1 = 0} or SS2 = {0 = 0} then
                VarsZeroSystem := OVZS;
                return Res;
            end if;
        end if;  
    end do;

    W := x^w*u;
    VarsZeroSystem := OVZS;
    return Res, make_return(W, SS1, SS2, params);
  
end proc:
################################################################################
q_DENOM := proc( an, a0, n, q, x ) 
local u, A, B, c, k, i, h_s;
      
    A := eval(an,  x=x*q^(-n));
    B := a0;
    h_s := sort([op(QDifferenceEquations:-QDispersion(B, A, q, x))]);
    u := 1 ;
    for k from nops(h_s) by -1 to 1 do
       c := evala(Gcdex(A, eval(B, x = x*q^(h_s[k])),x));
       u := u*mul(eval(c, x=x*q^(-i)), i = 0 .. h_s[k]);
       A := radnormal(A/c);
       B := radnormal(B/eval(c, x=x*q^(-h_s[k])));
    end do;
    return u;
end proc: 
################################################################################
AccurateQSummation := proc(L, Q::name, x::name, q::name) 
local A, params, opts, LL, ll, AdjA, L1, eq, y, sol, Res, r, R, i, OVZS; 

    OVZS := VarsZeroSystem;
    if nargs > 4 then
        if args[5]::'name' then
            params := {args[5]};
            opts := [args[6..-1]];
        else
            opts := [args[5..-1]];
        end if;
        hasoption(opts, 'parameters' = {'list','set'}('name'), 'params', 
                                                                        'opts');
        if opts <> [] then
            error "undetermined option %1", opts[1];
        end if;
    else
        params := {};
    end if;
    if Extended then
        VarsZeroSystem := VarsZeroSystem union {x} minus {q};
    else
        VarsZeroSystem := VarsZeroSystem union {x, q};
    end if;

    if not L::'polynom'('anything', Q) then
        error "the 1st argument are expected polynomial in %1", Q; 
    elif not L::'polynom'('polynom'('anything', {x, params[]}), Q) then
        error "only polynomial in %1 coefficients are handled", {x, params[]};
    elif not L::'polynom'('polynom'('ratpoly'('algnum', q), 
                               VarsZeroSystem minus {q} union params), Q) then
        error "only polynomial over ratpoly(algnum, %1) coefficients \
                                                                are handled", q;
    end if;

    A := OreTools:-SetOreRing([x, q], 'qshift');
    AdjA := OreTools:-AdjointRing(A); 

    LL := OreTools:-Converters:-FromPolyToOrePoly(L, Q); 
    L1 := OreTools:-AdjointOrePoly(LL, A);

    eq := OreTools:-Converters:-FromOrePolyToLinearEquation(L1, y, AdjA) = 1;
    sol := RationalSolution(eq, y(x), 'parameters' = params);
    if sol = NULL then
       VarsZeroSystem := OVZS;
       return NULL;
    end if;
    if nops([sol]) = 1 then
        sol := {}, sol;
    end if;

    Res := NULL;
    for i to nops([sol])/2 do
       r := sol[2*i];
       ll := eval_spec2(LL, sol[2*i-1], params);
       R := OreTools:-Quotient['left'](OreTools:-Add('OrePoly'(1), 
                        OreTools:-ScalarMultiply(-r, ll)), 'OrePoly'(-1, 1), A);
       Res := Res, sol[2*i-1], OreTools:-Converters:-FromOrePolyToPoly(R, Q);
    end do:
    VarsZeroSystem := OVZS;
    return remove( `=`, [Res], {})[];
end proc:
################################################################################
QZeilberger := proc(s1::ratpoly, s2::ratpoly, 
                    x1::name, x2::name, q::name, Q_n::name,
                    d::nonnegint)
local N, Cn, K, Ck, reqn, y, params, SS1, SS2, S, Sol, V, cnd, info, opts,
      p, r, s, t, a, i, j, k, P, vars, A, eq, l, psi, res, tmp,
      Rd, Eq, c, fvars, Ll, G, Res, SS, OVZS;

    OVZS := VarsZeroSystem;
    if nargs > 7 then
        if args[8]::'name' then
            params := {args[8]};
            opts := [args[9..-1]];
        else
            opts := [args[8..-1]];
        end if;
        hasoption(opts, 'parameters' = {'list','set'}('name'), 'params', 
                                                                        'opts');
        params := params;
    else
        params := {};
    end if;
    if Extended then
        VarsZeroSystem := VarsZeroSystem union {x1, x2} minus {q};
    else
        VarsZeroSystem := VarsZeroSystem union {x1, x2, q};
    end if;
    if not {s1, s2}::'set'('ratpoly'('polynom'('anything', params), {x1, x2})) 
    then
        error "rational fucntions in %1 are handled", {x1, x2, params[]};
    elif not {s1, s2}::'set'('ratpoly'('ratpoly'('algnum', q), 
                               VarsZeroSystem minus {q} union params)) then
        error "rational fucntions over ratpoly(algnum, %1) are handled", q;
    end if;

    K := x1;
    N := x2;
    Cn := radnormal(s2);
    Ck := radnormal(eval(s1, K = K/q));

    try reqn, y, K, q, params, SS1, SS2 := 
                  check_arguments(y(q*K)-eval(Ck,K=K*q)*y(K), y(K), args[8..-1])
    catch: error
    end try;

    if SS1 = {1 = 0} or SS2 = {0 = 0} then 
        VarsZeroSystem := OVZS;
        return NULL;
    end if;

    Cn := eval_spec2(Cn, SS1, params);
    Ck := eval_spec2(Ck, SS1, params);

    info := eval(op(3, reqn));
    SS := ZeroSystem(info['coeffs'][3], VarsZeroSystem, params);
    Res := NULL;
    if SS <> {1 = 0} then
#        Res := procname(s1, s2, x1, x2, q, Q_n, d, 'parameters' = params, 
#                                               'S1' = SS1 union SS, 'S2' = SS2);
        SS, SS2 := check_pair_of_systems(SS1, add_to_S2(SS, SS2), 
                                                                params)[1..2][];
        if SS = {1 = 0} or SS2 = {0 = 0} then
            VarsZeroSystem := OVZS;
            return Res;
        end if;
        if SS <> SS1 then 
            VarsZeroSystem := OVZS;
            return Res, procname(s1, s2, x1, x2, q, Q_n, d, 
                                  'parameters' = params, 'S1' = SS, 'S2' = SS2);
        end if;
    else
        Res := NULL;
    end if; 

    SS := ZeroSystem(info['coeffs'][2], VarsZeroSystem, params);
    if SS <> {1 = 0} then
#         Res := Res, procname(s1, s2, x1, x2, q, Q_n, d, 
#                        'parameters' = params, 'S1' = SS1 union SS, 'S2' = SS2);
        SS, SS2 := check_pair_of_systems(SS1, add_to_S2(SS, SS2), 
                                                                params)[1..2][];
        if SS = {1 = 0} or SS2 = {0 = 0} then
            VarsZeroSystem := OVZS;
            return Res;
        end if;
        if SS <> SS1 then 
            VarsZeroSystem := OVZS;
            return Res, procname(s1, s2, x1, x2, q, Q_n, d, 
                                  'parameters' = params, 'S1' = SS, 'S2' = SS2);
        end if;
    end if;

    p := numer(Cn); r := denom(Cn);
    s := eval(numer(Ck), K = K*q); 
    t := denom(Ck);

    P := a[0];
    vars := [seq(a[i], i = 0 .. d)];
    for i from 1 to d do
        Rd := eval_spec2(radnormal(lcoeff(P, vars[i])/eval(r, N = N*q^(i-1))),
                                                                   SS1, params);
        P := eval_spec2(P*denom(Rd)+a[i]*eval(p, N = N*q^(i-1))*numer(Rd),
                                                                   SS1, params);
    end do;
    Eq := Equation(P, lcoeff(P, a[0]), s, t, K, q, N, params, SS1, SS2);
    for i to nops([Eq])/2 do
        if Eq[i*2-1]::`=` then 
            Eq := Eq[1..i*2-2], {(lhs-rhs)(Eq[i*2-1]) = 0}, Eq[i*2..-1] 
        end if;
        if Eq[i*2-1] <> SS1 then
            try Res := Res, procname(s1, s2, x1, x2, q, Q_n, d, 
                 'parameters' = params, 'S1' = Eq[i*2-1] union SS1, 'S2' = SS2);
            catch:
            end try;
            SS, SS2 := check_pair_of_systems(SS1, add_to_S2(Eq[i*2-1], SS2), 
                                                               params)[1..2][];
            if SS = {1 = 0} or SS2 = {0 = 0} then
                VarsZeroSystem := OVZS;
                if Res <> NULL then 
                    return Res;
                else 
                    error "No Z-pair of order %1 was found", d; 
                end if;
            end if;
            if SS <> SS1 then 
                VarsZeroSystem := OVZS;
                try Res := Res, procname(s1, s2, x1, x2, q, Q_n, d, 
                                  'parameters' = params, 'S1' = SS, 'S2' = SS2);
                catch:
                end try;
                if Res <> NULL then 
                    return Res;
                else 
                    error "No Z-pair of order %1 was found", d; 
                end if;
            end if;
            next;
        end if;
        l := UpperBound(Eq[i*2][1], Eq[i*2][2], Eq[i*2][3], K, q);
        fvars := [seq(c[j], j = 0 .. l)];

        Sol := add(fvars[j+1]*K^j, j = 0 .. l); 
        V := Vector([fvars[],vars[], -1]); 
        if l >= 0 then
            res := Eq[i*2][1]*eval(Sol, K = K*q) - Eq[i*2][2]*Sol - Eq[i*2][3];
        else 
            res := Eq[i*2][3]; 
        end if;
        res := collect(expand(res), K); 
        S := [coeffs(res, K)]; 
        S := map(primpart, S, {params[], fvars[], vars[], N});
        S := eval_spec2(S, SS1, params);

        A := LinearAlgebra:-GenerateMatrix(S, [fvars[],vars[]], 'augmented');
        A := LinearAlgebra:-GaussianElimination(A);

        cnd := NULL;
        for k from LinearAlgebra:-Rank(A) by -1 to 1 do
            eq := radnormal(A[k].V);
            tmp := denom(eq);
            tmp := eval_spec2(tmp, SS1, params);
            if tmp = 0 then 
                tmp := ZeroSystem(denom(eq), VarsZeroSystem, params);
                SS := select(el -> lhs(el)::`*`, tmp);
                if SS <> {} then
                    tmp := op(1, lhs(SS[1]));
                    try Res := Res, procname(s1, s2, x1, x2, q, Q_n, d, 
                                        'parameters' = params, 
                                        'S1' = {tmp = 0} union SS1, 'S2' = SS2);
                        Res := Res, procname(s1, s2, x1, x2, q, Q_n, d, 
                                        'parameters' = params, 
                                        'S1' = SS1, 'S2' = SS2 union {tmp = 0});
                    catch:
                    end try;
                    VarsZeroSystem := OVZS;
                    if Res <> NULL then 
                        return Res;
                    else 
                        error "No Z-pair of order %1 was found", d; 
                    end if;
                end if;
                error "unimplemented" 
            end if;        
            SS := ZeroSystem(tmp, VarsZeroSystem, params);
            if SS <> {1 = 0} then
                try Res := Res, procname(s1, s2, x1, x2, q, Q_n, d, 
                                           'parameters' = params, 
                                               'S1' = SS1 union SS, 'S2' = SS2);
                catch:
                end try;
                SS, SS2 := check_pair_of_systems(SS1, add_to_S2(SS, SS2), 
                                                                params)[1..2][];
                if SS = {1 = 0} or SS2 = {0 = 0} then
                    VarsZeroSystem := OVZS;
                    if Res <> NULL then 
                        return Res;
                    else 
                        error "No Z-pair of order %1 was found", d; 
                    end if;
                end if;
                if SS <> SS1 then 
                    VarsZeroSystem := OVZS;
                    try Res := Res, procname(s1, s2, x1, x2, q, Q_n, d, 
                                                'parameters' = params, 
                                                         'S1' = SS, 'S2' = SS2);
                    catch:
                    end try;
                    if Res <> NULL then 
                         return Res;
                    else 
                        error "No Z-pair of order %1 was found", d; 
                    end if;
                end if;
            end if;
            eq := numer(eq);
            eq := eval_spec2(eq, SS1, params);
            if eq = 0 then next end if;
            j := min(map2(op, 1, indets(eq, c['integer'])));
            if j = infinity then
                j := min(map2(op, 1, indets(eq, a['integer'])));
                if j <> infinity then j := a[j] end if;
            else
                j := c[j];
            end if;
            if j = infinity then
                psi := eq;
            else
                psi := coeff(eq, j);
            end if;
            SS := ZeroSystem(psi, VarsZeroSystem, params);
            if SS <> {1 = 0} then
                try Res := Res, procname(s1, s2, x1, x2, q, Q_n, d, 
                                           'parameters' = params, 
                                               'S1' = SS1 union SS, 'S2' = SS2);
                catch:
                end try;
                SS, SS2 := check_pair_of_systems(SS1, add_to_S2(SS, SS2), 
                                                                params)[1..2][];
                if SS = {1 = 0} or SS2 = {0 = 0} then
                    VarsZeroSystem := OVZS;
                    if Res <> NULL then 
                        return Res;
                    else 
                        error "No Z-pair of order %1 was found", d; 
                    end if;
                end if;
                if SS <> SS1 then 
                    VarsZeroSystem := OVZS;
                    try Res := Res, procname(s1, s2, x1, x2, q, Q_n, d, 
                                                'parameters' = params, 
                                                         'S1' = SS, 'S2' = SS2);
                    catch:
                    end try;
                    if Res <> NULL then 
                        return Res;
                    else 
                        error "No Z-pair of order %1 was found", d; 
                    end if;
                end if;
            end if;
            if j = infinity then
                Sol := 0;
            else
                tmp := solve(eval(eq, {cnd}), j);
                tmp := eval_spec2(tmp, SS1, params);
                cnd := cnd, j = tmp;
            end if;
        end do;
        if eval(a[d], {cnd}) = 0 then 
            VarsZeroSystem := OVZS;
            if Res <> NULL then 
                return Res;
            else 
                error "No Z-pair of order %1 was found", d; 
            end if;
        end if;
        cnd := adjust({seq(j = eval(j, {cnd}), j = vars),
                       seq(j = eval(j, {cnd}), j = fvars)}, vars[d+1], params);
        Ll := add(vars[j+1]*Q_n^j, j = 0 .. d);
        Ll := radnormal(eval(Ll, cnd), 'expanded');
        G := eval(Sol*Eq[i*2][4]*denom(Ll), cnd);
        Ll := collect(numer(Ll), Q_n);
        Res := Res, make_return([Ll,radnormal(G)], SS1, SS2, params);

    end do;
    VarsZeroSystem := OVZS;
    if Res <> NULL then 
        return remove( `=`, [Res], {})[];;
    else 
        error "No Z-pair of order %1 was found", d; 
    end if;
end proc:
################################################################################
Equation := proc(P, R, s, t, y, q, N, params, SS1, SS2)
local C, g1, g0, den, v, u, u0, u1, w, i, u_i, f, Res, C_i, g0_i, g1_i; 
     
    g1 := R*s; 
    g0 := R*t; 
    C := t;
    g1 := radnormal(eval(g1, y = y/q)); 
    den := denom(g1);
    v := gcd(g0*den, numer(g1), 'g0', 'g1'); 
    g1 := eval(g1, y = y*q); 
    C := C*den/v;
    u := UniversalDenominator(g1*f(q*y) - g0*f(y), f(y),
                                 'parameters' = params, 'S1' = SS1, 'S2' = SS2);
    if nops([u]) = 1 then
        u := {}, u;
    end if;
    if nops([u]) = 0 then 
         error "unimplemented"
    end if;
    Res := NULL;
    for i to nops([u])/2 do
        u_i := numer(radnormal(u[2*i])); 
        C_i := eval(C, select(type, u[2*i-1], `=`))/u_i;
        gcd(u_i, eval(g0, select(type, u[2*i-1], `=`)), 'u0', 'g0_i'); 
        gcd(eval(u_i, y = y*q), eval(g1, select(type, u[2*i-1], `=`)),
                                                                  'u1', 'g1_i');
        w := lcm(u1, u0);
        Res := Res, u[i*2-1], 
                    [g1_i*quo(w, u1, y), g0_i*quo(w, u0, y), P*w, C_i];
    end do;
    return Res;
end proc:
################################################################################
UpperBound := proc(g1, g0, h, y, q)
local d, d1, d0, c, a, f, wq, l;
    d1 := degree(g1, y); 
    d0 := degree(g0, y); 
    d := degree(h, y);
    if d1 = d0 then
        c := radnormal(lcoeff(g0, y)/lcoeff(g1, y));
        if c::'polynom'('anything', q) then 
            if c = q^degree(c, q) then 
                return max(degree(c, q), d-d0); 
            end if;
        end if;
    end if;

    a := map(collect, [expand(g1), expand(g0)], {q, y}, 'distributed'); 
    f := collect(expand(h), {q, y}, 'distributed');

    wq := max(map(degree, a, q)); 
    l := max(wq, degree(f, y));
    return l;

end proc:
################################################################################
adjust := proc(sols, x, params)
local freevars, s, t, xr, j, ev;

    freevars := []; 
    s := sols;
    for t in s do
        if lhs(t) = rhs(t) and not member(lhs(t), params) then
            freevars := [freevars[],lhs(t)];
        end if; 
        if lhs(t) = x then 
            xr := rhs(t); 
        end if;
    end do;
    if Testzero(xr) then xr := 1 end if;
    for t in freevars do
        for j from 0 while Testzero(eval(xr, t = j)) do end do;
        ev := t = j;
        xr := eval(xr, ev);
        s := eval(s minus {t = t}, ev) union {ev};
    end do; 
    return map(y -> lhs(y) = radnormal(rhs(y)), s);
end proc:
################################################################################
ParametricGCRD2 := proc(A1, B1, x::name, p1, Par::name)
local Var, A, B, p, r, w, h, res, S, sr, Lc1, Lc2, Lc12, C1, C2,
      s1, s2, F1, F2, i, tmp;

    A := collect(numer(radnormal(A1)), x);
    B := collect(numer(radnormal(B1)), x);
    p := primpart(collect(p1, Par, radnormal), Par);

    if degree(p, Par) = 0 then
        return(Piecewise())
    end if;

    res := NULL;
    if p = 0 and degree(A, x) > 0 and degree(B, x) > 0 then
        p := UpperBoundOfCondition(A,B,x,Par);
        if degree(p, Par) = 0 then
            return(Piecewise(1))
        elif p <> 0 then
            res := 1;
        end if;
    end if;

    if p <> 0 then
        A := rem(A, p, Par);
        B := rem(B, p, Par);
    end if;

    Var := VarsZeroSystem minus {x};
    if A = 0 then
        return Piecewise(Filter(p, B, x, Var, Par))
    end if;

    if B = 0 then
        return Piecewise(Filter(p, A, x, Var, Par))
    end if;

    if degree(A, x) < degree(B, x) then
       i := A; A := B; B := i;
    end if;

    C1 := content(A, [x, Var[]]);
    C2 := content(B, [x, Var[]]);

    Lc1 := content(lcoeff(A, x), Var);
    Lc2 := content(lcoeff(B, x), Var);

    if p <> 0 then
        C1 := gcdex(C1, p, Par);
        C2 := gcdex(C2, p, Par);
        Lc1 := gcdex(Lc1, p, Par);
        Lc2 := gcdex(Lc2, p, Par);
    else
        C1 := Sqrfree(C1, Par);
        C2 := Sqrfree(C2, Par);
        Lc1 := Sqrfree(Lc1, Par);
        Lc2 := Sqrfree(Lc2, Par)
    end if;

    Lc12 := gcdex(Lc1, Lc2, Par);
    divide(C1, gcdex(C1, Lc12, Par), 's1');
    divide(C2, gcdex(C2, Lc12, Par), 's2');
    F1 := Filter(s1, B, x, Var, Par);
    F2 := Filter(s2, A, x, Var, Par);
    if degree(Lc12, Par) = 0 then
        res := F1, F2, res;
    else
        res := F1, F2, op(procname(A, B, x, Lc12, Par)), res;
    end if;

    if p <> 0 then
        divide( p, gcdex(p, Lc12*s1*s2, Par), 'p');
    end if;

    S := FractionFree(A, B, x);
    S := [seq(S[2][i], i = 3 .. S[1])];

    if S=[] then
        tmp := Filter(p, B, x, Var, Par);
        if nops([tmp]) = 1 then
            return Piecewise(res, tmp);
        else
            return Piecewise(tmp, res);
        end if;
    end if;

    r := content(lcoeff(S[-1],x), Var);
    if p = 0 then
        r := Sqrfree(r, Par);
        divide(r, Lc12*s1*s2, 'r');
    else
        r := gcdex(r, p, Par)
    end if;
    divide(p, r, 'p');
    tmp := Filter(p, S[-1], x, Var, Par);
    if nops([tmp]) = 1 then
        res := res, tmp;
    else
        res := tmp, res;
    end if;

    S := S[1..-2];

    while S <> [] do
        sr := S[-1];
        S := S[1..-2];
        h := gcdex(r, content(lcoeff(sr, x), Var), Par);
        divide(r, h, 'w');
        r := h;
        res := Filter(w, sr, x, Var, Par), res;
    end do;

    res := Filter(r, B, x, Var, Par), res;

    Piecewise(res)

end proc:
################################################################################
UpperBoundOfCondition := proc(A, B, x::name, Par::name)
local Var, d1, d2, max_degree, M, i, j, det_i,
      degree_i, tmp, X;

    Var := VarsZeroSystem minus {x};
    if Var = {} then Var := X else Var := Var[1] end if;

    d1 := degree(A, x);
    d2 := degree(B, x);
    max_degree := degree(A, Var)*(d2-1) +
                  degree(B, Var)*(d1-1);

    M := Matrix(
           [seq([seq(coeff(collect(x^j*B, x), x, i), i = 0 .. j+d2)],
                j = 0 .. d1-1),
            seq([seq(coeff(collect(x^j*A, x), x, i), i = 0 .. j+d1)],
                j = 0 .. d2-1)] );

    det_i := radnormal(LinearAlgebra['Determinant'](eval(M, Var=0)));
    for i to max_degree-1 while det_i = 0 do
       det_i := radnormal(LinearAlgebra['Determinant'](eval(M, Var=i)));
    end do;

    if det_i = 0 then return 0 end if;

    det_i := Sqrfree(det_i, Par);
    degree_i := degree(collect(det_i, Par, radnormal), Par);
    if degree_i = 0 then return 1 end if;

    j := 0;
    for i do
        det_i := gcdex(det_i, 
                       LinearAlgebra['Determinant'](eval(M, Var=i)), Par);
        tmp := degree(collect(det_i, Par, radnormal), Par);
        if tmp = 0 then
           return 1
        elif tmp = degree_i and j < 2 and tmp >= 10 then
           j := j+1;
        elif tmp = degree_i then
           return det_i
        else
           degree_i := tmp;
           j := 0;
        end if;
    end do;

    return det_i;
end proc:
################################################################################
FractionFree := proc(Ore1, Ore2, x::name)
local P1, P2, d1, d2, deg, S, d, l, a, b, i, e, R,  
      init, m, c1, c2, q;

    d1 := degree(Ore1, x); 
    d2 := degree(Ore2, x);
    if d1 >= d2 then
        P1 := Ore1; P2:=Ore2;
    else
        P1 := Ore2; P2 := Ore1;
        deg := d1; d1 := d2; d2 := deg; 
    end if;

    S := array(1 .. d2+2);
    d := array(1 .. d2+2);
    l := array(1 .. d2+2);
    a := array(1 .. d2+2);
    b := array(1 .. d2+2);
    S[1] := P1; S[2] := P2;
    d[1] := d1; d[2] := d2;
    l[2] := d[1] - d[2] + 1;

    m := 2; i := 2;
    while true do 
        i := i+1;
        if i = 3 then
            e := (-1)^l[2];
        else
            e := expand((-1)^l[i-1]*b[i-2]^(l[i-1]-1)*a[i-2]);
        end if;
        R := prem(S[i-2], S[i-1], x, 'init');
        if R = 0 then
            m := i-1; 
            break;
        end if;

        d[i] := degree(R, x);
        S[i] := radnormal(R/e);
        if d[i] = 0 then
            m := i; break;  
        else
            l[i] := d[i-1] - d[i] + 1;
            a[i-1] := lcoeff(S[i-1], x);
            if i = 3 then
                b[2] := expand(a[2]^(l[2]-1));            
            else
                if l[i-1] = 2 then
                    b[i-1] := a[i-1];
                else
                    c1 := expand(a[i-1]^(l[i-1]-1));
                    c2 := expand(b[i-2]^(l[i-1]-2));
                    divide(c1, c2, 'q');
                    b[i-1] := q;
                end if;
            end if;            
        end if;
    end do;
    return [m, S];
end proc:  
################################################################################
Piecewise := proc()
    if nargs > 1 then
        piecewise(args)
    elif nargs = 1 then 
        args
    else
        NULL
    end if;
end proc:
################################################################################
Sqrfree := proc(p, v::name)
local i;
    collect(mul(i[1], i = sqrfree(p, v)[2]), v)
end proc:
################################################################################
Filter := proc(p, d, x::name, Var::set(name), Par::name)
local q, con;

    if depends(p, Var) then return NULL end if;
    q := d;

    if degree(p,Par) = 0 then
        return NULL
    elif q <> 0 and p <> 0  then
        q := rem(q, p, Par);
        con := content(lcoeff(q, x), Var);
        q := q/con;
    end if;

    `if`(p <> 0, p = 0, NULL), map(collect, primpart(q, x), [Var[],Par])

end proc:
################################################################################
check_arguments := proc()
local  reqn, info, opts, params, SS2, SS1, y, x, q, QEq, a, slv, i, C;

    if args[1] <> (0=0) then
        if type(args[1], 'specfunc'('anything', 'QESol')) then
            reqn := args[1];
            opts := [args[2..-1]];
        else
            try reqn := QDifferenceEquations:-QECreate(args[1], args[2]);
            catch:
                error
            end try:
            opts := [args[3..-1]];
        end if;
        info := eval(op(3, reqn));
        q := info['q_par'];
        x := info['vars'];
        y := info['functions'];
        if Extended then
            VarsZeroSystem := VarsZeroSystem union {x} minus {q};
        else
            VarsZeroSystem := VarsZeroSystem union {x, q};
        end if;
    else
        reqn := 0;
        opts := [args[3..-1]];
    end if;

    if opts = [] then
        params := {};
    elif opts[1] :: 'name' then 
        params := {opts[1]};
        opts := opts[2..-1];
    else
        params := {};
    end if; 

    if not opts::'list'('name' = 'anything') then
        error "undetermined argument %1", 
                             remove(type, opts, 'list'('name' = 'anything'))[1];
    end if;
    hasoption(opts, 'parameters'={'list','set'}('name'), 'params', 'opts');
    params := params;
    VarsZeroSystem := VarsZeroSystem minus params;

    SS2 := {};
    hasoption(opts, 'S2'={'list','set'}(`=`), 'SS2', 'opts');
    SS2 := SS2;

    SS1 := {};
    hasoption(opts, 'S1'={'list','set'}(`=`), 'SS1', 'opts');
    SS1 := SS1;

    if opts <> [] then
        error "undetermined option %1", opts[1];
    end if;

    if reqn <> 0 then
        if info['linear'] = false then 
            error "only linear equations are handled"; 
        elif nops({y}) > 1 then
            error "only single variables are handled";
        elif nops({x}) > 1 then
            error "only equations in one variable are handled";
        elif not type(info['coeffs'],
                     'list'('polynom'('anything', {x,params[]}))) then
            error "only polynomial in %1 coefficients are handled", 
                                                                  {x, params[]};
        elif not type(info['coeffs'],  
                    'list'('polynom'('ratpoly'('algnum', q), 
                                indets(info['coeffs'], 'name') minus {q}))) then
            error "only polynomial over ratpoly(algnum,%1) coefficients \
                                                                are handled", q;
        end if;
    end if;

    try SS1, SS2, slv := check_pair_of_systems(SS1, SS2, params);
    catch: 
         error
    end try;
    if reqn <> 0  then
        a := eval_spec2(info['coeffs'], slv, params);
        QEq := eval(primpart(add(a[i]*C[i], i = 1 .. nops(a)), 
                                    {x, params[], seq(C[i], i = 1 .. nops(a))}),
                    {seq(C[i+2] = y(x*q^i), i = 0 .. nops(a)-1), C[1] = 1});
        if depends(QEq, y) then
            reqn := QDifferenceEquations:-QECreate(QEq, y(x));
            op(3,reqn)['coeffs'] := map(collect, eval(op(3, reqn)['coeffs']), 
                                                 VarsZeroSystem, 'distributed'); 
        else
            reqn := QEq;
        end if;
    end if;
    return reqn, y, x, q, params, SS1, SS2;
end proc:
################################################################################
ZeroSystem := proc(plnm, vars, vars1)
local SS;
    SS := collect(expand(plnm), vars, 'distributed'); 
    SS := {coeffs(SS, vars)}; 
    SS := map(`=`, SS, 0);
    check_system(SS, vars1);
end proc:
################################################################################
check_system := proc(SS, vars, name)
local S, t, S1, St;

    S := map(lhs-rhs, SS);
    S := map(collect, map( expand, map(primpart, S, vars)), vars);
    if S = {0} then
        return {0 = 0};
    end if;
    S := remove(Testzero, S);
    if not S::'set'('polynom'('anything', vars)) then
        error "only algebraic equations for %1 over radalgnum coefficients \
                                                 in %2 are handled", vars, name;
    end if;
    S := map(factors, S);
    S := map2(op, 2, S);
    S := map(el -> map2(op, 1, el), S);
    S := map(convert, S, `*`);
    S1 := S;
    for t in S do 
       St := select(divide, S minus {t}, t);
       if St <> {} then
           S1 := S1 minus St;
       end if;
    end do;
    return map(`=`, S1, 0);
end proc:
################################################################################
check_pair_of_systems := proc(ss1, ss2, params)
local tmp, SS1, SS2, slv;

    try SS1 := check_system(ss1, params, 'S1'); 
    catch: 
        error
    end try;
    try SS2 := check_system(ss2, params, 'S2');
    catch: 
        error
    end try;

    if member(1 = 0, SS1) or member(0 = 0, SS2 ) then 
        return {1 = 0}, {0 = 0}, {} 
    end if;

    if nops(SS2) = 1 then
        if (indets(SS1) intersect indets(SS2) intersect params) <> {} then
            tmp := map(el -> lhs(el)/gcd(lhs(el), lhs(SS2[1])) = 0, SS1);
            if SS1 <> tmp then
                return procname(tmp, SS2, params);
            end if;
        end if;
    end if;

    SS1 := remove(el -> el = (0 = 0),SS1); 
    SS2 := remove(el -> el = (1 = 0),SS2);

    tmp := select(el -> nops(indets(el) intersect params) = 1, SS1);
    if tmp <> {} then
        SS1, SS2, slv := procname( 
                               eval_spec2(SS1 minus {tmp[1]}, {tmp[1]}, params),
                                     eval_spec2(SS2, {tmp[1]}, params), params);
        if member(1 = 0, SS1) or member(0 = 0, SS2 ) then 
            return {1 = 0}, {0 = 0}, {} 
        end if;
        slv := slv union {tmp[1]};
        SS1 := SS1 union {tmp[1]};
    else
        slv := {};
    end if;

    return SS1, SS2, slv;
end proc:
################################################################################
make_return := proc(Res, S1, S2, params)
local i, ind, SS1, SS2, g, res;

    if S1 = {} and S2 = {} then return Res end if;
    if not OutputShort then
        SS1 := select(el -> lhs(el)::`*`, S1);
        if SS1 <> {} then
            g := op(1, lhs(SS1[1]));
            res := procname(Res, 
                      check_pair_of_systems({g = 0} union S1, S2, params)[1..2], 
                                                                        params);
            return res, procname(Res, 
                      check_pair_of_systems(S1, S2 union {g = 0}, params)[1..2],
                                                                        params);
        end if;
    end if; 
    SS1 := NULL;
    for i in S1 do
        ind := indets(i, 'name') intersect params;
        if nops(ind) = 1 then
            if not OutputShort or lhs(i)::'linear'(ind[1]) then
                SS1 := SS1, ind[1] = RootOf(i, ind[1]);
            else
                SS1 := SS1, i;
            end if;
        else
            SS1 := SS1, i;
        end if;
    end do;
    SS2 := NULL;
    for i in eval(S2, {SS1}) do
        ind := indets(i, 'name') intersect params;
        if nops(ind) = 1 then
            if lhs(i)::'linear'(ind[1]) then
                SS2 := SS2, ind[1] <> RootOf(i, ind[1]);
            else
                SS2 := SS2, lhs(i) <> 0;
            end if;
        else
            SS2 := SS2, lhs(i) <> 0;
        end if;
    end do;

    if nops( {SS1, SS2} ) = 1 then
        return SS1, SS2, eval(Res, {SS1});
    elif nops({SS2}) = 0 then
        return {SS1}, eval(Res, {SS1});
    else
        return seq([{i, SS1}, eval(Res, {SS1})][], i = {SS2});
    end if;
end proc:
################################################################################
eval_spec := proc(expr1, sbst)
local expr2, indts, x, y, i;
    expr2 := rhs(sbst);
    y := op(0, lhs(sbst));
    x := op(lhs(sbst));
    indts := indets(expr1, y('anything'));
    eval(expr1, {seq(i = eval(expr2, x = op(i)), i = indts)});
end proc:
################################################################################
eval_spec2 := proc(expr, sl, vars)
local t, t1, t2, slv;
    if sl::'set' then
        slv := remove(type, sl, `<>`);
    elif sl::`=` then
        slv := {sl};
    elif sl::`<>` then
        return expr;
    end if;
    if nops(slv) = 0 then return expr end if;
    slv := map(lhs-rhs, slv);
    slv := map(`=`, slv, 0);
    
    if expr::'polynom'('anything',vars) then
        t := indets(slv[1]) intersect vars;
        if nops(t) = 1 then 
            return procname(rem(expr, lhs(slv[1]), t[1]), slv[2..-1], vars);
        else
            return expr;
        end if;
    elif expr::`=` and rhs(expr) = 0 then
        t := procname(lhs(expr), slv, vars);
        t1 := indets(slv) intersect vars;
        t2 := indets(t) intersect vars;
        if t2 minus t1 = {} and t <> 0 then
            return 1 = 0;
        else
            return t = 0;
        end if;
    else
        return map(procname, expr, slv, vars);
    end if;
    
end proc:
################################################################################
add_to_S2 := proc(SS, SS2)
local i;
    if SS2 = {} then
        SS
    else  
       {seq(map(`=`, map(`*`, map(lhs, SS2), lhs(i)), 0)[], i = SS)}
    end if;
end proc:
end module:
################################################################################