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`;
