Am am sure it is based on the Euclidean algorithm for polynomials. (The modulus is in general a prime number).

Your system is numerically unstable.
Supposing that the constants are "exact", it is possible to make a mix of symbolic and numeric approaches.
The ideea is to convert the constants to rationals (for a precision of 15 digits) then solve the system symbolically wrt the variable ionic (because this is possible), then solve the univariate equation in ionic with very high precision.

Here are the results. The solution is unique [almost]. Some variables are negative.

solve( sum((1/(1+x+x^2))*y^n, n=1..infinity) = 1, parametric=true, real) assuming abs(y)<1;

    [{x = x, y = (x^2 + x + 1)/(x^2 + x + 2)}]

assuming is necessary, otherwise the series is divergent.

Download JInt-res.mw

restart;
`diff/f` := (u,x) -> f[op(procname)-1](u)*diff(u,x):

Examples

diff(f[4](x), x, x);
    f[2](x)

diff(x^2 + f[4](f[5](2*x)),  x,  x);
    2 + 4*f[2](f[5](2*x))*f[4](2*x)^2 + 4*f[3](f[5](2*x))*f[3](2*x)

f:=proc(x,y)
   local lm := op(procname);
   if not type(procname, 'indexed') or nops([lm])<>2 then error "Two indices needed" fi;
   x^lm[1] + y^lm[2]
end proc:

https://mapleprimes.com/questions/230075-Rewriting-General-Formula-#answer271028

Download quadric.mw

So, you must help Maple a bit.

Download ineqz.mw

Download cyl.mw

Here is a procedure for Conic classification (general case). It was extracted (& adapted) from a MathApp.
It will not be useful to you, unless you read about this subject.

P:=proc(a,b,c,d,e,f)  #  See ?MathApps,ConicSections
    local equ,disc,Delta, Type, Point, deg, equ2;
    equ := sort(a*x^2+b*x*y+c*y^2+d*x+e*y+f);
    Delta := 4*a*c*f-a*e^2+b*e*d-b^2*f-c*d^2;
    disc := b^2-4*a*c;
    if Delta <> 0 then

        if disc < 0 then
            if c*Delta > 0 then
                Type := "This equation has no real solutions.";
            elif b = 0 and a = c then
                Type := "A circle.";
            else
                Type := "An ellipse.";
            end if;
        elif disc = 0 then
            Type := "A parabola.";
        else
            Type := "A hyperbola.";
        end if;

    # From here on, Delta = 0 : we have the degenerate cases

    elif disc > 0 then #OK
        Type := "Two intersecting lines.";
    elif disc < 0 then #OK
        Type := "A single point.";
        Point := [(2*c*d-b*e)/disc, (2*a*e-b*d)/disc];
    elif equ = 0 then
        Type := "Every point is a solution.";
    else
        deg := degree(equ, [x,y]);
        if deg = 0 then
            Type := "This equation has no solutions.";
        elif deg = 1 then
            Type := "A line.";
            equ2 := equ;
        else
            # a,b,c not all 0
            if c = 0 then
                # c=b=e=0
                disc := d^2-4*a*f;
            else
                # if a = 0 then                
                    # a=b=d=0
                    disc := e^2-4*c*f;
                # else
                    # both discs above have same sign,
                    # so either one is sufficient
                # end if;
            end if;
            if disc > 0 then
                Type := "Two parallel lines.";
            elif disc = 0 then
                Type := "A line.";
                equ2 := select(has, factor(equ)*t, [x,y]);
                if type(equ2, `^`) and op(2, equ2) = 2 then
                    equ2 := op(1, equ2);
                end if;
            else
                Type := "This equation has no solutions.";
            end if;
        end if;
    end if;
    equ=0,'Type'=Type;
end proc:

Why do you say that the first ODE y' = p1/p2 is exact? It became exact, after you have multiplied it with the integrating factor p2.
Most ODEs are in this situation. Should all be declared exact?