################################################################ ## Copyright (c) August 2021, UCAS and JKU. All rights reserved. ## Author: Lixin Du ## 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 = be the free abelian group generated by shift operators # sigma_t, sigma_{x1}, ..., sigma_{xn}. Let G = 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) 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) rather than K(t, y). # 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) 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), 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) 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) 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 = and Gd = . ### 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) ## # iff L1 ( phi (a)) = 0 for some L1 in K(t) ################################ 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) 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), 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 = and Gd = . ### 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) ## # iff L1 ( phi (a)) = 0 for some L1 in K(t) ################################ 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 := ; A := ; 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); # 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), 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 ## 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 = 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 ## 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 = 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 = . # # 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 = , # for each i, sigma_t^ki1 (di), sigma_t^ki2 (di), ...are in distinct G-orbits with G = . ############################################################################################################################## ##################################################################################### # 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 ## 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 = , Gp/ 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)[] 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 = and G_{t, p}/ 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 ## 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: