(**************************************************************************) (* *) (* Algorithms Poly and Hyper *) (* *) (* Version of July 31, 1995 *) (* programmed by Marko Petkov"sek *) (* IMFM, Ljubljana, Slovenia *) (* marko.petkovsek@mat.uni-lj.si *) (* *) (* This file contains Mathematica implementations of algorithms Poly and *) (* Hyper which find polynomial resp. hypergeometric solutions of linear *) (* recurrences with polynomial coefficients. For details, see the paper *) (* *) (* M. Petkov"sek: Hypergeometric solutions of linear recurrences *) (* with polynomial coefficients, J. Symb. Comp. 14 (1992) 243 - 264. *) (**************************************************************************) BeginPackage["Hyper`"] Verbose::usage = "Verbose is an option of Poly and Hyper specifying the printout level. Default value: False. Other values: Automatic, True." Solutions::usage = "Solutions is an option of Hyper specifying the number of solutions returned. Default value: 1. Other values: All." Quadratics::usage = "Quadratics is an option of Hyper 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." Poly::usage = "Poly[eqn, y[n]] finds a basis for the (linear or affine) space of polynomial solutions of eqn which can be nonhomogeneous." Hyper::usage = "Hyper[eqn, y[n]] finds at least one hypergeometric solution of the homogeneous equation eqn over the field of rational numbers Q (provided any such solution exists). Hyper[eqn, y[n], Solutions -> All] finds a generating set (not necessarily linearly independent) for the space of solutions generated by hypergeometric terms over Q. Hyper[eqn, y[n], Quadratics -> True] finds solutions over quadratic extensions of Q. Solutions y[n] are given by their rational representations y[n+1]/y[n]. Warning: The worst-case time complexity of Hyper is exponential in the degrees of the leading and trailing coefficients of eqn." Begin["`Private`"] Unprotect[Power]; 0^0 = 1; Protect[Power]; Options[Hyper] = {Verbose -> False, Solutions -> 1, Quadratics -> False}; Options[Poly] = {Verbose -> False}; Poly[eqn_, a_[n_], opts___] := Module[{min, max, dd, polyCoeffs, i, eq, lhs, rhs, verbose}, verbose = Verbose /. {opts} /. Options[Poly]; eq = Numerator[Together[ If[Head[eqn] === Equal, eqn[[1]]-eqn[[2]], eqn] ]]; {min, max} = {Min[#], Max[#]}& @ (First /@ Cases[{eq}, a[_], Infinity]/.n->0); dd = max - min; eq = Collect[Expand[eq /. n -> n - min], Table[a[n+i], {i,0,dd}]]; rhs = - Select[eq, FreeQ[#, a]&]; lhs = eq + rhs; polyCoeffs = Module[{i}, Coefficient[lhs, #]& /@ Table[a[n+i], {i,0,dd}]]; Return[PolyK[polyCoeffs, rhs, a[n], verbose]] ]; Coeff[_, _, _?Negative] = 0 Coeff[p_, n_, i_] := Coefficient[p, n, i] 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}]]] PolyK[polyCoeffs_, rhs_, (a_)[n_], verbose_:False] := Module[{dd = Length[polyCoeffs] - 1, m, degreePoly, s = 0, d, degrees, deg, unknowns, lhs, i, sol, poly, sub, k, cc, genh, part, verb = MemberQ[{True, Automatic}, verbose]}, m = Max[(Exponent[#1, n] & ) /@ polyCoeffs]; For[degreePoly = 0, degreePoly === 0, s++, degreePoly = Module[{i, j}, Factor[Sum[Binomial[d, j]* Sum[i^j*Coef[polyCoeffs, i + 1, n, m - s + j], {i, 0, dd}], {j, 0, s}]]]]; s--; If[verb, Print[""]; Print["Degree polynomial: ", degreePoly /. d -> n]]; degrees = Flatten[Union[Factor /@ (d /. Solve[degreePoly == 0, d]) /. d -> {}]]; If[verb, Print["Roots found: ", degrees]]; degrees = Append[degrees, Exponent[rhs, n] - m + s]; 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[n^k, {k, 0, deg}]; lhs = polyCoeffs . Table[a[n + i], {i, 0, dd}]; sub = lhs /. a[x_] :> (poly /. n -> x); sub = Reverse[CoefficientList[Expand[sub - rhs], n]]; 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, #] & ) /@ l] MakeList[expr_, head_:List] := If[Head[expr] === head, Apply[List, expr], {expr}] Hyper[eqn_, a_[n_], opts___] := Module[{min, max, dd, 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, verbose, solutions, quadratics}, verbose = Verbose /. {opts} /. Options[Hyper]; solutions = Solutions /. {opts} /. Options[Hyper]; quadratics = Quadratics /. {opts} /. Options[Hyper]; eq = Numerator[Together[If[Head[eqn] === Equal, eqn[[1]] - eqn[[2]], eqn]]]; verb = MemberQ[{True, Automatic}, verbose]; {min, max} = ({Min[#1], Max[#1]} & )[First /@ Cases[eq, a[_], Infinity] /. n -> 0]; dd = max - min; unknowns = Table[a[n + i], {i, 0, dd}]; eq = Collect[eq, unknowns]; eq = Factor /@ (eq /. n -> n - min); rhs = -Select[eq, FreeQ[#1, a] & ]; lhs = eq + rhs; polyCoeffs = (Expand[Coefficient[lhs, #1]] & ) /@ unknowns; m = Max[(Exponent[#1, n] & ) /@ polyCoeffs]; lc = Last[polyCoeffs]; tc = First[polyCoeffs]; degseq = (Exponent[#1, n] & ) /@ polyCoeffs; ld = Last[degseq]; td = First[degseq]; powers = Table[z^i, {i, 0, dd}]; lcoefs = (Coeff[#1[[1]], n, #1[[2]]] & ) /@ Thread[{polyCoeffs, degseq}]; llc = Last[lcoefs]; ltc = First[lcoefs]; If[verbose, Print[""]; Print["order = ", dd]; 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, dd}]; mm = Max[seq]; pos = First /@ Position[seq, mm]; If[verb, Print[""]; Print["d = ", d, "; seq = ", seq, "; pos = ", pos]]; If[Length[pos] <= 1, Null, zeqn = lcoefs[[pos]] . powers[[pos]] == 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[sub === {}, {}, z /. sub]]; Zlist = Complement[Zlist, {0}]; If[verbose, Print[""]; Print["Zlist = ", Zlist]]; For[ii = 1, ii <= Length[Zlist], ii++, Z = Zlist[[ii]]; If[!FreeQ[Blist[], Blist], bl = If[quadratics === True, FactorList2[lc /. n -> n - dd + 1, n], FactorList[lc /. n -> n - dd + 1]]; Blist[] = Select[bl, !FreeQ[#1, n] & ]; If[Blist[] === {}, Blist[] = {{1, 0}}]; birr = First /@ Blist[]; degsb = (Exponent[#1, n] & ) /@ 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]], n, #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, n], FactorList[tc]]; Alist[] = Select[al, !FreeQ[#1, n] & ]; If[Alist[] === {}, Alist[] = {{1, 0}}]; airr = First /@ Alist[]; degsa = (Exponent[#1, n] & ) /@ 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]], n, #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]]]; bc = Apply[Times, birr^(expb - bt[[bb]])] /. n -> n + dd - 1; For[aa = 1, aa <= Length[at], aa++, A = Apply[Times, airr^at[[aa]]]; ac = Apply[Times, airr^(expa - at[[aa]])]; If[verb, Print[""]; Print[" A = ", A, " B = ", B, " Z = ", Z]]; polyeq = Module[{i, j, k}, Join[{ltc*ac*Product[B /. n -> n + k, {k, 0, dd - 2}]}, Table[polyCoeffs[[i + 1]]* Product[A /. n -> n + j, {j, i - 1}]* Product[B /. n -> n + k, {k, i, dd - 2}]*Z^i, {i, dd - 1}], {llc*bc*Product[A /. n -> n + j, {j, dd - 1}]*Z^dd}]]; polyeq = Expand /@ polyeq; If[verbose, Print[""]; Print[" ", polyeq]]; verbp = Switch[verbose, True, Automatic, _, False]; sol = PolyK[polyeq, 0, a[n], verbp]; If[sol === {}, If[verbose, Print[" No non-zero polynomial \ solution."]; Print[""]], If[verb, Print[""]; Print[" ===> C = ", sol]]; sol = Module[{i}, Table[Coefficient[sol, C[i]], {i, Count[{sol}, C[_], Infinity]}]]; sol = (sol /. n -> n + 1)/sol; sol = Factor /@ ((Z*A*sol)/B); If[solutions === 1, Return[sol]]; found = Union[found, sol]]]]]]]]; Return[found]] 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])) F2[p_, x_] := Apply[Times, x - (x /. Solve[p == 0, x])] End[] EndPackage[]