################################################################
## Copyright (c) August 2021, UCAS and JKU. All rights reserved.
## Author: Lixin Du <dulixin17@mails.ucas.ac.cn>
## Load the modules:
##         ShiftOrbitDecomposition.mm,
##         OrbitalDecomposition.mm,
##         RationalReduction.mm,
##         RationalSummation.mm.
## Modified on December 15, 2021
#################################################################


RationalTelescoping := module()
      description "Decide the Existencee of Telescopers for Multivariate Rational Functions in Shift Case";
      export
                ExistenceAnnihilator,
                IsTelescoperable,
                IsTelescoperable2,
                AdditiveDecomposition2,
                VerifyTelescoper,
                VerifyTelescoper2,
                VerifyDecomposition2,
                ShiftApply,
                myLCLM,
                TransitionMatrix2,
                RandomRational,
                RandomRational2;
      local
                GenerateSubgroup,
                Transformation2,
                difference_case;
      option
                package;


############################# Data  Description ###########################################
# 0. Setting: the field K is of characteristic zero. In this algorithm, K can be taken
# as a field of multivariate rational functions Q(u), where Q is the rational number
# field and u = [u1, ..., um] are parameters.
#
# 1. Shift operators and difference operators:
# Let f be a multivariate rational function in K(t, x) with x = [x1, ..., xn]. For each variable v
# in {t,x1, ..., xn}, the shift operator sigma_v is defined as an F-automorphism of K(t, x) such that
#
#        sigma_v (v) = v + 1,     and     sigma_v (w) = w,  for all w in {t,x1, ..., xn}, and w <> v
#
# The difference operator is defined as
#
#                                Delta_v := sigma_v - 1,
#
# where 1 is the identity map of K(t, x).
# Remark: In the following algorithms, we shall use the list representation of shifts,
#         (see Data Description in OrbitalDecomposition module).
#
# 2. Group actions under shifts
# Let Gt = <sigma_t, sigma_{x1}, ..., sigma_{xn}> be the free abelian group generated by shift operators
# sigma_t, sigma_{x1}, ..., sigma_{xn}. Let G = <sigma_{x1}, ..., sigma_{xn}> be the subgroup of Gt.
# Let Gtp and Gp be the isotropy group of a polynomial p in Gt and G, respectively.
##############################################################################################

################################ Task  Description ######################################################
# 1. Existence Problem of Telescopers:
# 1.1 Given a multivariate rational function f in K(t, x, y) where x = [x1, ..., xn] and y = [y1, ..., ym]
# decide whether there exist a nonzero linear recurrence operator L in K(t)<St> such that
#
#           L (f) = Delta_{x1} (g1) + ... + Delta_{xn} (gn).
#
# for some g1,....,gn in K(t, x, y), where St is the shift operator w.r.t. t. If such an operator
# L exists, we say L is a telescoper of type (sigma_t; sigma_{x1}, ..., sigma_{xn}) for f and
# g1, ..., gn are the certificates of L.
#
# Remark: 1. An operator L = a0 + a1 St + ... + al St^l is represented by a polynomial in St instead of
#         an Ore polynomial in St.
#         2. y = [y1, ..., ym] are parameters and L belongs to K(t)<St> rather than K(t, y)<St>.
# 1.2 Let tau_0 \in G_t \setminus G. Let tau_1, ..., tau_k \in G be linearly independent over Z.
# Given a multivariate rational function f in K(t, x, y), decide whether there exist a nonzero linear
# recurrence operator L in K(t)<T> such that
#
#           L (f) = Delta_{tau_1} (g1) + ... + Delta_{tau_k} (gk).
#
# for some g1,....,gk in K(t, x, y), where T(h) = tau_0 (h) for any h \in K(t, x, y). If such an
# operator L exists, we say L is a telescoper of type (tau_0; tau_1, ..., tau_r) for f and
# g1, ..., gn are the certificates of L.
# 2.  Additive Decomposition for creative telescoping
# Given a multivariate rational function f in K(t, x, y), compute a triple (L, [g1, ..., gn], r)
# with L in K(t)<St>, g1, ..., gn, r in K(t, x, y) such that
#
#             L (f) = Delta_{x1} (g1) + ... + Delta_{xn} (gn) + r
#
# satisfying the property
#
#        f has a telescoper of type (sigma_t; sigma_{x1}, ..., sigma_{xn})   iff   r = 0.
########################################################################################################


##############################################################################################################
# Export
# Name: ExistenceAnnihilator
# Calling sequence:
#                    ExistenceAnnihilator(f, x, t, 'St')
# Input:  f,  a multivariate rational function in K(t, x);
#         x,  a variable name or a list of variable names;
#         t,  a variable name;
#         St,  an operator name representing the shift operator w.r.t. t.
# Output: [true, L],  if there exists a nonzero L in K(t)<St> such that
#                                   L (f) = 0,
#                     where St (f) := sigma_t (f) and L is represented by a polynomial in St.
#         [false, f],  for otherwise.
# NOTE: This procedure requires that St does not appear in f_ as a variable or parameter.
####################################I###########################################################################


ExistenceAnnihilator := proc(f_, x, t, St)
    local f, a, b, b1, b2, xt, L, flag, T, n, pair, a1, b11, termlist, p, pc, pp, L1, A;
    f := normal(expand(f_));
    a := expand(numer(f));
    b := expand(denom(f));
    b1 := content(b, x);
    b2 := primpart(b, x);
    xt := [t,op(x)];
    L := [];

#### decide the existence of Annihilator ###
    if t in indets(b2) then
        printf("%a contains a variable %a \n", b2, t);
        flag := false;
    else
        flag := true;

        ## obtain an annihilator for each term (which is separable) in a ##
        termlist := op(a);
        for p in termlist do
            pc := content(p, x);
            pp := primpart(p, x);
            St := 'St';
            L1 := pc * St - ShiftEquivalenceTesting:-Shift(pc, t, 1);
            L1 := primpart(L1, St);
            L := [op(L), L1];
        end do;
        n := nops(L);
        L := myLCLM(L, t, St)[1];

        ## L (a) = 0, so (L * b1) (a/ b1) = 0. Compute L * b1 in Ore algebra##
        A := OreTools:-SetOreRing(t, 'shift');
        AddConversionRule('shift', difference_case);
        L := OreTools[Converters]:-FromPolyToOrePoly(L, St);
        L := OreTools:-Multiply(L, OrePoly(b1), A);
        L := OreTools[Converters]:-FromOrePolyToPoly(L, St);
    end if;

    if flag = false then
        [flag, f];
    else
        [flag, L];
    fi;
end proc:


##############################################################################################################
# Export
# Name: IsTelescoperable
# Calling sequence:
#                   IsTelescoperable(f, x, t, 'St', y)
# Input:  f,  a multivariate rational function in K(t, x, y);
#         x,  a variable name or a list of variable names;
#         t,  a variable name;
#         St, an operator name;
#         y,  (optional) a parameter name or a list of parameter names.
# Output: true, L, g,  if f has a telescoper L in  K(t)<St> such that
#                                 L (f) = Delta_{x1} (g1) + .... + Delta_{xn} (gn)
#                      where g = [g1, ..., gn] with gi in K(t, x, y) and x = [x1, ..., xn].
#         false,  if f does not have any telescoper of tye (sigma_t; simga_{x1}, ..., simga_{xn}).
###############################################################################################################


IsTelescoperable := proc(f, x_, t, St, y_)
    local n, N, x, xt, x1, tuple, g, ListRemainder, Gtp, L, G, flag, i, orbit, tau0, Gd, g1, g2, d, f1,\
    k0, T, k, fracs, numers, s, a, a1, e, d_, b1, L0, r1, y, z, u, T1, LL, h, result, Gtd, j, pair, L1, b;
    x := convert(x_, list);
    if nargs = 5 then
        y := convert(y_, list);
        z := [op(x), op(y)];
    else
        z := x;
    end if;
    xt := [t, op(x)];
    n := nops(x);
    x1 := x[1];

    tuple := RationalReduction:-OrbitalReduction(f, x1, x, t);
    g := tuple[1];
    ListRemainder := tuple[2]; # list representation of remainder
    Gtp := tuple[3];

#### f = Delta_{x1} (g1) + ... + Delta_{xn} (gn) + r ###
    LL := [1]; # list of Telescopers
    G := [g]; # list of certificates

#### the remainder after orbital reduction is 0 ###
    if ListRemainder = [] then
        return true, 1, g;
    end if;
    i := 1;
    while i <= nops(ListRemainder) do
        orbit := ListRemainder[i];
        d := orbit[1];
        orbit := orbit[2];
        Gtd := Gtp[i];
        tau0 := Gtd[1];
        Gd := GenerateSubgroup(Gtd);

        #### when rank (Gtd/ Gd) = 0, a/ d^e is has a telescoper ##
        # iff  a/ d^e is (sigma_{x1}, ..., sigma_{xn})-summable ######
        if ArrayTools:-IsZero(Array(tau0)) then
            for fracs in orbit do
                s := op(fracs[1]);
                d_ := ShiftEquivalenceTesting:-Shift(d, t, s);
                numers := subsop(1 = NULL, fracs);
                tuple := RationalSummation:-ListSummable([[d_,numers]], [Gd], x);
                printf("case 1: \n");
                if nops([tuple]) = 1 then
                    f1 := normal(RationalReduction:-FromListToRemainder([[d_, numers]], x));
                    printf("case 1: %a is not summable w.r.t. %a\n", f1, x);
                    return false;
                end if;
                g1 := normal(tuple[2]);
                LL := [op(LL), 1];
                G := [op(G), g1];
            end do;
        else

        #### when rank (Gtd/Gd) = 1, let Gtd/ Gd = <tau0 + Gd> and Gd = <tau_1, ..., tau_r>.  ###
            k0 := tau0[1];

            T := TransitionMatrix2(Gtd);
            ## calculate the rank of Gd ##
            if n = 1 or ArrayTools:-IsZero(Array(Gtd[2][1])) then
                k := 0;
            else
                k := nops(Gtd[2]);
            end if;
            printf("rank Gtd = %a, Transition Matrix = %v\n ", k,T);

            ## decide the existence of telescoper for each fraction a/ sigma_t^l (d)^e #############
            ## a/ sigma_t^s (d)^e has a telescoper of type (sigma_t; sigma_{x1}, ..., sigma_{xn}) ##
            ## iff a has a telescoper of type (tau_0; tau_1, ..., tau_r) ###########################
            for fracs in orbit do
                s := fracs[1];
                numers := subsop(1 = NULL, fracs);
                j := 1;
                while j <= nops(numers) do
                    a := numers[j][1];
                    # printf("a = %a\n",a);
                    e := numers[j][2];
                    d_ := ShiftEquivalenceTesting:-Shift(d^e, t, s);
                    a1 := op(RationalSummation:-TransitionMap([a], T, xt));
                    # printf("a1 = %a\n",a1);
                    if k = 0 then
                        printf("case 2: ");
                        ## when rank (Gd) = 0, a/ sigma_t^s (d)^e has a telescoper of type #############
                        # (sigma_t, sigma_{x1}, ..., sigma_{xn}) iff L0 (a) = 0 for some L0 in K(t)<tau0> ##
                        # iff L1 ( phi (a)) = 0 for some L1 in K(t)<St> ################################
                        pair := ExistenceAnnihilator(a1, z, t, St);
                        if pair[1] = false then
                            printf("case 2: %a does not have annihilators w.r.t. %a\n", a1, t);
                            return false;
                        else
                            ## L1 (t, St) (a1) = 0 ##
                            L1 := pair[2];
                            # print(L1);
                            b1 := [0];
                        end if;

                    else
                        printf("case 3: \n");
                        if nargs = 5 then
                            u :=  [op(x[k + 1 .. nops(x)]), op(y)];
                        else
                            u :=  x[k + 1 .. nops(x)];
                        end if;
                        tuple := IsTelescoperable(a1, x[1 .. k], t, 'St', u);
                        if nops([tuple]) = 1  then
                            u := x_[1..k];
                            printf("case 3: %a does not have a telescoper w.r.t. (%a, %a)\n", a1, t, u);
                            return false;
                        end if;

                        ## L1 (t, St) (a1) = Delta_{x1} (b_1) + ... + Delta_{x_k}) (b_k) + c1 ##
                        L1 := tuple[2];
                        b1 := tuple[3];
                    end if;

                        ## L0 (t, tau_0) (a) = Delta_{tau_1} (phi^{-1} (b_1)) + ... + Delta_{tau_k} (phi^{-1} (b_k)) + phi^{-1} (c1) ##
                        L0 := subs([t = t/ k0], L1);
                        #printf("L0=%a\n", L0);
                        ## L (t, St) (a/ d_) = Delta_{x1} (g1) + ... + Delta_{xn} (gn) + L0(t, tau_0)(a)/ d_, ##

                        pair := Transformation2(L0, a, d_, tau0, x, t, St);
                        L := pair[1];
                        g1 := pair[2];
                        LL := [op(LL), L];

                    if not ArrayTools:-IsZero(Array(b1))  then
                        T1 := LinearAlgebra:-MatrixInverse(T);
                        b := RationalSummation:-TransitionMap(b1, T1, xt);
                        #print(b);
                        ## L0(t, tau_0)(a)/ d_ = Delta_{tau_1} (b1/ d_) + ... + Delta_{tau_k} (bk/d_) + c/ d_ ###################
                        ## transformation: Delta_{tau_1} (h1) + ... + Delta_{tau_k} (hk) = Delta_{x1} (g1) + ... + Delta_{xn} (gn) ###
                        h := map( x -> x/ d_, b);
                        g2 := RationalSummation:-Transformation(h, Gd, x);
                    else
                        g2 := convert(LinearAlgebra:-ZeroVector(n), list);
                    end if;
                    g := zip((x, y) -> x + y, g1, g2);
                    G := [op(G), g];

                    j := j + 1;
                end do;
            end do;
        end if;
        i := i + 1;
    end do;

    if nops(LL) = 1 then
        return true, op(LL), op(G);
    else
        result := myLCLM(LL, t, St, G);
        return true, op(result);
    fi;
end proc:



################################################################################################
# Export
# Name: IsTelescoperable2
# Calling sequence:
#                   IsTelescoperable2(f, x, t, 'T', tau, y)
# Input: f,  a multivariate rational function;
#        x,  a list of n variable names;
#        t,  a variable name;
#        T, an operator name;
#        tau,  a list tau = [[a00, ..., a0n], [[a10, ..., a1n], [a20, ..., a2n], ..., [ak0, ..., akn]]]
#              with aij being integers such that a00 is nonzero and a10 = ... = ak0 = 0;
#        y,  (optional) a parameter name or a list of parameter names.
# Output: true, L, g,  if f has a telescoper L in  K(t)<T> such that
#                              L(t, tau_0) (f) = Delta_{tau_1} (g1) + .... + Delta_{tau_k} (gk)
#                      where g = [g1, ..., gk] with gi in K(t, x, y),
#                      for each i, tau_i is the shift operator represented by [ai1, ..., ain].
#         false,  if f does not have any telescoper of type (tau_0; tau_1, ..., tau_k).
# Remark: We require that the a00 is nonzero, a10 = a20 = ... = ak0 = 0, and
#         tau_1, ..., tau_k are linearly independent over Z.
#################################################################################################


IsTelescoperable2 := proc(f, x, t, T, tau, y)
    local A, f1, k, tuple, h, L, B, u, g, a00;
    A := TransitionMatrix2(tau);
    f1 := op(RationalSummation:-TransitionMap([f], A, [t, op(x)]));
    # print(f1);
    k := nops(tau[2]);
    if nargs = 6 then
        u :=  [op(x[k + 1 .. nops(x)]), op(y)];
    else
        u :=  x[k + 1 .. nops(x)];
    end if;
    tuple := IsTelescoperable(f1, x[1 .. k], t, T, u);
    if nops([tuple]) = 1 then
        return false;
    else
        L := tuple[2];
        h := tuple[3];

        #### now L(t, St) phi(f) = Delta_{x1} (h1) + ... + Delta_{xk} (hk) #############################
        ## so L(t/a00, T) f = Delta_{tau_1} (phi^{-1}(h1)) + ... + Delta_{tau_k} (phi^{-1}(hk)) ####
        B := LinearAlgebra:-MatrixInverse(A);
        a00 := tau[1][1];
        L := subs([t = t/a00], L);
        g := RationalSummation:-TransitionMap(h, B, [t, op(x)]);
        return true, L, g;
    end if;
end proc:


##############################################################################################################
# Export
# Name: AdditiveDecomposition2
# Calling sequence:
#                   AdditiveDecomposition2(f, x, t, 'St', y)
# Input:  f,  a multivariate rational function;
#         x,  a variable name or a list of variable names;
#         t,  a variable name;
#         St, an operator name;
#         y,  (optional) a parameter name or a list of parameter names.
# Output: [L, [g1, ..., gn], r],  a list such that
#                                       L (f) = Delta_{x1} (g1) + .... + Delta_{xn} (gn) + r
#                              satisfying
#                                  f has a telescoper of type (St; sigma_{x1}, ..., sigma_{xn})   iff    r = 0.
#                              where L in K(t)<St>, gi in K(t, x, y) and x = [x1, ..., xn].
###############################################################################################################


AdditiveDecomposition2 := proc(f, x_, t, St, y_)
    local n, N, x, xt, x1, tuple, g, e, ListRemainder, Gtp, L, G, R, i, orbit, tau0, Gd, g1, g2, r1, d, f1,\
    k0, T, k, fracs, numers, s, a, a1, d_, b1, c1, L0, y, z, u, T1, LL, r, h, result, Gtd, j, pair, L1, b, c;
    x := convert(x_, list);
    if nargs = 5 then
        y := convert(y_, list);
        z := [op(x), op(y)];
    else
        z := x;
    end if;
    xt := [t, op(x)];
    n := nops(x);
    x1 := x[1];

    tuple := RationalReduction:-OrbitalReduction(f, x1, x, t);
    g := tuple[1];
    ListRemainder := tuple[2]; # list representation of remainder
    Gtp := tuple[3];

#### f = Delta_{x1} (g1) + ... + Delta_{xn} (gn) + r ###
    LL := [1]; # list of Telescopers
    G := [normal(g)]; # list of certificates
    R := [0]; # list of remainders

#### the remainder after orbital reduction is 0 ###
    if ListRemainder = [] then
        return [1, g, 0];
    end if;
    i := 1;
    while i <= nops(ListRemainder) do
        orbit := ListRemainder[i];
        d := orbit[1];
        orbit := orbit[2];
        Gtd := Gtp[i];
        tau0 := Gtd[1];
        Gd := GenerateSubgroup(Gtd);

        #### when rank (Gtd/ Gd) = 0, a/ d^j is has a telescoper ##
        # iff  a/ d^j is (sigma_{x1}, ..., sigma_{xn})-summable ######
        if ArrayTools:-IsZero(Array(tau0)) then
            for fracs in orbit do
                s := fracs[1];
                d_ := ShiftEquivalenceTesting:-Shift(d, t, s);
                numers := subsop(1 = NULL, fracs);
                tuple := RationalSummation:-ListAdditiveDecomposition([[d_,numers]], [Gd], x);
                g1 := normal(tuple[1]);
                r1 := normal (tuple[2]);
                printf("case 1: \n");
                if nops([tuple]) = 1 then
                    printf("case 1: %a is not summable w.r.t. %a\n", r1, x);
                end if;
                LL := [op(LL), 1];
                G := [op(G), g1];
                R := [op(R), r1];
            end do;
        else

        #### when rank (Gtd/Gd) = 1, let Gtd/ Gd = <tau0 + Gd> and Gd = <tau_1, ..., tau_r>.  ###
            k0 := tau0[1];

            T := TransitionMatrix2(Gtd);
            ## calculate the rank of Gd ##
            if n = 1 or ArrayTools:-IsZero(Array(Gtd[2][1])) then
                k := 0;
            else
                k := nops(Gtd[2]);
            end if;
            printf("rank Gd = %a, Transition Matrix = %v\n ", k,T);

            ## decide the existence of telescoper for each fraction a/ sigma_t^l (d)^j #############
            ## a/ sigma_t^s (d)^j has a telescoper of type (sigma_t; sigma_{x1}, ..., sigma_{xn}) ##
            ## iff a has a telescoper of type (tau_0; tau_1, ..., tau_r) ###########################
            for fracs in orbit do
                s := op(fracs[1]);
                numers := subsop(1 = NULL, fracs);
                j := 1;
                while j <= nops(numers) do
                    a := numers[j][1];
                    e := numers[j][2];

                    d_ := ShiftEquivalenceTesting:-Shift(d^e, t, s);
                    a1 := op(RationalSummation:-TransitionMap([a], T, xt));
                    if k = 0 then
                        printf("case 2: ");
                        ## when rank (Gd) = 0, a/ sigma_t^s (d)^e has a telescoper of type #############
                        # (sigma_t, sigma_{x1}, ..., sigma_{xn}) iff L0 (a) = 0 for some L0 in K(t)<tau0> ##
                        # iff L1 ( phi (a)) = 0 for some L1 in K(t)<St> ################################
                        pair := ExistenceAnnihilator(a1, z, t, St);
                        if pair[1] = false then
                            g1 := convert(LinearAlgebra:-ZeroVector(n), list);
                            r1 := a/ d_;
                            printf("case 2: %a does not have annihilators w.r.t. %a\n", a1, t);
                            L1 := 1;
                            b1 := [0];
                            c1 := a1;
                        else
                            ## L1 (t, St) (a1) = 0 ##
                            L1 := pair[2];
                            b1 := [0];
                            c1 := 0;
                        end if;

                    else
                        printf("case 3: \n");
                        if nargs = 5 then
                            u :=  [op(x[k + 1 .. nops(x)]), op(y)];
                        else
                            u :=  x[k + 1 .. nops(x)];
                        end if;
                        tuple := AdditiveDecomposition2(a1, x[1 .. k], t, 'St', u);
                        if normal(tuple[3]) <> 0 then
                            u := x_[1..k];
                            printf("case 3: %a does not have a telescoper w.r.t. (%a, %a)\n", a1, t, u);
                        end if;

                        ## L1 (t, St) (a1) = Delta_{x1} (b_1) + ... + Delta_{x_k}) (b_k) + c1 ##
                        L1 := tuple[1];
                        b1 := tuple[2];
                        c1 := tuple[3];
                    end if;

                        ## L0 (t, tau_0) (a) = Delta_{tau_1} (phi^{-1} (b_1)) + ... + Delta_{tau_k} (phi^{-1} (b_k)) + phi^{-1} (c1) ####
                        L0 := subs([t = t/ k0], L1);

                        ## L (t, St) (a/ d_) = Delta_{x1} (g1) + ... + Delta_{xn} (gn) + L0(t, tau_0)(a)/ d_, ##
                        pair := Transformation2(L0, a, d_, tau0, x, t, St);
                        L := pair[1];
                        g1 := pair[2];
                        LL := [op(LL), L];

                    if c1 <> 0 or not ArrayTools:-IsZero(Array(b1))  then
                        T1 := LinearAlgebra:-MatrixInverse(T);
                        pair := RationalSummation:-TransitionMap([op(b1),c1], T1, xt);
                        b := subsop(nops(pair) = NULL, pair);
                        c := pair[nops(pair)];
                        r := c/d_;

                        ## L0(t, tau_0)(a)/ d_ = Delta_{tau_1} (b1/ d_) + ... + Delta_{tau_k} (bk/d_) + c/ d_ ###################
                        ## transformation: Delta_{tau_1} (h1) + ... + Delta_{tau_k} (hk) = Delta_{x1} (g1) + ... + Delta_{xn} (gn) ###
                        h := map( x -> x/ d_, b);
                        g2 := RationalSummation:-Transformation(h, Gd, x);
                    else
                        g2 := convert(LinearAlgebra:-ZeroVector(n), list);
                        r := 0;
                    end if;
                    g := zip((x, y) -> x + y, g1, g2);
                    G := [op(G), g];
                    R := [op(R), r];

                    j := j + 1;
                end do;
            end do;
        end if;
        i := i + 1;
    end do;
    if nops(LL) = 1 then
        [op(LL), op(G), op(R)];
    else
        result := myLCLM(LL, t, St, G, R);
        result;
    end if;
end proc:


##############################################################################################################
# Export
# Name: ShiftApply
# Calling sequence:
#                    ShiftApply(L, f, t, 'St')
# Input:  L,  a polynomial in St;
#         f,  a multivariate rational functions;
#         t,  a variable name;
#         St,  an operator name representing the shift operator w.r.t. t.
# Output: g, where g = L(t, St) (f).
###############################################################################################################


ShiftApply := proc(L_, f, t, St)
    local A, L, f1;
    A := OreTools:-SetOreRing(t, 'shift');
    AddConversionRule('shift', difference_case);
    L := OreTools[Converters]:-FromPolyToOrePoly(L_, St);
    f1 := normal(expand(OreTools:-Apply(L, f, A)));
end proc:


##############################################################################################################
# Export
# Name: VerifyDecomposition2
# Calling sequence:
#                    VerifyDecomposition2(L, f, g, r, x, t, 'St')
# Input:  L,  a polynomial in St;
#         f, r,  two multivariate rational functions;
#         g,  a list of multivariate rational functions;
#         x,  a variable name or a list of variable names;
#         t,  a variable name;
#         St,  an operator name representing the shift operator w.r.t. t.
# Output: true,  if L (f) = Delta_{x1} (g1) + ... + Delta_{xn} (gn) + r.
#         false,  otherwise.
###############################################################################################################


VerifyDecomposition2 := proc(L_, f, g, r, x, t, St)
    local f1;
    f1 := ShiftApply(L_, f, t, St);
    RationalReduction:-VerifyDecomposition(f1, g, r, x);
end proc:



##############################################################################################################
# Export
# Name: VerifyTelescoper
# Calling sequence:
#                    VerifyTelescoper(L, f, g, x, t, 'St')
# Input:  L,  a polynomial in St;
#         f,  a multivariate rational functions;
#         g,  a list of multivariate rational functions;
#         x,  a variable name or a list of variable names;
#         t,  a variable name;
#         St,  an operator name representing the shift operator w.r.t. t.
# Output: true,  if L (f) = Delta_{x1} (g1) + ... + Delta_{xn} (gn).
#         false,  otherwise.
###############################################################################################################


VerifyTelescoper := proc(L, f, g, x, t, St)
    local flag;
    flag := VerifyDecomposition2(L, f, g, 0, x, t, St);
    flag;
end proc:


################################################################################################
# Export
# Name: VerifyTelescoper2
# Calling sequence:
#                   VerifyTelescoper2(L, f, g, x, t, 'T', tau)
# Input: L,  a polynomial in St;
#        f,  a multivariate rational functions;
#        g,  a list of multivariate rational functions;
#        x,  a list of variable names;
#        u,  a parameter name or a list of parameter names;
#        tau,  a list tau = [[a11, ..., a1n], [[a21, ..., a2n], ..., [ak1, ..., akn]]]
#              with aij being integers.
# Output: true  if f = Delta_tau1 ( g1 ) + ... + Delta_tau_k ( gk )
#                    where x = [x1, ..., xn], g = [g1, ..., gk],
#                    for each i, tau_i is the shift operator represented by [ai1, ..., ain].
#         false, otherwise.
#################################################################################################


VerifyTelescoper2 := proc(L, f, g, x, t, T, tau)
    local row, col, val, c, f1, h, i, flag, u, k;
    c := PolynomialTools:-CoefficientVector(L, T);
    row, col, val := ArrayTools:-SearchArray(convert(c, Array));
    f1 := add(val[i] * ShiftEquivalenceTesting:-Shift(f, [t, op(x)], (row[i] - 1) * tau[1]), i = 1 .. ArrayTools:-Size(val, 1));
    h := 0;
    i := 1;
    k := nops(tau[2]);
    for i to k do
        h := h + ShiftEquivalenceTesting:-Shift(g[i], [t, op(x)], tau[2][i]) - g[i];
    end do;
    flag := Testzero(f1 - h);
end proc:


##########################################################################################
# Export
# Name: RandomRational
# Calling sequence:
#                    RandomRational(x)
# Input: x,  a variable name or a list of variable names.
# Output: f,  a multivariate rational function in x whose denominator is a product
#             linear polynomials.
##########################################################################################


RandomRational := proc(x)
    local n, f, i, p, q, j, g;
    n := nops(x);
    f := 0;
    for i to 1 do
        p := randpoly(x, coeffs = rand(-5 .. 5), degree = 1);
        q := 1;
        for j to n do
            q := q * randpoly(x[1..j], coeffs = rand(-2 .. 2), degree = 1);
            if q = 0 then
                q := 1;
            end if;
        end do;
        g := p/ q;
        f := f + g;
    end do;
    f;
end proc:


##########################################################################################
# Export
# Name: RandomRational2
# Calling sequence:
#                    RandomRational2(x)
# Input: x,  a variable name or a list of variable names.
# Output: f,  a multivariate rational function in x whose denominator is a product
#             some polynomials of degree 1 or 2.
##########################################################################################


RandomRational2 := proc(x)
    local n, f, i, p, q, j, g;
    n := nops(x);
    f := 0;
    for i to 1 do
        p := randpoly(x, coeffs = rand(-5 .. 5), degree = 1);
        q := 1;
        for j to n do
            q := q * randpoly(x[j..n], coeffs = rand(-2 .. 2), degree = RandomTools[LinearCongruence]:-GenerateInteger(range =1 .. 2));
            if q = 0 then
                q := 1;
            end if;
        end do;
        g := p/ q;
        f := f + g;
    end do;
    f;
end proc:


##############################################################################################################
# Local
# Name: GenerateSubgroup
# Calling sequence:
#                    GenerateSubgroup(Gpt)
# Input:  Gpt,  a list [tau_0, [tau_1, ..., tau_r]] representing a basis of the isotropy group Gpt
#               w.r.t. variables t and x.
# Output: Gp,  a list  [sigma_1, ..., sigma_r] representing a basis of the isotropy group Gp w.r.t. x,
#              this means [0, op(sigma_1)] = tau_1.
###############################################################################################################


GenerateSubgroup := proc(Gpt)
    local Gp, F;
    Gp := Gpt[2];
    F := tau -> subsop(1 = NULL, tau);
    Gp := map(F, Gp);
end proc:


####################################################################################################
# Local
# Name: TransitionMatrix2
# Calling sequence:
#                   TransitionMatrix2(Gtp)
# Input: Gtp,  a list [tau_0, [tau_1, ..., tau_r] representing a 'normalized' basis of Gtp with
#              rank (Gtp/Gp) = 1, i.e, tau_0[1] <> 0, and tau_i[1] = 0, for i = 1, ..., r;
# Output: A,  an invertible matrix A satisfying the following property:
#         for any rational function f(x) in K(x) with x = [x1, ..., xn],
#
#               f(x) has a telescoper of type (tau_0; tau_1, ..., tau_r)
#         iff
#               f(xA) has a telescoper of type (sigma_t; sigma_{x1}, ..., sigma_{xr}).
######################################################################################################


TransitionMatrix2 := proc(Gtp)
    local A, tau0, N, Gp, B, n, v, A1;
    tau0 := convert(Gtp[1], Matrix);
    if ArrayTools:-IsZero(Array(Gtp[2][1])) then
        N := LinearAlgebra:-IdentityMatrix(nops(Gtp[1]));
        N := LinearAlgebra:-DeleteRow(N, 1);
        A := Matrix([[tau0], [N]]);
    else
        Gp := GenerateSubgroup(Gtp);
        B := RationalSummation:-TransitionMatrix(Gp);
        n := nops(Gtp[1])-1;
        v := Matrix(1 .. n, 1, 0);
        A1 := <v|B>;
        A := <tau0,A1>;
    end if;
end proc:


######################################################################################################
# Local
# Name: Transformation2
# Calling sequence:
#                    Transformation2(L0, a, b, tau0, x, t, 'St')
# Input: L0,  a polynomial in St;
#        a,  a multivariate rational function;
#        b,  a multivariate rational function;
#        tau0,  a list [k0, k1, ...] representing a shift operator w.r.t. t and x;
#        x,  a variable name or a list of variable names;
#        t,  a variable name.
# Output: [L, g],  a list such that L is a polynomial in St:
#
#                             L = sum_l e_l * St^{lk0},
#
##                 if subs(St = Tau, L0) = sum_i e_l Tau^l := L0 (t, Tau), and g = [g1, ..., gn] is a list of
#                  multivariate rational functions satisfying the property:
#
#                          L(t, St) (a/ b) = Delta_{x1} (g1) + .... + Delta_{xn} (gn) + L0(t, Tau) (a) / b
#
##                 if tau0 (b) = b and x = [x1, ..., xn].
# Remark: In particular, if tau0(b) = b, L0 (a) = 0, then this procedure returns a telescoper L for a/b.
#######################################################################################################


Transformation2 := proc(L0, a, b, tau0, x, t, St)
    local n, g, k0, c, v, s, i, j, g1, L, aj, tau, pair;
    n := nops(x);
    g := convert(LinearAlgebra:-ZeroVector(n), list);
    k0 := tau0[1];
    c := [coeffs(L0, St, 'v')];
    v := [v];
    s := map(degree, v, St);
    i := 1;
    while i <= nops(c) do
        j := s[i];
        aj := c[i]*ShiftEquivalenceTesting:-Shift(a, t, j * k0);
        tau := subsop(1 = NULL, tau0);
        tau := -j * tau;
        pair := RationalReduction:-AbramovReduction(aj, b, x, tau);
        g1 := pair[1];
        g := zip((x, y) -> x + y, g, g1);
        i := i + 1;
    end do;

#### L = L (t, St^k0)
    L := subs(St = St^k0, L0);
    [L, g];
end proc:


#################################################################################################################
# Local
# Name: myLCLM
# Calling sequence:
#                   myLCLM(P, t, St, H, R)
# Input:  P,  a list [P1, ..., Pk] of nonzero polynomials representing a list of shift operators in K(t)<St>;
#         t,  a variable name;
#         St,  a shift operator name;
#         H,  (optional) a list [h1, ..., hk] of certificates, where hj = [hj1, ..., hjn] are certificates for Pi;
#         R,  (optional) a list [r1, ..., rk] of multivariate rational functions.
# Output: [L, g, r],  a list such that L is the LCLM of P1, ..., Pk, and g = [g1, ..., gn] are the certificates
#                   for L, r is a remainder, that means
#                                 L = R1 * P1 = ... = Rk * Pk,
#                   for some Ri in K(t)<St>, gj = R1 (h1j) + ... + Rk (hkj) for j = 1, ..., n, and
#                                 r = R1 (r1) + ... + Rk (rk).
# Remark: 1. If there is no give H, then return g as [].
#         2. Let x = [x1, ..., xn]. If Pi (fi) = Delta_{x1} (hi1) + ... + Delta_{xn} (hin) + ri for i = 1, ..., k, then
#
#                   L (f1 + ... + fk) = R1 * P1 (f1) + ... + Rk * Pk (fk) + L(r1 + ... + rk)
#                                     = R1 * (Delta_{x1} (h11) + ... + Delta_xn (h1n)) + ... + Rk * (Delta_{x1} (hk1) + ... + Delta_{xn} (hkn)) + r
#                                     = Delta_{x1} (g1) + ... + Delta_{xn} (gn) + r.
###################################################################################################################


myLCLM := proc(P, t, St, H, R )
    local A, f, Q, F, L, Rfactor, gj, j, h, g, r;
    A := OreTools:-SetOreRing(t, 'shift');

#### convert polynomials to Ore polynomials
    OreTools[Converters]:-AddConversionRule('shift', difference_case);
    f := p -> OreTools[Converters]:-FromPolyToOrePoly(p, St);
    Q := map(f, P);

    F := OreTools:-LCM(op(Q), A, 'cofactors = true');
    L := F[1];

#### Apply cofactors on the certificates H ####
    Rfactor := F[2];
    g := [];


    if nargs > 3  then
        for j from 1 to nops(H[1]) do
            gj := add(OreTools:-Apply(Rfactor[i], H[i][j], A), i = 1 .. nops(Rfactor));
            g := [op(g), gj];
        end do;
        if nargs = 5 then
            r := normal(add(OreTools:-Apply(Rfactor[i], R[i], A), i = 1 .. nops(R)));
        end if;
    end if;

    L := OreTools[Converters]:-FromOrePolyToPoly(L, St);

    if nargs = 3 then
        return [L];
    elif nargs = 4 then
        return [L, g];
    else
        return [L, g, r];
    end if;
end proc:


##########################################################################
# Local
# This code is from OreTools[Converters]_FromOrePolyToPoly.
# It is used to define the rule for an Ore algebra in the difference case.
##########################################################################


difference_case := proc(eq, func, A)
    local reqn, info, fcn, E, receq, req, func_set, n, place, i, var;
    var := OreTools[Properties][GetVariable](A);

    # Do some argument checking and processing
    reqn := LREtools['REcreate'](eq, func(var), {});
    info := eval(op(4, reqn));

    if info['linear']=false then
        error "only linear equations are handled"
    end if;
    n := info['vars'];
    if not (nops([n])=1 and type(n, name)) then
        error "only single equations are handled"
    end if;
    fcn := op(2, reqn)[1];
    # locate n in arguments to fcn
    for place to nops(fcn) while op(place, fcn) <> n do end do;

    # Transform the equation from a difference equation to a
    # difference operator
    req := args[1];
    func_set := indets(req, specfunc(anything, info['functions']));
    receq := collect( subs( map( (x, y, z, i) -> x = (z+1)^(op(i, x)-y),
        func_set, n, 'E', place), req), 'E', Normalizer);
    'OrePoly'(seq(coeff(receq, 'E', i), i=0..degree(receq, 'E')))
end proc:


end module:


#################################################################
## Copyright (c) August 2021, UCAS and JKU. All rights reserved.
## Author: Lixin Du <dulixin17@mails.ucas.ac.cn>
## Load the modules:
##         ShiftOrbitDecomposition.mm,
##         OrbitalDecomposition.mm,
##         RationalReduction.mm.
## Modified on December 15, 2021
#################################################################


RationalSummation := module()
      description "Decide the Summability of Multivariate Rational Functions";
      export
                TransitionMatrix,
                TransitionMap,
                Transformation,
                IsSummable,
                IsSummable2,
                AdditiveDecomposition,
                VerifySummation,
                VerifySummation2,
                ListSummable,
                ListAdditiveDecomposition,
                RandomSummableRational;
      local
                VerifyTransformation;
      option
                package;


################################ Notation  Description #################################################
# 0. Setting: the field K is of characteristic zero. In this algorithm, K can be taken
# as a field of multivariate rational functions Q(u), where Q is the rational number
# field and u = [u1, ..., um] are parameters.
#
# 1. Shift operators and difference operators:
# Let f be a multivariate rational function in K(t, x) with x = [x1, ..., xn]. For each variable v
# in {t,x1, ..., xn}, the shift operator sigma_v is defined as an F-automorphism of K(t, x) such that
#
#        sigma_v (v) = v + 1,     and     sigma_v (w) = w,  for all w in {t,x1, ..., xn}, and w <> v
#
# The difference operator is defined as
#
#                                Delta_v := sigma_v - 1,
#
# where 1 is the identity map of K(t, x).
# Remark: In the following algorithms, we shall use the list representation of shifts,
#         (see Data Description in OrbitalDecomposition module).
#
# 2. Group actions under shifts
# Let G = <sigma_{x1}, ..., sigma_{xn}> be the free abelian group generated by shift operators
# sigma_{x1}, ..., sigma_{xn}. Let Gp be the isotropy group of a polynomial p in G.
########################################################################################################


################################ Task  Description ######################################################
# 1. Summability Testing:
# 1.1 Given a multivariate rational function f in K(x) with x = [x1, ..., xn], decide whether
# there exist g1,....,gn in K(x) such that
#
#             f = Delta_{x1} (g1) + ... + Delta_{xn} (gn).
#
# If such gi exist, we say f is (sigma_x1,...,sigma_xn)-summable.
#
# 1.2  Let tau_1, ..., tau_k \in G be linearly independent over Z.  Given a multivariate rational
# function f in K(x) with x = [x1, ..., xn], decide whether there exist g1,....,gn in K(x) such that
#
#             f = Delta_{tau_1} (g1) + ... + Delta_{tau_k} (gk).
#
# If such gi exist, we say f is (tau_1, ..., tau_k)-summable.
#
# 2. Additive Decomposition for the summability problem
# Given a multivariate rational function f in K(x) with x = [x1, ..., xn], compute rational functions
# g1, ..., gn and r in K(x) such that
#
#             f = Delta_{x1} (g1) + ... + Delta_{xn} (gn) + r
#
# satisfying the property
#
#            f is (sigma_x1, ..., sigma_xn)-summable    iff   r = 0.
########################################################################################################


#########################################################################################
# Export
# Name: TransitionMatrix
# Calling sequence:
#                   TransitionMatrix(k)
# Input: k,  a list [[a11, ..., a1n], [a21, ..., a2n], ..., [ar1, ..., arn]]
#            with aij being integers or rational numbers.
# Output: A,  an invertible matrix A = (aij)_{1 <= i, j <= n}.
# Remark: the matrix A in the output satisfies the following property:
#         for any rational function f(x) in K(x),
#
#               f(x) is (tau_1, ..., tau_r)-summable
#         iff
#               f(xA) is (sigma_x1, ..., sigma_xr)-summable
#
#         where x = [x1, ..., xn] and tau_i = (sigma_x1)^(ai1) * ... * (sigma_xn)^(ain)
#         is the shift operator represented by [ai1, ..., ain] for i = 1, ..., r.
#########################################################################################


TransitionMatrix := proc(k)
    local A, G, B, place, i, j, n, number_set, b;

#### perform Gaussian elimination ####
    A := convert(k, Matrix);
    G := LinearAlgebra:-GaussianElimination(A);
    B := convert(G, Array);

#### find the place of the first nonzero element in each row of the matrix ####
    place := [];
    for i to nops(k) do
        j := ArrayTools:-SearchArray(B[i], 1, location = first);
        place := [op(place), j[1]];
    end do;
    n := nops(k[1]);
    number_set := {seq(i, i = 1 .. n)};
    place := number_set minus convert(place, set);

#### put 1 on appropriate places ####
    for i in place do
        b := LinearAlgebra:-Transpose(convert(Vector(n, shape = unit[i]), Matrix));
        A := Matrix([[A], [b]]);
    end do;
    A;
end proc:


########################################################################################
# Export
# Name: TransitionMap
# Calling sequence:
#                   TransitionMap(f, A, x)
# Input: f,  a list of rational functions;
#        A,  a matrix with entries being rational numbers;
#        x,  a variable name or a list of variable names.
# Output: g,  a list [g1, ..., gr] of rational functions such that
#                             gi(x) = fi(xA)
#             for i = 1, ..., r, where f = [f1, ..., fr].
# Remark: the input f must be a list. The reason is for a single rational function, such
#         as f = x1 + x2, the command convert(f, list) gives [x1,x2], which cannot be
#         recognized as [f].
#########################################################################################


TransitionMap := proc(f::list, A, x_)
    local F, G, B, vx, vy, y, h1, h, g, x;
    G := [];
    F := convert(f, list);
    x := convert(x_, list);
    for h in F do
        vx := convert(x, Matrix);
        vy := LinearAlgebra:-Multiply(vx, A);
        y := convert(vy, list);
        y := zip((x, y) -> x = y, x, y);
        h1 := subs(y, h);
        G := [op(G), h1];
    end do;
    G;
end proc:


######################################################################################################
# Export
# Name: Transformation
# Calling sequence:
#                    Transformation(h, tau, x)
# Input: h,  a list of multivariate functions;
#        tau,  a shift operator or a list of shift operators;
#        x,  a variable name or a list of variable names.
# Output: g,  a list of multivariate rational functions such that
#
#                  Delta_{tau_1} (h1) + ... + Delta_{tau_r} (hr)
#                = Delta_{x1} (g1) + .... + Delta_{xn} (gn),
#
#             where  tau = [tau_1, ..., tau_r], x = [x1, ..., xn],
#             h = [h1, ..., hr] and g = [g1, ..., gn].
# REMAEK: THERE IS AN ALTERNATIVE WAY TO ACHIEVE Transformation BY THE FOLLOWING FORMULA
#            tau1 - 1 = (sigma_x1 - 1) * (...) + ... + (sigma_xn - 1) * (...).
#######################################################################################################


Transformation := proc(h, tau_, x_)
    local r, n, g, i, f, L, g1, x, tau;
    r := nops(h);
    x := convert(x_, list);
    tau := convert(tau_, list);
    n := nops(x);
    g := convert(LinearAlgebra:-ZeroVector(n), list);
    i := 'i';
    for i to r do
       f := ShiftEquivalenceTesting:-Shift(h[i], x, tau[i]);

       # tau1 (f) - f = Delta_{x1} (...) + ... + Delta_{xn} (...) + r with r = 0 #
       L := RationalReduction:-AbramovReduction(f, 1, x, tau[i]);
       g1 := L[1];
       g := zip((x, y) -> x + y, g, g1);
    end do;
    g;
end proc:


####################################################################################################
# Export
# Name: ListSummable
# Calling sequence:
#                   ListSummable(L, Gp, x)
# Input: L,  a list [[d1, [a11, a12, ...]], [d2, [a21, a22, ...]], ...] representing a remainder
#            in orbital reduction for the summability problem;
#        Gp,  a list [Gd1, Gd2, ...] consisting of the isotropy group of d1, d2, ...;
#        x,  a variable name or a list of variable names.
# Output: true, g  if f = Delta_x1 (g1) + ... + Delta_xn (gn) is (sigma_x1, ..., sigma_xn)-summable,
#                  where x = [x1, ..., xn].
#         false,  if f is not (sigma_x1, ..., sigma_xn)-summable.
####################################################################################################


ListSummable := proc(L, Gp, x)
    local h, i, n, flag, orbit, d, e, r1, numers, flag1, Gd, Gd1, pair, j, a, T, k, a1, L1, Gp1, T1, b1, b2, g, b, g1;

    n := nops(x);
    flag := true;
    g := convert(LinearAlgebra:-ZeroVector(n), list);
    if L = [] then
        return true, g;
    end if;

#### run over all the orbits in L ###
    i := 1;
    while i <= nops(L) do
        orbit := L[i];
        d := orbit[1];
        numers := orbit[2];
        Gd := Gp[i];

        if n = 1 or ArrayTools:-IsZero(Array(Gd[1]))  then

        #### when Gd = [[0,0,..]] is a trivial group, a/ d^j is (sigma_{x1}, ..., sigma_{xn})-summable iff a = 0 ####
            r1 := RationalReduction:-FromListToRemainder([orbit], x);
            printf("%a is not summable w.r.t. %a\n", r1, x[1 .. n]);
            return false;
        else

        ##### when Gd = [tau_1, ..., tau_k] is nontrivial, a/ d^j is (sigma_{x1}, ..., sigma_{xn})-summable
        # iff a is (tau_1, ..., tau_k)-summable iff phi(a) is (sigma_{x1}, ..., sigma_{x_k})-summable for some
        # difference K-automorphism phi determined by a matrix T ###
            j := 1;
            while j <= nops(numers) do
                a := numers[j][1];
                e := numers[j][2];
                T := TransitionMatrix(Gd);
                k := nops(Gd);
                printf("Transition Matrix = %v \n ", T);
                a1 := op(TransitionMap([a], T, x));
                # print(a1);
                b1, L1, Gp1 := op(RationalReduction:-OrbitalReduction(a1, x[1], x[1 .. k]));

                if L1 <> [] then
                    pair := ListSummable(L1, Gp1, x[1 .. k]);
                    if nops([pair]) = 1 then
                       return false;
                    else
                       b2 := pair[2];
                       b1 := zip((x, y) -> x + y, b1, b2);
                   end if;
                end if;

                if not ArrayTools:-IsZero(Array(b1))  then
                    T1 := LinearAlgebra:-MatrixInverse(T);

                    #### now a = Delta_{tau_1} (phi^{-1} (b_1)) + ... + Delta_{tau_k} (phi^{-1} (b_k)) + phi^{-1} (c1) ####
                    # set b =[phi^{-1} (b_1), ..., phi^{-1} (b_k)] ########################################
                    b := TransitionMap(b1, T1, x);
                    # print(b);
                    #### transformation: Delta_{tau_1} (h1) + ... + Delta_{tau_k} (hk) = Delta_x1 (g1) + ... + Delta_xn (gn) ###
                    h := map( x -> x/ d^e, b);
                    g1 := Transformation(h, Gd, x);
                    g := zip((x, y) -> x + y, g, g1);
                end if;

                j := j + 1;
            end do;
        end if;
        i := i + 1;
    end do;
    flag, g;
end proc:


####################################################################################################
# Export
# Name: ListAdditiveDecomposition
# Calling sequence:
#                   ListAdditiveDecomposition(L, Gp, x)
# Input: L,  a list [[d1, [a11, a12, ...]], [d2, [a21, a22, ...]], ...] representing a remainder
#            in orbital reduction for the summability problem;
#        Gp,  a list [Gd1, Gd2, ...] consisting of the isotropy group of d1, d2, ...;
#        x,  a variable name or a list of variable name.
# Output: [[g1, ..., gn], r],  a list such that
#                      f = Delta_x1 (g1) + ... + Delta_xn (gn) + r,
#             satisfying
#                    f is (sigma_x1, ..., sigma_xn)-summable   iff   r = 0,
#             where x = [x1, ..., xn].
####################################################################################################


ListAdditiveDecomposition := proc(L, Gp, x, g_)
    local h, i, n, g, r,  orbit, d, e, numers, Gd, r1, Gd1, pair, g1, j, a, T, k, a1, b1, L1, Gp1, b2, c1, T1, b, c;

    n := nops(x);
    r := 0;
    g := convert(LinearAlgebra:-ZeroVector(n), list);

    if L = [] then
        return [g, 0];
    end if;

#### run over all the orbits in L ###
    i := 1;
    while i<= nops(L) do
        orbit := L[i];
        d := orbit[1];
        numers := orbit[2];
        Gd := Gp[i];

        if n = 1 or ArrayTools:-IsZero(Array(Gd[1]))  then

        #### when Gd = [[0,0,..]] is a trivial group, a/ d^e is (sigma_{x1}, ..., sigma_{xn})-summable iff a = 0 ####
            r1 := RationalReduction:-FromListToRemainder([orbit], x);
            r := r + r1;
        else

        ##### when Gd = [tau_1, ..., tau_k] is nontrivial, a/ d^e is (sigma_{x1}, ..., sigma_{xn})-summable
        # iff a is (tau_1, ..., tau_k)-summable iff phi(a) is (sigma_{x1}, ..., sigma_{x_k})-summable for some
        # difference K-automorphism phi determined by a matrix T ###
            j := 1;
            while j <= nops(numers) do
                a := numers[j][1];
                e := numers[j][2];

                T := TransitionMatrix(Gd);

                k := nops(Gd);
                printf("Transition Matrix = %v \n ", T);
                #### compute b1 = [b_1, ..., b_k] and c such that ####################
                ## a1(x) = phi(a) = Delta_{x1} (b_1), ..., Delta_{x_k}) (b_k) + c1 ###
                a1 := op(TransitionMap([a], T, x));
                # print(a1);
                # tt := time();
                # print(begin_AddDecomp);

                b1, L1, Gp1 := op(RationalReduction:-OrbitalReduction(a1, x[1], x[1 .. k]));

                if L1 <> [] then
                    b2, c1 := op(ListAdditiveDecomposition(L1, Gp1, x[1 .. k]));
                    b1 := zip((x, y) -> x + y, b1, b2);
                else
                    c1 := 0;
                end if;

                #### now a = Delta_{tau_1} (phi^{-1} (b_1)) + ... + Delta_{tau_k} (phi^{-1} (b_k)) + phi^{-1} (c1) ####
                if c1 <> 0 or not ArrayTools:-IsZero(Array(b1)) then
                    T1 := LinearAlgebra:-MatrixInverse(T);
                    if c1 <> 0 then
                        c := op(TransitionMap([c1], T1, x));
                        r := r + c/ d^e;
                    end if;
                    if not ArrayTools:-IsZero(Array(b1)) then
                        b := TransitionMap(b1, T1, x);
                        # print(b);
                        #### transformation: Delta_{tau_1} (h1) + ... + Delta_{tau_k} (hk) = Delta_{x1} (g1) + ... + Delta_{xn} (gn) ###
                        h := map( x -> x/ d^e, b);
                        g1 := Transformation(h, Gd, x);
                        g := zip((x, y) -> x + y, g, g1);
                    end if;
                end if;
                j := j + 1;
            end do;
        end if;
        i := i + 1;
    end do;
    [g, r];
end proc:


################################################################################################
# Export
# Name: IsSummable
# Calling sequence:
#                     IsSummable(f, x)
# Input: f,  a multivariate rational function;
#        x,  a variable name or a list of variable names.
# Output: true, g,  if f = Delta_{x1} (g1) + ... + Delta_{xn} (gn)  is (sigma_{x1}, ..., sigma_{xn})-summable,
#                  where x = [x1, ..., xn], g = [g1, ..., gn] and gi are rational funcitons.
#         false,  if is not (sigma_x1, ..., sigma_xn)-summable.
#################################################################################################


IsSummable := proc(f, x_)
    local x, x1, R, L, Gp, g, g1, flag;
    x := convert(x_, list);
    x1 := x[1];

    R := RationalReduction:-OrbitalReduction(f, x1, x);
    g := R[1];
    L := R[2];
    Gp := R[3];

    R := ListSummable(L, Gp, x);

    if nops([R]) = 1 then
        return false;
    else
        flag := R[1];
        g1 := R[2];
    end if;
    g := zip((x,y) -> x + y, g, g1);

    flag, g;
end proc:


################################################################################################
# Export
# Name: VerifySummation
# Calling sequence:
#                VerifySummation(f, g, x)
# Input: f,  a multivariate rational function;
#        g,  a multivariate rational function or a list of multivariate rational functions;
#        x,  a variable name or a list of variable names.
# Output: true,  if f = Delta_{x1} (g1) + ... + Delta_{xn} (gn),
#                    where g = [g1, ..., gn] and x = [x1, ..., xn].
#         false, otherwise.
#################################################################################################


VerifySummation := proc(f, g, x)
    local flag;
    flag := RationalReduction:-VerifyDecomposition(f, g, 0, x);
    flag;
end proc:



