###########################################################################
######                   Generalized Exponents                       ######
###########################################################################

# Mark van Hoeij, 2020.

# Input:
#
#  L must be in C[x, x^(-1), tau, tau^(-1)], use primpart(L,tau) to end up there.
#
#  x, tau are names. Multiplication in C[x,tau] is tau * x = (x+1) * tau
#
#  ext (optional) is a set that describes the field of constants, if not
#      specified then ext = the field generated by the coefficients in L.
#
GenExp := proc(L, x, tau::name, ext::set)  local l,d,i,c,v,RES,s,r,m,P;

    if nargs=1 or type(x,set) then
	return thisproc(L, _Env_LRE_x, _Env_LRE_tau, _passed[2..-1])
    elif not has(L,tau) then
	error "expecting an operator in", tau
    elif not type(L, ratpoly(anything,x)) then
	error "expecting coeffients to be rational functions in", x
    elif not type(L, polynom(anything,x)) then
	return thisproc(primpart(L,tau), _passed[2..-1])
    elif _npassed=3 then
	return thisproc(_passed, indets(L, {'RootOf','radical','nonreal'}))
    fi;

    # Compute degrees of the coefficients and use evala to make sure
    # the degrees are correct. Example 1: degree(1+(1+1/t-(t+1)/t)*x, x);
    #
    l, d := ldegree(L,tau), degree(L,tau);
    for i from l to d do
	c := coeff(L, tau, i);
	v[i] := degree(c, x);
	if c<>0 and evala(coeff(c,x,v[i])) = 0 then
		return thisproc(collect(L, tau, i -> evalaLcoeff(i,x)), x, tau, ext)
	fi
    od;

    # Now use the degrees to compute the tau-NewtonPolygon + Polynomials.
    # Factor these polynomials over Q(ext) and then compute the rest of
    # the generalized exponents using GenExp_sc:
    #
    RES := NULL; m := d;
    while m > l do
	# Next slope in the tau-polygon:
	s := min(seq((v[m]-v[i])/(m-i), i = l .. m-1));
	# Points on slope s:
	r := {seq(`if`( v[m]-v[i] = s*(m-i), i, NULL), i = l .. m)};
	# First point on slope s:
	m := min(op(r));
	# Newton polynomial:
	P := add(coeff(coeff(L,tau,i),x,v[i])*x^((i-m)/denom(s)), i=r);
	P := select(has, evala(Factors(P, ext))[2], x);
	RES := RES, GenExp_sc(L,x,tau,ext,l,d, s,P,max(seq(i[2],i=P)),VarTR())
    od;
    RES := {seq(ShortForm(op(i),evalb(_passed[-1]='LongForm')), i=[RES])};
    userinfo(3,'LREtools',"Generalized Exponents", RES);
    RES
end:

# Produce variables to represent Generalized Exponents.
#
VarTR := proc() local N,T0,T1,T2,T3,T4,T5,T6,T7,T8; options remember;
	N,T0,T1,T2,T3,T4,T5,T6,T7,T8
end:
VarT := proc(n) local T; options remember;
	if n<9 then VarTR()[n+2] else T[n] fi
end:

# Use evala to make sure that lcoeff(nonzero poly) really is not zero.
# See Example 1 above why this can be necessary.
evalaLcoeff:=proc(a,x) local d,i;
	d := degree(a,x);
	if a = 0 or evala(coeff(a,x,d)) <> 0 then return a fi;
	thisproc(add(coeff(a,x,i)*x^i, i=ldegree(a,x)..d-1),x)
end:

# If K is a difference field and r is in K then tau -> r*tau is an
# endomorphism of K[tau]. The image of L is denoted as Symprod(L, r*tau-1).
# The reason for this name is because the solutions of that operator
# are products of solutions of L and the solution of r*tau-1.
#
# How to compute a Generalized Exponent?  The idea is: Repeatedly make use
# of SymProd until the operator has a generalized exponent e = 1, which
# is equivalent to the indicial equation having 0 as a root.
# For that to happen, the first necessary condition is that the tau-polygon has
# a slope equal to 0, and its polynomial has a root 1.
# This first necessary condition is equivalent to saying that the Delta polygon
# has a slope > 0. As long as that slope is < 1 we can make another SymProd to
# make that slope higher, until we have a slope >= 1. That corresponds to a
# non-constant indicial equation, whose roots then give whats left of e.

# Send L = sum(a_k tau^k) to Ls = SymProd(L, T1*tau - 1) = sum(a_k (T1 tau)^k).
# This sends slope s to 0. Then compute c = RootOf(poly for slope s) and
# compute SymProd(Ls, c*tau-1) which sends c to 1. Then call GenExpDelta.
#
GenExp_sc := proc(L,x,tau,ext,l,d, s,P,Mult, N,T0,T1)
    local m, mm, k, a_k, dg, v, c, ORD, i, Ls, RES, ds, ns, j, e;

    m := min(seq(s*k-degree(coeff(L,tau,k),x), k=l..d)); # min(val(a_k T1^k))
    mm := `if`(s=0, 0, mods(m/s, denom(s))); # Makes v-(k-mm)*s below in Z.
    for k from l to d do
	a_k := coeff(L,tau,k);	   # L = sum(a_k tau^k, k=l..d)
	dg := degree(a_k, x);	   # degree(a_k,x) = -val(a_k)
	v := s*k-dg-m;		   # v = val(a_k T^s) - min(val). Hence v >= 0.
	if v > Mult then	   # v is used to trim a_k (v > Mult means trim(a_k) = 0)
	    c[k] := NULL
	else
	    ORD := Mult-ceil(v)+1; # Use only terms from a_k T1^k/T0^v with val <= Mult
	    # This represents the Symmetric Product sum(a_k (T1 tau)^k)) applied to x^N:
	    c[k] := [convert(series(
		  # Rewrite a_k in terms of T0 = 1/x and scale it to T0-valuation 0:
		series(add(coeff(a_k,x,i)*T0^(dg-i), i = dg-Mult+ceil(v) .. dg), T0=0,ORD)*
		  # (T1 * tau)^k is the following series times T1^k * tau^k :
		series(`if`(k<0,mul(1+i*T0, i=k..-1),1/mul(1+i*T0, i=1..k-1))^s, T0=0,ORD)*
		  # tau^k(x^N)/x^N where N is a variable is:
		series((1+k*T0)^N, T0=0, ORD),    T0=0,ORD),polynom),
		  # The factor T0^v = T0^(v-(k-const)*s) * T1^(k-const) belongs to tau^k:
		T0^(v-(k-mm)*s) * T1^(k-mm), v, k]
	fi
    od;
    # Symmetric product of L that sent the slope s to 0 and scaled m = min(val) to 0:
    Ls := [seq(c[k], k=l..d)];  RES := NULL;
    for i in P do
	c := RootOf(i[1],x);          # Goal: SymProd(Ls, c*tau - 1) to send c to 1
	ds, ns := denom(s), numer(s); #       Insert c so that T1^ds = c * T0^ns
	e := [[T0, T0    = 1/x      , 1      ,   1    ,   1           ],
	      [T1, T1^ds = c * T0^ns, s      ,   ds   , degree(i[1],x)]];
	#     [Tn, equation(Tn)     , val(Tn), Ram(Tn), AlgDeg(Tm)    ]
	v := [seq(`if`(j[3]>i[2],NULL, [
		TrimT0N(j[1], i[2], i[2]-j[3], N,T0), # In C[T0,N]. Valuation zero.
		EvalMonomial(j[2],e), j[3], # Monomial and its valuation j[3] >= 0
		j[4] # power of tau         # Call to EvalMonomial is optional here.
	]), j=Ls)];
	RES := RES, GenExpDelta(ext,v,x,e,N, [i[2],i[2]], 1)
    od;
    RES
end:
# f is C[N, T0]. Remove terms with N-degree > Ncut or valuation > T0cut
#
TrimT0N := proc(f, Ncut, T0cut, N,T0) local i;
    add(TrimN(coeff(f, T0, i), Ncut, N) * T0^i, i=0..floor(T0cut))
end:
TrimN := proc(f, Ncut, N) local i;
    add(coeff(f, N, i) * N^i, i=0..Ncut)
end:
EvalMonomial := proc(M, e) local p,q,i,T,k,d,m;
    # e = [ [T0,    T0 = 1/x     ,   1    ,  1,       1   ],
    #       [T1, T1^r1 = c1*T0^n1, val(T1), r1, algDeg(c1)],  ...]
    p := M; q := 1;
    for i from nops(e) to 2 by -1 do T := e[i,1]; if has(p, T) then
	k := degree(p, T); d := e[i,4]; m := modp(k,d); # use mods?
	q := q * T^m;
	p := coeff(p, T, k) * EvalMonRem(T, k-m, e[1..i], d)
    fi od;
    p * q
end:
# Compute T^k where k is a multiple of d and T in {T1, T2, ...}
#
EvalMonRem := proc(T,k,e,d) options remember;
	if k = 0 then return 1 fi;
	evala(EvalMonomial(thisproc(T,k-sign(k)*d,e,d),e) * rhs(e[-1,2])^sign(k))
end:

# The Delta Polygon is the Newton Polynomial of L written as an element of
#	C(( T0^(1/r) ))[Delta]
# where Delta = tau - 1.  We do not explicitly write L this way, instead,
# we apply L to x^N where N is a variable. That produces an element Npoly in
#	C(( T0^(1/r) ))[[ N ]] 
# and we compute the Newton Polygon of that. To ensure that we only compute
# slopes of interest, we cut Npoly down to N-degree Top[1] (we throw away
# everything of degree > Top[1]). We also throw away terms with valuation
# greater than Top[2] since those can not affect the Newton Polygon/Polynomials.
#
GenExpDelta := proc(ext,Ls,x, e, N, Top,SlopeBound,Npoly)
	local s,P,   k,cf,v,r,m,R,oldR,   i,d,T,Tn;
    if _npassed=7 then return thisproc(_passed,
	collect(add(i[1]*i[2], i=Ls),{N,seq(i[1],i=e)},'distributed'))
    fi:
    for k from Top[1] to 0 by -1 do
	# Evaluate only terms needed for the maximal slope.
	# This way in recursion, no terms will be recomputed.     max(val):
	cf[k], v[k] := LowestValuationTerm(coeff(Npoly, N, k), e, Top[2]-(Top[1]-k)*s,
		Top[2]-(Top[1]-k)*SlopeBound -`if`(k=Top[1],1/1000,0) ); # < min(val)
	s := `if`(k=Top[1], 0, max(s, (v[Top[1]]-v[k])/(Top[1]-k))) # max slope found so far
    od:	                                                            # 0 <= s < SlopeBound <= 1
	r := {seq(`if`(v[Top[1]]-v[k] = s*(Top[1]-k), k, NULL), k=0..Top[1])};
	if s=0 then m := 0 else m := min(op(r)) fi;
	oldR := mul(k[4],k=e); R := denom(s*oldR); # old and new ramification

	# Newton Polynomial for the maximal slope:
	P := add(cf[k]*x^((k-m)/R), k=r);
	d := degree(P,x);
	T := TermWithVal(R*s, e);
	P := add(EvalMonomial(coeff(P,x,i)*numer(T)^(d-i)*denom(T)^i,e)*x^i, i=0..d);
	P := select(has,evala(Factors(evala(Primpart(P,x)),ext union indets(e,RootOf)))[2],x);
	Tn := VarT(nops(e));
	if s = 0 then
	    P := IntDiff([seq(collect(i[1]/lcoeff(i[1],x),x,evala)$i[2],i=P)],x,oldR);
	    return seq([e,i,x,N], i=P)
	fi;                                                      # val, ram, algdeg
	seq(SymprodTn([op(e),[Tn,Tn^R = RootOf(i[1],x)*e[1,1]^R/T, 1-s,  R , degree(i[1],x)]],
	    [i[2], Top[2]-s*(Top[1]-i[2])], s, ext,Ls,x,N), i=P),
	# Treated slope s with SymprodTn. Use recursion for slopes < s
	`if`(m=0, NULL, thisproc(_passed[1..5], [m, v[m]], s, Npoly))
end:
ShortForm := proc(e,f,x,N,LongForm) local i;	# Convert to a shorter format.
	if _Env_RightFactors_SymCubeCase = true and nops(f[2])=1 then
		NULL	# Used only in ReduceToOrderTwo
	elif LongForm then
		[op(e),[op(f), 1, 1, degree(f[1],x)]]
	else
		[e[2,1]*mul(1+e[i,1],i=3..nops(e))*(1+N/x),
		 seq(eval(e[i,2],e[1,2]),i=2..nops(e)), N = [seq(RootOf(f[1],x)+i, i=f[2])]]
	fi
end:
ValTerm := proc(t, e) local i;  # Each entry of e:  [Tn, equation, val(Tn), ram, algdeg]
	add(degree(t,i[1])*i[3], i=e)
end:
TermWithVal := proc(v,e) # Product of T0,T1,... with valuation v
	local w, i, q, f, m;
	w := v; q := 1;
	for i from nops(e) to 1 by -1 do f := e[i]; if f[4] > 1 then
		m := mods(w/f[3], f[4]);
		q := q * f[1]^m;  w := w - m*f[3]
	fi od;
	q * f[1]^w
end:
SymprodTn := proc(e, Top,s, ext,Ls,x,N)
	local L,i,j,A;
	L := NULL;     # Compute a symmetric product:
	for j in Ls do if j[3] <= Top[2] then
		A := Tntauk(j[4], e[1,1], e[-1,1], e[-1,3], Top[2]-j[3]);
		A := collect(A*j[2],{seq(i[1],i=e)},'distributed');
		A := [`if`(type(A,`+`),op(A),A)];
		A := add(`if`(ValTerm(i,e) <= Top[2], EvalMonomial(i,e), 0), i=A);
		L := L,[TrimT0N(j[1], Top[1], Top[2]-j[3], N, e[1,1]), A, j[3], j[4]]
	fi od;
	GenExpDelta(ext,[L],x, e, N, Top, s)
end:
# Compute ((1+Tn)*tau)^k = output * tau^k. Compute terms up to valuation V.
# 0 < val(Tn) = v < 1
#
Tntauk := proc(k, T0, Tn, v, V)
	local p,s,i;
	if   V <   v or k=0 then 1
	elif V < 2*v or k=1 then 1 + k*Tn
	else
		p := Tn^floor(1/v);
		s := `if`(k<0,1/mul(1+Tn/(1+i*T0*p)^v, i=k..-1),
				mul(1+Tn/(1+i*T0*p)^v, i=0..k-1));
		s := convert(series(s, Tn=0, floor(V/v)+1), polynom);
		collect(eval(s, T0=T0/p),{T0,Tn},'distributed')
	fi
end:
# Output: only the term of lowest valuation, if its
# valuation is in (BelowMinVal, MaxVal].
#
LowestValuationTerm := proc(F, e, MaxVal, BelowMinVal)
	local B,v,m,M,i,w;
	B := BelowMinVal;
	if F = 0 then return 0, infinity fi:
    do
	v := MaxVal; m := 0; M := 0;
	for i in [`if`(type(F,`+`),op(F),F)] do
		w := ValTerm(i, e);
		if  w <= B then M := M + i
		elif w = v then m := m + i
		elif w < v then m :=     i; v := w
		fi;
	od; # Uncomment this for bugtesting:
	    # if evala(M) <> 0 then error "bug detected" fi;
	# Either apply EvalMonomial to terms of m or to 2nd argument for GenExpDelta
	m := evala(m);
	if m <> 0     then return m, v        fi;
	if v = MaxVal then return 0, infinity fi;
	B := v;	    # lprint("Terms cancelling out")
    od:
end:
# v is a list of irreducible polynomials. Compute which of them
# are integer/r shifts of each other, determine the "minimal" ones, and
# the non-negative integers you have to shift them with to get the other ones.
#
IntDiff := proc(v::list, x::name, r::posint)
    local i,j,d,e,large,dif,n;
    large := {}; n := nops(v); for i to n do dif[i] := 0; d[i] := degree(v[i],x) od;
    for i from 1 to n-1 do for j from i+1 to n do if d[i]=d[j] then
	e := evala(coeff(v[i],x,d[i]-1)-coeff(v[j],x,d[i]-1));
	if type(r*e/d[i],integer) and evala(v[i]-subs(x=x+e/d[i],v[j]))=0 then
	    if e >= 0 then
		dif[i] := dif[i], e; if e > 0 then large := large union {j} fi
	    fi;
	    if e <= 0 then
		dif[j] := dif[j], -e; if e < 0 then large := large union {i} fi
	    fi
	fi
    fi od od;
	[seq(`if`(member(i,large),NULL,[v[i],sort([dif[i]])]),i=1..n)]
end:


###############################################################################
###     To get a degree bound we have to multiply generalized exponents     ###
###     (taking a Norm when defined over an extension), and an AdjustTerm   ###
###     which is the sum of the valuations of the differences.              ###
###############################################################################

eNorm := proc(e, x, tau, ext) local ext2, p, S, i, s, R, r, a;
	ext2 := ext union indets(e[2], RootOf);
	# p := evala(Norm(tau-rhs(e[2,2]), ext2, ext))^mul(i[4]*i[5], i=e[3..-1])
	# has NewtonPolygon terms, but this p has more terms at non-integer slopes:
	p := eTopCoeffs(e, x, tau, ext);
	S := 0;
	for i in e[3..-1] do 
		ext2 := ext2 union indets(i,RootOf);
		S := S * i[5];
		if type(i[2],list) then  # This is the last entry of e
			s := -coeff(i[1],x,i[5]-1)/coeff(i[1],x,i[5])
		elif numer(i[3])=1 then
			s := EvalMonomial( -(-i[1])^denom(i[3]) / e[1,1], e) / denom(i[3])
		else
			next
		fi;
		S := S + evala(Trace(s, ext2, ext));
	od;
	r := mul(i[4],i=e);     # total ramification in e
	R := mul(i[5],i=e) * r; # total degree [C((t))[e]:C((t))]
	a := ( AdjustmentTerm(e,e) - R*(e[2,3]+1) )/2;
	S := collect(coeff(p,tau,0) * (1 + (r * S - a)/ x), x, evala);
	(-1)^R * S, R, R * e[-1, 2], e[2, 3], p
end:
AdjustmentTerm := proc(e1, e2) local e, N, S, d, i, m, mn, n, v1, v2, vm;
	v1 := e1[2, 3];
	v2 := e2[2, 3];
	vm := min(v1, v2);
	d := 2;
	while d < nops(e1) and d < nops(e2) and e1[d] = e2[d] do
		d := d+1
	od;
	# Now d is the first entry where e1 and e2 disagree.

	S := vm; N := 1; m := 0;
	for i from 2 to d-1 do
		n := e1[i,4] * e1[i,5];
		N := n * N;
		mn := min( e1[i+1,3],e2[i+1,3]);
		S := n^2 * S + N * (mn-m);
		m := mn
	od;
	S * mul(mul(i[4]*i[5], i=e[d..-1]), e=[e1,e2])
end:

# In addition to computing the "Norm + AdjustmentTerm", which gives the dominant
# part of the determinant, we also compute the top coefficients. That includes
# the Newton polygon, Newton polynomial, as well as (in case of non-integer slope)
# the coefficient of the highest integer point in the polygon.
#
eTopCoeffs := proc(e, x, tau, ext)
	local S, dS, ext2, i, id, t;
	S := tau - expand(mul(1+e[i,1], i=3..nops(e)-1));
	id := {tau, seq(e[i,1],i=1..nops(e)-1)};
	S := add(`if`(ValTerm(t,e)<1, EvalMonomial(t,e), 0), t=S);
	for i from nops(e)-1 to 2 by -1 do if i=2 or has(S, e[i,1]) then
		if i=2 then S := subs(tau=tau/e[2,1], S) fi;
		if lhs(e[i,2])=e[i,1] then
			S := subs(e[i,2],S)
		else
			S := resultant(numer(S), numer(lhs(e[i,2])-rhs(e[i,2])), e[i,1])
		fi;
		S := collect(S/lcoeff(S,tau), id, 'distributed');
		if not type(S,`+`) then error "unexpected format" fi;
		dS := degree(S,tau);
		S := add(`if`( ValTerm(t,e)+`if`(i=2,e[2,3]*(degree(t,tau)-dS),0)<1,
			EvalMonomial(t,e), 0), t=S)
	fi od;
	S := eval(S, e[1,2]);
	ext2 := ext union indets(e,RootOf);
	if ext2 <> ext or e[-1,5]>1 then
		S := collect(evala(Norm(S, ext2, ext))^e[-1,5], {x,tau}, 'distributed');
		dS := degree(S,tau);
		S := add(`if`(-degree(t,x)+e[2,3]*(degree(t,tau)-dS)<1, t, 0), t=S)
	fi;
	collect(S,tau,Normalizer)
end:

# Faster version, does the same.
#
eTopCoeffs := proc(e, x, tau, ext)
	local S, dS, ext2, i,j, id, t, z,pz,S0;
	S := tau - mul(1+EvalMonomial(e[i,1],e), i=3..nops(e)-1);
	id := {tau, seq(e[i,1],i=1..nops(e)-1)};
	for i from nops(e)-1 to 2 by -1 do
		if i > 2 and not has(S, e[i,1]) then next fi;
		dS := degree(S,tau);
		if i=2 then S := e[2,1]^dS * subs(tau=tau/e[2,1], S) fi;
		S := collect(S, id, 'distributed', Normalizer);
		S := add(`if`(ValTerm(t,e)+`if`(i=2,e[2,3]*(degree(t,tau)-dS),0)<1,
			EvalMonomial(t,e), 0), t=S);
		S0 := S;
		pz := NumberTheory['CyclotomicPolynomial'](e[i,4],z);
		for j to e[i,4]-1 do
			S := rem(S * rem(subs(e[i,1] = z^j * e[i,1], S0),pz,z), pz,z);
			S := collect(S, id, 'distributed', Normalizer);
			dS := degree(S,tau);
			S := add(`if`(ValTerm(t,e)+`if`(i<3,e[2,3]*(degree(t,tau)-dS),0)<1,
				EvalMonomial(t,e), 0), t=S);
		od;
		S := collect(S, id, 'distributed', Normalizer);
		if has(S, {z, e[i,1]}) then error "unexpected" fi
	od;
	S := eval(S, e[1,2]);
	ext2 := ext union indets(e,RootOf);
	if ext2 <> ext or e[-1,5]>1 then
		S := collect(evala(Norm(S, ext2, ext))^e[-1,5], {x,tau}, 'distributed');
		dS := degree(S,tau);
		S := add(`if`(-degree(t,x)+e[2,3]*(degree(t,tau)-dS)<1, t, 0), t=S)
	fi;
	collect(S,tau,Normalizer)
end:



###############################################################################
###     Now form all possible products, which gives the dominant part of    ###
###     the determinant and top coefficients for all possible RFactors.     ###
###############################################################################

DetFactors := proc(L, x, tau, ext) local i,j,k,S;
	S := SubLists(GenExp(_passed,'LongForm'), x, tau, ext);
	eval({seq([i[1], i[2], {seq( [k, mul( `if`(j[2] = k, j[1]^j[3], 1), j=i[3])],
		k = map(l -> l[2], i[3]))}], i = S)}, VarT(1)=tau)
end:
SubLists := proc(E::set, x, tau, ext)
	local j, k, ANS1, ANS2, L, P, R, a, d, le, i, p, ve;
	if E={} then return E fi;
	P[0] := 1;
	le := E[1];	# A generalized exponent in LongForm
	d, R, L, ve, p := eNorm(le, x, tau, ext);
	if nops(L)>1 then
		a := AdjustmentTerm(le, le)
	fi;
	for i from 1 to nops(L) do
		# Maybe these can be better written as cnd lists:
		P[i] := MultiplyLocalDeterminants(P[i-1], d * (1+(L[-i] - (i-1)*a)/x), x)
	od;
	# OUTPUT IF nops(E) = 1:   { [P[i], R*i, {[p^i, ve]}, {[i, le]} ] , i > 0 }
	# Otherwise, we have to combine these with MultiplyLocalDeterminant, +, union
	# and now use the last one for further adjustment-terms.
	#
	ANS1 := {seq( [P[i], R*i, {[p,ve,i,R,le]}], i=1..nops(L))};
	if nops(E)>1 then
		ANS2 := thisproc(E minus {le}, x, tau, ext);
		ANS1 := ANS1 union ANS2 union {seq(seq( [MultiplyLocalDeterminants(
		  P[i], j[1]*(1 - add(i*k[3]*AdjustmentTerm(le, k[5]), k=j[3])/x), x),
		  R*i + j[2], j[3] union {[p,ve,i,R,le]}], i=1..nops(L)), j=ANS2)}
	fi;
	ANS1;
end:
# Local determinants are of the form (...)*x^d + (...)*x^(d-1)
# where lower-degree terms are irrelevant. So to multiply them,
# we can simply discard the lower degree terms.
#
MultiplyLocalDeterminants := proc(d1, d2, x) local p,d;
	p := d1*d2; d := degree(p,x);
	Normalizer(coeff(p,x,d))*x^d + Normalizer(coeff(p,x,d-1))*x^(d-1)
end:

