Theorem51 := proc(L, x, tau) local i, c, dA, b, r, g;
	if degree(L,tau) <> 3 then
		error "expecting an operator of order 3"
	fi;
	for i from 0 to 2 do
		c[i] := normal( coeff(L,tau,i)/coeff(L,tau,3) )
	od:
	if c[1] = 0 and c[2] = 0 then
		return "Case 3(a)"
	elif c[1] = 0 or c[2] = 0 then
		return false  # (I in Theorem 5.1(1) is not zero)
	elif normal(c[0] - subs(x=x-1, c[1]) * c[2]) = 0 then
		return "Case 3(b)"
	fi;
	dA := 1 - subs(x=x-1,c[1]) * c[2]/c[0];
	b := normal( (1 - subs(x=x-1,c[2]) * c[2]/c[1])/dA );
	b := normal(b);
	r := normal( subs(x=x-2, c[2])/ (subs(x=x-1, b)-1) );
	g := normal( ( 1-c[1]/c[0] * subs(x=x-1, c[1]/c[2]) )/dA );
	if normal(subs(x=x+1, b) - g) <> 0 then
		return false
	fi;
	[b, r]
end;


# The extended Euclidean Algorithm for D = C(x)[tau]
#
GCDEX := proc(a, b) 
	local r1, q1, s, t; 
	if b = 0 then 
		1/content(a, _Env_LRE_tau), 0
	else
		q1, r1 := LREtools[RightDivision](a, b);
		t, s := procname(b, r1);
		s, collect(t-s*q1, _Env_LRE_tau, normal)
	fi
end:

# Compute the m'th SYMMetric power of L1:
symmpower := proc(L1, m)  global x,tau ;
	local c, n1, sb1, i, P, var, j, IsZero, EQN, s, NumberFree, k, N;
	n1 := degree(L1, tau); 
	sb1 := u(x + n1) = -add(coeff(L1, tau, i)*u(x + i)/lcoeff(L1, tau), i = 0 .. n1 - 1); 
	P[0] := u(x)^m; 
	N := binomial(m + n1 - 1, n1 - 1); 
	for i to N do 
		P[i] := subs(sb1, subs(sb1, subs(x = x + 1, P[i - 1])))
	od;
	var := {seq(u(x + i), i = 0 .. n1 - 1), seq(u(x + j), j = 0 .. n1 - 1)}; 
	IsZero := collect(add(c[i]*P[i], i = 0 .. N), var, distributed); 
	EQN := {coeffs(IsZero, var)}; 
	var := {seq(c[i], i = 0 .. N)}; 
	s := solve(EQN, var); 
	NumberFree := add(`if`(lhs(i) = rhs(i), 1, 0), i = s); 
	if 1 < NumberFree then
		s := solve(s union {seq(c[k] = 0, k = N + 2 - NumberFree .. N)}, var)
	fi;
	sort(collect(primpart(subs(s, add(c[i]*tau^i, i = 0 .. N)), tau), tau, factor), tau)
end:

ReduceOrder := proc(L3) 
            local x,tau, L6, L1, ux, i, G2, G, S, vars, rel, rel1, rel2, C, pt, LG, rel3, X, a, b,i1,i2,new_ANS,ANS, bvals;
            tau := _Env_LRE_tau;
            x := _Env_LRE_x;
            # Step 1
            L6 := symmpower(L3,2); 
            # Step 2
         if degree(L6,tau) < 6 then
            i := Theorem51(L3, x, tau);
            if not type(i,list) then return i fi;
            return 1, 1, tau^2+tau+i[1], i[2]
	 fi;

            # Step 3
            L1 := LREtools[RightFactors](L6, 1);
            if nops(L1)<>1 then
		"not 2-solvable (Assuming absolute irreducibility was checked)"
            fi;
            L1 := L1[1];
            # Step 4
            ux[3] := solve(LREtools[OperatorToRecurrence](L3, u(x)), u(x + 3)); 
            ux[0] := u(x); 
           for i to 5 do ux[i] := normal(subs(x = x + 1, u(x + 3) = ux[3], ux[i - 1])); end do; 
            # Step 5
            G2 := add(a[i]*ux[i]^2, i = 0 .. 5); 
            G := add(b[i]*ux[i], i = 0 .. 2); 
            # Step 6
            S := {coeffs(collect(G2 - G^2, {seq(ux[i], i = 0 .. 2)}, distributed), {seq(ux[i], i = 0 .. 2)})}; 
            G2 := add(a[i]*tau^i, i = 0 .. 5); 
            vars := indets(G2) minus {tau}; 
            # Step 7 and 8 
            LREtools[RightDivision](subs(solve(S, vars), G2), L1)[2]; 
            rel := collect(%, {b[0], b[1], b[2]}); 
            rel1 := b[0] = X[1] - 1/2*coeff(rel, b[0], 1); 
            rel:= collect(subs(rel1, rel), {X[1], b[1], b[2]},Normalizer);
         if degree(rel,b[1]) = 2 then
            bvals := NULL;
            rel2 := b[1] = X[2] - 1/2*coeff(rel, b[1], 1)/coeff(rel, b[1], 2); 
            rel3:= b[2] = X[3];
            # Step 9
            C := collect(primpart(subs(rel2,rel3, rel), {X[1], X[2], X[3]}), {X[1], X[2], X[3]}); 
            # Step 10, Running the code in https://www.math.fsu.edu/~hboukaed/Implementations/Conic Code.txt
            pt := Conic(expand(coeff(C, X[1], 2)), expand(coeff(C, X[2], 2)), expand(coeff(C, X[3], 2))); 
            if pt = FAIL then return FAIL fi;
          for i1 in [1,-1] do for i2 in [1,-1] do
            bvals := bvals, solve({pt[1]*i1 - X[1], pt[2]*i2 - X[2], pt[3] - X[3], rel1, rel2, rel3}, {X[1], X[2], X[3], b[0], b[1], b[2]}); 
          od od;
            else
              bvals := solve({rel1, X[1], b[1]-1, b[2]}, {X[1], b[0], b[1], b[2]});
           fi;
        for S in [bvals] do
            # Step 11
            G := subs(S, add(b[i]*tau^i, i = 0 .. 2)); 
            # Step 12
            LG := LREtools[RightDivision](LREtools[LCLM](L3, G), G)[1]; 
            LG := collect(LG/lcoeff(LG, tau), tau, factor); 
            i := Theorem51(LG, x, tau);
            if type(i, list) then
		new_ANS := G, GCDEX(G,L3)[1], tau^2+tau+i[1], i[2]
            else
		new_ANS := i
            fi;
            if not assigned(ANS) or length([ANS]) > length([new_ANS]) then
		ANS := new_ANS
	    fi;
        od;
	ANS
end;