###########################################################################
######         Compute Right-Factors of Recurrence Operators         ######
###########################################################################

# The algorithm starts by computing generalized exponents, defined in
#      www.math.fsu.edu/~hoeij/papers/issac10/J.pdf
# These are used to efficiently compute hypergeometric solutions
# of the d'th exterior power system described in
#      www.math.fsu.edu/~hoeij/slides/bronstein/slides.pdf
# Solutions of this system that satisfy Plucker relations correspond to
# factors of L.
#
# Mark van Hoeij, 2020.

# Input: A recurrence operator L in C(x)[tau].
#        An integer d.
#        (optional): a set ext giving the algebraic extensions in C.
#
# Output: If d > 0:  { Right-Factors of order d }.
#         If d = 0:  [ Degree-bounds for Right-Factors of order 1..n ].
#
RFactors := proc(L, d::nonnegint, ext::set)
local n, i, M, v, x, tau, pL, P;
	x, tau := _Env_LRE_x, _Env_LRE_tau;
	n := degree(L, tau);
	if L = 0 or n = FAIL then
		error "First argument should be a non-zero polynomial in the shift operator", tau
	elif _npassed < 3 then
		return thisproc(_passed, indets(L, {'RootOf','radical','nonreal'}))
	elif ldegree(L,tau)<>0 then
		userinfo(1,{'LRETools','rsolve'},cat("Operators are factored in C(",x,")[",tau,",",
		  tau,"^(-1)]. The factor ",tau^ldegree(L,tau)," is ignored as it is a unit in this ring."));
		return thisproc(`LREtools/mult`(tau^(-ldegree(L,tau)),L), _passed[2..-1])
	elif not type(L, polynom(ratpoly(anything,x), tau)) then
		error "First argument should be polynomial in", tau, "and rational in", x
	elif d > n then
		return {}
	elif d = n then
		return `if`(d=0, [], {SimplifyR(L, x, tau)})
	elif ext = {} then
		Normalizer := normal
	else
		Normalizer := evala
	fi;
	pL := collect(primpart(L,tau), tau);
	P := PairsAB(pL, x, tau, ext);
	if d=0 then
		# compute degree-bounds only but no actual factors
		M := {seq([v[3,2], degreeBound(pL, x, n, v[3,2], v, v[1])], v=P)};
		return [seq(max(seq(`if`(v[1]=i, v[2], NULL), v=M)), i=1..n)]
	fi;
	P := [seq(`if`(v[3,2]=d, v, NULL), v=P)];
	if P = [] then return {} fi:
	# We computed all pairs, then discarded pairs of order <> d. We could have
	# avoided computing pairs of order > d but it makes little CPU difference.
	M := collect(pL/factor(lcoeff(pL,tau)),tau,Normalizer);
	M := ExteriorPowerSystem(M, tau, n, d);
	i := evalb(ext union indets(L) minus {x,tau} = {}); # input over Q?
	userinfo(3,'LREtools', "Number of cases", nops(P));
	M := {seq(PolySols(M, x, tau, n, d, v, v[1] * lcoeff(v[3,1],x), i), v=P)};
	forget(coeffRemember); forget(ValueTopCoeff);
	map(SimplifyR, {seq(PluckerRelations(n, d, x, tau, i, ext), i=M)}, x, tau)
end:
SimplifyR := proc(R, x, tau)
	if type(R,list) then return [thisproc(R[1], x, tau),op(2..-1,R)] fi:
	# sort(map(sort, collect(primpart(R*`tools/sign`(R,[tau,x]),tau),tau,factor),x), tau)
	# Newer version:
	`LREtools/PrimPart`(R, x, tau)
end:

degreeBound := proc(L, x, n, d, v, r)
local i;
	if v[2] - degree(CommonFactor(n-d, x, r), x) < 0 then # implies lcoeff(R,tau)=0
		-infinity
	elif d < n then
		v[2] - degree(CommonFactor(n-d-1, x, r), x) # Bounds deg(leading coeff)
		     + add(degree(i[2],x), i=v[3,3])        # adds NewtonPolygon
	else
		degree(L,x) # Handle d=n separately because CommonFactor
	fi                  # does not work for n-d-1 = -1
end:


###############################################################################
######                       Exterior Power System                       ######
###############################################################################

# Input: L in C(x)[tau]
# Output: matrix M such that order-d factors of L are hypergeometric solutions of MY = tau(Y).
ExteriorPowerSystem := proc(L, tau, n, d)
local B, b;
	B := Subsets(n,d);  # if b = [1,3,4] in B then view it as tau^0 wedge tau^2 wedge tau^3
	Matrix([seq(WriteInTermsOfBasis(L, tau, b, B, n), b=B)]) # = the action of tau on B.
end:

# Subsets of {1..n} given as ordered lists
Subsets := proc(n::nonnegint, d::nonnegint)
options remember;
local i;
	if d=0 then
		[ [ ] ]  # the only list with 0 elements is [ ]
	elif n=0 then
		[ ]
	else
		[ op(thisproc(n-1, d)), seq( [op(i),n], i=thisproc(n-1, d-1)) ]
	fi
end:

WriteInTermsOfBasis := proc(L, tau, b, B, n)
local i, tau_b, j;
	tau_b := [seq(i+1, i=b)];	# b = tau^i1 wedge ... So tau_b = tau^(i1+1) wedge ...
	if tau_b[-1] <= n then
		WriteInTermsOf(tau_b, B, 1)
	else
		tau_b := tau_b[1..-2];
		add(`if`(member(i, tau_b), 0, WriteInTermsOf(sort([op(tau_b), i]), B,
		 -coeff(L, tau, i-1)/coeff(L, tau, n) * mul(`if`(j>i,-1,1), j=tau_b))), i=1..n)
	fi
end:

WriteInTermsOf := proc(b, B, c)
local i;
	for i do if B[i]=b then return [0$(i-1), c, 0$(nops(B)-i)] fi od
end:


###############################################################################
######    Solving M Y = r tau(Y) where M is the Exterior Power System    ######
###############################################################################


CommonFactor := proc(s, x, r)
options remember, system;
	denom( ProductShifts(eval(r, x=x-s), x, s) )
	# lcm(%, thisproc(s-1,x,r)) if s>0 is also a valid bound,
	# but it turns out to be the same.
end:

ProductShifts := proc(r, x, n)
options remember, system;
local i;
	Normalizer( mul(eval(r, x=x+i), i=0..n-1) )
end:

ValueTopCoeffNotSorted := proc(S, x, tau, d, G, g)
local k;
	if member(S[-1],S[1..-2]) then return 0 fi;
	k := 1;
	while S[k] < S[-1] do k := k+1 od;
	(-1)^(d-k) * ValueTopCoeff([op(S[1..k-1]),S[-1],op(S[k..-2])], x, tau, d, G, g)
end:
ValueTopCoeff := proc(S, x, tau, d, G, g)
options remember;
local i,j,DB;
	if max(S) = d then
		return 1
	fi;
	DB := floor( -add((S[j]-j)*g[j], j=1..d) );
	coeff(-add(coeffRemember(G,x,tau,i) * ValueTopCoeffNotSorted([op(S[1..-2]), S[-1]-d+i],
		x, tau, d, G, g), i=0..d-1), x, DB) * x^DB
end:

coeffRemember := proc(G,x,tau,i)
options remember; # (table cleared at the end of RFactors)
local j,L;
	if _npassed=3 then
		L := mul(j[2], j=G);
		collect(L/lcoeff(L,tau), tau, normal)
	else
		L := coeff(thisproc(G,x,tau),tau,i);
		lcoeff(L,x) * x^degree(L,x)
	fi
end:

# Solve M Y = r tau(Y) where M is an exterior power matrix.
PolySols := proc(M, x, tau, n, d, v, r, overQ)
local a::nothing, j, Ansatz, CF, DB, G, S, S1, Sn, vars, g, i, isz, l, s, sb,
	 top1, SB, SBN, ANS, vl, Y, MAT, e, c, varY, npiv, piv, sol;

	G := sort([op(v[3,3])], (i1,i2) -> evalb(i1[1] > i2[1]));
	# Example: [[-1/2, tau^2-x], [-2/3, tau^3-x^2]]
	g := [seq(i[1]$degree(i[2],tau), i=G)];

	i := 0; SB := {};
	Sn := Subsets(n,d);
	for S in Sn do if member(1,S) then
		i := i+1;
		CF[S] := CommonFactor(n-max(S), x, r);
		DB[i] := v[2] - add((S[j]-j)*g[j], j=1..d) - degree(CF[S],x);
		sb := {};
		if i = 1 then
			if DB[1] < 0 then
				return NULL # DegreeBound DB[1] is negative
			fi;
			userinfo(4,'LREtools', "Building equations for type", tau-r,
			  "with top coefficients", mul(j[2],j=G),
			  "and component-wise Content Bound", [seq(CommonFactor(n-j, x, r), j=d..n)]);
			top1 := a[i, DB[i]] * lcoeff(CF[S], x);
		else
			DB[i] := floor(DB[i]);
			sb := a[i,DB[i]] = lcoeff(ValueTopCoeff(S,x,tau,d,G,g), x) * top1 / lcoeff(CF[S],x);
			SB := SB union {lhs(sb)-rhs(sb)}
		fi;
		Ansatz[S] := subs(sb, add(a[i,j] * x^j, j=0..DB[i])) # Part of Ansatz without common factor
	fi od;
	for S in Subsets(n,d) do if not member(1,S) then
		s := min(S)-1;
		S1 := [seq(j-s,j=S)];	# Ansatz[S] is a shift of Ansatz[S1]
		CF[S] := Normalizer(ProductShifts(r, x, s) * eval(CF[S1], x = x+s));
		Ansatz[S] := subs(x=x+s, Ansatz[S1]) # don't Normalize this
	fi od;
	Ansatz := [seq(CF[S] * Ansatz[S], S = Sn)];
	# Further improvements to "content bound" CF[S] may not be needed as it
	# tends to be close to content( ConvertToOperator(Ansatz,tau,d), tau).

	# Put primary variables a[1..d, *] at the end (ordering affects CPU time)
	g := [seq(i-j, j=0..i-1)];
	vars := [seq(seq(a[i,DB[i]-j], j=`if`(i=1,0,1)..DB[i]),i=g)];
	c := nops(vars);

	# Replacing many variables a[i,j] by 1 variable Y speeds polynomial arithmetic:
	varY := {seq(vars[i] = Y^i, i=1..c)}:
	l := Vector( eval(Ansatz, varY) );
	isz := M . l - r * subs(x=x+1, l); # should be zero

	g := nops(Sn); # Reverse ordering of isz and take coefficients:
	isz := map(coeffs,[seq(`if`(member(n, Sn[g-i]),
		sort(collect(numer(isz[g-i]),x),x),
		NULL), # Discard entries that are zero by construction
	i=0..g-1)], x);

    if overQ then
	MAT := Matrix(nops(isz),c);
	i := 0;
	for g in isz do e := expand(g); if e<>0 then
		i := i+1;
		for j in `if`(type(e,`+`),e,[e]) do
			MAT[i,degree(j,Y)] := eval(j, Y=1)
		od;
	fi od:
	MAT, piv, npiv := SolveTools:-LinearSolvers:-RationalDenseNullSpace(MAT, ':-skiptest'=true);
	if npiv = [] then return NULL fi;
	sol := {seq(vars[i] = add(MAT[i,j] * vars[npiv[j]], j=1..nops(npiv)), i = piv)};
	l := eval(Ansatz, sol);
	i := convert(l, Vector);
	# Skipped test in NullSpaceOverQ, to prove correctness check that this is zero:
	MAT := M . i - r * subs(x=x+1, i);
	MAT := {seq(`if`(member(n, Sn[i]), normal(MAT[i]), NULL), i=1..nops(Sn))} minus {0};
    fi;
	if MAT <> {} then
		userinfo(2,'LREtools',`if`(overQ, "Check failed, redoing solve", "Input is not over Q"));
		sol := subs([seq(Y^(c-i) = vars[c-i], i=0..c-1)], map(collect,{op(isz)},Y));
		sol := SolveTools:-Linear(sol, vars);
		l := eval(Ansatz, sol);
	fi;
	ANS := NULL;
    do
	while v[2] - degree(l[1],x) > 0 do
		if l[1] = 0 then
			return ANS # stop because order(l) < d
		fi;
		SBN := eval(SB, {seq(i = a[op(1,i),op(2,i)-v[2]+degree(l[1],x)],
		  i = indets(SB) minus indets([_passed]))});  # Not indets intersect vars
		if Normalizer(eval(SBN, sol)) minus {0} = {} then
			break # SBN holds
		fi;
		# Require SBN to hold so that top coeffs match GenExp data.
		# That reduces dim(output), reducing EQ in PluckerRelations.
		sol := SolveTools:-Linear(sol union SBN, vars);
		l := eval(l, sol)
	od;
	vl := {op(vars)} intersect indets(l);
	# lprint("Number of solutions", nops(vl));
	if vl = {} then
		return ANS
	else
		ANS := ANS, [seq( map(coeff, l, i), i=vl)];
		if nops(vl)=1 then return ANS fi:
		sol := SolveTools:-Linear(sol union {lcoeff(l[1],x)}, vars);
		l := eval(l, sol);
		# In case the leading coefficient is not normalized:
		while l[1]<>0 and Normalizer(lcoeff(l[1],x))=0 do
			l := subs( lcoeff(l[1],x) = 0, l)
		od:
	fi
    od
end:


################################################################################
######    Use Plucker Relations to select solutions that give RFactors    ######
################################################################################


# Converts a solution v of the exterior power system to an operator.
ConvertToOperator := proc(v, tau, d)
local i;
	add( (-1)^i * v[i+1] * tau^(d-i), i=0..d)
end:

# Input: R is a linear combination of vars.
# Output: change-of-basis that typically makes R smaller.
RowReduceR := proc(R, x, tau, d, vars)
local cc,i,treated,sb,c,cf,v,tmp::nothing;
	if has(R, tau) then
		cc := [lcoeff(R,tau), seq(coeff(R,tau,i), i=0..d-1)]
	else
		cc := R
	fi;
	treated := {};
	sb := NULL;
	for c in cc do for i from degree(c,x) to 0 by -1 while treated<>vars do
		cf := subs(sb, coeff(c,x,i));
		v := (indets(cf) intersect vars) minus treated;
		if v={} then next fi;
		sb := sb, subs(tmp=v[1], SolveTools:-Linear({cf = tmp}, {v[1]}));
		treated := treated union {v[1]}
	od od;
	sb
end:

# Input: n, d, V = list-of-poly-sols
PluckerRelations := proc(n, d, x, tau, V, ext)
local R, i, vars, P, B, EQ, l, s, Y::nothing, sb, N, RI;
	N := nops(V);
	R := add(_Z[i] * ConvertToOperator(V[i],tau,d), i=1..N);
	R := primpart(R, {tau,seq(_Z[i],i=1..N)});
	EQ := {};
	if d > 1 and d < n-1 then
		# Convert lists to polynomials so we can divide by the gcd:
		P := add(collect(add(_Z[i] * V[i][l], i=1..N), x) * Y^l, l=1..nops(V[1]));
		P := primpart(P, Y); # In ex2 this is faster
		B := Subsets(n,d);
		# Now construct binomial(n-1, d)-1 quadratic equations for the _Z[i]
		EQ := {seq(seq(add(
			  eval(coeff(R, tau, l-i), x=x+i) * CoeffOfP(s, l, B, P, Y),
			  l=i..d+i), s = Subsets(i+d-1, d-2) )	# s wedge tau^i R = 0
			,i=1..n-1-d)};
		EQ := {seq(coeffs(collect(i,x),x), i=EQ)}
	fi;
	vars := {seq(_Z[i],i=1..N)};
	sb := RowReduceR(R, x, tau, d, vars); # Change-of-basis to make R smaller.
	l := indets(subs(sb,lcoeff(lcoeff(R,tau),x))) intersect vars;
	if nops(l) <> 1 then error "Due to RowReduceR this should have only one variable" fi;
	sb := sb, {l[1] = _Z[N], _Z[N] = l[1]};
	R := collect(subs(sb,R), tau, Normalizer);
	RI := subs(_Z[N] = 1, R);    # Inhomogeneous R and
	vars := vars minus {_Z[N]};  # its variables
	EQ := map(numer@evala, subs(sb, _Z[N]=1, EQ)) minus {0};
	if EQ = {} then
		return `if`(N=1, RI, [RI, "a family of dimension", N-1, "with equations", {} ])
	fi;
	EQ := Groebner['Basis'](EQ, tdeg(op(vars)));
	if member(1,EQ) then
		NULL
	else
		l := {seq(`if`(degree(i,vars) = 1, i, NULL), i=EQ)};
		if l <> {} then
			l := SolveTools:-Linear(l, vars);
			RI, EQ := subs(l, RI), subs(l, EQ);
			vars := indets(map(rhs,l)) intersect vars;
			EQ := Groebner['Basis'](EQ, tdeg(op(vars)))
		fi;
		N := nops(vars);
		RI, EQ, vars := op(subs( {seq(vars[i] = _Z[i], i=1..N)}, [RI, EQ, vars]));
		if N = 0 then
			RI
		elif N = 1 and has(EQ,_Z[1]) then # Solving a univariate quadratic equation:
			seq(eval(RI, i), i=solve(evala(map(Factor,EQ,ext)), vars))
		else
			[RI, "a family of dimension", Groebner:-HilbertDimension(EQ,vars),
			 "with equations", {op(EQ)}]
		fi
	fi
