##########################################################################
#
# Module: Risch							   
#                                     					   
###########################################################################
Risch :=module()
   option package;
   description "algorithms for computing Risch pair in transcendental Liouvillian towers";
   export
	NoNewConst,
	ReductionPair,
	TowerDerivative;

   local
	RischPair0,
	RischPair,
	RischPairHyperExp,
 	AuxiliaryHyperExp,
	ComplementHyperExp,
	RischPairPrimi,
 	AuxiliaryPrimi,
	ComplementPrimi;



##############################################
# Name: RischPair0
# Calling sequence:   RischPair0(x, f, g)
# Input: x, a variable;
#            f, a normalizd rational function in C(x);
#            g, a rational function in C(x);
# Output: (s, r), the f-residual pair of g
##############################################
RischPair0 := proc(x, f, g) 
	local fs, gs, s, r, a, b, L; 
	fs, gs := normal(f), normal(g); 
	a, b := Rational:-NumerAndMonicDenom(x, f); 
	L := Reduction:-ResidualForm0(g, a, b, x); 
	s, r := L[1], normal(L[2] + L[3]/b); 
	return s, r; 
end proc:



###################################################
# Name: RischPair
# Calling sequence:   RischPair(x, T, dT, f, g)
# Input: x, a variable;
#          T, a list of indeterminates, each T[i] is a Liouvillian monomial;
#          dT, a list of derivatives corresponding to T;
#           f, an element in C(x,T);
#           g, an element in C(x,T);
# Output: (s, r), the f-residual pair of g
###################################################
RischPair := proc(x, T, dT, f, g) 
	local z,eta,s,r,fs;
	fs:=normal(f);
	z,eta:=op(Reduction:-GCanonicalForm(x,T,dT,fs));

	if T = [] then 
		s,r:=RischPair0(x,z,eta*g);
	else
		if degree(dT[-1], T[-1]) = 1 then 
			s,r:= RischPairHyperExp(x, T, dT, z, eta*g); 
		else 
			s,r:=RischPairPrimi(x, T, dT, z, eta*g); 
		end if; 
	end if;
	
	return normal(s/eta), normal(r/eta);
end proc:




########################################################
# Name: RischPairHyperExp
# Calling sequence:     RischPairHyperExp(x, T, dT, f, g)
# Input: x, a variable;
#            T, a list of indeterminates, with T[-1] is a hyperexponential monomial;
#            dT, a list of derivatives corresponding to T;
#            f, a normalized element in C(x,T);
#            g, a rational function in C(x,T);
# Output: (s, r), the f-risch pair of g
########################################################
RischPairHyperExp := proc(x, T, dT, f, g) 
	local n, t, fs, gs, a, b, s, h, p, u, q, W, v1, w1, b1, c1, v2, w2, b2, c2, c, ut, qt, ct; 
	n := nops(T); 

	if n = 0 then 
		return RischPair0(x, f, g); 
	end if; 

	t := T[n]; 
	fs, gs := normal(f), normal(g); 
	if gs = 0 then 
		return 0, 0; 
	end if; 

	a, b := Rational:-NumerAndMonicDenom(t, fs); 
	s, h, p := op(Reduction:-GKernelAndShellReduction(x, T, dT, gs, a, b)); 
	if p = 0 then 
		return s, h; 
	end if; 

	u, q := AuxiliaryHyperExp(x, T, dT, fs, p); 
	if q = 0 then 
		return normal(s + u), h; 
	end if; 

	W := ComplementHyperExp(x, T, dT, fs); 
	if W = [] then 
		return normal(s + u), normal(h + q/b); 
	end if;
 
	if nops(W) = 1 then 
		v1, w1, b1, c1 := op(W[1]); 
		c := Rational:-Coefficient([x, op(T)], q, b1)/c1; 
		ut, qt := normal(u + c*v1), normal(q - c*w1); 
	end if; 

	if nops(W) = 2 then 
		v1, w1, b1, c1 := op(W[1]); 
		v2, w2, b2, c2 := op(W[2]); 
		c := normal(Rational:-Coefficient([x, op(T)], q, b1)/c1); 
		ut, qt := normal(u + c*v1), normal(q - c*w1); 
		ct := normal(Rational:-Coefficient([x, op(T)], qt, b2)/c2); 
		ut, qt := normal(ut + ct*v2), normal(qt - ct*w2); 
	end if; 

	return normal(s + ut), normal(h + qt/b); 

	end proc:


############################################################
# Name:  AuxiliaryHyperExp
# Calling sequence:     AuxiliaryHyperExp(x, T, dT, f,  p)
# Input: x, a variable;
#            T, a list of indeterminates, with T[-1] is a hyperexponential monomial;
#            dT, a list of derivatives corresponding to T;
#            f, a normalized element in C(x,T);
#            p, a Laurent polynomial in T[-1], `
# Output: (u, q) is an auxiliary pair of p with respect to f
#############################################################
AuxiliaryHyperExp := proc(x, T, dT, f, p) 
	local n, t, fs, ps, ld, a, b, m, am, bm, a0, b0, hd, hc, td, tc, u, q, lambda, s, r, mu, ut, qt, A; 
	n := nops(T); 
	t := T[n]; 
	fs, ps := normal(f), normal(p); 
	ld := normal(dT[n]/t); 
	if ps = 0 then 
		return 0, 0; 
	end if; 

	a, b := Rational:-NumerAndMonicDenom(t, fs); 
	m := max(degree(a, t), degree(b, t)); 
	am, bm := coeff(a, t, m), coeff(b, t, m); 
	a0, b0 := coeff(a, t, 0), coeff(b, t, 0); 
	hd, hc := Rational:-HdegAndHcoeff(t, ps); 
	td, tc := Rational:-TdegAndTcoeff(t, ps); 
	if hd < m and 0 <= td then 
		return 0, ps; 
	end if; 

	if m <= hd and bm = 0 then 
		u, q := normal(hc/am)*t^(hd - m), 0; 
		ps := ps - b*DField:-Derivative(x, T, dT, u) - a*u; 
	end if;

	if hd < m and td < 0 and b0 = 0 then 
		u, q := normal(tc/a0)*t^td, 0; 
		ps := ps - b*DField:-Derivative(x, T, dT, u) - a*u; 
	end if; 

	if m <= hd and bm <> 0 then 
		lambda := am/bm + (hd - m)*ld; 
		#A := Reduction:-GCanonicalForm(x, [op(1 .. n - 1, T)], [op(1 .. n - 1, dT)], lambda); 
		s, r := RischPair(x, [op(1 .. n - 1, T)], [op(1 .. n - 1, dT)], lambda, hc); 
		u, q := s*t^(hd - m), bm*r*t^hd; 
		ps := ps - b*DField:-Derivative(x, T, dT, u) - a*u - q; 
	end if; 

	if hd < m and td < 0 and b0 <> 0 then 
		mu := a0/b0 + td*ld; 
		#A := Reduction:-GCanonicalForm(x, [op(1 .. n - 1, T)], [op(1 .. n - 1, dT)], mu); 
		s, r := RischPair(x, [op(1 .. n - 1, T)], [op(1 .. n - 1, dT)], mu, tc/b0); 
		u, q := s*t^td, b0*r*t^td;
		ps := ps - b*DField:-Derivative(x, T, dT, u) - a*u - q; 
	end if; 

	ut, qt := AuxiliaryHyperExp(x, T, dT, fs, ps); 

	return normal(u + ut), normal(q + qt); 
end proc:


############################################################
# Name: ComplementHyperExp
# Calling sequence:      ComplementHyperExp(x, T, dT, f)
# Input: x, a variable;
#            T, a list of indeterminates, with T[-1] is a hyperexponential monomial;
#            dT, a list of derivatives corresponding to T;
#            f, a normalized element in C(x,T);
# Output: a finite presentation of an f-complement
#############################################################
ComplementHyperExp := proc(x, T, dT, f) 
	local n, t, fs, L, a, b, m, am, bm, a0, b0, M, N, Ts, 
		  dTs, dt, J, k, u, l, v, r, p, s, q, bs, c, st, qt, bt, ct, cs, CR, temp, subsflag, i, j; 
	n := nops(T); 
	t := T[n]; 
	fs := normal(f); 

	CR := DataList:-Pairsearch(x, T, dT, fs); 
	if CR <> false then 
		subsflag := [_x=x,seq(_t[i] = T[i], i = 1 .. n)]; 
		CR := subs(subsflag, CR); 
		return CR; 
	end if; 

	a, b := Rational:-NumerAndMonicDenom(t, fs); 
	m := max(degree(a, t), degree(b, t)); 
	am, bm := coeff(a, t, m), coeff(b, t, m); 
	a0, b0 := coeff(a, t, 0), coeff(b, t, 0); 

	if m = 0 then 
		DataList:-UpdateList(x, T, dT, fs, []); 
		return []; 
	end if; 

	M, N := [], []; 
	if bm <> 0 then M := ParametricLogDer:-LiouvillianPLD(x, T, dT, am/bm); end if; 
	if M <> [] and 0 < M[1] then M := []; end if; 
	if b0 <> 0 then N := ParametricLogDer:-LiouvillianPLD(x, T, dT, a0/b0); end if; 
	if N <> [] and N[1] <= 0 then N := []; end if; 

	if M = [] and N = [] then 
		DataList:-UpdateList(x, T, dT, fs, []); 
		return []; 
	end if; 

	J := []; 
	if M <> [] and N = [] then 
		k, u := -M[1], 1/M[2]; 
		r, p := AuxiliaryHyperExp(x, T, dT, fs, b*DField:-Derivative(x, T, dT, u*t^k) + a*u*t^k); 
		b, c := Rational:-BasisElement([x, op(T)], p); 
		J := [[u*t^k - r, p, b, c]]; 
	end if; 

	if M = [] and N <> [] then 
		l, v := -N[1], 1/N[2]; 
		s, q := AuxiliaryHyperExp(x, T, dT, fs, b*DField:-Derivative(x, T, dT, v*t^l) + a*v*t^l); 
		b, c := Rational:-BasisElement([x, op(T)], q); 
		J := [[v*t^l - s, q, b, c]]; 
	end if; 

	if M <> [] and N <> [] then 
		k, u := -M[1], 1/M[2]; 
		l, v := -N[1], 1/N[2]; 
		r, p := AuxiliaryHyperExp(x, T, dT, fs, b*DField:-Derivative(x, T, dT, u*t^k) + a*u*t^k); 
		s, q := AuxiliaryHyperExp(x, T, dT, fs, b*DField:-Derivative(x, T, dT, v*t^l) + a*v*t^l); 
		bs, cs := Rational:-BasisElement([x, op(T)], p); 
		c := Rational:-Coefficient([x, op(T)], q, bs)/cs; 
		st, qt := normal(v*t^l - s - c*(u*t^k - r)), normal(q - c*p); 
		bt, ct := Rational:-BasisElement([x, op(T)], qt); 
		J := [[u*t^k - r, p, bs, cs], [st, qt, bt, ct]]; 
	end if; 

	DataList:-UpdateList(x, T, dT, fs, J); 
	return J; 
end proc:

########################################################
# Name: RischPairPrimi
# Calling sequence:     RischPairPrimi(x, T, dT, f, g)
# Input: x, a variable;
#            T, a list of indeterminates, with T[-1] is a primitive monomial;
#            dT, a list of derivatives corresponding to T;
#            f, a normalized element in C(x,T);
#            g, a rational function in C(x,T);
# Output: (s, r), the f-risch pair of g
########################################################
RischPairPrimi := proc(x, T, dT, f, g) 
	local n, t, fs, gs, a, b, s, h, p, au, q, W, w1, w2, w3, w4, w5, eta, r, l, b1, c1, c,
		  d, m, am, bm, i, uu, v, u, ra, j, b2, c2, ss, su, sv, bb, cc, dd, R,v0,b0,c0,d0; 
	n := numelems(T); 
	if n = 0 then return RischPair0(x, f, g); end if; 

	t := T[n]; 
	fs, gs := normal(f), normal(g); 
	if gs = 0 then return 0, 0; end if; 

	a, b := Rational:-NumerAndMonicDenom(t, fs); 
	s, h, p := op(Reduction:-GKernelAndShellReduction(x, T, dT, gs, a, b)); 
	if p = 0 then return s, h; end if; 

	au, q := AuxiliaryPrimi(x, T, dT, fs, p); 
	if q = 0 then return normal(s + au), h; end if; 

	W := ComplementPrimi(x, T, dT, fs); 
	if W = [] then return normal(s + au), normal(h + q/b); end if; 

	w1, w2, w3, w4, w5:= op(W); 
	eta, r, l := op(w1); 
	b1, c1, c := op(w3); 
	v0, b0, c0:=op(w5);
	d := degree(q, t); 
	m := max(degree(a, t), degree(b, t)); 
	am, bm := coeff(a, t, m), coeff(b, t, m); 
	if w2 = [] and w4 = [] then 
		for i from d by -1 to m do 
			R := b*DField:-Derivative(x, T, dT, eta*t^(i - m + 1)) + a*eta*t^(i - m + 1); 
			uu, v := AuxiliaryPrimi(x, T, dT, fs, R); 
			u := eta*t^(i - m + 1) - uu; 
			ra := Rational:-Coefficient([x, op(T[1 .. n - 1])], coeff(q, t, i), b1)/(c1*(i - m + 1) + c); 
			q := q - ra*v; 
			au := au + ra*u; 
		end do; 

		# reduction with v0
		if v0<>0 then
			d0:=degree(v0,t);
			ra := Rational:-Coefficient([x, op(T[1 .. n - 1])], coeff(q, t, d0), b0)/c0;
			q := q - ra*v0; 
			au := au + ra*eta; 
		end if;
	end if; 

	if w2 = [] and w4 <> [] then 
		j, b2, c2 := op(w4); 
		for i from d by -1 to m  do 
			R := b*DField:-Derivative(x, T, dT, eta*t^(i - m + 1)) + a*eta*t^(i - m + 1); 
			uu, v := AuxiliaryPrimi(x, T, dT, fs, R); 
			u := eta*t^(i - m + 1) - uu; 
			if i = j then 
				ra := Rational:-Coefficient([x, op(T[1 .. n - 1])], coeff(q, t, i), b2)/c2; 
			else 
				ra := Rational:-Coefficient([x, op(T[1 .. n - 1])], coeff(q, t, i), b1)/(c1*(i - m + 1) + c); 
			end if; 
			q := q - ra*v; 
			au := au + ra*u; 
		end do; 

		# reduction with v0
		if v0<>0 then
			d0:=degree(v0,t);
			ra := Rational:-Coefficient([x, op(T[1 .. n - 1])], coeff(q, t, d0), b0)/c0;
			q := q - ra*v0; 
			au := au + ra*eta; 
		end if;
	end if; 

	if w2 <> [] then 
		ss, su, sv, bb, cc := op(w2); 
		dd := degree(sv, t); 
		for i from d by -1 to m - 1 do 
			if i <> ss + m - 1 then 
				if i<>m-1 then
					R := b*DField:-Derivative(x, T, dT, eta*t^(i - m + 1)) + a*eta*t^(i - m + 1); 
					uu, v := AuxiliaryPrimi(x, T, dT, fs, R); u := eta*t^(i - m + 1) - uu; 
					ra := Rational:-Coefficient([x, op(T[1 .. n - 1])], coeff(q, t, i), b1)/(c1*(i - m + 1) + c); 
					q := q - ra*v; 
					au := au + ra*u; 
				else 
					# reduction with v0
					if v0<>0 then 
						d0:=degree(v0,t);
						ra := Rational:-Coefficient([x, op(T[1 .. n - 1])], coeff(q, t, d0), b0)/c0;
						q := q - ra*v0; 
						au := au + ra*eta; 
					end if;
				end if;
				if i = dd then 
					ra := Rational:-Coefficient([x, op(T[1 .. n - 1])], coeff(q, t, i), bb)/cc; 
					q := q - ra*sv; 
					au := au + ra*su; 
				end if; 
			end if; 
		end do; 
		if dd < m - 1 and bb <> 0 then 
			ra := Rational:-Coefficient([x, op(T[1 .. n - 1])], coeff(q, t, dd), bb)/cc; 
			q := q - ra*sv; 
			au := au + ra*su; 
		end if; 
	end if; 

	return normal(s + au), normal(h + q/b); 
end proc:


########################################################
# Name:  AuxiliaryPrimi
# Calling sequence:    AuxiliaryPrimi(x, T, dT, f, g)
# Input: x, a variable;
#            T,   a nonempty list of indeterminates, with T[-1] is a primitive monomial;
#            dT, a list of derivatives corresponding to T;
#            f, a normalized element in C(x,T);
#            p, a polynomial in T[-1], 
# Output: (u, q) is an auxiliary pair of p with respect to f
########################################################
AuxiliaryPrimi := proc(x, T, dT, f, p) 
	local t, n, am, bm, a, b, d, c, s, r, ps, fs, gs, u, q, xi, eta, m, pd; 
	n := numelems(T); 
	t := T[n]; 
	fs, ps := normal(f), normal(p); 
	a, b := Rational:-NumerAndMonicDenom(t, fs); 
	m := max(degree(a, t), degree(b, t)); d := degree(ps, t); 
	am, bm := coeff(a, t, m), coeff(b, t, m); 
	u, q := 0, 0; 
	if bm = 0 then 
		while m <= d do 
			pd := coeff(ps, t, d); 
			u := u + pd*t^(d - m)/am; 
			ps := normal(ps - b*DField:-Derivative(x, T, dT, pd*t^(d - m)/am) - a*pd*t^(d - m)/am); 
			d := degree(ps, t); 
		end do; 
	else 
		#xi, eta := op(Reduction:-GCanonicalForm(x, T[1 .. n - 1], dT[1 .. n - 1], am/bm)); 
		while m <= d do 
			pd := coeff(ps, t, d); 
			s, r := RischPair(x, T[1 .. n - 1], dT[1 .. n - 1], am/bm, pd/bm); 
			u := u + s*t^(d - m); 
			q := q + r*t^d; 
			ps := normal(ps - b*DField:-Derivative(x, T, dT, s*t^(d - m)) - a*s*t^(d - m) - r*bm*t^d); 
			d := degree(ps, t); 
		end do; 
	end if; 
	return normal(u), normal(q + ps); 
end proc:



############################################################
# Name: ComplementPrimi
# Calling sequence:     ComplementPrimi(x, T, dT, f)
# Input: x, a variable;
#            T,   a nonempty list of indeterminates, with T[-1] is a primitive monomial;
#            dT, a list of derivatives corresponding to T;
#            f, a normalized element in C(x,T);
# Output: a finite presentation of an f-complement
#############################################################
ComplementPrimi := proc(x, T, dT, f) 
	local n, t, fs, a, b, am, bm, m, W,w5, xi, eta, r, l, ra, u, v, 
		 uu, tempu, tempv, tempuu, v0, s, j, b0,c0, bb, cc, b1, c1, b2, c2, c, d, cv, CR, subsflag, i,mu; 
	n := numelems(T); 
	t := T[n]; 
	fs := normal(f); 
	a, b := Rational:-NumerAndMonicDenom(t, fs); 
	m := max(degree(a, t), degree(b, t)); 
	am, bm := coeff(a, t, m), coeff(b, t, m); 

	CR := DataList:-Pairsearch(x, T, dT, fs); 
	if CR <> false then 
		subsflag := [_x = x, seq(_t[i] = T[i], i = 1 .. n)]; 
		CR := subs(subsflag, CR); 
		return CR; 
	end if; 

	if bm = 0 then 
		W := []; 
		DataList:-UpdateList(x, T, dT, fs, W); 
		return W; 
	end if; 

	#xi, eta := op(Reduction:-GCanonicalForm(x, T[1 .. n - 1], dT[1 .. n - 1], am/bm)); 
	eta:=ParametricLogDer:-LiouvillianLD(x, T[1 .. n - 1], dT[1 .. n - 1],am/bm);
	if eta=[] then 
		W := []; 
		DataList:-UpdateList(x, T, dT, fs, W); 
		return W; 
	end if; 
	
	eta:=1/eta;
	r := RischPair(x, T[1 .. n - 1], dT[1 .. n - 1], am/bm, eta*dT[n])[2]; 
	l := RischPair(x, T[1 .. n - 1], dT[1 .. n - 1], am/bm, DField:-Derivative(x,T,dT,eta)*coeff(b, t, m - 1)+eta*coeff(a, t, m - 1))[2]; 
	s := -normal(l/r); 
	b1, c1 := Rational:-BasisElement([x, op(T[1 .. n - 1])], r); 
	c := Rational:-Coefficient([x, op(T[1 .. n - 1])], l, b1); 

	v0 :=normal( b*DField:-Derivative(x, T, dT, eta) + a*eta); #v0=P_{fs}(eta) 
	if v0<>0 then 
		b0,c0:= Rational:-BasisElement([x, op(T[1 .. n - 1])], lcoeff(v0,t));
		w5:=[v0,b0,c0];
	else 
		w5:=[0,0,1];
	end if;


	if type(s, integer) = false or type(s, positive) = false then 
		j := -normal(c/c1); 
		if type(j, integer) = false or type(j, positive) = false then 
			W := [[eta, r, l], [], [b1, c1,c], [],w5]; 
		else 
			b2, c2 := Rational:-BasisElement([x, op(T[1 .. n - 1])], j*r + l); 
			W := [[eta, r, l], [], [b1, c1, c], [j, b2, c2],w5]; 
		end if; 
		DataList:-UpdateList(x, T, dT, fs, W); 
		return W; 
	else #degree of v0 is m-1
		uu, v := AuxiliaryPrimi(x, T, dT, fs, b*DField:-Derivative(x, T, dT, eta*t^s) + a*eta*t^s); 
		u := eta*t^s - uu; 
		d := degree(v, t); 
		if d < m - 1 then 
			bb, cc := Rational:-BasisElement([x, op(T[1 .. n - 1])], lcoeff(v, t)); 
			W := [[eta, r, l], [s, u, v, bb, cc], [b1, c1, c], [],w5]; 
			DataList:-UpdateList(x, T, dT, fs, W); 
			return W; 
		end if; 
		do 
			d := degree(v, t); 
			if d <> m-1 then 
				cv := c1*(d - m + 1) + c; 
				tempuu, tempv := AuxiliaryPrimi(x, T, dT, fs, b*DField:-Derivative(x, T, dT, eta*t^(d - m + 1)) + a*eta*t^(d - m + 1)); 
				tempu := eta*t^(d - m + 1) - tempuu; 
				ra := Rational:-Coefficient([x, op(T[1 .. n - 1])], lcoeff(v, t), b1)/cv;
				v := normal(v - ra*tempv); u := u - ra*tempu; 
			else 
				ra := Rational:-Coefficient([x, op(T[1 .. n - 1])], lcoeff(v, t), b0)/c0;
				v := normal(v - ra*v0); u := u - ra*eta; 
			end if; 
		until degree(v, t) = d or degree(v, t) < m - 1; 
		bb, cc := Rational:-BasisElement([x, op(T[1 .. n - 1])], lcoeff(v, t)); 
		W := [[eta, r, l], [s, u, v, bb, cc], [b1, c1, c], [],w5]; 
		DataList:-UpdateList(x, T, dT, fs, W);
		 return W; 
	end if; 
end proc:





############################################################
# Name: NoNewConst
# Calling sequence:   NoNewConst(x, T, dT)
# Input: x, a variable;
#            T,   a nonempty list of indeterminates;
#            dT, a list of derivatives corresponding to T.
# Output:  true if tower is regular, false otherwise
#############################################################
NoNewConst:=proc(x, T, dT)
	local i;
	for i from 1 to numelems(T) do
		if degree( normal(dT[i]), T[i] )=0 then 
			if normal( Risch(x, T[1..i-1], dT[1..i-1], 0, dT[i])[2] )=0 then return false; end if;
		else 
			if  ParametricLogDer:-RationalCombination(x, T[1..i-1], dT[1..i-1], normal(dT[i]/T[i]))<>false then return false; end if;
		end if;
	end do;

	return true;

end proc:
	


######################################################################
# Name: ReductionPair
# Calling sequence:   ReductionPair(f, g)   #need  to set tower using command SetTower(x,T,dT)
# Input: f, an element in C(x,T);
#           g, an element in C(x,T);
# Output: the f-reduction pair of g
######################################################################
ReductionPair:=proc(f, g)
	local K;

	K:=DataList:-GetTower();
	return RischPair(K, f, g);

end proc:


######################################################################
# Name: TowerDerivative
# Calling sequence:   TowerDerivative(f)   #need  to set tower using command SetTower(x,T,dT)
# Input: f, an element in C(x,T);
# Output: the derivative of f
######################################################################
TowerDerivative:=proc(f)
	local K;

	K:=DataList:-GetTower();
	return DField:-Derivative(K, f);

end proc:



end module:
