# Isomorphism from D = C(x)[tau] to Dm = C(x)[tau^m]
isom_D_Dm := proc(L, m) 
    subs(_Env_LRE_tau = _Env_LRE_tau^m, _Env_LRE_x = _Env_LRE_x / m, L); 
end:

# Inverse isomorphism (from Dm to D)
isom_Dm_D := proc(L, m) 
    subs(_Env_LRE_tau = _Env_LRE_tau^(1/m), _Env_LRE_x = _Env_LRE_x  * m, L); 
end:

# section_op(L, m, 1) returns "L downarrow m", the smallest operator in Dm that is divisible by L.
# section_op(L, m)    returns L^(m) which is the image of "L downarrow m" under isom_Dm_D
#
section_op := proc(L, m)
	local a, ord, Lm, i, eqns, solutions, nf, tau;
	if nargs = 2 then return procname(L, m, m) fi;
	tau := _Env_LRE_tau;
  	ord := degree(L, tau); 
	Lm := add(a[i]*tau^(m*i), i = 0 .. ord); 
	eqns := {coeffs(RightDivision(Lm, L)[2], tau)}; 
	solutions := solve(eqns, {seq(a[i], i = 0 .. ord)}); 
	nf := add(`if`(lhs(i) = rhs(i), 1, 0), i = solutions); 
        if nf > 1 then
		eqns := {op(eqns), seq(a[i] = 0, i = ord - nf + 2 .. ord)}; 
		solutions := solve(eqns, {seq(a[i], i = 0 .. ord)})
	fi;
	Simpl_op(isom_Dm_D(eval(Lm, solutions), m), tau)
end:

Simpl_op := proc(L, tau)
	sort(collect(primpart(L,tau),tau,factor),tau)
end:

# Input: irreducible L.  This assumption is not checked here (use RightFactors to check).
#
# Output:  true if L is absolutely irredubile, otherwise [p, {R1,...}] where p is a prime
# and where isom_D_Dm(R.i) generate a non-trivial submodule(s) of (D/DL)-downarrow-p,
# proving explicitly that M is not absolutely irreducible.
#
AbsFactorization := proc(L) 
	local i,ord,p,c,G,L_tilde,L_p,S,x,tau;
	x, tau := _Env_LRE_x, _Env_LRE_tau;
	ord := degree(L,tau);
	for p in  ifactors(ord)[2] do
	    if p[1]=2 and args[-1] = `Hom method` then
		# Step 1(a) in Section 3.3:
		L_tilde := subs(tau = -tau, L);
		# Step 1(b) in Section 3.3: (this step needs the file "Hom.txt")
		G := Hom(L, L_tilde);
		if G <> [] then
			if nops(G) > 1 then error "Input was not irreducible" fi;
			G := G[1];
			c := LREtools[MultiplyOperators](subs(tau=-tau,G), G);
			c := LREtools[RightDivision](c, L)[2];
			if has(c, {x,tau}) then error "should not happen" fi;
			G := collect(G/sqrt(c), tau, normal);
			return [2, {seq(Simpl_op(isom_Dm_D(LREtools[RightDivision](LREtools[LCLM](L, 1 + i*G), 1 + i*G)[1], 2), tau), i = {1,-1})}]
		fi
	    else
		L_p := section_op(L, p[1]); 
		if degree(L_p, tau) < ord then 
			return [p[1], {1}]
		fi;
		S := LREtools[RightFactors](L_p, ord/p[1]);
		if S <> {} then
			return [p[1],S]
		fi
	    fi
	od;
	# Step 2
	true
end: