
# Authors: Yisen Wang <wangyisen@amss.ac.cn>
######################################################################################
read "MultivariateILD_test";
with(LREtools);
AddOreSato := module()
  description "additive decompositions of WZ-forms";
  export
        PartialFraction,
        SortFraction,
        CompletePartialFraction,
        AbramovReduction,
        VerifyWZForm,
        VerifyWZFormDecomp,
        WZFormDecomp,
        IntegerLinearTypePoly,
        TransUniWZForms;
  local
        ListAdd,
        ShiftSum,
        CheckZero,
        SortType,
        ChooseMin;

  option package;



###############################################################################
# This text consists of two parts:
# Abramov's Reduction of Rational Functions, Additive Decomposition of WZ-forms.

############################################################################

#Task 1: Abramov's Reduction_Written by Yisen Wang

##########################Task Description#####################################

# Given a rational function f (x), write it as an addition delta(g, x)+sum(sum(A_ij/B_i^j, j=1..w_i), i=1..n) where B_i's are not <sigma_x>-equivalent.


##########################################################################
# Name: Partial Fraction
# Calling sequence:
#				PartialFraction(a, d, x)
# Given a field K.
# Input:   a	a polynomial in K[x];
#          d	a list of polynomials in K[x];
#          x	a variable name.
#
# Output: Error	if d[i]'s are not coprime;
# else	   a0	a polynomial in K[x];
#           A	a list of polynomials in K[x],
#	such that  a/product(di, i=1..n)=a0+sum(A[i]/d[i], i=1..n), where n is the length of d.
##########################################################################
PartialFraction := proc(a, d, x)
    local n, i, j, q, t, r, A, a0, b, w, flag;
    flag := 0;
    n := nops(d);
    w := product(d[i], i = 1 .. n);
    a0 := quo(a, w, x); r := rem(a, w, x);
    A := [seq(b[i], i = 1 .. n)];
    if n > 1 then 
    for i from 1 to n - 1 do
        w := simplify(w/d[i]);
        if 0 < degree(gcdex(w, d[i], x), x) then
            flag := 1;
            break;
        else
            gcdex(w, d[i], r, x, 'A[i]', 't');
            q := quo(t, w, x); r := rem(t, w, x);
            a0 := simplify(a0 + q);
        end if;
    end do;
    A[n] := r;
    else A :=[r];
    end if;
    if flag = 1 then
        return Error;
    else
        return a0, A;
    end if;
end proc;

#################################################################################################
# SortFraction
# Calling sequence:
#				SortFraction(f)
# Input: 
#	f  a multivariate rational function over the rational number field;
# Output: 
#	a  a multivariate polynomial which is the numerator of f;
#	d  a list of the irreducible factors of the denominator of f;
#	e  a list of the integers with the same length of d such that denom(f)=mul(d[i]^e[i], i=1..n) where n = nops(d).
#################################################################################################
SortFraction:=proc(f)
  local c,  L, n, a, b, d, e, i;
  a := numer(f); b := denom(f);
  c := factors(b)[1];
  a := simplify(a/c);
  L := factors(b)[2];
  n := nops(L);
  d := []; e := [];
  for i from 1 to n do 
	d := [op(d), L[i][1]];
	e := [op(e), L[i][2]];
  end do;
  return (a, d, e);
end proc;
	
################################################################################################
# Name: CompletePartialFraction
# Calling sequence:
#				CompletePartialFraction(a, d, e, x)
# Given a field K.
# Input:   a	a polynomial in K[x];
#          d	a list of polynomials in K[x];
#          e	a list of positive integers with the same length of d;
#          x	a variable name.
#
# Output: Error	if d[i]'s are not coprime;
# else	   a0	a polynomial in K[x];
#        output	a list of length n  whose ith component is [di,A[i]],
#	such that  a/product(d[i]^e[i], i=1..n)=a0+sum(sum(A[i][j]/d[i]^j, j=1..e[i]), i=1..n), where n is the length of d and e.
################################################################################################
CompletePartialFraction := proc(a, d, e, x)
    local n, a0, i, j, q, A, b, c, B, D, output;
    n := nops(d);
    output := [seq(b[i], i = 1 .. n)];
    A := [seq(b[i], i = 1 .. n)];
    for i from 1 to n do
        A[i] := [seq(c[j], j = 1 .. e[i])];
    end do;
    D := [seq(b[i], i = 1 .. n)];
    for i from 1 to n do
        D[i] := d[i]^e[i];
    end do;
    if PartialFraction(a, D, x)=Error then return Error; 
    else 
        a0, B := PartialFraction(a, D, x);
   end if;
    for i from 1 to n do
        for j from e[i] by -1 to 1 do
            q := quo(B[i], d[i], x);  A[i][j] := rem(B[i], d[i], x);
            B[i] := q;
        end do;
        a0 := simplify(a0 + B[i]);
    end do;
    for i from 1 to n do
        output[i] := [d[i], A[i]];
    end do;
    return a0, output;
end proc;

#####################################################################################################
# Name: ListAdd
# Calling Sequence: 
#			ListAdd(a, b)
# Input :
#	a a list of length n;
#	b a list of length m.
#Output:
#	c a list of length max(m,n) such that c[i] = a[i] + b[i], 1<=i <= max(m,n).
#	Note that if i > n then a[i] := 0 and if i > m then b[i] := 0. 
#####################################################################################################
ListAdd := proc(a,b)
  local n, m, N, i, a0, b0, c;
  n := nops(a);
  m := nops(b);
  N := max(m,n);
  a0 := a; b0 := b;
  if N > n then
 	 for i from n+1 to N do a0 := [op(a0), 0]; end do;
  end if;
  if N > m then 
	for i from m+1 to N do b0 := [op(b0), 0]; end do;
  end if;
  c := [];
  for i from 1 to N do c := [op(c), simplify(a0[i]+b0[i])]; end do;
  return (c);
end proc;

#####################################################################################################
# Name: ShiftSum
# Calling Sequence:
#			ShiftSum(f, x, n1, n2)
# Input: 
#	f    an expression w.r.t. x;
#	x    a variable name;
#	n1  the lower shift bound;
#	n2  the upper shift bound.
# Output: 
#	g  an expression w.r.t. x s.t. g = sum(shift(f, x, t), t=n1..n2).
#####################################################################################################
ShiftSum:=proc(f, x, n1, n2)
  local i, g;
  g := 0;
  for i from n1 by 1 to n2 do
	g := g + shift(f, x, i);
  end do;
  return g;
end proc;

######################################################################################################
# Name: CheckZero
# Calling sequence:
#			CheckZero(a)
# Input: 
#	a  a list of expressions;
# Output: 
#	0 if all the elements of a are zero;
#	1 if some element is not zero.
######################################################################################################
CheckZero := proc(a)
  local flag, n, i;
  n := nops(a);
  flag := 0;
  for i from 1 to n do 
	if a[i] <> 0 then flag := 1; break; end if;
  end do;
  return flag;
end proc;
  
#####################################################################################################
# Name: AbramovReduction
# Calling Sequence:
#			AbramovReduction(f, x)
# Given a field K.
# Input:
#          f    a rational function w.r.t. x;
#          x	a variable name.
#
# Output: 
# 	a0	    a polynomial in K[x];
#  Output2	a list of length m  whose ith component is B[i];
#  Output3	a list of length m whose ith component is A[i] of length w[i];
#  such that  f = \De(a0)+sum(sum(A[i][j]/B[i]^j, j=1..w[i]), i=1..m), where m is the length of B and w.
# Furthermore, B[i] and B[j] are not <sigma_x>-shift equivalent.
####################################################################################################
AbramovReduction:=proc(a, d, e, x)
  with(LREtools);
  local a0, L, i,  j,  n,  B,  b,  A, s, k, l, Output2, Output3, t, w;
  n := nops(d);
  B := d; w := e;
  a0, L := CompletePartialFraction(a, d, e, x);
  A := [seq(b[i], i = 1 .. n)];
  for i from 1 to n do
	A[i] := L[i][2];
  end do;
  a0 := sum(a0, x);
  for i from 1 to n do 
	if B[i] <> 0 then 
	for j from (i+1) to n do
		if B[j] <> 0 then 
			s := dispersion(B[i], B[j], x); t := dispersion(B[j], B[i], x);
			if s <> FAIL then k := s[1]; w[i] := max(w[i], w[j]);
			   for l from 1 to w[j] do A[j][l] := shift(A[j][l], x, -k);
			   a0 := a0 + ShiftSum(A[j][l]/(B[i]^l), x, 0, k-1); end do;
			   A[i] := ListAdd(A[i], A[j]); 
			   B[j] := 0; w[j]=0;
			 end if;
			if t<> FAIL then k := t[1];w[i] := max(w[i], w[j]);
			   for l from 1 to w[j] do A[j][l] := shift(A[j][l], x, k);
			   a0 := a0 - ShiftSum(A[j][l]/(B[i]^l), x, -k, -1); end do;
			   A[i] := ListAdd(A[i], A[j]); 
			   B[j] := 0; w[j]=0;
			 end if;
		end if;
	end do;
	end if;		
  end do;
  Output2 := []; Output3 := [];
  for i from 1 to n do 
    if B[i] <> 0 and CheckZero(A[i])<>0 then Output2 := [op(Output2), B[i]]; Output3 := [op(Output3), A[i]]; end if;
  end do;
  return (a0, Output2, Output3);