end:
# s wedge l corresponds to an element of basis B.
# Find that element and return the corresponding coefficient of P.
#
CoeffOfP := proc(s, l , B, P, Y)
local i,b,sg;
	if member(l, s) then
		0
	else
		sg := mul(`if`(i > l, -1, 1), i=s); # = sign of sort-reordering:
		b := sort( [op(s),l] );
		b := [1, seq(i+1, i=b)]; # Using only wedge products that include tau^0
		for i do if B[i]=b then return sg * coeff(P, Y, i) fi od
	fi
end:


##################################################################################
###     Construct r = c * A/B for which M Y = r tau(Y) could have PolySols.    ###
###     Pairs A,B come from trailing/leading coefficient while c and degree    ###
###     bound come from Generalized Exponents.                                 ###
##################################################################################


# Example:
#	A = [[x+1/2,1], [x,2], [x-1,1]];
#	P = x
# Then Output = [[x + 1/2, 1]],
#		[[x, 0], [x, 0], [x - 1, 1]]
#
FindIntShifts := proc(A, P, x)
local res, resA, i, d;
	resA := NULL; res := NULL;
	for i in A do
		if IntShift(P, i[1], x, 'd') then
			res := res, [i[1],d]$i[2]
		else
			resA := resA, i
		fi
	od;
	[resA], sort([res],proc(i,j) evalb(i[2]<j[2]) end)
end:

# Check if two monic irreducible poly's are an integer shift from each other.
IntShift := proc(A, P, x, d)
local n,e;
	n := degree(A,x);
	if n=degree(P,x) then
		e := Normalizer(coeff(A,x,n-1)-coeff(P,x,n-1))/n;
		if type(e,integer) and Normalizer(A-subs(x=x+e,P))=0 then
			if _npassed=4 then
				d := e
			fi;
			return true
		fi
	fi;
	false
end:

# Leading and Trailing Singularities after some desingularization.
#
LTSing := proc(L, x, tau)
local i, l, d, M;
	l, d := ldegree(L,tau), degree(L,tau);
	i := iquo(l+d, 2);
	M := coeff(L, tau, i+1) * subs(x=x+1, L) * tau - subs(x=x+1, coeff(L, tau, i)) * L;
	M := primpart(M, tau);
	[gcd(coeff(L,tau,d), subs(x=x-1, coeff(M,tau,d+1))), gcd(coeff(L,tau,l), coeff(M,tau,l))];
end:

# Faster version
#
LTSing := proc(L, x, tau, ext)
local M,l,d,m,c,A,i,r, c1,c0;
	l, d := ldegree(L,tau), degree(L,tau);
	A := seq(`LREtools/g_p/factors`(coeff(L,tau,c),x,ext), c=[d,l]);
	if ext union (indets(L) minus {x,tau}) <> {} then return A fi:
	m := iquo(l+d, 2);
	c1 := coeff(L, tau, m+1);  c0 := subs(x=x+1, coeff(L, tau, m));
	M := c1 * subs(x=x+1, collect(L*tau,tau)) - c0 * L;	# (c1 tau - c0) * L
	for c to 2 do for i to nops(A[c]) do
	    m[c,i] := A[c,i,2];
	    if degree(A[c,i,1],x)=1 then
		r := solve(A[c,i,1],x) - `if`(c=1,1,0);
		if eval(M,x=r) = 0 and Normalizer(eval([c1,c0][c], x=r))<>0 then
		    m[c,i] := `if`(m[c,i]=1, 0, m[c,i] - ldegree(convert(subs(x=x+r,
				series(M, x=r, m[c,i])),polynom), x));
		    userinfo(3,'LREtools',cat("Reduced ",
			`if`(c=1,"lead","trail"),"ing singularity"),A[c,i,1])
		fi
	    fi
	od od:
	seq([seq(`if`(m[c,i]>0, [A[c,i,1], m[c,i]], NULL), i=1..nops(A[c]))], c=1..2)
end:


# Look for pairs A,B such that we have to find solutions of M the form
#  	Sol(tau - c * A/B) * Polynomial
#
PairsAB := proc(L, x, tau, ext)
options remember, system;
local Af, Bf, i, Pseq, P, vA, vB, w, dP, resP, s;
	Pseq := NULL;
	Bf, Af := LTSing(_passed);	# Example:
					#	Af = [[x+1/2,1], [x,2], [x-1,1]]
	do				#	Bf = [[x-5/2,2], [x,1], [x-2,1], [x+2,1]]
		# Select the next factor P:
		i := [op(Af), op(Bf)];
		if i = [] then break fi;
		P := i[1,1];
		Af, vA := FindIntShifts(Af, P, x);
		Bf, vB := FindIntShifts(Bf, P, x);

		# Example: If P was x, then we now have
		# 	Af = [[x + 1/2, 1]]
		#	Bf = [[x - 5/2, 2]]
		#	vA = [[x, 0], [x, 0], [x - 1, 1]]
		#	vB = [[x + 2, -2], [x, 0], [x - 2, 2]]
		#
		# For A, use the last ones from vA in resP.
		# For B, use the first ones from vB in resP (these are moved from vB to w)

		w := []; resP := [vA,w];

		# Now loop over the nops(vA)+nops(vB)+1 combinations:
		do
			if [vA,vB]=[[],[]] then
				break
			elif vA = [] or (vB<>[] and vA[1,2] > vB[1,2]) then
				w := [op(w), vB[1]];
				vB := vB[2..-1]
			else
				vA := vA[2..-1]
			fi;
			resP := resP, [vA,w]
		od;

		dP := degree(P,x);
		# Give list of all [A,B,s,d] and add this list to Pseq
		# Here s corresponds to a slope of the NewtonPolygon
		# So s = -v where v is as in AAECC, Vol. 17, 83-115 (2006).
		# The d is the same as in that paper (e.g. a polynomial solution
		# of degree 5 would correspond to some [c,s,d] with c=1, s=0, d=5.
		# Here [c,s,d] corresponds to a factor tau - c*x^s*(1+d*t+O(t^2))
		# where t = 1/x.  In other words:   tau - c*t^(-s)*(1+d*t+O(t^2))
		Pseq := Pseq, [seq([mul(i[1],i=w[1]), mul(i[1],i=w[2]),
			(nops(w[1])-nops(w[2])) * dP, Normalizer(
				 add(coeff(i[1], x, dP-1), i=w[1])
				-add(coeff(i[1], x, dP-1), i=w[2]))], w=[resP])]
	od;
	# Check p-curvature? It'd have to be fast (p = 2 or 3) for it to make sense.
	w := DetFactors(L,x,tau,ext); forget(EvalMonRem);
	if assigned(_Env_RightFactors_data) then
		w := {seq(`if`(DetFactorsSelect(i, x, tau), i, NULL), i=w)}
	fi;
	s := {seq(degree(i[1],x),i=w)};
	{seq(`if`(member(P[3],s), seq(MatchAB(op(P),i,x),i=w), NULL), P=CombinationsAB(Pseq))}
end:
# Check if an A/B can be matched.
#
MatchAB := proc(A,B,s,d,L,x)
local db;
	if s=degree(L[1],x) then
		db := Normalizer( coeff(L[1],x,s-1)/lcoeff(L[1],x)-d);
		if type(db, nonnegint) then
		    if not assigned(_Env_RightFactors_data) or DeterminantSelect(A/B, x) then
			return [A/B, db, L]  # db = degree-bound for lcoeff(R, tau)
		    fi
		fi
	fi;
	NULL
end:

CombinationsAB := proc(L)
local i,j;
	if _npassed=0 then
		[ [1,1,0,0] ]
	elif _npassed=1 then
		L
	else
		[seq(seq([i[1]*j[1], i[2]*j[2], i[3]+j[3], i[4]+j[4]],
		 j=L),i=thisproc(_passed[2..-1]))]
	fi
end:

