(*The package "qDesingularization.m" needs two packages "HolonomicFunctions.m" and "Singular.m"*)
(*"HolonomicFunctions.m" can be downloaded from: 
http://www.risc.jku.at/research/combinat/software/ergosum/installation.html#download 
*)
(*"Singular.m" can be downloaded from: 
http://www.risc.jku.at/research/combinat/software/Singular/index.html
You also need to put "Singular.m" in the same directory as "qDesingularization.m" *)

<<qDesingularization.m;

HolonomicFunctions Package version 1.7.3 (21-Mar-2017)
written by Christoph Koutschan
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

rec = q^(1+n)*(1+q^(1+n))*(-1+q^(1+2*n))*(-1+q^(3+2*n))-q^(2+2*n)*(-1+q^n)*(-1+q^(3+2*n))*f[n]+(-1+q^(1+n))^2*(1+q^(1+n))*(1-q^(1+n)-q^(1+2*n)-q^(3+2*n)-q^(3+3*n)+q^(4+4*n))*f[1+n]-q^(2+2*n)*(-1+q^(2+n))*(-1+q^(1+2*n))*f[2+n]; (* rec is the q-difference equation in "rec.twist.knot.-1.m" *)

rec1 = HomNormRe[rec, f, q] (*take the homogeneous part of rec and make a normalization f[a_] -> f[a]/(q^a - 1) *)

P = ReToOp[rec1, f, n, QS[x, q^n]]

k=OrderBound[P,q,x,QS[x,q^n]] (* compute an order bound for a leading (trailing) desingulaized operator of P*)

M=ContSubmodule[P,k,x,QS[x,q^n]]; (* compute a generating set of the k-th submodul of Cont(P), which is the q-Weyl closure of P *)

L=LeadDesingularizedOperator[P,q,x,QS[x,q^n]] (* compute a leading desingularized operator of P*)

T=TrailDesingularizedOperator[P,q,x,QS[x,q^n]] (* compute a trailing desingularized operator of P *)

B=BiDesingularizedOperator[P,q,x,QS[x,q^n]] (* compute a bi-desingularized operator of P*)

g=InhomPart[rec,f] (*take the inhomogeneous part of rec*)

rec1=OpToRe[P,L,f,n,g] (* tranform L into a inhomogeneous q-difference equation, whose homogeneous part is L(f(n)), inhomogeneous part is (L/P)(g)*)

rec2=OpToRe[P,T,f,n,g] (* tranform T into a inhomogeneous q-difference equation, whose homogeneous part is T(f(n)), inhomogeneous part is (T/P)(g)*)

rec3=OpToRe[P,B,f,n,g] (* tranform B into a inhomogeneous q-difference equation, whose homogeneous part is B(f(n)), inhomogeneous part is (B/P)(g)*)