end proc;

##############################################################################################################



# Task 2:  Additive Decomposition of WZ-Forms

#############################################Task Description####################################################
						
#To decompose a WZ-form into an exact form plus several uniform WZ-forms.
#A WZ-form (f[1],...,f[n]) with each f[i] as a n-variable rational function satifies that delta(f_j, vars[i])=delta(f_i, vars[j]), 1<=i, j<=n.
#An exact form is a WZ-form that there exists a term g such that f[i]=delta(g, vars[i]).
#A uniform WZ-form is a WZ-form that f[i]'s are of the same integer-linear type. 
# For any uniform WZ-form (f_1,...,f_n) of integer-linear type v, there exists a univariate rational function r(z), such that f_i=sum(r(v_1x_1+...+v_nx_n+k), k=0..v_i)

################################################################################################################
# Name: IntegerLinearTypePoly
# Calling sequence:
#		IntegerLinearTypePoly(p, vars)
# Input:
#	p	an irreducible polynomial;
#	vars	a list of variables.
# Output:
#	false	if p is not an integer linear polynomial 
#	L	the list of the intgerlinear type of p;
###############################################################################################################
IntegerLinearTypePoly := proc(p, vars)
  local S, L, Z;
  if IsIrredIntegerLinear(p, vars) = false then return false;
     else S := IntegerLinearDecomp(p, vars, Z, false, false); L := indices(S[2]); return op(L);
  end if;
end proc:

##############################################################################################################
# Name: VerifyWZForm
# Calling sequence:
#		VerifyWZForm(f, vars)
# Input:
#	f	a n-tuple of rational functions;
#	vars	a list of variables with the same length of L;
# Output:
#	true	if f is a WZ-form;
#	false	otherwise.
##############################################################################################################
VerifyWZForm := proc(f, vars)
  local n, i, j, flag;
  n := nops(f);
  flag := 0;
  for i from 1 to n do
	for j from (i+1) to n do
	if simplify((shift(f[i], vars[j], 1)- f[i])-(shift(f[j], vars[i], 1)- f[j]))<> 0 then flag := 1; break; end if;
	end do;
	if flag = 1 then break; end if;
  end do;
  if flag = 0 then return true; else return false; end if;
end proc:



##############################################################################################################
# Name: SortType
# Calling sequence:
#		SortType(B, vars)
# Input: 
#	B	a list of multivariate irreducible integer-linear polynomials;
#	vars	a list of variables.
# Output:
#	C	a sort of the index such that each corresponding element in B is a list of the polynomials in the same type. 
#	exactly classify the elements of B depending on the integer-linear type.
###############################################################################################################
SortType:=proc(B, vars)
  local m, i, j, T, C, S, L;
  T := [];
  m := nops(B);
  for i from 1 to m do 
	T := [op(T), IntegerLinearTypePoly(B[i], vars)];
  end do;
  S := {};
  C := [];
  for i from 1 to m do
	if evalb(i in S)= false then 
	L := [i]; S := S union {i};
	for j from (i+1) to m do
		if evalb(j in S) = false then 
		if T[i] = T[j] then
		   L := [op(L), j]; S := S union {j};
		end if;
		end if;
	end do;
	C := [op(C), L];
	end if;
  end do;
  return C;
end proc:
 
   
##############################################################################################################
# Name:  WZFormDecomp
# Calling sequence:
#		WZFormDecomp(f, vars)
#  Input:  
#	f	a n-tuple (f[1],...,f[n]) which is a WZ-form;
#	vars	a list of variables with the same length of L;
#  Output:
#	E	a rational function;
#	U	a set of uniform WZ-forms such that L=(delta(E, vars[i]) , i=1..n)+U, where delta is the difference operator;
#################################################################################################################
WZFormDecomp := proc(L, vars)
  local  n, i, j, E, S, U, m, s, B, A, a0, index, g, e, w, l, F, k, c;
  F := [];
  n := nops(L);
  U := [];
  E := 0;
  if L[1] <>0  and n >1 then 
  a0, B, A := AbramovReduction(SortFraction(L[1]), vars[1]);
  E := a0;
  index := SortType(B, vars);
  m := nops(index);
  for i from 1 to m do
	g := 0;
	for e in index[i] do 
	w := nops(A[e]);
	for j from 1 to w do 
		g := g + simplify(A[e][j]/(B[e]^j));
	end do;
	end do;
	l := [g];
	for j from 2 to n do 
		l := [op(l), expand(sum(shift(g, vars[j], 1) - g, vars[1]))];
	end do;
	U:= [op(U), l];
  end do;
  for i from 2 to n do  
	F := [op(F), simplify(L[i]- shift(E, vars[i], 1)+E-add(U[j][i], j=1..m))];
  end do;
  elif L[1] = 0 and n >1 then
  F:= subsop(1 = NULL, L);   
  elif n=1 and L[1] <> 0 then 
  E := AbramovReduction(SortFraction(L[1]), vars[1])[1];
  c := simplify(L[1]-delta(E, vars[1]));
  if c <> 0 then U:= [[c]]; end if;
  F:=[]; 
  end if;
  if F <>[] then 
  l := WZFormDecomp(F, subsop(1 = NULL, vars))[1];
  s := WZFormDecomp(F, subsop(1 = NULL, vars))[2];
  k := nops(s);
  for i from 1 to k do 
	s[i] := [0, op(s[i])];
  end do;
  return (simplify(E+l), [op(U), op(s)]);
  else return (E, U);
  end if;
end proc;


#####################################################################################################
# Name: VerifyWZFormDecomp 
# Callsing sequence :
#		VerifyWZFormDecomp(L, E, U, vars)
# Input:
#	L	an n-tuple which is a WZ- form;
#	E	representing the exact form (delta(E,  vars[i]), i=1..n);
#	U	the uniform WZ-form;
#	vars	a list of variables.
# Output:
#	true	if E+U is really a decomposition of L;
#	false	otherwise.
########################################################################################################
VerifyWZFormDecomp := proc(L, E, U, vars)
  local n, m, i, j, flag;
  flag := 0;
  n := nops(L);
  m := nops(U);
  for i from 1 to n do 
	if simplify(L[i]-delta(E, vars[i])-add(U[j][i], j=1..m)) <> 0 then flag := 1; break; end if;
  end do;
  if flag = 0 then return true; else return false; end if;
end proc;



######################################################################################################
# Name: ChooseMin
# Calling sequence:
#		ChooseMin(L)
# Input:
#	L	a list of integers;
# Output:
#	i	an index such that the absolute value of L[i] is nonzero and minimal.
####################################################################################################
ChooseMin := proc(L)
  local n, i, j, k;
  n := nops(L);
  for i from 1 to n do
	if L[i]<>0 then k := i; 
	for j from (i+1) to n do if abs(L[j])< abs(L[i]) then k := j; end if; end do;
	break; end if;
  end do;
  return k;
end proc:

######################################################################################################
# Name: TransUniWZForms
# Calling sequence:
#		TransUniWZForm(U, vars, Z)
# Input:
#	U	a list of uniform WZ-form;
#	vars	a list of variables;
#	Z	a new variable;
# Output:
#	T	a table [L[i]=R[i], i=1..m] where m = nops(U);	
# such that for each i, U[i][j]=add(subs(z=add(L[i][j]vars[j], j=1..n)+k, R[i][z]), k=0..L[i][j]-1).
######################################################################################################
TransUniWZForms := proc(U, vars, Z)
  local n, m, i, j, L, S, w, k, t, R, rec, y, T;
  L := []; S := U; w := []; R:=[]; T:=[];
  n := nops(vars);
  m := nops(U);
  for i from 1 to m do
	for j from 1 to n do 
	if U[i][j] <> 0 then w := [op(w), j]; L := [op(L), IntegerLinearTypePoly(denom(U[i][j]), vars)]; break; end if;
	end do;
  end do;
  for i from 1 to m do 
       	for j from 1 to n do
	S[i][j] := simplify(subs(vars[w[i]]=1/L[i][w[i]]*(Z-add(L[i][k]*vars[k], k=w[i]+1..n)), U[i][j]));
      	end do;
  end do;
  for i from 1 to m do
	t := ChooseMin(L[i]);
	if L[i][t] > 0 then rec := add(y(Z+k), k=0..L[i][t]-1)-S[i][t]; R:=[op(R), ratpolysols(rec, y(Z), {}, output = onesol)]; end if;
	if L[i][t] < 0 then rec := -add(y(Z+k), k=L[i][t]..-1)-S[i][t]; R:=[op(R), ratpolysols(rec, y(Z), {}, output = onesol)]; end if;
	T := [op(T), L[i]=R[i]];
  end do;
  return table(T);
end proc;

end module;