#######################################################################################################################
# Export
# Name: IsSummable2
# Calling sequence:
#                   IsSummable2(f, x, tau)
# Input: f,  a multivariate rational function;
#        x,  a list of n variable names;
#        tau,  a list tau = [[a11, ..., a1n], [a21, ..., a2n], ..., [ak1, ..., akn]]
#              with aij being integers.
# Output: true, [g1, ..., gk],  if  f = Delta_{tau_1} (g1) + .... + Delta_{tau_k} (gk) is (tau_1, ..., tau_k)-summable
#                               where x = [x1, ..., xn],
#                               for each i, tau_i is the shift operator represented by [ai1, ..., ain].
#         false, if f is not (tau_1, ..., tau_k)-summable.
# Remark: We require that tau_1, ..., tau_k are linearly independent over Z.
#######################################################################################################################


IsSummable2 := proc(f, x, tau)
    local T, f1, k, tuple, h, T1, g;
    T := TransitionMatrix(tau);
    f1 := op(TransitionMap([f], T, x));
    k := nops(tau);
    tuple := IsSummable(f1, x[1 .. k]);
    if nops([tuple]) = 1 then
        return false;
    else
        h := tuple[2];

        #### now phi(f) = Delta_{x1} (h1) + ... + Delta_{xk} (hk) ##########
        ## so f = Delta_{tau_1} (phi^{-1}(h1)) + ... + Delta_{tau_k} (phi^{-1}(hk)) ####
        T1 := LinearAlgebra:-MatrixInverse(T);
        g := TransitionMap(h, T1, x);
        return true, g;
    end if;
end proc:


################################################################################################
# Export
# Name: VerifySummation2
# Calling sequence:
#                   VerifySummation2(f, g, x, tau)
# Input: f,  a multivariate rational functions;
#        g,  a list of multivariate rational functions;
#        x,  a list of variable names;
#        tau,  a list tau = [[a11, ..., a1n], [a21, ..., a2n], ..., [ak1, ..., akn]]
#              with aij being integers.
# Output: true,  if f = Delta_{tau_1} ( g1 ) + ... + Delta_{tau_k} ( gk )
#                where x = [x1, ..., xn], g = [g1, ..., gk],
#                for each i, tau_i is the shift operator represented by [ai1, ..., ain].
#         false, otherwise.
#################################################################################################


VerifySummation2 := proc(f, g, x, tau)
    local h, i, flag, k;
    h := 0;
    i := 1;
    k := nops(tau);
    for i to k do
        h := h + ShiftEquivalenceTesting:-Shift(g[i], x, tau[i]) - g[i];
    end do;
    flag := Testzero(f-h);
end proc:



################################################################################################
# Export
# Name: AdditiveDecomposition
# Calling sequence:
#                   AdditiveDecomposition(f, x)
# Input: f,  a multivariate rational function;
#        x,  a variable name or a list of variable names.
# Output: [[g1, ..., gn], r],  a list such that
#                      f = Delta_{x1} (g1) + ... + Delta_{xn} (gn) + r,
#             satisfying
#                    f is (sigma_{x1}, ..., sigma_{xn})-summable   iff   r = 0,
#             where x = [x1, ..., xn].
#################################################################################################


AdditiveDecomposition := proc(f, x_)
    local x, x1, R, g, L, Gp, flag, g1, r;
    x := convert(x_, list);
    x1 := x[1];

    R := RationalReduction:-OrbitalReduction(f, x1, x);
    g := R[1];
    L := R[2];
    Gp := R[3];

    R := ListAdditiveDecomposition(L, Gp, x);
    g1 := R[1];
    r := R[2];

    g := zip((x,y) -> x + y, g, g1);

    [g, r];
end proc:


##########################################################################################
# Export
# Name: RandomSummableRational
# Calling sequence:
#                    RandomSummableRational(x)
# Input: x,  a variable name or a list of variable names.
# Output: f,  a (sigma_x1, ..., sigma_xn)-summable rational function if x = [x1, ..., xn].
##########################################################################################


RandomSummableRational := proc(x)
    local n, f, i, j, p, q, g;
    n := nops(x);
    f := 0;
    for i to n do
        p := randpoly(x, coeffs = rand(-5 .. 5), degree = 2);
        q := 1;
        for j to n do

            q := q * randpoly(x[1..j], coeffs = rand(1 .. 5), degree = 1)^rand(1..3)();
        end do;
        g := p/q;
        f := f + ShiftEquivalenceTesting:-Shift(g, x[i], 1) - g;
    end do;
    f;
end proc:


##############################################################################################################
# Local
# Name: VerifyTransformation
# Calling sequence:
#                    VerifyTransformation(h, g, tau, x)
# Input: h,  a list of multivariate functions;
#        g,  a list of multivariate functions;
#        tau,  a list of shift operators;
#        x,  a variable name or a list of variable names.
# Output: true,  if Delta_{tau_1} (h1) + ... + Delta_{tau_r} (hr) = Delta_{x1} (g1) + .... + Delta_{xn} (gn),
#                 where x = [x1, ..., xn], tau = [tau_1, ..., tau_r], h = [h1, ..., hr] and g = [g1, ..., gn].
#         false,  for otherwise.
###############################################################################################################


VerifyTransformation := proc(h, g, tau, x)
    local i, f1, f2, y;
    f2 := 0;
    i := 1;
    for i to nops(h) do
        f2 := f2 + ShiftEquivalenceTesting:-Shift(h[i], x, tau[i]) - h[i];
    end do;
    f1 := 0;
    i := 1;
    y := convert(x, list);
    for i to nops(y) do
        f1 := f1 + ShiftEquivalenceTesting:-Shift(g[i], y[i], 1) - g[i];
    end do;
    Testzero(f1 - f2);
end proc:


end module:


################################################################
## Copyright (c) August 2021, UCAS and JKU. All rights reserved.
## Author: Lixin Du <dulixin17@mails.ucas.ac.cn>
## Load the modules:
##         OrbitDecomposition.mm,
##         OrbitalDecomposition.mm.
## Modified on December 15, 2021
#################################################################


RationalReduction := module()
      description "Compute an Orbital Reduction for Multivariate Rational Functions";
      export
              AbramovReduction,
              OrbitalReduction,
              VerifyReduction,
              VerifyDecomposition,
              FromListToRemainder,
              CollectMultiplicity;
      local
              AbramovUniReduction,
              ReduceOneOrbit,
              ReduceOneOrbit2;
      option
               package;


############################# Notation  Description ###########################################
# 0. Setting: the fields F and K are of characteristic zero.
# 1. Shift operators and difference operators:
# Let f be a multivariate rational function in F(x) with x = [x1, ..., xn]. For each variable v
# in x, the shift operator sigma_v is defined as an F-automorphism of F(x) such that
#
#             sigma_v (v) = v + 1,     and     sigma_v (w) = w,  for all w in x and w <> v
#
# The difference operator is defined as
#
#                                Delta_v := sigma_v - 1,
#
# where 1 is the identity map of F(x).
#
# 2. Group action under shifts
# Let G = <sigma_x1, ..., sigma_xn> be the free abelian group generated by shift operators
# sigma_x1, ..., sigma_xn. Let p be a non-constant polynomial in F[x] and H be a subgroup of G.
# The set
#
#             [p]_H := { tau(p) | p in H }
#
# is called the H-orbit of p. Then for two polynomials p, q in F[x],
#
#      [p]_H = [q]_H      iff      p = tau(q) for some tau in H.
#
# In this case, we say p and q are H-equivalent. The isotropy group of p in H is
#
#           Hp := { tau in H | tau(p) = p }.
#
# Remark: In the following algorithms, we shall use the list representation of shifts,
#         (see Data Description in OrbitalDecomposition module).
##############################################################################################


################################ Task  Description ############################################################################
# 1. Orbital reduction for the summability problem
# Let f be a rational function in F(x) with x = [x1, ..., xn]. An orbital reduction of f
# for the summability problem is a list
#
#           [[g1, ..., gn], r],   with r = [[d1, [a11, e11], [a12, e12]...], [d2, [a21, e21], [a22, e22], ...], ...]
#
# satisfying the following properties:
#
#                          f = Delta_x1 (g1) + ... + Delta_xn (gn) + r
# and
#                r = sum_i sum_j aij/ di^j = a11/ d1^e11 + a12/ d1^e12 + ... + a21/ d2^e21 + a22/ d2^e22 + ...,
#
# where d1, ... are irreducible polynomials in x,
# a11, a12, ..., a21, a22, ... are polynomials in x1,
# e11, e12, ..., e21, e22, ... are positive integers,
# deg_x1 (aij) < deg_x1 (di) for all i, j,
# d1, ... are in distinct G-orbits with G = <sigma_x1, ..., sigma_xn>.
#
# 2. Orbital reduction for the existence problem of telescopers
# Let f be a rational function in K(t, x) with x = [x1, ..., xn]. An orbital reduction of f
# for the existence problem of telescopers is a list
#
#        [[g1, ..., gn], r], with r = [[d1, [[k11, [a111, e111], [a112, e112], ...], [k12, [a121, e121], ...], ...]], [d2, ...], ...]
#
# satisfying the following properties:
#
#                          f = Delta_x1 (g1)  + ... +  Delta_xn (gn) + r
# and
#             a111                        a112                                    a121                            a211
#   r = ----------------------- + -------------------------   + ... +  ------------------------ +  ... +  --------------------- + ...
#      sigma_t^{k11} (d1)^e111    (sigma_t^{k11} (d1))^e112             sigma_t^{k12} (d1)^e121         sigma_t^{k21} (d2)^e211
#
# where d1, ... are irreducible polynomials in x,
# a111, a121, ..., a112, ..., a211, ... are polynomials in x1,
# e111, e121, ..., e112, ..., e211, ... are positive integers,
# deg_x1 (aij) < deg_x1 (di) for all i, j,
# d1, ... are in distinct Gt-orbits with Gt = <sigma_t, sigma_x1, ..., sigma_xn>,
# for each i, sigma_t^ki1 (di), sigma_t^ki2 (di), ...are in distinct G-orbits with G = <sigma_x1, ..., sigma_xn>.
##############################################################################################################################



#####################################################################################
# Export
# Name: AbramovReduction
# Calling sequence:
#                    AbramovReduction(a, b, x, tau)
# Input: a,  a multivariate polynomial or rational function;
#        b,  a multivariate polynomial or rational function;
#        x,  a variable name or a list of variable names;
#        tau,  an integer or a list of integers representing a shift operator.
# Output: [[g1, ..., gn], a1],  a list such that
#                           a/ tau(b) = Delta_x1 (g1) + ... + Delta_xn (gn) + a1/ b,
#                       where a1 = tau^{-1} (a) and x = [x1, ..., xn].
#####################################################################################


AbramovReduction := proc(a, b, x_, tau_)
    local n, a1, b1, i, j, g, L, tau, x, r;
    n := nops(x_);
    x := convert(x_, list);
    tau := convert(tau_, list);
    a1 := a;
    b1 := b;
    j := tau;
    g := [];
    i := 1;
    for i from 1 to n do
        j[i] := 0;
        b1 := ShiftEquivalenceTesting:-Shift(b, x, j);
        L := AbramovUniReduction(a1, b1, x[i], tau[i]);
        g := [op(g), L[1]];
        a1 := L[2];
    end do;
    r := a1;
    [g, r];
end proc:


#########################################################################################
# Export
# Name: OrbitalReduction
# Calling sequence:
#                    OrbitalReduction(f, x1, x)
#                    OrbitalReduction(f, x1, x, t)
# Input: f,  a multivariate rational function;
#        x1,  a main variable name in x;
#        x,  a variable name or a list of variable names;
#        t,  (optional) a variable name.
# Output: [[g1, ..., gn], r, Gp], a list such that
#                       f = Delta_x1 ( g1 ) + .... + Delta_xn ( gn ) + r,
#             where g1, ..., gn ara rational functions,
#             r is a list representing the remainder in orbital reduction,
#             Gp is a list consisting of the isotropy group of the denominator of
#             all simple fractions in r, more precisely:
##            if there is no variable t, then r is a list
#                     [[d1, [a11, a12, a13, ...]], [d2, [a21, ...]], ...]
#                representing a remainder of f in the orbital reduction
#                for the summability problem, and
#                Gp = [Gd1, Gd2, ...], where Gdi is a list representing a basis of
#                the isotropy group of di w.r.t. x
##            if there exists a variable t, then r is a list
#                     [[d1, [[[k11], a11, a12, ...], [[k12], ...], ...]], [d2, ...], ...]
#                represents a remainder of f in the orbital reduction
#                for the existence problem of telescopers, and
#                Gp = [Gtd1, Gtd2, ...], where Gdi is a list representing a 'normarlized'
#                basis of the isotropy group of di w.r.t. t and x.
#########################################################################################

OrbitalReduction := proc(f, x1, x, t)
    local F, y, Gp, Gp_, g1, r1, g, r, n, f0, i, orbit, Gd, pair;
    #tt := time();
    n := nops(x);
    if _params['t'] <> NULL then
        y := [t, op(x)];
    else
        y := x;
    end if;


    F := OrbitalDecomposition:-OrbitalPartialFraction(f, x1, y, 'Gp_');
    Gp_;
#### reduce the polynomial part of f which is always summable ###
    f0 := F[1];
    g := convert(LinearAlgebra:-ZeroVector(n), list);
    if f0 <> 0 then
        g1 := sumtools:-gosper(f0, x1);
        g[1] := g1;
    end if;

#### reduce each orbit in orbital partial fractions of f ###
    F := F[2];
    r := [];
    Gp := [];
    i := 1;
    while i <= nops(F) do
        orbit := F[i];
        Gd := Gp_[i];
        if _params['t'] = NULL then
            pair := ReduceOneOrbit(orbit, Gd, x);
        else
            pair := ReduceOneOrbit2(orbit, Gd, x, t);
        end if;
        g1 := pair[1];
        r1 := pair[2];
        Gd := pair[3];
        g := zip((x, y) -> x + y, g, g1);
        if r1 <> [] then
            r := [op(r), r1];
            Gp := [op(Gp), Gd];
        end if;
        i := i + 1;
    end do;
    #tt := time() - tt;
    #print(time_for_orbital_reduction);
    #print(tt);
    [g, r, Gp];

end proc:


#################################################################################################
# Export
# Name: VerifyReduction
# Calling sequence:
#                VerifyReduction(f, g, L, x)
#                VerifyReduction(f, g, L, x, t)
# Input: f,  a multivariate rational function;
#        g,  a multivariate rational function or a list of multivariate rational functions;
#        L,  a list representing a remainder in the orbital reduction;
#        x,  a variable name or a list of variable names;
#        t,  (optional) a variable name.
# Output: true,  if f = Delta_{x1} ( g1 ) + ... + Delta_{xn} ( gn ) + r,
#                where r is a rational function represented by L.
#         false, otherwise.
# Remark: If there is no variable t, then return true if and only if
#                     L = [[d1, [a11, a12, a13, ...]], [d2, [a21, ...]], ...]
#               represents a remainder of f in the orbital reduction
#               for the summability problem.
#         If there exists a variable t, then return true if and only if
#                     L = [[d1, [[[k11], a11, a12, ...], [[k12], ...], ...]], [d2, ...], ...]
#               represents a remainder of f in the orbital reduction
#               for the existence problem of telescopers.
##################################################################################################


VerifyReduction := proc(f, g, L, x, t)
    local r, d1, d2, y, i, j, orbit, sol;

    if _params['t'] = NULL then
        y := x;
        r := FromListToRemainder(L, x);
    else
        y := [t, op(x)];
        r := FromListToRemainder(L, x, t);
    end if;

#### verify orbital conditions ###
    i := 2;
    for orbit in L do
        d1 := orbit[1];
        j := i;
        while j <= nops(L) do
            d2 := L[j][1];
            sol := OrbitalDecomposition:-ShiftEquivalent(d1, d2, y);
            if sol <> [] then
                printf("%a and %a are shift equivalent w.r.t. %a\n", d1, d2, y);
                return false;
            fi;
            j := j + 1;
        end do;
        i := i + 1;
    end do;

    VerifyDecomposition(f, g, r, x);
end proc:


################################################################################################
# Export
# Name: VerifyDecomposition
# Calling sequence:
#                VerifyDecomposition(f, g, r, x)
# Input: f, r,  two multivariate rational functions;
#        g,  a multivariate rational function or a list of multivariate rational functions;
#        x,  a variable name or a list of variable names.
# Output: true  if f = Delta_x1 ( g1 ) + ... + Delta_xn ( gn ) + r,
#                    where g = [g1, ..., gn] and x = [x1, ..., xn].
#         false, otherwise.
#################################################################################################


VerifyDecomposition := proc(f, g, r, x_)
    local n, h, i, flag, x;
    n := nops(x_);

#### if x is single variable and not in a list, write it as a list ###
    x := convert(x_, list);
    h := 0;
    i := 1;
    for i to n do
        h := h + ShiftEquivalenceTesting:-Shift(g[i], x[i], 1) - g[i];
    end do;
    flag := Testzero(f - h - r);
end proc:


##########################################################################################
# Export
# Name: FromListToRemainder
# Calling sequence:
#                FromListToRemainder(L, x)
#                FromListToRemainder(L, x, t)
# Input: L,  a list representing a remainder in orbital reduction;
#        x,  a variable name or a list of variable names;
#        t,  (optional) a variable name.
# Output: r,  the remainder represented by L.
# Remark: If there is no variable t, then
#                     L = [[d1, [a11, a12, a13, ...]], [d2, [a21, ...]], ...]
#               represents a remainder of f in the orbital reduction
#               for the summability problem.
#         If there exists a variable t, then
#                     L = [[d1, [[[k11], a11, a12, ...], [[k12], ...], ...]], [d2, ...], ...]
#               represents a remainder of f in the orbital reduction
#               for the existence problem of telescopers.
#########################################################################################


FromListToRemainder := proc(L, x, t)
    local r;
    if _params['t'] = NULL then
        r := OrbitalDecomposition:-FromListToPartialFraction([0, L], x, input = 'remainder');
    else
        r := OrbitalDecomposition:-FromListToPartialFraction([0, L], [t, op(x)], input = 'CTremainder');
    end if;
    r;
end proc:


#################################################################
# Local
# Name: AbramovUniReduction
# Calling sequence:
#                    AbramovUniReduction(a, b, x, k)
# Input: a,  a univariate polynomial or rational function;
#        b,  a univariate polynomial or rational function;
#        x,  a variable name;
#        k,  an integer.
# Output: [g, a1],  a list of rational functions such that
#                              a/ b(x + k) = Delta_x (g) + a1/ b,
#                   where a1 = a(x - k).
#################################################################


AbramovUniReduction := proc(a, b, x, k)
    local g, a1, i;
    if a=0 then
        g := 0;
        a1 := 0;
        return [g, a1];
    fi;
    if 0 <= k then
        g := add(ShiftEquivalenceTesting:-Shift(a, x, i - k)/ShiftEquivalenceTesting:-Shift(b, x, i), i = 0 .. k - 1);
    else
        g := -add(ShiftEquivalenceTesting:-Shift(a, x, i)/ShiftEquivalenceTesting:-Shift(b, x, k + i), i = 0 .. -k - 1);
    fi;
    a1 := ShiftEquivalenceTesting:-Shift(a, x, -k);
    [g, a1];
end proc:


#########################################################################
# Local
# Name: ReduceOneOrbit
# Calling sequence:
#                ReduceOneOrbit(L, Gd_, x)
# Input: L,  a list [d, [[tau1, a11, a12, ...], [tau2, ...], ...]] representing
#            an orbit in the orbital partial fraction decomposition;
#        Gd_,  a list representing the isotropy group of d w.r.t. x;
#        x,  a variable name or a list of variable names.
# Output: [g, r, Gd], a list such that
##                        [g, r] with r = [d, [c1, c2, ...]]
#                 represents the orbital reduction of L
#                 for the summability problem, and
#                 Gd is the isotropy group of d w.r.t. x.
# Remark: If the remainder of L is zero, then we assign r = [] and Gd = [].
##########################################################################


ReduceOneOrbit := proc(L, Gd_, x)
    local d, F, g, c, tau, i, a1, g1, b, pair, numers, r, Gd, j, e, u, a;
    d := L[1];
    F := L[2];
    g := convert(LinearAlgebra:-ZeroVector(nops(x)), list);
    c := [];

#### there is no reduction for the first tuple, because tau = [0, 0, ...] ###
    c := subsop(1 = NULL, F[1]);
    if nops(F) = 1 then
        return [g, [d, c], Gd_];
    end if;

#### compute a reduction for each simple fraction ###
    i := 2;
    for i from 2 to nops(F) do
        tau := F[i][1];
        numers := subsop(1 = NULL, F[i]);
        j := 1;
        while j<= nops(numers) do
            a := numers[j][1];
            e := numers[j][2];
            b := d^e;
            pair := AbramovReduction(a, b, x, tau);
            g1 := pair[1];
            a1 := pair[2];
            g :=  zip((x, y) -> x + y, g, g1);
            c := [op(c),[a1,e]];

            j := j + 1;
        end do;
    end do;

####  collect all numerators with the same multiplicity ####
    c := CollectMultiplicity(c);

#### prepare the output when c <> [] ###
    if c <> [] then
        Gd := Gd_;
        r := [d, c];
    else
        Gd := [];
        r := [];
    fi;

    [g, r, Gd];
end proc:



##########################################################################
# Local
# Name: ReduceOneOrbit2
# Calling sequence:
#                ReduceOneOrbit2(L, Gtd_, x, t)
# Input: L,  a list [d, [[tau1, a11, a12, ...], [tau2, ...], ...]] representing
#            an orbit in the orbital partial fraction decomposition;
#        Gtd_,  a list representing the isotropy group of d w.r.t. x and t;
#        t,  a main variable name;
#        x,  a variable name or a list of variable names.
# Output: [g, r, Gtd], a list such that
#                       [g, r] with r = [d, [[[k1], a11, a12, ...], [[k2], ...], ...]]
#                 represents the orbital reduction of L for
#                 the existence problem of telescopers, and
#                 Gtd is the isotropy group of d w.r.t. t and x.
# Remark: If the remainder of L is zero, then we assign r = [] and Gtd = [].
##########################################################################


ReduceOneOrbit2 := proc(L, Gtd_, x, t)
    local d, F, g, e, c, tau, a1, g1, b, pair, numers, r, Gd, j, a, m, Gtd, tau0, s, i, c1, k, u, k0;
    u := [u1];
    d := L[1];
    F := L[2];
    g := convert(LinearAlgebra:-ZeroVector(nops(x)), list);
    c := [];

    Gtd := OrbitalDecomposition:-BasisNormalization(Gtd_);
    tau0 := Gtd[1];

#### there is no reduction for the first tuple, because tau = [0, 0, ...] ###
    c := subsop(1 = NULL, F[1]);
    k := 0;
    if nops(F) = 1 then
        return [g, [d, [[k, op(c)]]], Gtd];
    end if;
    c := [[k, op(c)]];

#### compute a reduction for other simple fractions ###
    i := 2;

    for i from 2 to nops(F) do
        tau := F[i][1];
        tau := OrbitalDecomposition:-FromShiftRelationToMin(tau, tau0);
        k0 := tau[1];
        tau := subsop(1 = NULL, tau);
        numers := subsop(1 = NULL, F[i]);
        c1 := [];
        j := 1;
        while j<= nops(numers) do
            a := numers[j][1];
            e := numers[j][2];
            b := ShiftEquivalenceTesting:-Shift(d^e, t, k0);
            pair := AbramovReduction(a, b, x, tau);
            g1 := pair[1];
            a1 := pair[2];
            g :=  zip((x, y) -> x + y, g, g1);
            c1 := [op(c1), [a1, e]];
            j := j + 1;
        end do;
        c := [op(c),[k0, op(c1)]];
    end do;

####  collect all numerators with the same G-orbit and the same multiplicity ####
   c := [ListTools:-Categorize((x,y)-> x[1] = y[1], c)];
   i := 1;
    while i <= nops(c) do
        u := c[i];
        k0 := u[1][1];
        if nops(u) > 1 then
            c1 := [seq(op(subsop(1 = NULL, u[i])), i = 1..nops(u))];
            c1 := CollectMultiplicity(c1);
            if c1 <> [] then
                c := subsop(i = [k0, op(c1)], c);
            else
                c := subsop(i = NULL, c);
                i := i - 1;
            end if;
        else
            c := subsop(i = op(u), c);
        end if;
        i := i + 1;
    end do;


#### test whether the remainder is zero ###
    if c <> [] then
        r := [d, c];
    else
        r := [];
        Gtd := [];
    fi;
    [g, r, Gtd];
end proc:


#######################################################################################
# Local
# Name: CollectMultiplicity
# Calling sequence:
#                CollectMultiplicity(L)
# Input: L,  a list [[a1, e1], ..., [an, en]],
#            where ai are rational functions and ei are positive integers.
# Output: F,  a list [[c1, m1], ..., [ck, mk]] such that
#             {m1 < ... < mk} is a subset of {e1, ..., en}
#             with ci = \sum_{j : ej = mi} aj  and  ci <> 0.
# Remark: this procedure is used to collect all numerators with the same multiplicity.
########################################################################################


CollectMultiplicity := proc(L)
    local F, i, u, a, e;
    F := [ListTools:-Categorize((x,y)-> x[2] = y[2], L)];
    i := 1;
    while i <= nops(F) do
        u := F[i];
        if nops(u) > 1 then
            a := add(u[j][1], j = 1..nops(u));
            e := u[1][2];
            if Testzero(a) = false then
                F := subsop(i = [a,e], F);
            else
                F := subsop(i = NULL, F);
                i := i - 1;
            end if;
        else
            F := subsop(i = op(u),F);
        end if;
        i := i + 1;
    end do;
    F;
end proc:

end module:


################################################################
## Copyright (c) August 2021, UCAS and JKU. All rights reserved.
## Author: Lixin Du <dulixin17@mails.ucas.ac.cn>
## Load the module:
##         ShiftEquivalenceTesting.mm,
## Modified on December 15, 2021
#################################################################


OrbitalDecomposition := module ()
      description "Orbital Irreducible Partial Fraction Decomposition of Multivariate Rational Functions";
      export
              ShiftEquivalent,
              OrbitalFactorization,
              OrbitalPartialFraction,
              VerifyPartialFraction,
              IsotropyGroup,
              IsotropyGroup2,
              FromListToPolynomial,
              FromShiftRelationToMin,
              FromListToPartialFraction,
              BasisNormalization;
      option
              package;


######################## Data  Description ####################################################
# Notation:
# The base field K is a field of characteristic 0. In this algorithm, K can be taken
# as a field of multivariate rational functions Q(u), where Q is the rational number
# field and u = [u1, ..., um] are parameters.
#
# List representation:
# 1. A multivariate rational function f(x1, ..., xn) over K is represented by the
# following list
#
#          [f, x],
#
# where f is a rational function and x = [x1, ..., xn] are variables.
#
# 2. A shift operator tau in the field of multivariate rational functions K(x)
# where x = [x1, ..., xn] is represented by the following list
#
#                  tau = [a1, ..., an],
#
# this means tau(f(x1, ..., xn)) = f(x1 + a1, ..., xn + an) for any f in K(x).
################################################################################################


############################### Task  Description ###################################################
# 1. Shift equivalence testing:
# Given two multivariate polynomials p, q in K[x], decide whether
# there exist integers k1, ..., kn, such that
#
#                       q(x1, ..., xn) = p(x1 + k1, ..., xn + kn).
#
# If exist such integers, find all of them. In this case, we say p and q are
# shift equivalent w.r.t. x.
#
# 2. Isotropy group:
# Let p in K[x]. The set of all integers k = [k1, ..., kn] such that
#
#                      p(x1, ...., xn) = p(x1 + k1, ..., xn + kn)
#
# forms a free abelian group, called the isotropy group w.r.t. x, denoted by Gp.
# Let Hp denote the following subgroup of Gp:
#
#                    Hp = { k = [k1, ..., kn] in Gp | k1 = 0 }.
#
# Then rank (Gp/ Hp) = 0 or 1. Compute rank (Gp/ Hp) and if rank (Gp/ Hp) = 1,
# find a basis [tau_1, ... , tau_r] of Gp such that
#
#                     Hp = <tau_2, ..., taut_r>,   Gp/ Hp = <tau_1 + Hp>.
#
# 3. Orbital factorization of a polynomial:
# An orbital irreducible factorization of a multivariate polynomial p in K[x]
# is a list
#
#   [c, [[d1, [[d11, tau_11, e11], [d12, tau_12, e12], ...]], [d2, [[d21, tau_21, e21], [d22, tau_22, e22], ...]], ...]],
#
# satisfying the following properties:
#
#      p = c * (tau_11 (d1))^e11 * (tau_12 (d1))^e12 * ... * (tau_21 (d2))^e21 * (tau_22 (d2))^e22 * ...
#
# where d1, ... are irreducible polynomials and in distinct shift orbits, that is, for all 1 <= i, j,
#
#           di and dj are shift equivalent w.r.t. x  if and only if  i = j,
#
# c in K, is a constant,
# tau_11, tau_12, ..., tau_21, tau_22, ... are the list representation of shift operators w.r.t. x,
# e11, e12, ..., e21, e22, ... are positive integers,
# d1, .. are "unit normal" irreducible polynomials (see Procedure factors in Maple)
# dij = tau_ij (di) for all 1 <=i, j.
#
# 4. Orbital partial fraction decomposition of a rational function:
# Let f be a multivariate rational function in K(x). An orbital irreducible partial
# fraction decomposition of f w.r.t. one variable x1 is a list
#
# [f0, [[d1, [[tau_11, [a111, e111], [a112, e112], ...], [tau_12, [a121, e121], [a122, e122], ...], ...]], [d2, [[tau_21, [a211, e211], ...], ...]], ...]]
#
# such that
#
#                 a111                  a112                       a121                  a122                       a211
# f = f0 + ------------------- + ------------------  + ... + ------------------ + ------------------ + ... + ------------------ + ...
#          (tau_11 (d1))^e111    (tau_11 (d1))^e112          (tau_12 (d1))^e121   (tau_12 (d1))^e122         (tau_21 (d2))^e211
#
# where di are irreducible polynomials and in distinct shift orbits, that is, for all 1 <= i, j,
#
#     di and dj are shift equivalent w.r.t. x  if and only if  i = j,
#
# f0 is a polynomial in x1,
# ailj are polynomials in x1 and deg_x1 (ailj) < deg_x1(di) for each 1 <= i, j, l,
# eilj are positive integers,
# tau_il are the list representation of shift operators w.r.t. x.
#######################################################################################################


#################################################################
# Export
# Name: ShiftEquivalent
# Calling sequence:
#                     ShiftEquivalent(p, q, x)
#                     ShiftEquivalent(p, q, x, output = 'type')
# Input: p, q,  two multivariate polynomials;
#        x,  a variable name or a list of variable names;
#        output,  (optional) should be a type of solutions in the form:
#                         output = 'particular' or 'basis'.
# Output: L,  a list representing the shift equivalence relation of p and q w.r.t. x;
#         [],  if p and q are not shift equivalent.
#
##             If there is no condition on 'type', then return L as a general solution of
#                                     q(x) = p(x + a).
##             If 'type' = 'particular', then return L as a particular solution.
##             If 'type' = 'basis', then return L as a list
#                                       [a0, [a1, ..., ar]],
#                        such that a general solution is in the form
#                                   a0 + Z1 * a1 + ... + Zr * ar,
#                        where a0, a1, ..., ar are lists of integers and
#                        Z1, ..., Zr are arbitrary integers.
#                  In this case, if there exists only one solution a = a0 s.t. p (x) = q(x + a),
#                  then return  L as a list
#                                          [a0, [[0,0,...]]].
# Remark: 1. If p = q, then there is at least one trivial solution a = [0, 0, ..., 0],
#            such that p (x) = p(x + a).
#################################################################


ShiftEquivalent := proc(p, q, x, {output := 'general'})
    local sol, L, n, L0, L1, i, var, value, y, F;

#### compute a general solution ###
    sol := ShiftEquivalenceTesting:-SET_ADP_expansion(p, q, x);
    L := sol;

####  type of solutions ###
    if L <> [] then
        var := convert(indets(sol), list);
        n := nops(var);
        if output = 'particular' or output = 'basis' then
            # compute a particular solution #
            value := convert(LinearAlgebra:-ZeroVector(n), list);
            y := zip((x, y)-> x = y, var, value);
            L0 := subs(y, sol);
            if output = 'particular' then
                L := L0;
            elif output = 'basis' then
                L1 := [];
                sol := zip((x, y) -> y - x, sol, L0);
                if n <> 0 then
                    for i from 1 to n do
                        value := convert(Vector(n, shape = unit[i]), list);
                        y := zip((x, y) -> x = y, var, value);
                        L1 := [op(L1), subs(y, sol)];
                    end do;
                else
                    L1 := [convert(LinearAlgebra:-ZeroVector(nops(x)), list)];
                end if;
                L := [L0, L1];
            end if;
        end if;
    end if;
    L;
end proc:


################################################################################################
# Export
# Name: OrbitalFactorization
# Calling sequence:
#                     OrbitalFactorization(p, x)
#                     OrbitalFactorization(p, x, 'Gp')
# Input: p,  a multivariate polynomial;
#        x,  a variable name or a list of variable names.
# Output: L, a list
#                          [c, [[d1, [[d11, tau11, e11], [d12, tau12, e12], ...]], [d2, ...], ..]],
#            representing an orbital irreducible factorization of p w.r.t. x.
# Optional: the third argument if presented is assigned
#                           Gp = [Gd1, Gd2, ...]
#           with Gdi being the isotropy group of di, for i = 1, 2, ...
# Remark: c is a constant,  d11 = tau11(d1), d12 = tau12(d1), ...
################################################################################################


OrbitalFactorization := proc(p, x, Gp)
    local irrlist, n, pc, pp, d1, e1, tau1, flag, equivclass, d2, pair, tau2, e2, s, sol, i, tau, L, R, Gp_, Gd1;
    #tt := time();
    pc := content(p, x);
    pp := primpart(p, x);
    n := nops(x);
    irrlist := factors(pp);
    pc := pc * irrlist[1]; ## pc may not be 1
    irrlist := irrlist[2];
    L := [];
    if nargs = 3 then
        Gp_ := [];
    end if;
    while irrlist<>[] do
        d1 := irrlist[1][1];
        e1 := irrlist[1][2];
        tau1 := convert(LinearAlgebra:-ZeroVector(n), list);
        equivclass := [[d1, tau1, e1]];
        irrlist := subsop(1 = NULL, irrlist);
        s := nops(irrlist);

    #### find all factors that are shift equivalent to d1 ###
        i := 1;
        flag := false;
        while i<= s do
            pair := irrlist[i];
            d2 := pair[1];

            # check whether there exists tau s.t. d2 = tau(d1)
            # note: d2 and d1 have the same isotropy group
            if nargs = 2 or flag = true then
                sol := ShiftEquivalent(d1, d2, x, output = 'particular');
            else
                sol := ShiftEquivalent(d1, d2, x, output = 'basis');
            end if;
            if sol <>[] then
                if nargs = 2 or flag = true then
                    tau := sol;
                else
                    tau := sol[1];
                    Gd1 := sol[2];
                end if;
                e2 := pair[2];
                tau2 := tau;
                equivclass := [op(equivclass), [d2, tau2, e2]];
                irrlist := subsop(i = NULL, irrlist);
                flag := true;
                s := s - 1;
                i := i - 1;
            end if;
            i := i + 1;
        end do;
        if nargs = 3 then
            if flag = false then
                Gd1 := IsotropyGroup(d1, x);
            end if;
            Gp_ := [op(Gp_), Gd1];
        end if;
        L := [op(L), [d1, equivclass]];
    end do;
    #tt := time() - tt;
    #print(time_for_orbital_factorization);
    #print(tt);
    if nargs = 3 then
        Gp := Gp_;
    end if;
    L := [pc, L];
end proc:


#################################################################
# Export
# Name: OrbitalPartialFraction
# Calling sequence:
#                     OrbitalPartialFraction(f, x1, x)
#                     OrbitalPartialFraction(f, x1, x, 'Gp')
# Input: f,  a multivariate rational function;
#        x1,  a main variable name;
#        x,   a variable name or a list of variable names.
# Output: L,  a list [f0, [[d1, [[tau_11, [a111, e111], [a121, e121], ...], [tau_12, ...], ...]], [d2, ...], ...]] representing
#             an orbital irreducible partial fraction decomposition of f w.r.t. x1 over Q(u).
# Optional: the third argument if presented is assigned
#                           Gp = [Gd1, Gd2, ...]
#                  with Gdi being the isotropy group of di, for i = 1, 2, ...
# Remark: 1. tau_ij are shift operators w.r.t. x.
#         2. x contains x1.
#################################################################


OrbitalPartialFraction := proc(f_, x1, x, Gp)
    local f,  p, q, L, tt, polyorbits, Gp_, flist, e, col, factorclass, i, equivclass, j, F, f0, d, orbit, Numer, a, tau;

#### f = p/ q ###
    ##tt := time();
    f := normal(f_);
    p := numer(f):
    q := denom(f):
    q := primpart(q, x1):
    p := normal(expand(f * q));
    L := OrbitalFactorization(q, x,'Gp_');
    p := p/ L[1]: ## L[1] may not be 1
    polyorbits := L[2]:

#### prepare a factor list for partial fraction decomposition ###
    flist := []:
    i := 1:
    while i <= nops(polyorbits) do
        equivclass := polyorbits[i][2]:

        # delete tau in equivclass #
        # slower # factorclass := map(L -> subsop(2 = NULL, L), equivclass);
        j := 1;
        factorclass :=[]:
        while j <= nops(equivclass) do
            factorclass := [op(factorclass), subsop(2 = NULL, equivclass[j])]:
            j := j + 1:
        end do:
        flist := [op(flist), op(factorclass)]:
        i := i + 1;
    end do:
    flist:

#### obtain a partial fraction decomposition ###
   # printf("p = %a, flist = %a, x1 = %a \n", p, flist, x1);

    F := convert([p, op(flist)], parfrac, x1);
    f0 := F[1];
    F := subsop(1 = NULL, F);

#### insert tau_ij and obtain an orbital partial fraction decomposition ###
    L := [];
    i := 1;
    while i<= nops(polyorbits) do
        equivclass := polyorbits[i];
        d := equivclass[1];
        orbit := equivclass[2];
        j := 1;

    #### collect partial fractions which are in the same orbit as d ###
        Numer := [];
        while j<= nops(orbit) do
            tau := orbit[j][2];

        #### collect numerators with denominater tau(d) ###
        #### obtain a nested list representation ##########
            a := subsop(1 = NULL,F[1]);
            if a <> 0 then
                e, col, a :=ArrayTools:-SearchArray(Array(a));
                e := convert(e, list);
                a := convert(a, list);
                a := ListTools:-Interleave(a,e);
                a := ListTools:-LengthSplit(a,2);
                F := subsop(1 = NULL, F);
                a := [tau, a];
                Numer := [op(Numer), a];
            else
                if nargs = 4 then
                    Gp_ := subsop(j = NULL, Gp_);
                    orbit := subsop(j = NULL, orbit);
                    j := j - 1;
                end if;
            end if;
            j := j + 1;
        end do;
        L := [op(L), [d, Numer]];
        i := i + 1;
    end do;

    if nargs = 4 then
        Gp := Gp_;
    end if;
    #tt := time() - tt;
    #print(time_for_orbital_partial_fraction);
    #print(tt);
    L := [f0, L];
end proc:


#################################################################################
# Export
# Name: FromListToPolynomial
# Calling sequence:
#                     FromListToPolynomial(L, x)
# Input: L,  a list  [c, [[d1, [[d11, tau11, e11], [d12, tau12, e12], ...]], [d2, ...], ..]]
#            representing an orbital factorization of a multivariate polynomial w.r.t. x;
#        x,   a variable name or a list of variable names.
# Output: p,  a multivariate polynomial represented by L.
################################################################################


FromListToPolynomial := proc(L, x)
    local F, p, u, c, d, d1, orbit, tau;
    c := L[1];
    F := L[2];
    p := c;
    for orbit in F do
        d := orbit[1];
        for u in orbit[2] do
            tau := u[2];
            d1 := ShiftEquivalenceTesting:-Shift(d, x, tau);
            p := p * d1^u[3];
        end do;
    end do;
    p;
end proc:


##############################################################################################
# Export
# Name: FromListToPartialFraction
# Calling sequence:
#                     FromListToPartialFraction(L, x)
#                     FromListToPartialFraction(L, x, input = 'type')
# Input: L,  a list representing a multivariate rational function w.r.t. x;
#        x,   a variable name or a list of variable names;
#        input,  (optional) should be a type of L in the form:
#                         input = 'remainder' or 'CTremainder'.
# Output: f, a multivariate rational function represented by L.
#
#         If there is no condition on 'tpye', then L should be a list in the form
#                 [f0, [[d1, [[tau_11, a111, a121, ...], [tau_12, ...], ...]], [d2, ...], ...]]
#             representing an orbital partial fraction decomposition w.r.t. x
#         If 'type' = 'remainder', then L should be a list in the form
#                 [0, [[d1, [a11, a12, ...]], [d2, [a21, a22...]], ...]]
#             representing a remainder in additive decomposition for the summability problem,
#             (see Task Description in RationalReduction module).
#         If 'type' = 'CTremainder', then L should be a list in the form
#                  [0, [[d1, [[[k1], a11, a12, ...], [[k2], ...], ...]], [d2, ...], ...]]
#             representing a remainder in additive decomposition for the existence problem
#             of telescopers, where ki represents the shift operator w.r.t. the first variable
#             of x, (see Task Description in RationalReduction module).
##############################################################################################


FromListToPartialFraction := proc(L, x_, {input := 'orbital'})
    local f0, F, f, numers, x, orbit, a, u, e, j, d, d1, tau;
    f0 := L[1];
    F := L[2];
    f := f0;
    x := convert(x_, list);
    for orbit in F do
        d := orbit[1];
        orbit := orbit[2];
        if input='remainder' then
            orbit := [orbit];
        end if;
        for numers in orbit do
            if input = 'remainder' then
                d1 := d;
            elif input='CTremainder' then
                tau := numers[1];
                d1 := ShiftEquivalenceTesting:-Shift(d, x[1], tau);
                numers := subsop(1 = NULL, numers);
            else
                tau := numers[1];
                d1 := ShiftEquivalenceTesting:-Shift(d, x, tau);
                numers := subsop(1 = NULL, numers);
            end if;
            j := 1;
            for u in numers do
                a := u[1];
                e := u[2];
                f := f + a/ d1^e;
                j := j + 1;
            end do;
        end do;
    end do;
    f;
end proc:


##################################################################################################
# Export
# Name: VerifyPartialFraction
# Calling sequence:
#                       VerifyPartialFraction(L, f, x1, x)
# Input: L,  a list [f0, [[d1, [[tau_11, a111, a121, ...], [tau_12, ...], ...]], [d2, ...], ...]];
#        f,  a multivariate rational function;
#        x1,  a main variable name;
#        x,  a variable name or a list of variable names.
# Output: true,  if L is a list representing orbital irreducible partial fraction
#                      decomposition of f w.r.t. x1 over K
#         false,  otherwise.
# Remark: 1. verify degree conditions: deg_x1 (aijl) < deg_x1(di)
#         2. verify orbital conditions:
#                         di and dj are not shift equivalent w.r.t. x if i <>j
#         3. verify  f - f0 - sum_i sum_j sum_l aijl/ tau_ij (di)^(j) = 0
##################################################################################################


VerifyPartialFraction := proc(L, f, x1, x)
    local f0, F, g, numers, orbit, a, d, d1, tau, u, e, i, j, d2, sol;
    f0 := L[1];
    F := L[2];
    if degree(f0, x1) < 0 and f0 <> 0 then
        printf("%a is not a polynomial in %a \n", f0, x1);
        return false;
    fi;
    #g := f0;

#### verify degree conditions ###
    for orbit in F do
        d := orbit[1];
        orbit := orbit[2];
        for numers in orbit do
            tau := numers[1];
            d1 := ShiftEquivalenceTesting:-Shift(d, x, tau);
            numers := subsop(1 = NULL, numers);
            #j := 1;
            for u in numers do
                a := u[1];
                e := u[2];
                if degree(a, x1) >= degree(d1, x1) then
                    printf("deg_{%a}(%a)<deg_{%a}(%a) is not true \n", x1, a, x1, d1);
                    return false;
                end if;
                # g := g + a/ (d1^(e));
                # j := j + 1;
            end do;
        end do;
    end do;

#### verify orbital conditions ###
    i := 2;
    for orbit in F do
        d1 := orbit[1];
        F := subsop(1 = NULL, F);
        j := i;
        while j <= nops(F) do
            d2 := F[j][1];
            sol := ShiftEquivalent(d1, d2, x);
            if sol<>[] then
                printf("%a and %a are shift equivalent w.r.t. %a\n", d1, d2, x);
                return false;
            end if;
            j := j + 1;
        end do;
        i := i + 1;
    end do;
    g := FromListToPartialFraction(L, x);
    Testzero(f - g);
end proc:


#####################################################################################################
# Export
# Name: IsotropyGroup
# Calling sequence:
#                      IsotropyGroup(p, x)
# Input: p,  a multivariate polynomial;
#        x,  a variable name or a list of variable names.
# Output: L,  a list [tau_1, ... , tau_r] representing a basis of the isotropy group Gp,
##             if Gp is a nontrivial group,
#         [[0, 0, ..., 0]],  if Gp is trivial.
####################################################################################################


IsotropyGroup := proc(p, x)
    local sol, n, tau, L, v;
    sol := ShiftEquivalent(p, p, x, output = 'basis');
    sol := sol[2];
    if sol = [] then
        n := nops(x);
        tau := LinearAlgebra:-ZeroVector(n);
        L := [convert(tau, list)];
    else
        L := sol;
    end if;
    L;
end proc:


###################################################################################################
# Export
# Name: IsotropyGroup2
# Calling sequence:
#                       IsotropyGroup2(p, t, x)
# Input: p, a multivariate polynomial in t and x;
#        t, a variable name;
#        x, a variable name or a list of variable names.
# Output: L, a list [tau_0, [tau_1, ..., tau_r]] representing a basis of the isotropy
#            group G_{t, p} w.r.t t and x such that
#                              G_p = <tau_1, ...taut_r>   and   G_{t, p}/ G_p = <tau_0 + G_p>,
#            where G_{t, p} is the isotropy of p with respect to t and x,
#            G_p is the isotropy group of p with respect to x.
# Remark: If G_{t, p}/ G_p is a trivial group, i.e., rank (G_{t, p}/ G_p) = 0,
#         then assign tau_0 as [0, 0, ..., 0].
###################################################################################################


IsotropyGroup2 := proc(p, t, x)
    local Gpt, L;
    Gpt := IsotropyGroup(p, [t, op(x)]);
    L := BasisNormalization(Gpt);
end proc:


#################################################################
# Export
# Name: FromShiftRelationToMin
# Calling sequence:
#                      FromShiftRelationToMin(phi, tau0)
# Input: phi,  a list representing a shift operator;
#        tau0,  a list representing an element in isotropy group with minimal index
#               for the first variable.
# Output: phi0,  a list representing a shift operator with minimal index
#                for the first variable (modulo tau 0).
# Remark: the input tau0 is allowed to be a zero list.
#################################################################


FromShiftRelationToMin := proc(phi, tau0)
    local a, b, q, T, r, phi0;
    if tau0[1] = 0 then
        phi0 := phi;
    else
        a := phi[1];
        b := tau0[1];
        r:=a mod b;
        q := iquo(a - r, b);
        T := (x, y) -> x - q * y;
        phi0 := zip(T, phi, tau0);
    end if;
    phi0;
end proc:


#######################################################################################
# Export
# Name: BasisNormalization
# Calling sequence:
#                    BasisNormalization(Gp)
# Input: Gp,  a list [sigma_1, sigma2, ...] representing a basis of the isotropy group.
# Output: L,  a list [tau_0, [tau_1, ...]] representing a 'normalized' basis of the
#                 isotropy group Gp w.r.t. the first variable.
########################################################################################


BasisNormalization := proc(Gp)
    local A, n, tau_0, Hp, L;
    A := convert(Gp, Matrix);
    n := nops(Gp[1]);

#### In fact, it normalizes Gp w.r.t. x1, x2, ...respectively ###
    A := LinearAlgebra:-HermiteForm(A);
    if A[1, 1] = 0 then
        tau_0 := convert(LinearAlgebra:-ZeroVector(n), list);
    Hp := convert(A, list, nested);
    else
        tau_0 := convert(LinearAlgebra:-Row(A, 1), list);
        Hp := subsop(1 = NULL, convert(A, list, nested));
        if Hp = [] then
            Hp := [convert(LinearAlgebra:-ZeroVector(n), list)];
        end if;
    end if;
    L := [tau_0, Hp];
end proc:


end module:


############################################################################
## Copyright (c) September 2021, Beihang University. All rights reserved.
## Author: Hanqian Fang <hqfang_math@163.com>
## Modified on September, 2021
############################################################################


ShiftEquivalenceTesting := module()
     description "Test whether two polynomials are shift equivalent";
     export
             Shift,
             SET_ADP_expansion;
     local
             DegClassify,
             HomogeneousNoMoreOne;
     option
             package;

######################################################################
# Export
# Name: SET_ADP_expansion
# Calling sequence:
#                   SET_ADP_expansion(p, q, x)
# Input: p, q,  two multivariate polynomials;
#        x,  a variable name or a list of variable names.
# Output: k,  a list representation of a general integer solution of
#                          q(x) = p(x + k).
#         [],  if there is no integer solutions.
# Remark: 1. A general solution may contains some free variables.
#         2. This is an integer version of SET_ADP_expansion
######################################################################


SET_ADP_expansion := proc(p, q, x)
    local n, y, f, g, dl_sub, ht_sub, d_sub, hl, i, j, a, sol_particular, eqs, expr, exprs, tail, eq, para, k, vars, globalvars, parameter;
    n := nops(x);
    y := convert(x, list);
    f := expand(p);
    g := expand(q);
    i := 'i';
    a := [seq(i[j], j = 1 .. n)];
    if f = 0 and g = 0 then
        return a;
    end if;
    if degree(f, convert(y, set)) <= degree(expand(f - g), convert(y, set)) then
        return [];
    end if;
   
#### extract and classify the coefficients of f(y + i) - g(y) w.r.t. y, which are polynomials in i ####
    dl_sub, ht_sub := DegClassify([coeffs(expand(Shift(f, y, a) - g), y)], a);
    sol_particular := [seq(i[j] = 1, j = 1 .. n)];
    eqs := {};
    for d_sub in dl_sub do

    #### for each polynomial f in classification, substitute one particular solution into f - H^1 (f) - H^0 (f)) in i ####
        hl := ht_sub[d_sub];
        exprs := {};
        for expr in hl do
            tail := HomogeneousNoMoreOne(expr, a);
            expr := subs(sol_particular, expr - tail) + tail;
            exprs := {op(exprs), expr};
        end do;

        #### extract the coefficients of every polynomial of G w.r.t parameters, if exist ###
        vars := indets(exprs);
        globalvars := convert(a, set);
        parameter := vars minus globalvars;
        exprs := map(coeffs, exprs, parameter);

        #### solve the linear system ####
        eq := solve(exprs, convert(a, set));
        if eq = () or isolve(eq) = () then
            return [];
        end if;
        eqs := solve(eqs union eq, convert(a, set));
        if eqs = () then
            return [];
        else
            eq := isolve(eqs);
        end if;
        if eq = () then
            return [];
        end if;
        eq := subs(eq, a);

        #### update the particular solution ####
        para := convert(indets(eq), list);
        k := nops(para);
        sol_particular := subs(zip((x, y) -> x = y, para, convert(LinearAlgebra:-ZeroVector(k), list)), eq);
        sol_particular := [seq(i[j] = sol_particular[j], j = 1 .. n)];
    end do;
    return eq;
end proc:


#################################################################
# Export
# Name: Shift
# Calling sequence:
#                     Shift(f, x, k)
# Input: f,  a multivariate polynomial or rational function;
#        x,  a variable name or a list of variable names;
#        k,  an integer or a list of integers.
# Output: g,   a rational function such that
#
#                             g(x) = f(x + a).
#
# Remark: If x = [x1, ..., xn] and a = [a1, ..., an] are two lists,
#              then the notation g(x) = f(x + a) means
#
#                    g(x1, ..., xn) = f(x1 + a1, ..., xn + an).
#################################################################


Shift := proc(f, x, a)
    local y, g;
    y := zip((x, y) -> x = x + y, x, a);
    g := subs(y, f);
    g;
end proc:


######################################################################
# Local
# Name: DegClassify
# Calling sequence:
#                     DegClassify(l,x)
# Input: l,  a list of multivariate polynomials;
#        x,  a variable name or a list of variable names.
# Output: ht,   a table such that
#              ht = ([d_1 = H^{d_1},d_2 = H^{d_2},...,d_k =H^{d_k}]);
#         dl,    a list such that
#              dl = [d_1,d_2, ..., d_k];
#         where d_1 < d_2 < ... < d_k;
#         H^{d_i} is a sublist of l where its element f satsfying
#               totaldegree(f) = d_i.
#####################################################################


DegClassify := proc(l, x)
    local g, y, lc, lm, tl, d, i, ht, dl, hl;
    y := convert(x, set); ht := table();
    dl := [];
    for i in l do
        i := expand(i);
        d := degree(i, y);
        if d in dl then
            ht[d] := [op(ht[d]), i];
        else
            ht[d] := [i];
            dl := [op(dl), d];
        end if;
    end do;
    dl := sort(dl, `<`);
    dl, ht;
end proc:


##################################################################
# Local
# Name: HomogeneousNoMoreOne
# Calling sequence:
#                      HomogeneousNoMoreOne(f,x)
# Input: f,  a multivariate polynomial;
#        x,  a variable name or a list of variable names.
# Output: expr,   an expression such that
#                  expr = H^1+ H^0;
#         where f = H^d+ H^{d-1} + ... + H^0;
#         H^i is homogeneous of degree i.
##################################################################


HomogeneousNoMoreOne := proc(f, x)
    local g, y, lc, lm, tl, d, i, expr;
    g := expand(f);
    y := convert(x, set);
    lc, lm := Groebner:-LeadingTerm(g, plex(op(y)));
    if lc * lm - g = 0 then
        tl := [g];
    else
        tl := [op(g)];
    end if;
    expr := 0;
    for i in tl do
        d := degree(i, y);
        if d = 0 or d = 1 then
            expr := expr + i;
        end if;
    end do;
    expr;
end proc:
end module: