(***************************************************************************) (* *) (* Algorithms qPoly and qHyper *) (* *) (* Version of July 13, 1995 *) (* programmed by Marko Petkov"sek *) (* IMFM, Ljubljana, Slovenia *) (* marko.petkovsek@mat.uni-lj.si *) (* *) (* This file contains Mathematica implementations of algorithms qPoly and *) (* qHyper which find polynomial resp. q-hypergeometric solutions of *) (* linear q-difference equations with polynomial coefficients. *) (* For details, see the paper: *) (* S.A. Abramov, M. Petkov"sek, Finding all q-hypergeometric solutions *) (* of q-difference equations, Proc. FPSAC '95, B. Leclerc and J.-Y. Thibon,*) (* Eds., Univ. Marne-la-Vall'ee, 1995, pp. 1 - 10. *) (***************************************************************************) BeginPackage["qHyper`"] q::usage = "q can be a symbol, or a number different from zero and roots of unity." Verbose::usage = "Verbose is an option of qPoly and qHyper specifying the printout level. Default value: False. Other values: Automatic, True." Solutions::usage = "Solutions is an option of qHyper specifying the number of solutions returned. Default value: 1. Other values: All." Quadratics::usage = "Quadratics is an option of qHyper specifying whether quadratic factors of the leading and trailing coeficients of the recurrence irreducible over Q should be completely factored. Default value: False. Other values: True." qHyper::usage = "qHyper[eqn, y[x]] finds at least one q-hypergeometric solution of the homogeneous equation eqn over the field of rational numbers Q (provided any such solution exists). qHyper[eqn, y[x], Solutions -> All] finds a generating set (not necessarily linearly independent) for the space of solutions generated by q-hypergeometric terms over Q. qHyper[eqn, y[x], Quadratics -> True] finds solutions over quadratic extensions of Q. Solutions y[x] are given by their rational representations y[q x]/y[x]. Warning: The worst-case time complexity of qHyper is exponential in the degrees of the leading and trailing coefficients of eqn." qPoly::usage = "qPoly[eqn, y[x]] finds a basis for the (linear or affine) space of polynomial solutions of eqn which can be nonhomogeneous." AssocEqn::usage = "AssocEqn[eqn, y[x]] constructs the equation satisfied by the power series coefficients of solutions of eqn." Begin["`Private`"] Unprotect[Power]; 0^0 = 1; Protect[Power]; Options[qHyper] = {Verbose -> False, Solutions -> 1, Quadratics -> False}; Options[qPoly] = {Verbose -> False}; qHyper[eqn_, a_[x_], opts___] := Module[{min, max, rho, polyCoeffs, m, degreePoly, s = 0, d, degrees, deg, sol, poly, sub, k, eq, degseq, lc, tc, ld, td, Alist, Blist, Zlist, seq, mm, pos, z, exps, Atempl, Btempl, lcoefs, powers, lhs, rhs, unknowns, i, degsa, degsb, expa, expb, airr, birr, aa, bb, found = {}, llc, ltc, at, bt, Z, A, B, polyeq, lca, lcb, verb, verbp, tcoefs, zeqn, mindeg, be, blc, verbose, solutions, quadratics}, verbose = Verbose /. {opts} /. Options[qHyper]; solutions = Solutions /. {opts} /. Options[qHyper]; quadratics = Quadratics /. {opts} /. Options[qHyper]; eq = Numerator[Together[If[Head[eqn] === Equal, eqn[[1]] - eqn[[2]], eqn]]]; verb = MemberQ[{True, Automatic}, verbose]; {min, max} = ({Min[#1], Max[#1]} & )[(Log[q, #1] & ) /@ (First /@ Cases[{eq}, a[_], Infinity] /. x -> 1) /. logrule]; rho = max - min; unknowns = Table[a[x*q^i], {i, 0, rho}]; eq = Factor /@ Collect[Expand[eq /. x -> x/q^min], unknowns]; rhs = -Select[eq, FreeQ[#1, a] & ]; lhs = eq + rhs; polyCoeffs = (Expand[Coefficient[lhs, #1]] & ) /@ unknowns; m = Max[(Exponent[#1, x] & ) /@ polyCoeffs]; lc = Last[polyCoeffs]; tc = First[polyCoeffs]; degseq = (Exponent[#1, x] & ) /@ polyCoeffs; ld = Last[degseq]; td = First[degseq]; powers = Table[z^i, {i, 0, rho}]; lcoefs = (Coeff[#1[[1]], x, #1[[2]]] & ) /@ Thread[{polyCoeffs, degseq}]; llc = Last[lcoefs]; ltc = First[lcoefs]; If[verbose, Print[""]; Print["order = ", rho]; Print["lc = ", lc]; Print["tc = ", tc]; Print["ld = ", ld]; Print["td = ", td]; Print["degree sequence = ", degseq]; Print["leading coeffs = ", lcoefs] ]; For[d = -ld, d <= td, d++, seq = degseq + Table[i*d, {i, 0, rho}]; mm = Max[seq]; If[!FreeQ[Blist[], Blist], bl = If[quadratics === True, FactorList2[lc, x], FactorList[lc]]; Blist[] = Select[bl, !FreeQ[#1, x] & ]; If[Blist[] === {}, Blist[] = {{1, 0}}]; birr = First /@ Blist[]; degsb = (Exponent[#1, x] & ) /@ birr; If[Length[Complement[Union[degsb], {0, 1}]] > 0, Print["\n Warning: irreducible factors of degree > 1 in leading \ coefficient;\n some solutions may not be found"]]; expb = Last /@ Blist[]; lcb = (Coefficient[#1[[1]], x, #1[[2]]] & ) /@ Thread[{birr, degsb}]; birr = Expand[birr/lcb]; Blist[] = Thread[{birr, expb}]; Btempl = Bags[expb, degsb]; If[verbose, Print[""]; Print[" Blist = ", Blist[]]]]; If[!FreeQ[Alist[], Alist], al = If[quadratics === True, FactorList2[tc, x], FactorList[tc]]; Alist[] = Select[al, !FreeQ[#1, x] & ]; If[Alist[] === {}, Alist[] = {{1, 0}}]; airr = First /@ Alist[]; degsa = (Exponent[#1, x] & ) /@ airr; If[Length[Complement[Union[degsa], {0, 1}]] > 0, Print["\n Warning: irreducible factors of degree > 1 in trailing \ coefficient;\n some solutions may not be found"]]; expa = Last /@ Alist[]; lca = (Coefficient[#1[[1]], x, #1[[2]]] & ) /@ Thread[{airr, degsa}]; airr = Expand[airr/lca]; Alist[] = Thread[{airr, expa}]; Atempl = Bags[expa, degsa]; If[verbose, Print[" Alist = ", Alist[]]]]; For[da = Max[0, d], da <= Min[td, ld + d], da++, db = da - d; If[verb, Print[""]; Print[" da = ", da, ", ", "db = ", db]]; at = Atempl[[da + 1]]; bt = Btempl[[db + 1]]; For[bb = 1, bb <= Length[bt], bb++, B = Apply[Times, birr^bt[[bb]]] /. x -> x*q^(1 - rho); be = Collect[Expand[B], x]; blc = Coefficient[be, x, Exponent[be, x]]; B = Expand[B/blc]; bc = blc*Apply[Times, birr^(expb - bt[[bb]])]; For[aa = 1, aa <= Length[at], aa++, A = Apply[Times, airr^at[[aa]]]; ac = Apply[Times, airr^(expa - at[[aa]])]; polyeq = Module[{i, j, k}, Join[{ltc*ac*Product[B /. x -> x*q^k, {k, 0, rho - 2}]}, Table[polyCoeffs[[i + 1]]* Product[A /. x -> x*q^j, {j, i - 1}]* Product[B /. x -> x*q^k, {k, i, rho - 2}]*z^i, {i, rho - 1}], {llc*bc*Product[A /. x -> x*q^j, {j, rho - 1}]*z^rho}]]; polyeq = Expand /@ polyeq; mindeg = Min[(Exponent[#1, x, Min] & ) /@ Collect[polyeq, x]]; tcoefs = (Coeff[#1, x, mindeg] & ) /@ polyeq; polyeq = (Expand[#1/x^mindeg] & ) /@ polyeq; zeqn = Apply[Plus, tcoefs] == 0; sub = Solve[zeqn, z]; If[!FreeQ[sub, ToRules], Print["\n Warning: algebraic roots of degree > 4 in characteristic \ polynomial;\n some solutions may not be found"]]; sub = Select[sub, FreeQ[#1, ToRules] & ]; Zlist = Union[If[MemberQ[{{},{{}}}, sub], {}, Together /@ (z /. sub)]]; Zlist = Complement[Zlist, {0}]; If[verbose, Print[""]; Print["Zlist = ", Zlist]]; For[ii = 1, ii <= Length[Zlist], ii++, Z = Zlist[[ii]]; If[verb, Print[""]; Print[" A = ", A, " B = ", B, " Z = ", Z]]; verbp = Switch[verbose, True, Automatic, _, False]; sol = qPolyK[polyeq /. z -> Z, 0, a[x], verbp]; If[sol === {}, If[verbose, Print[""]; Print[" No non-zero polynomial solution."]], If[verb, Print[""]; Print[" ===> C = ", sol]]; sol = Module[{i}, Table[Coefficient[sol, C[i]], {i, Count[{sol}, C[_], Infinity]}]]; sol = (sol /. x -> q*x)/sol; sol = Factor /@ ((Z*A*sol)/B); If[solutions === 1, Return[sol]]; found = Union[found, sol]]]]]]]; Return[found]] logrule = {Log[q^(n_)]/Log[q] -> n} Coeff[_, _, _?Negative] = 0 Coeff[p_, n_, i_] := Coefficient[p, n, i] FactorList2[p_, x_] := Module[{pp, lc, lcoefs, kvadr}, pp = FactorList[p]; lc = Select[pp, FreeQ[#1, x] & ]; If[lc === {}, lc = {{1, 1}}]; pp = Complement[pp, lc]; lcoefs = ({LC[First[#1], x], Last[#1]} & ) /@ pp; lc = {Apply[Times, (First[#1]^Last[#1] & ) /@ Join[lc, lcoefs]], 1}; pp = Thread[{Expand[First /@ pp/First /@ lcoefs], Last /@ pp}]; kvadr = Select[pp, Exponent[First[#1], x] === 2 & ]; pp = Complement[pp, kvadr]; kvadr = Apply[Join, (F2[#1, x] & ) /@ kvadr]; Return[Join[{lc}, pp, kvadr]]] LC[p_, x_] := Coefficient[p, x, Exponent[p, x]] F2[{p_, d_}, x_] := ({#1, d} & ) /@ (x - (x /. Solve[p == 0, x])) Bags[{}, {}] = Bags[{0}, {0}] Bags[exps_, degs_] := Module[{tuples, values, i}, tuples = Flatten[Apply[Outer, Join[{List}, (Range[0, #1] & ) /@ exps]], Length[exps] - 1]; values = (Apply[Plus, #1 . degs] & ) /@ tuples; Return[Table[tuples[[Apply[Join, Position[values, i]]]], {i, 0, exps . degs}]]] qPolyK[polyCoeffs_, rhs_, (a_)[x_], verbose_] := Module[{rho = Length[polyCoeffs] - 1, d, degreePoly, s = 0, dd, degrees, deg, unknowns, lhs, i, sol, poly, sub, k, cc, genh, part, verb = MemberQ[{True, Automatic}, verbose]}, d = Max[(Exponent[#1, x] & ) /@ polyCoeffs]; degreePoly = Sum[Coef[polyCoeffs, i + 1, x, d]*x^i, {i, 0, rho}]; If[verb, Print[""]; Print["Degree polynomial: ", degreePoly]]; degrees = Select[(Log[q, #1] & ) /@ (* Flatten[Union[Factor /@ (x /. Solve[degreePoly == 0, x]) /. *) Flatten[Union[Factor /@ (x /. Select[{ToRules[Roots[degreePoly == 0, x, Cubics -> False, Quartics -> False]]}, FreeQ[#, Roots]&]) /. x -> {}]] /. logrule, FreeQ[#1, Log] & ]; If[verb, Print["Roots found: ", degrees]]; degrees = Append[degrees, Exponent[rhs, x] - d]; degrees = Select[degrees, IntegerQ[#1] && NonNegative[#1] & ]; If[verb, Print["Possible degrees: ", degrees]]; deg = Max[degrees]; If[deg === -Infinity, Return[{}]]; unknowns = Array[C, deg + 1, 0]; poly = unknowns . Table[x^k, {k, 0, deg}]; lhs = polyCoeffs . Table[a[x*q^i], {i, 0, rho}]; sub = lhs /. a[u_] :> (poly /. x -> u); sub = Reverse[CoefficientList[Expand[sub - rhs], x]]; If[verbose, Print["Triangular system: "]; Print[""]; Print[MatrixForm[Thread[sub == 0]]]]; sub = TriSolve[sub, unknowns]; If[sub === Fail, Return[{}], sol = Collect[poly /. sub, unknowns]]; If[sol === 0, Return[{}]]; sol = MakeList[sol, Plus]; part = Select[sol, FreeQ[#1, C] & ]; genh = Complement[sol, part]; sol = Factor[Apply[Plus, part]] + Apply[Plus, (Factor[Numerator[Together[#1]]] & ) /@ genh]; cc = Cases[{sol}, C[_], Infinity]; Return[sol /. Thread[cc -> Array[C, Length[cc]]]]] Coef[_, _, _, _?Negative] = 0 Coef[p_, i_, n_, j_] := Coefficient[p[[i]], n, j] TriSolve[eqns_, unknowns_] := Module[{sub = {}, sub1, x, eq, eqn, term, un}, sub1 = Scan[(eq = Collect[Numerator[Together[#1 /. sub]], unknowns]; If[eq == 0, Null, Return[Fail], If[FreeL[eq, unknowns], Return[Fail], eqn = MakeList[eq, Plus]; term = Scan[If[Apply[Or, Function[x, !FreeQ[#1, x]] /@ unknowns], Return[#1]] & , eqn]; un = First[Select[MakeList[term, Times], MemberQ[unknowns, #1] & ]]; sub1 = {un -> (term - eq)/(term/un)}; sub = Join[sub /. sub1, sub1]]]) & , eqns]; Return[If[sub1 === Fail, Fail, sub]]] FreeL[expr_, l_] := Apply[And, (FreeQ[expr, #1] & ) /@ l] MakeList[expr_, head_:List] := If[Head[expr] === head, Apply[List, expr], {expr}] qPoly[eqn_, a_[x_], opts___] := Module[{min, max, rho, polyCoeffs, i, eq, lhs, rhs, verbose}, verbose = Verbose /. {opts} /. Options[qPoly]; eq = Numerator[Together[If[Head[eqn] === Equal, eqn[[1]] - eqn[[2]], eqn]]]; {min, max} = ({Min[#1], Max[#1]} & )[(Log[q, #1] & ) /@ (First /@ Cases[{eq}, a[_], Infinity] /. x -> 1) /. logrule]; rho = max - min; eq = Factor /@ Collect[Expand[eq /. x -> x/q^min], Table[a[x*q^i], {i, 0, rho}]]; rhs = -Select[eq, FreeQ[#1, a] & ]; lhs = eq + rhs; polyCoeffs = Module[{i}, (Coefficient[lhs, #1] & ) /@ Table[a[x*q^i], {i, 0, rho}]]; Return[qPolyK[polyCoeffs, rhs, a[x], verbose]]] AssocEqn[eqn_, (a_)[x_]] := Module[{min, max, rho, d, polyCoeffs, i, eq, lhs, rhs}, eq = Numerator[Together[If[Head[eqn] === Equal, eqn[[1]] - eqn[[2]], eqn]]]; {min, max} = ({Min[#1], Max[#1]} & )[(Log[q, #1] & ) /@ (First /@ Cases[{eq}, a[_], Infinity] /. x -> 1) /. logrule]; rho = max - min; eq = Factor /@ Collect[Expand[eq /. x -> x/q^min], Table[a[x*q^i], {i, 0, rho}]]; rhs = -Select[eq, FreeQ[#1, a] & ]; lhs = eq + rhs; polyCoeffs = Module[{i}, (Coefficient[lhs, #1] & ) /@ Table[a[x*q^i], {i, 0, rho}]]; d = Max[(Exponent[#1, x] & ) /@ polyCoeffs]; Return[Sum[a[x*q^s]*Sum[Coef[polyCoeffs, i + 1, x, d - s]*q^(i*s)*x^i, {i, 0, rho}], {s, 0, d}]]]; End[] EndPackage[]