BeginPackage["qDesingularization`", {"Singular`"}]
`InhomPart::usage = "InhomPart[rec, f] gives the inhomogeneous part of rec, where rec is a q-difference equation with unknown funcion f";
`HomPart::usage = "HomPart[rec, f] gives the homogeneous part of rec, where rec is a q-difference equation with unknown funcion f";
`HomNormRe::usage = "HomNormRe[rec, f, q] takes the homogeneous part of a q-difference rec and make a normalization f[a_] -> f[a]/(q^a - 1), where f[n] is an unknown function appearing in rec";
`ReToOp::usage = "ReToOp[rec, f, n, S] transforms a homogeneous q-difference equation rec into a q-difference operator in QQ(q, x)[S], where f[n] is an unknown function appearing in rec";
`OpToRe::usage = "OpToRe[P, B, f, n, g] transforms a q-difference operator B into an inhomogeneous q-difference equation with unknown function f[n], inhomogeneous part (B/P)(g), where B/P is the left quotient of B w.r.t. a q-difference operator P, g is the inhomogeneous part of another q-difference equation";
`qDispersionIrr::usage = "qDispersionIrr[u, v, q, x] gives the dispersion of u and v, where u(x) and v(x) are two primitive irreducible polynomials in QQ[q][x] with positive degree in x, u(0) is not equal to 0";
`qDispersion::usage = "qDispersion[u, v, q, x] gives the dispersion of u and v using irreducible factorization over QQ[q][x], where u(x) and v(x) are two polynomials in QQ[q][x], and u(0) is not equal to 0";
`qDispersionRes::usage = "qDispersionRes[u, v, q, x] gives the dispersion of u and v using resultant-based algorithm, where u(x) and v(x) are two polynomials in QQ[q][x], and u(0) is not equal to 0";
`OrderBound::usage = "OrderBound[P, q, x, S] gives an order bound for a desingularized operator of P, which is a q-difference operator in QQ[q][x][S], but regarded as an operator in QQ(q)[x][S]";
`ContSubmodule::usage = "ContSubmodule[P, k, x, S] gives the k-th submodule of the q-Weyl closure of P, where P is a q-difference operator in QQ[q][x][S], but regarded as an operator in QQ(q)[x][S]";
`qWeylClosure::usage = "qWeylClosure[P, k, x, S] gives a generating set of the q-Weyl closure of P, where P is a q-difference operator in QQ[q][x][S], but regarded as an operator in QQ(q)[x][S]";
`MinCoeffIdealGB::usage = "MinCoeffIdealGB[Id, q, x] gives {g, c} using Groebner bases over QQ[q][x], where g is an element in the ideal generated by Id in QQ(q)[x] with minimial degree in x, and c is a finite set of polynomials in QQ[q, x] such that g = c[1] Id[1] + ... + c[k] Id[k], where k is the length of Id";
`MinCoeffIdealEuclidean::usage = "MinCoeffIdealEuclidean[Id, x] gives {g, c} using Euclidean algorithm over QQ(q)[x], where g is an element in the ideal generated by Id in QQ(q)[x] with minimial degree in x, and c is a finite set of polynomials in QQ(q)x] such that g = c[1] Id[1] + ... + c[k] Id[k], where k is the length of Id";
`LeadDesingularizedOperator::usage = "LeadDesingularizedOperator[P, q, x, S] gives a leading desingularized operator of P, where P is a q-difference operator in QQ[q][x][S], but regarded as an operator in QQ(q)[x][S] ";
`TrailDesingularizedOperator::usage = "TrailDesingularizedOperator[P, q, x, S] gives a trailing desingularized operator of P, where P is a q-difference operator in QQ[q][x][S], but regarded as an operator in QQ(q)[x][S]";
`BiDesingularizedOperator::usage = "BiDesingularizedOperator[P, q, x, S] gives a bi-desingularized operator of P, where P is a q-difference operator in QQ[q][x][S], but regarded as an operator in QQ(q)[x][S]";
Begin["Private`"]
(*
This is a hack for being able to use HolonomicFunctions commands inside the Private context.
We cannot load HolonomicFunctions in BeginPackage. We have to load it at the very end.
*)
`OreAlgebra = RISC`HolonomicFunctions`OreAlgebra;
`ToOrePolynomial = RISC`HolonomicFunctions`ToOrePolynomial;
`ApplyOreOperator = RISC`HolonomicFunctions`ApplyOreOperator;
`OreReduce = RISC`HolonomicFunctions`OreReduce;
`LeadingCoefficient = RISC`HolonomicFunctions`LeadingCoefficient;
`LeadingExponent = RISC`HolonomicFunctions`LeadingExponent;
`OrePolynomialListCoefficients = RISC`HolonomicFunctions`OrePolynomialListCoefficients;
`OreTimes = RISC`HolonomicFunctions`OreTimes;
`Extended = RISC`HolonomicFunctions`Extended;
(*
Name: InhomPart
Input: rec - a univariate inhomogeneous linear q-difference equation
f - a formal function appearing in rec
Output: g - the inhomogeneous part of rec
*)
InhomPart[rec_, f_] := (rec /. f[_] -> 0);
(*
Name: HomPart
Input: rec - a univariate inhomogeneous linear q-difference equation
f - a formal function appearing in rec
Output: rec1 - the homogeneous part of rec
*)
HomPart[rec_, f_] := Collect[rec - InhomPart[rec, f], f[_]];
(*
Name: HomNormRe
Input: rec - a univariate inhomogeneous linear q-difference equation
f - a formal function appearing in rec
q - a variable which is not a root of unity of QQ
Output: rec1 - the homogeneous part of rec and make a nomalization: f[a_] -> f[a]/(q^a - 1)
*)
HomNormRe[rec_, f_, q_] :=
Module[{a},
Return[Together[HomPart[rec, f] /. f[a_] -> f[a]/(q^a - 1)]];
];
(*
Name: ReToOp
Input: rec - a univariate homogeneous linear q-difference equation
f - a formal function appearing in rec
n - a variable appearing in rec
S - a q-dilation operator
Output: P - a nonzero q-difference operator in QQ[q][x][S], which represents rec
*)
ReToOp[rec_, f_, n_, S_] := ToOrePolynomial[rec, f[n], OreAlgebra[S]];
(*
Name: OpToRe
Input: P - a nonzero q-dilation operator in QQ[q][x][S], which represents
the homogeneous part of some q-difference equation rec_old
B - a desingularized operator (three kinds) of P, which belongs to QQ[q][x][S], but regarded as an element in QQ(q)[x][S]
f - a formal function appearing in rec_old
n - a variable appearing in rec_old
g - the inhomogeneous part of rec_old
Output: rec - a new univariate inhomogeneous linear q-difference equation whose homogeneous part is B(f) and
inhomogeneous part is (B/P)(g)
*)
OpToRe[P_, B_, f_, n_, g_] := ApplyOreOperator[B, f[n]] + Together[ApplyOreOperator[OreReduce[B, {P}, Extended -> True][[3, 1]], g]];
(*
Name: qDispersionIrr
Input: u - a primitive irreducible polynomial in QQ[q][x] of positive degree in x, where u is not divisible by x
v - a primitive irreducible polynomial in QQ[q][x] of positive degree in x
q - a variable which is not a root of unity of QQ
x - an indeterminate, which is equal to q^n, n is also a variable
Output: m - a nonnegative integer which is equal to
max{ r, 0 | r in NN, deg gcd(u(q^r x), v(x)) > 0 }
*)
qDispersionIrr[u_, v_, q_, x_] :=
Module[{d, p, e},
If[(u /. x -> 0) === 0, Throw["the first argument is divisible by x, please enter a new one"]];
d = Exponent[u, x]; (*take the degree of u w.r.t. x *)
If[Not[d === Exponent[v, x]] || ((v /. x -> 0) === 0), Return[0], (* if the degree of u and v w.r.t. x is different, then the q-dispersion is 0 *)
p = Together[(Coefficient[u, x, 0] Coefficient[v, x, d])/(Coefficient[v, x, 0] Coefficient[u, x, d])];
If[Not[PolynomialQ[p, q]], Return[0], (* if p is not a polynomial w.r.t. q, then the q-dispersion is 0 *)
e = Exponent[p, q];
If[Not[Expand[p - q^e] === 0] , Return[0], (* if p is not equal to q^e, then the q-dispersion is 0 *)
If[Not[Mod[e, d] === 0], Return[0], (* if e is not divisible by d, then the q-dispersion is 0 *)
If[Not[Expand[(u /. x -> q^(e/d) x)/(q^e Coefficient[u, x, d]) - v/Coefficient[v, x, d]] === 0], Return[0], Return[e/d]];
];
];
];
];
];
(*
Name: qDispersion (* compute the q-dispersion by the factorization-based algorithm *)
Input: u - a polynomial in QQ[q][x], where u is not divisible by x
v - a polynomial in QQ[q][x]
q - a variable which is not a root of unity
x - an indeterminate, which is equal to q^n, n is also a variable
Output: m - a nonnegative integer which is equal to
max{ r, 0 | r in NN, deg gcd(u(q^r x), v(x)) > 0 }
*)
qDispersion[u_, v_, q_, x_] :=
Module[{u1, v1, r1, r2, m, a},
If[(u /. x -> 0) === 0, Throw["the first argument is divisible by x, please enter a new one"]];
If[Exponent[u, x] === 0 || Exponent[v, x] === 0, Return[0]]; (* if u or v is a polynomial in QQ[q], then the dispersion of u and v is equal to 0 *)
u1 = Last[FactorTermsList[u, x]]; (* pull out the content of u w.r.t. x, u1 is the primitive part of u *)
v1 = Last[FactorTermsList[v, x]]; (* pull out the content of v w.r.t. x, v1 is the primitive part of v *)
r1 = First /@ Delete[FactorList[u1], 1]; (* compute a list of irreducible factors of u1 without content and exponents *)
r2 = First /@ Delete[FactorList[v1], 1]; (* compute a list of irreducible factors of v1 without content and exponents *)
m = 0; (* initialize the q-dispersion to be 0 *)
For[i = 1, i <= Length[r1], i++,
For[j = 1, j <= Length[r2], j++,
a = qDispersionIrr[r1[[i]], r2[[j]], q, x]; (* compute the q-dispersion of r1[[i]] and r2[[j]] *)
If[a > m, m = a];
];
];
Return[m];
];
(*
Name: qDispersionRes (* compute the q-dispersion by the resultant-based algorithm *)
Input: u - a polynomial in QQ[q][x], where u is not divisible by x
v - a polynomial in QQ[q][x]
q - a variable which is not a root of unity
x - an indeterminate, which is equal to q^n, n is also a variable
Output: m - a nonnegative integer which is equal to
max{ r, 0 | r in NN, deg gcd(u(q^r x), v(x)) > 0 }
*)
qDispersionRes[u_, v_, q_, x_] :=
Module[{uu, r, w, s, t, m},
If[(u /. x -> 0) === 0, Throw["the first argument is divisible by x, please enter a new one"]];
r = Resultant[u /. x -> w x, v, x];
s = Exponent[r, w, Min];
t = Exponent[Coefficient[r, w, s], q, Min];
For[m = t, m >= 0, m--,
If[Expand[r /. w -> q^m] === 0, Return[m]];
];
Return[0];
];
(*
Name: OrderBound
Input: P - a nonzero q-difference operator in QQ[q][x][S] but regarded as an element in QQ(q)[x][S], where S x = q x S
q - a variable which is not a root of unity
x - an indeterminate, which is equal to q^n, n is also a variable
S - a q-difference generator
Output: k - an order bound for a desingularized operator of P
*)
OrderBound[P_, q_, x_, S_] :=
Module[{l, t, e},
l = LeadingCoefficient[P]; (* the leading coefficient of P *)
t = Coefficient[P, S, 0]; (* the trailing coefficient of P *)
e = Exponent[l, x, Min]; (* the lowest exponent of x in l *)
If[e >= 1, l = Together[l/x^e]]; (* divide out the factor x^e in l *)
Return[qDispersion[l, t, q, x] + LeadingExponent[P][[1]]] (* return the order bound for a desingularized operator of P*)
];
(*
Name: ContSubmodule
Input: P - a nonzero q-difference operator in QQ[q][x][S] but regarded as an element in QQ(q)[x][S], where S x = q x S
k - a nonnegative integer which is greater than or equal to the order of L
x - an indeterminate
S - a q-difference generator
Output: M - a generating set of k-th submodule of Cont(P), which is equal to
{T in Cont(P) | the order of P is less than or equal to k }, where
Cont(P) = QQ(q, x)[S] P \cap QQ(q)[x][S]
*)
ContSubmodule[P_, k_, x_, S_] :=
Module[{A, a, coeff, dim, syz, M},
A = Sum[a[i]**S^i, {i, 0, k}]; (* make an ansatz of order k with unknow coefficients a[i]*)
coeff = OreReduce[A, {Together[P]}]; (*reduce A w.r.t. P *)
coeff = Numerator[OrePolynomialListCoefficients[coeff]]; (* take numerators of the nonzero coefficients in coeff *)
coeff = Table[Coefficient[coeff, a[i]], {i, 0, k}]; (* for each i, take the coefficient of a[i] *)
syz = SingularSyz[coeff, {x}]; (* compute a generating set of the syzygy module generated by coeff as a QQ(q)[x]-module *)
M = Table[ToOrePolynomial[Sum[syz[[i, j + 1]]**S^j, {j, 0, k}], OreAlgebra[S]], {i, Length[syz]}];(* compute a generating set of the k-th submodule of Cont(P) *)
Return[M];
];
(*
Name: qWeylClosure
Input: P - a nonzero q-difference operator in QQ[q][x][S] but regarded as an element in QQ(q)[x][S], where S x = q x S
q - a variable which is not a root of unity
x - an indeterminate
S - a q-dilation operator
Output: M - a generating set of the q-Weyl closure of P, which is equal to QQ(q, x)[S] P \cap QQ(q)[x][S]
*)
qWeylClosure[P_, q_, x_, S_] := ContSubmodule[P, OrderBound[P, q, x, S], x, S];
(*
Name: MinCoeffIdealGB (* compute the element in some coefficient ideal with minimal degree in x by GB *)
Input: Id - a finite set of polynomials in QQ[q, x], but regarded as elements in QQ(q)[x]
q - a variable which is not a root of unity
x - an indeterminate
Output: {g, c} - g is an element in the ideal generated by Id in QQ(q)[x] with minimial degree in x,
c is a finite set of polynomials in QQ[q, x] such that
g = c[1] Id[1] + ... + c[k] Id[k], where k is the length of Id
*)
MinCoeffIdealGB[Id_, q_, x_] :=
Module[{g, m, i},
g = SingularGroebner[Id, {x, q}, MonomialOrder -> Lexicographic]; (* compute the Groebner basis of Id w.r.t. lexicographic order x > q *)
m = MonomialList[#, {x, q}][[1]]& /@ g; (* take the leading monomials of g *)
i = 1;
For[j = 2, j <= Length[m], j++,
If[Exponent[m[[j]], x] < Exponent[m[[i]], x], i = j]; (* take the monomial in m with minimal exponent in x*)
];
c = SingularLift[Id, {g[[i]]}, {x, q}][[1]]; (* compute the cofactors such that g[[i]] = c[1] Id[1] + ... + c[k] Id[k], where k is the length of Id *)
Return[{g[[i]], c}];
];
(*
Name: MinCoeffIdealEuclidean (* compute the element in some coefficient ideal with minimal degree in x by Euclidean algorithm *)
Input: Id - a finite set of polynomials in QQ[q, x], but regarded as elements in QQ(q)[x]
x - an indeterminate
Output: {g, c} - g is an element in the ideal generated by Id in QQ(q)[x] with minimial degree in x,
c is a finite set of polynomials in QQ(q)[x] such that
g = c[1] Id[1] + ... + c[k] Id[k], where k is the length of Id
*)
MinCoeffIdealEuclidean[Id_, x_] :=
Module[{k = Length[Id], g, c, a, b},
If[k === 1, Return[{Id[[1]], 1}]];
g = Id[[1]];
c = {1};
For[i = 2, i <= k, i++,
{g, {a, b}} = PolynomialExtendedGCD[g, Id[[i]], x];
c = Append[a c, b];
];
Return[{g, c}];
];
(*
Name: LeadDesingularizedOperator
Input: P - a nonzero q-difference operator in QQ[q][x][S] but regarded as an element in QQ(q, x)[S], where S x = q S
q - a variable which is not a root of unity
x - an indeterminate
S - a q-difference generator
Output: L - a leading desingularized operator of P, which belongs to QQ[q][x][S], but regarded as an element in QQ(q)[x][S].
*)
LeadDesingularizedOperator[P_, q_, x_, S_] :=
Module[{k, M, l, c},
k = OrderBound[P, q, x, S]; (* compute an order bound for a leading desingularized operator of P *)
M = ContSubmodule[P, k, x, S]; (* compute a generating set of the k-th submodule of Cont(P) *)
l = Coefficient[#, S, k]& /@ M; (* take the coefficient of S^k for each element in M *)
c = MinCoeffIdealGB[l, q, x][[2]]; (* take the cofactors in MinCoeffIdealGB of l*)
Return[Factor[Sum[OreTimes[c[[i]], M[[i]]], {i, Length[l]}]]]; (* return a leading desingularized operator of P *)
];
(*
Name: TrailDesingularizedOperator
Input: P - a nonzero q-difference operator in QQ[q][x][S] but regarded as an element in QQ(q, x)[S], where S x = q S
q - a variable which is not a root of unity
x - an indeterminate
S - a q-difference generator
Output: T - a leading desingularized operator of P, which belongs to QQ[q][x][S], but regarded as an element in QQ(q)[x][S].
*)
TrailDesingularizedOperator[P_, q_, x_, S_] :=
Module[{k, M, t, g, m, i, c},
k = OrderBound[P, q, x, S]; (* compute an order bound for a trailing desingularized operator of P *)
M = ContSubmodule[P, k, x, S]; (* compute a generating set of the k-th submodule of Cont(P) *)
t = Coefficient[#, S, 0]& /@ M; (* take the coefficient of S^k for each element in M *)
c = MinCoeffIdealGB[t, q, x][[2]]; (* take the cofactors in MinCoeffIdealGB of t *)
Return[Factor[Sum[OreTimes[c[[i]], M[[i]]], {i, Length[t]}]]]; (* return a trailing desingularized operator of P *)
];
(*
Name: BiDesingularizedOperator
Input: P - a nonzero q-difference operator in QQ[q][x][S] but regarded as an element in QQ(q, x)[S], where S x = q S
q - a variable which is not a root of unity
x - an indeterminate
S - a q-difference generator
Output: T - a leading desingularized operator of P, which belongs to QQ[q][x][S], but regarded as an element in QQ(q)[x][S].
*)
BiDesingularizedOperator[P_, q_, x_, S_] :=
Module[{L, T, m},
L = LeadDesingularizedOperator[P, q, x, S];
T = TrailDesingularizedOperator[P, q, x, S];
If[LeadingExponent[L][[1]] > LeadingExponent[T][[1]], Return[L + T]];
m = LeadingExponent[T][[1]] - LeadingExponent[L][[1]] + 1;
Return[Factor[OreTimes[S^m, L] + T]];
];
End[]
EndPackage[]
<< RISC`HolonomicFunctions`;