LFS := module() export EG, HypergeometricSolution, RationalSolution, Sdetermine, IsQHypergeometricTerm, RNFQShift, qSolve, Solve, dualnormal; local proper_expr, proper_name, SetQDifferenceRing, SetQIsHypergeometricTerm, SetQRNF, SetQInitPoint, InitPointQShift, InhomogeneousToHomogeneous, UniversalDenominatorForHomogeneous, UniversalDenominatorForHomogeneousSystemInNormalForm, ParticularRationalSolution, PartRatSolForSystemInNormalForm, RationalSolutionForInhomogeneous, RationalSolutionForHomogeneous, DegreeBoundForHomogeneous, PartSolForSystemInNormalForm, ParticularHypergeometricSolution, HypergeometricSubstitution, HypergeometricSubstitutionToSystemInNormalForm, InducedRecurrenceForHomogeneous, qDENOM, EG_, Shift_Column, Shift_Row, Row_Length, cases, HGNormal, HGNormal_step1, check_unknown, what_is_case, qratio, qratio_prod, SubsPower, default_dualnormal, SolveRec, qdifferenceToOrePoly, _c, terminate_search_for_rational; cases := ['differential', 'difference', 'qdifference']; proper_name := proc(x, excld := {}) local exclude, tmp; exclude := {'false', infinity, 'true', FAIL, 'undefined', excld[]}; tmp := select(type, exclude, 'indexed'); while tmp <> {} do tmp := map2(op, 0, tmp); exclude := exclude union tmp; tmp := select(type, tmp, 'indexed'); end do; x::'scalar' and not has(x, exclude); end proc; proper_expr := proc(x, excld := {}) local exclude, tmp; exclude := excld; tmp := select(type, exclude, 'indexed'); while tmp <> {} do tmp := map2(op, 0, tmp); exclude := exclude union tmp; tmp := select(type, tmp, 'indexed'); end do; x::'scalar' and not has(x, exclude); end proc; SetQDifferenceRing := proc(X, q) local p, x, qq; qq := convert(q, 'string'); if qq = "x" then x := "xx"; elif qq = "p" then p := "pp"; end if; OreTools:-SetOreRing(X, 'qdifference', 'sigma' = # (p,x)->eval(p,x=x*q) parse(sprintf("(%s,%s)->eval(%s,%s=%s*(%a))", p, x, p, x, x, q)), 'sigma_inverse' = # (p,x)->eval(p,x=x/q) parse(sprintf("(%s,%s)->eval(%s,%s=%s/(%a))", p, x, p, x, x, q)), 'delta' = proc (p, x) 0 end proc, 'theta1' = 1) end proc: SetQIsHypergeometricTerm := proc(q, n) local t, x, R, strs; if not q::{'name','posint'} then return parse(sprintf( "proc() \\ error \"The case 'qdifference' with base %a is not determined\"\\ end proc", q)); end if; strs := {convert(q, 'string'), convert(n, 'string')}; if "x" in strs then if "x1" in strs then x := "x2"; else x := "x1"; end if; end if; if "t" in strs then if "t1" in strs then t := "t2"; else t := "t1"; end if; end if; if "R" in strs then if "R1" in strs then R := "R2"; else R := "R1"; end if; end if; # (t, x1, R1)->LFS:-IsQHypergeometricTerm( # eval(t, x = q^n), n, q^n = x, 'R') parse(sprintf( "(%s, %s, %s)->\\ LFS:-IsQHypergeometricTerm(\\ eval(%s, %s = %a^%a), %a, %a^%a = %s, %s)", t, x, R, t, x, q, n, n, q, n, x, R)); end proc: SetQRNF := proc(q) local R, x, strs; if not q::'name' then return parse(sprintf( "proc(R, x) local y, sol; \\ sol := QDifferenceEquations:-RationalSolution(\\ y(%a*x)-R*y(x), y(x), output = basis);\\ if sol = NULL then return 1, R; \\ else return sol[1], 1 end if \\ end proc", q)); end if; strs := convert(q, 'string'); if "x" = strs then x := "x1"; elif "R" = strs then R := "R1"; end if; # (R, x) -> LFS:-RNFQShift(R, x, q) parse(sprintf( "(%s, %s) -> LFS:-RNFQShift(%s, %s, %a)", R, x, R, x, q)); end proc: SetQInitPoint := proc(q, n) local x, R, strs; strs := {convert(q, 'string'), convert(n, 'string')}; if "x" in strs then if "x1" in strs then x := "x2"; else x := "x1"; end if; end if; if "R" in strs then if "R1" in strs then R := "R2"; else R := "R1"; end if; end if; # (R, x) -> LFS:-InitPointQShift(R, x, q, n) parse(sprintf( "(%s, %s)-> (%a)(%s, %s, %a, %a)", R, x, eval(LFS:-InitPointQShift), R, x, q, n)); end proc: ################################################################### # Unified normal both for rational and algebraic numbers ################################################################### default_dualnormal := e->`if`(has(e,RootOf),evala(Normal(e)),normal(e)); dualnormal := default_dualnormal; ############################################################################### # Input list/set form: # Sys -- q-difference system # vars -- list of unknown functions, [y1(x),y2(x)] for example # Input (matrix-form): # Sys -- linear equation with matrix coefficients # A_n . y(q^n x) + ... A_1 . y(q x) + A_0 . y(x) = 0 # y(x) -- unknown function ############################################################################### ############################################################################### EG := proc(lt) local pp, number_of_equations, number_of_unknowns, EM, k, sys, cfs, i; if nargs = 0 or not lt in {'lead', 'trail'} then error "The %-1 argument ('lead' or 'trail') is expected", 1; end if; try pp := Sdetermine(args[2..-1]) catch: if nops([lastexception]) > 2 and searchtext("%-1", lastexception[2]) > 0 then error lastexception[2], lastexception[3]+1, lastexception[4..-1]; else error; end if; end try; if not pp['homogeneous'] then error "A homogeneous system is expected"; end if; cfs := [entries(pp['explicit_form'], 'nolist')]; cfs := {map(op@convert, cfs, 'set')[]}; if not cfs::'set'('ratpoly'('algebraic', pp['var'])) then error "The system's coefficients of type ratpoly in %1 are expected", pp['var']; end if; dualnormal := default_dualnormal; number_of_equations := LinearAlgebra:-RowDimension(pp['explicit_form'][0]); number_of_unknowns := LinearAlgebra:-ColumnDimension(pp['explicit_form'][0]); EM := Matrix([seq(pp['explicit_form'][pp['order']-i], i=0..pp['order'])]); EG_(lt, EM, pp['order'], number_of_unknowns, number_of_equations, pp['var'], pp['Ring']); if number_of_unknowns = 1 and number_of_equations = 1 then sys := EM[1, pp['order']+1]*pp['functions'][1] + add(EM[1, pp['order']-k+1]* OreTools:-Apply(OrePoly(0$k, 1), pp['functions'][1], pp['Ring']), k = 1..pp['order']); if convert(EM[1], 'set') <> {0} then userinfo(3, LFS, "The system is of full rank"); return sys=0, true; else userinfo(3, LFS, "The system is not of full rank"); return sys=0, false; end if; elif number_of_equations >= number_of_unknowns and convert(EM[number_of_unknowns], 'set') <> {0} then userinfo(3, LFS, "The system is of full rank"); sys := add(EM[1..number_of_unknowns, k*number_of_unknowns+1..(k+1)*number_of_unknowns]. OreTools:-Apply(OrePoly(0$(pp['order']-k), 1), op([0,0],pp['functions'][1])(pp['var']), pp['Ring']), k=0..pp['order']-1) + EM[1..number_of_unknowns, -number_of_unknowns .. -1]. op([0,0],pp['functions'][1])(pp['var']); return sys=0, true; else sys := add(EM[1..-1, k*number_of_unknowns+1..(k+1)*number_of_unknowns]. OreTools:-Apply(OrePoly(0$(pp['order']-k), 1), op([0,0],pp['functions'][1])(pp['var']), pp['Ring']), k=0..pp['order']-1) + EM[1..-1, -number_of_unknowns .. -1]. op([0,0],pp['functions'][1])(pp['var']); userinfo(3, LFS, "The system is not of full rank"); return sys=0, false; end if; end proc: ###################################################################### # HypergeometricSolution # Output: hypergeometric solutions for the given system ###################################################################### HypergeometricSolution := proc() local opts, Args, Field, define_field, outpt, pp, p1, j, Res, Res1, cfs, indts, dummy, x, norm; if membertype(`=`({identical('rhs_normal'), identical('Normal'), identical('output'), identical('constant_field'), identical('earlyterminate')}, 'anything'), [args], 'j') then Args := [[args][1..j-1][], 'Dummy' = true, [args][j..-1][]]; else Args := [args]; end if; opts, Args := selectremove(el -> el::`=`({identical('rhs_normal'), identical('Normal'), identical('output'), identical('constant_field'), identical('earlyterminate')}, 'anything'), Args); #terminate_search_for_rational := false; terminate_search_for_rational := true; hasoption(opts, 'earlyterminate'=truefalse, 'terminate_search_for_rational', 'opts'); try pp := Sdetermine(Args[]); catch: error; end try; x := pp['var']; cfs := [entries(pp['explicit_form'], 'nolist')]; cfs := {map(op@convert, cfs, 'set')[]}; indts := remove(type, cfs, 'ratpoly'('anything', x)); if indts <> {} then error "The system's coefficients of type ratpoly in %1 are expected but '%2' is", x, indts[1]; end if; define_field := false; indts := indets(cfs, 'name') minus {x, constants}; if assigned(pp['q_par']) and pp['q_par']::'name' then indts := select(proper_name, indts, {x, pp['q_par']}); indts := indts union {pp['q_par']}; else indts := select(proper_name, indts, {x}); end if; if not hasoption(opts, 'constant_field', 'Field', 'opts') then define_field := true; if indts = {} then Field := 'rational'; else Field := 'ratpoly'('rational', indts); end if; end if; try cfs := remove(type, cfs, 'ratpoly'(Field, x)); catch: error cat("Impossible to check system's coefficients: ", [lastexception][2]), [lastexception][3..-1][]; end try; if cfs <> {} then error "The system's coefficients of type ratpoly in %1 over the field '%2' are expected but '%3' is", x, Field, cfs[1]; end if; if pp['order'] = 0 then try Res := LinearAlgebra:-LinearSolve(pp['explicit_form'][0], pp['RHS']); catch: userinfo(3, 'LFS', StringTools[FormatMessage](lastexception[2..-1])); userinfo(3, 'LFS', "There is no rational solution"); return NULL; end try; userinfo(3, 'LFS', "The solution is found"); return Res; end if; if not pp['homogeneous'] then Args, opts := selectremove( el -> el::`=`({identical('Normal')}, 'anything'), opts); try Res := HGNormal(pp['RHS'], pp['Ring'], Args[], 'define_cf' = define_field, 'constant_field'=Field, 'impossible_exponent'=indts); catch "The variable n", "The case": error; catch: j := cat("The right-hand side normalizing: ", StringTools[FormatMessage](lastexception[2..-1])); error j; end try; if nops([Res]) = 2 then if define_field and op(0, Res[2]) = 'ratpoly' then indts := indts union op(2, Res[2]); end if; pp['RHS_normal'] := Res[1]; else pp['RHS_normal'] := Res; end if; end if; if not pp['homogeneous'] then if pp['RHS_normal']::'Vector[column]' and convert(pp['RHS_normal'], 'set') = {0} then pp['homogeneous'] := true; userinfo(3, LFS, "The right-hand side is normalized"); userinfo(3, LFS, "The system is homogeneous"); elif nops(pp[RHS_normal]) = 1 then userinfo(3, 'LFS', "The right-hand side has one hypergeometric term"); else userinfo(3, 'LFS', "The right-hand side has %1 non-similar hypergeometric terms", nops(pp[RHS_normal])); end if; end if; if define_field then if indts = {} then userinfo(3, 'LFS', "The system's constant field is the rational number field"); else userinfo(3, 'LFS', "The system's constant field is the rational number field extended by %1", indts); end if; end if; if pp['homogeneous'] then error "The homogeneous case is not implemented yet, sorry"; end if; hasoption(opts, 'output', 'outpt', 'opts'); if outpt <> 'partsol' then error "Only the particular solution case is implemented now, sorry"; elif opts <> [] then error "Invalid optional arguments; first one is '%1'", opts[1]; end if; Res := 0; for j to nops(pp['RHS_normal']) do if convert(pp['RHS_normal'][j][3], 'set') = {0} then next; end if; if pp['source_system']::'Matrix' then try Res1 := PartSolForSystemInNormalForm(pp, j); catch: error; end try; else try Res1 := ParticularHypergeometricSolution(pp, j); catch "The system is inconsistent or not full rank": userinfo(3, LFS, "The system is inconsistent or not full rank"); return NULL; catch: error; end try; end if; if Res1 = NULL then userinfo(3, 'LFS', "There is no particular solution"); return NULL; end if; Res1 := Res1, selectremove(type, pp['RHS_normal'][j][1]*dummy[1](x)*dummy[2](x), 'ratpoly'('anything', x)); if Res1[2] <> 1 then Res := Res + map(`*`, map(dualnormal, map(`*`, Res1[1], Res1[2])), Res1[3]/dummy[1](x)/dummy[2](x)); else Res := Res + map(`*`, Res1[1], Res1[3]/dummy[1](x)/dummy[2](x)); end if; end do; if convert(Res, 'set') = {0} then error "The homogeneous case is not implemented yet, sorry"; end if; userinfo(3, 'LFS', "The particular hypergeometric solution is found"); return Res; end proc: ###################################################################### # ParticularHypergeometricSolution # Output: particular solution for a system with any order which # right hand side is a hypergeometric term ###################################################################### ParticularHypergeometricSolution := proc(pp, j) local S, R::nothing, y::nothing; if pp['case'] = 'differential' and pp['RHS_normal'][j][2] = 0 or pp['case'] in {'difference', 'qdifference'} and pp['RHS_normal'][j][2] = 1 then userinfo(3, LFS, "%1. The case of a rational right hand side", j); S := copy(pp, 'deep'); elif pp['case'] in {'differential', 'difference', 'qdifference'} then S := HypergeometricSubstitution(pp, 1, pp['RHS_normal'][j][2]); S := copy(S, 'deep'); userinfo(3, LFS, sprintf("%a. Substitute %a to the system", j, y(pp['var'])=pp['RHS_normal'][j][1]*R(pp['var']))); else error "The case '%1' is not implemented", pp['case']; end if; S['RHS_normal'] := pp['RHS_normal'][j][3]; try return ParticularRationalSolution(S); catch "The system is inconsistent or not full rank": error; catch: return NULL; end try; end proc: ###################################################################### # HypergeometricSubstitutionToSystemInNormalForm ###################################################################### HypergeometricSubstitutionToSystemInNormalForm := proc(pp, u, Cert) local p2, R::nothing, P::nothing, s; p2 := table(['input' = 'matrices', 'matsize' = pp['matsize'], 'var' = pp['var'], 'case' = pp['case'], 'homogeneous' = true, 'Ring' = pp['Ring']]); if pp['case'] = 'qdifference' then p2['q_par'] := pp['q_par']; end if; userinfo(3, LFS, sprintf("Substitute %a", R(pp['var'])=u*P(pp['var']))); s := OreTools:-Properties:-GetSigma(pp['Ring'])(u, pp['var']); p2['source_system'] := LinearAlgebra:-MatrixAdd(pp['source_system'], -OreTools:-Properties:-Getdelta(pp['Ring'])( u, pp['var'])/s, u/s, 1); return p2; end; ###################################################################### # HypergeometricSubstitution ###################################################################### HypergeometricSubstitution := proc(pp, Cfs, Cert) local S, p1, k, j1, R::nothing, dy, dy1, h, x, sbst; option remember; if pp['input'] = 'matrices' then return HypergeometricSubstitutionToSystemInNormalForm(args); end if; x := pp['var']; sbst := OreTools:-Apply(OrePoly(0, 1), h(x), pp['Ring']) = Cert*h(x); S := table('sparse', [0=Cfs*pp['explicit_form'][0]]); dy := OrePoly(Cfs*h(x)); for k to pp['order'] do dy := OreTools:-Multiply(OrePoly(0, 1), dy, pp['Ring']); dy := eval(dy, sbst); dy1 := eval(dy, h(x)=1); for j1 from 0 to k do S[j1] := S[j1] + op(j1+1, dy1)*pp['explicit_form'][k]; end do; end do; p1 := table([ 'Ring' = pp['Ring'], 'var' = pp['var'], 'homogeneous' = true, 'functions' = pp['functions'], 'case' = pp['case'], 'order' = pp['order'], 'explicit_form' = copy(S, 'deep')]); if p1['case'] = 'qdifference' then p1['q_par'] := pp['q_par']; end if; return p1; end proc: PartRatSolForSystemInNormalForm := proc(pp) local m, tmp, x, k, l, u, M, a1, a2, Sol, c, v, indts; eval(pp); m := pp['matsize']; x := pp['var']; if pp['canonical'] <> true then pp['canonical'] := true; if pp['case'] in {'qdifference', 'difference'} then tmp := LinearFunctionalSystems:-CanonicalSystem( 'sh', pp['source_system'], pp['RHS'], pp['var'], pp['case']); if tmp = NULL then pp['canonical'] := "Incompatible system"; pp['part_rational_solution'] := NULL; return NULL; elif LinearAlgebra:-ColumnDimension(tmp[-1]) = 1 then userinfo(3, LFS, "The rational solution is found"); pp['canonical'] := tmp[-1]; pp['part_rational_solution'] := tmp[-1][1..-1,1]; return tmp[-1][1..-1,1]; elif LinearAlgebra:-ColumnDimension(tmp[1]) < m then pp['canonical'] := false; pp['source_system'] := tmp[1]; pp['RHS'] := tmp[2]; m := LinearAlgebra:-ColumnDimension(tmp[1]); pp['matsize'] := m; pp['canonical'] := tmp[-1]; end if; end if; end if; if pp['case'] in {'qdifference', 'difference'} then u := Vector(m, i -> primpart(lcm(denom(pp['RHS'][i]), seq(denom(pp['source_system'][i, k]), k = 1 .. m)), x)); M := LinearAlgebra:-DiagonalMatrix(u).pp['source_system']; M := map(dualnormal, LinearAlgebra:-MatrixInverse(M)); a1 := lcm(map(denom, convert(M, 'set'))[]); a2 := lcm(convert(u, 'set')[]); userinfo(3, LFS, "V = %1, W = %2", factor(a2), factor(a1)); if pp['case'] = 'qdifference' then a1 := dualnormal(a1/(x^ldegree(collect(a1, x), x))); a2 := dualnormal(a2/(x^ldegree(collect(a2, x), x))); u := 1/qDENOM(a2, a1, 1, pp['q_par'], x); else u := 1/`LREtools/DENOM`(a2, a1, 1, x); end if; end if; pp['homogeneous_system'] := InhomogeneousToHomogeneous(pp); if pp['case'] in {'qdifference', 'difference'} then pp['homogeneous_system']['part_of_denom'] := u; end if; try Sol := RationalSolutionForHomogeneous(pp['homogeneous_system']); catch: userinfo(3, 'LFS', StringTools[FormatMessage](lastexception[2..-1])); error; end try; if Sol = NULL then return NULL; end if; c := Sol[-1] - 1; indts := indets(c, _c['integer']); if nops(indts) > 0 then v := indts[1]; c := collect(c, v); Sol := eval(Sol, v=-coeff(c,v,0)/coeff(c,v,1)); elif c<>0 then userinfo(3, LFS, "There is no rational solution for inhomogeneous equation"); return NULL; end if; indts := indets(Sol, _c['integer']); Sol := ; Sol := map(dualnormal, Sol); pp['part_rational_solution'] := Sol; return Sol; end proc: ###################################################################### # PartSolForSystemInNormalForm # Output: particular solution for a system in the normal form which # right hand side is a hypergeometric term ###################################################################### PartSolForSystemInNormalForm := proc(p1, j) local pp, Sol, y::nothing, R::nothing, P::nothing; pp := table(['input' = 'matrices', 'matsize' = p1['matsize'], 'var' = p1['var'], 'case' = p1['case'], 'Ring' = p1['Ring']]); if pp['case'] = 'qdifference' then pp['q_par'] := p1['q_par']; end if; if OreTools:-Properties:-GetTheta1(p1['Ring']) = p1['RHS_normal'][j][2] then userinfo(3, LFS, "%1. The case of a rational right hand side", j); pp['source_system'] := copy(p1['source_system']); pp['RHS'] := copy(p1['RHS_normal'][j][3]); elif pp['case'] = 'differential' then userinfo(3, LFS, sprintf("%a. Substitute %a to the system", j, y(p1['var'])=p1['RHS_normal'][j][1]*R(p1['var']))); pp['source_system'] := LinearAlgebra:-MatrixAdd(p1['source_system'], -p1['RHS_normal'][j][2]); pp['RHS'] := p1['RHS_normal'][j][3]; elif pp['case'] in {'difference', 'qdifference'} then userinfo(3, LFS, sprintf("%a. Substitute %a to the system", j, y(p1['var'])=p1['RHS_normal'][j][1]*R(p1['var']))); pp['source_system'] := p1['source_system']/p1['RHS_normal'][j][2]; pp['RHS'] := p1['RHS_normal'][j][3]/p1['RHS_normal'][j][2]; else error "The case '%1' is not implemented", pp['case']; end if; p1[j] := copy(pp, deep); if p1['canonical'] = true then p1[j]['canonical'] := true; end if; Sol := PartRatSolForSystemInNormalForm(p1[j]); if j = 1 then if p1[j]['canonical'] = true then p1['canonical'] := true; elif p1[j]['canonical'] = "Incompatible system" then userinfo(3, 'LFS', "Incompatible system"); p1['canonical'] := "Incompatible system"; return NULL; else userinfo(3, LFS, "The system is not canonical"); p1['canonical'] := false; end if; end if; return Sol; end proc: ###################################################################### # RationalSolution # Output: rational solutions for the given system RationalSolution := proc() local Args, j, opts, pp, cfs, Res, outpt, indts, c, Field; if membertype(`=`({identical('rhs_normal'), identical('output'), identical('Normal')}, 'anything'), [args], 'j') then Args := [[args][1..j-1][], 'Dummy' = true, [args][j..-1][]]; else Args := [args]; end if; opts, Args := selectremove( el -> evalb(el::`=` and lhs(el) in {'output', 'earlyterminate'}), Args); outpt := 'basis'; hasoption(opts, 'output', 'outpt', 'opts'); # terminate_search_for_rational := false; terminate_search_for_rational := true; hasoption(opts, 'earlyterminate'=truefalse, 'terminate_search_for_rational', 'opts'); if select(el -> evalb(el::`=` and lhs(el) = 'rhs_normal'), Args) <> [] then try pp := Sdetermine(Args[]) catch: error; end try; else if not hasoption(opts, 'Normal', 'dualnormal', 'opts') then dualnormal := default_dualnormal; end if; try pp := Sdetermine(Args[]); catch: error; end try; if pp['homogeneous'] <> true then try pp['RHS_normal'] := map(dualnormal, pp['RHS']); catch: error cat("The right-hand side normalizing: ", StringTools[FormatMessage](lastexception[2..-1])); end try; if pp['RHS_normal']::'Vector[column]' and LinearAlgebra:-RowDimension(pp['RHS_normal']) = m then pp['homogeneous'] := LinearAlgebra:-Equal( pp['RHS_normal'], Vector(m)); else pp['homogeneous'] := false; end if; userinfo(3, LFS, "The right-hand side is normalized"); if pp['homogeneous'] then userinfo(3, LFS, "The system is homogeneous"); else userinfo(3, LFS, "The system is inhomogeneous"); end if; else pp['RHS_normal'] := pp['RHS']; end if; end if; if not pp['homogeneous'] and pp['order'] > 0 and not (pp['RHS_normal']::'Vector[column]' and LinearAlgebra:-RowDimension(pp['RHS_normal']) = LinearAlgebra:-RowDimension(pp['RHS'])) then error "The 'rhs_normal' option is wrong"; end if; cfs := [entries(pp['explicit_form'], 'nolist'), pp['RHS_normal']]; cfs := {map(op@convert, cfs, 'set')[]}; if not cfs::'set'('ratpoly'('algebraic', pp['var'])) then error "The system's coefficients of type ratpoly in %1 are expected", pp['var']; end if; indts := indets(cfs, 'name') minus {x, constants}; if assigned(pp['q_par']) and pp['q_par']::'name' then indts := select(proper_name, indts, {x, pp['q_par']}); indts := indts union {pp['q_par']}; else indts := select(proper_name, indts, {x}); end if; Field := 'ratpoly'('rational', indts); try cfs := remove(type, cfs, 'ratpoly'(Field, x)); catch: error cat("Impossible to check system's coefficients: ", [lastexception][2]), [lastexception][3..-1][]; end try; if cfs <> {} then error "The system's coefficients of type ratpoly in %1 over the field '%2' are expected but '%3' is", x, Field, cfs[1]; end if; if outpt = 'partsol' and pp['order'] = 0 then try return pp['explicit_form'][0]^(-1).pp['RHS']; catch: userinfo(3, 'LFS', "There is no particular solution"); return NULL; end try; elif outpt = 'partsol' and not pp['homogeneous'] then try Res := ParticularRationalSolution(pp); catch "The system is inconsistent or not full rank": userinfo(3, LFS, "The system is inconsistent or not full rank"); return NULL; catch: error; end try; if Res = NULL then userinfo(3, 'LFS', "There is no particular solution"); return NULL; elif not pp['source_system']::'Matrix' and LinearAlgebra:-RowDimension(pp['functions']) = 1 then return Res[1]; else return Res; end if; elif outpt = 'basis' and not pp['homogeneous'] then try Res := RationalSolutionForInhomogeneous(pp); catch "The system is inconsistent or not full rank": userinfo(3, LFS, "The system is inconsistent or not full rank"); return NULL; catch: error; end try; if Res = NULL then userinfo(3, 'LFS', "There is no particular solution"); return NULL; elif not pp['source_system']::'Matrix' and LinearAlgebra:-RowDimension(pp['functions']) = 1 then return [map(el -> el[1], Res[1]), Res[2][1]]; else return Res; end if; elif outpt = 'basis' and pp['homogeneous'] then try Res := RationalSolutionForHomogeneous(pp); catch "The system is inconsistent or not full rank": userinfo(3, LFS, "The system is not full rank"); return NULL; catch: error; end try; if Res = NULL then return []; end if; indts := indets(Res, _c['integer']); Res := [seq(map(dualnormal, eval(eval(, c=1), map(`=`, indts, 0))), c = indts)]; if Res = [] then userinfo(3, LFS, "Zero solution is found"); elif not pp['source_system']::'Matrix' and LinearAlgebra:-RowDimension(pp['functions']) = 1 then return map(el -> el[1], Res); else return Res; end if; else error "Only two cases are implemented: output=partsol (for inhomogeneous systems) and output=basis."; end if; end proc; ###################################################################### # InhomogeneousToHomogeneous # Output: an homogeneous system which is equivalent the given # inhomogeneous system ###################################################################### InhomogeneousToHomogeneous := proc(p1) local number_of_unknowns, number_of_equations, S, j, y; # if p1['input'] = matrices then if not assigned(p1['explicit_form']) then eval(p1); S := table(['Ring' = p1['Ring'], 'var' = p1['var'], 'homogeneous' = true, 'case' = p1['case'], 'input' = 'matrices', 'matsize' = p1['matsize']+1]); if S['case'] = 'qdifference' then S['q_par'] := p1['q_par']; end if; S['source_system'] := <, Vector['row'](S['matsize'], {(S['matsize']) = OreTools:-Properties:-GetTheta1(S['Ring'])})>; userinfo(3, LFS, "The homogeneous equation has been constructed"); return S; end if; number_of_equations := LinearAlgebra:-RowDimension( p1['explicit_form'][0]); number_of_unknowns := LinearAlgebra:-ColumnDimension( p1['explicit_form'][0]); S := table([ 'Ring' = p1['Ring'], 'var' = p1['var'], 'homogeneous' = true, 'functions' = [1 .. -1, 1], 'case' = p1['case'], 'order' = p1['order']]); if p1['case'] = 'qdifference' then S['q_par'] := p1['q_par']; end if; S['explicit_form'][0] := <, Vector['row'](number_of_unknowns+1)>; S['explicit_form'][1] := <| <(0$number_of_equations), 1>>; for j from 2 to p1['order'] do S['explicit_form'][j] := <, Vector['row'](number_of_unknowns+1)>; end do; if p1['case'] in {'qdifference', 'difference'} then S['explicit_form'][0][-1, -1] := -1; end if; userinfo(3, LFS, "The homogeneous equation has been constructed"); return S; end proc: ###################################################################### # ParticularRationalSolution # Output: a particular rational solution for an inhomogeneous system ###################################################################### ParticularRationalSolution := proc(p1) local S, c, v, Sol, indts; S := InhomogeneousToHomogeneous(p1); try Sol := RationalSolutionForHomogeneous(S); catch: userinfo(3, 'LFS', StringTools[FormatMessage](lastexception[2..-1])); error; end try; if Sol = NULL then return NULL; end if; c := Sol[-1] - 1; indts := indets(c, _c['integer']); if nops(indts) > 0 then v := indts[1]; c := collect(c, v); Sol := eval(Sol, v=-coeff(c,v,0)/coeff(c,v,1)); elif c<>0 then userinfo(3, LFS, "There is no rational solution for inhomogeneous equation"); return NULL; end if; indts := indets(Sol, _c['integer']); Sol := ; Sol := map(dualnormal, Sol); return Sol; end proc: ###################################################################### # RationalSolutionForInhomogeneous # Output: a general rational solution for an inhomogeneous system ###################################################################### RationalSolutionForInhomogeneous := proc(p1) local S, Sol, c, indts, v; S := InhomogeneousToHomogeneous(p1); try Sol := RationalSolutionForHomogeneous(S); catch: userinfo(3, 'LFS', StringTools[FormatMessage](lastexception[2..-1])); error; end try; if Sol = NULL then return NULL; end if; c := Sol[-1] - 1; indts := indets(c, _c['integer']); if nops(indts) > 0 then v := indts[1]; c := collect(c, v); Sol := eval(Sol, v=-coeff(c,v,0)/coeff(c,v,1)); elif c<>0 then userinfo(3, LFS, "There is no rational solution for inhomogeneous equation"); return NULL; end if; indts := indets(Sol, _c['integer']); Sol := [Sol[1..-2], ]; Sol := [-Sol[2], Sol[2]]; Sol := [[seq(map(dualnormal, eval(eval(Sol[1], c=1), map(`=`, indts, 0))), c = indts)], map(dualnormal, Sol[2])]; return Sol; end proc: ###################################################################### DegreeBoundForHomogeneous := proc(p1, bound) local deg, number_of_unknowns, number_of_equations, cns, vert_shft, rM, n, ncol, nlin, MM, p, rr, i, j, k, TM, maxn, t; if assigned(p1['bound_of_polynomial_solution_degree']) then return p1['bound_of_polynomial_solution_degree'][1]; end if; InducedRecurrenceForHomogeneous(p1); deg := -infinity; if p1['input'] = 'matrices' then number_of_unknowns := p1['matsize']; number_of_equations := p1['matsize']; else number_of_unknowns := LinearAlgebra:-RowDimension(p1['functions']); number_of_equations := LinearAlgebra:-RowDimension( p1['explicit_form'][0]); end if; rM := copy(p1['induced_recurrence'][1]), p1['induced_recurrence'][2..3]; n := p1['induced_recurrence'][4]; ncol := rM[2]-rM[3]+1; nlin := LinearAlgebra:-RowDimension(rM[1]); # Triangualize the trailing matrix and bound the degree of # the polynomial solution (by methods) cns := table(); vert_shft := Vector(number_of_unknowns); userinfo(3, LFS, "EG_trail for the induced recurrence"); EG_('trail', rM[1], ncol-1, number_of_unknowns, nlin, n, p1['induced_recurrence'][5], table({(1)=cns, (2)=vert_shft})); if number_of_equations < number_of_unknowns or convert(rM[1][number_of_unknowns], 'set') = {0} then error "The system is not full rank"; end if; MM := LinearAlgebra:-SubMatrix(rM[1], 1..number_of_unknowns, number_of_unknowns*(ncol-1)+1..number_of_unknowns*ncol); p := dualnormal( LinearAlgebra:-LUDecomposition(MM, 'output'='determinant')); userinfo(3, LFS, "The determinant of the trailing matrix is %1", factor(p)); t := rM[3] + max(convert(vert_shft, 'set')); userinfo(3, LFS, "The indicial equation is %1=0", factor(eval(p, n = n - t))); rr := Solve(p, n); # Check constraints if nargs > 1 and bound::'integer' then deg := bound; elif nops(rr)>0 then userinfo(3, LFS, "The maximal integer root of the indicial equation is %1", rr[-1] + t); for i from nops(rr) by -1 to 1 do maxn := rr[i]; if assigned(cns[maxn]) then TM := Matrix([[map(dualnormal,eval(MM,n=maxn))], seq([Vector['row'](number_of_unknowns, [seq(cns[maxn][j][k], k = number_of_unknowns*(ncol-1)+1 .. number_of_unknowns*ncol)])], j = 1 .. nops(cns[maxn])) ]); if LinearAlgebra:-Rank(TM) < number_of_unknowns then break; else maxn := -infinity; end if; else break; end if; end do; if maxn > -infinity then deg := max(deg, maxn + t); end if; else userinfo(3, LFS, "The indicial equation has no integer roots"); end if; p1['bound_of_polynomial_solution_degree'] := deg, rr; p1['EGt_induced_recurrence'] := rM[1], cns, vert_shft, {}, {}; return deg; end proc; ###################################################################### # RationalSolutionForHomogeneous # Output: a general rational solution for an homogeneous system ###################################################################### RationalSolutionForHomogeneous := proc(p1) local u, Syst, Sol, R::nothing, P::nothing, rM, bound, i, syst, slv, s1, indts, v::nothing, st, st1, number_of_unknowns, pu, y::nothing; if p1['input'] = 'matrices' then number_of_unknowns := p1['matsize']; else number_of_unknowns := LinearAlgebra:-RowDimension(p1['functions']) end if; st := time(): if terminate_search_for_rational = true then st1 := time(): bound := DegreeBoundForHomogeneous(p1); userinfo(2, LFS, "Time of degree bound computing", time()-st1); if bound = -infinity then userinfo(2, LFS, "The first control point. There is no rational solution: %1", time()-st); return NULL; end if; end if; st1 := time(); try u := UniversalDenominatorForHomogeneous(p1); catch: error; end try; userinfo(2, LFS, "Time of universal denominator constructing", time()-st1); if u = FAIL then return NULL; end if; if terminate_search_for_rational = true then bound := bound + degree(denom(u), p1['var']); userinfo(3, LFS, "The maximal integer root of the indicial equation plus the degree of the universal denominator is %1", bound); if bound < 0 then userinfo(2, LFS, "The second control point. There is no rational solution:", time()-st); return NULL; end if; end if; if (terminate_search_for_rational = false or p1['case'] in {'differential', 'difference'}) and u <> 1 then userinfo(3, LFS, sprintf("Substitute %a", R(p1['var'])=u*P(p1['var']))); st1:=time(); Syst := HypergeometricSubstitution(p1, u, OreTools:-Properties:-GetTheta1(p1['Ring'])); userinfo(2, LFS, sprintf("Time of substitution %a", R(p1['var'])='U'(p1['var'])*P(p1['var'])), time()-st1); else Syst := p1; end if; st1 := time(); DegreeBoundForHomogeneous(Syst, bound); userinfo(2, LFS, "Time of construction of the induced recurrence and EG doing", time()-st1); if terminate_search_for_rational = false and Syst['bound_of_polynomial_solution_degree'][1] < 0 then if Syst['bound_of_polynomial_solution_degree'][1] = - infinity then userinfo(2, LFS, "There is no integer root of the indicial equation"); else userinfo(2, LFS, "The maximum integer root of the indicial equation is %1<0", Syst['bound_of_polynomial_solution_degree'][1]); end if; userinfo(2, LFS, "There is no rational solution:", time()-st); return NULL; end if; rM := Syst['induced_recurrence'][1..3]; #print(add(map(factor,rM[1][1 .. -1, 1 + number_of_unknowns*(i - rM[3]) .. number_of_unknowns + number_of_unknowns*(i - rM[3])]) . (v(n - i + rM[3] + rM[2])), i = rM[3] .. rM[2])); if terminate_search_for_rational = true and not p1['case'] in {'differential', 'difference'} then bound := degree(denom(u), p1['var']); # pu := Sdetermine(denom(u)*y(p1['var']), y(p1['var'])); # pu['case'] := p1['case']; # pu['Ring'] := p1['Ring']; # userinfo(3, LFS, "Construct the induced recurrence for the universal denominator:"); # InducedRecurrenceForHomogeneous(pu); if p1['case'] = 'difference' then userinfo(3, LFS, "Construct a truncation of Laurent solution at %1: %2 + ... + %3 + %4", Syst['var']=infinity, v(Syst['bound_of_polynomial_solution_degree'][1])* product(Syst['var'] - k, k = 0 .. Syst['bound_of_polynomial_solution_degree'][1] - 1), v(-bound-rM[2]+rM[3])*product(Syst['var'] - k, k = 0 .. -bound-rM[2]+rM[3] - 1), O(Syst['var']^(-bound-rM[2]+rM[3]-1))); else userinfo(3, LFS, "Construct a truncation of Laurent solution at %1: %2 + ... + %3 + %4", Syst['var']=infinity, v(Syst['bound_of_polynomial_solution_degree'][1])* Syst['var']^Syst['bound_of_polynomial_solution_degree'][1], v(-bound+rM[2]-rM[3])*Syst['var']^(-bound-rM[2]+rM[3]), O(Syst['var']^(-bound-rM[2]+rM[3]-1))); end if; elif p1['case'] = 'difference' then userinfo(3, LFS, "Construct a polynomial solution: %1 + ... + %2", v(Syst['bound_of_polynomial_solution_degree'][1])* product(Syst['var'] - k, k = 0 .. Syst['bound_of_polynomial_solution_degree'][1] - 1), v(0)); bound := 0; else userinfo(3, LFS, "Construct a polynomial solution: %1 + ... + %2", v(Syst['bound_of_polynomial_solution_degree'][1])* Syst['var']^Syst['bound_of_polynomial_solution_degree'][1], v(0)); bound := 0; end if; st1 := time(); try Sol := SolveRec:-ConstructSolutionByTransformed( [Syst['bound_of_polynomial_solution_degree']], Syst['case'], [Syst['EGt_induced_recurrence']], number_of_unknowns, rM[2] - rM[3] + 1, LinearAlgebra:-RowDimension(rM[1]), rM[2], rM[3], Syst['induced_recurrence'][4], Syst['var'], -bound); catch: userinfo(3, LFS, lastexception); Sol := FAIL; end try; if terminate_search_for_rational = true and not p1['case'] in {'differential', 'difference'} then else userinfo(2, LFS, sprintf("Time of %a computing", P(p1['var'])), time() - st1); end if; if Sol = FAIL or Sol = NULL or convert(Sol, 'set') = {0} or Sol = [] then userinfo(2, LFS, "There is no rational solution: %1", time()-st); return NULL; end if; userinfo(5, LFS, "The truncation has been computed:", ); if terminate_search_for_rational = true and not p1['case'] in {'differential', 'difference'} then if bound > 0 then st1 := time(); Sol := map(el -> series((el + O(p1['var']^(-(bound+rM[2]-rM[3]+1))))* denom(u), p1['var']=infinity, bound+rM[2]-rM[3]+1), Sol); userinfo(2, LFS, "Time of multiplying the truncation by the universal denominator", time() - st1); userinfo(5, LFS, "The result of multiplying the truncation by the universal denominator:", ); Sol := map(el -> convert(el, 'polynom'), Sol); bound := min(map(ldegree, Sol, p1['var'])); if bound < 0 then indts := indets(Sol, _c['integer']); syst := {seq(map(coeff, Sol, p1['var'], i)[], i = bound .. -1)} minus {0}; try slv := solve(syst, indts); catch: slv := NULL; while syst <> {} do s1 := solve(syst[1], indts); slv := remove(evalb, eval(slv, s1) union s1); syst := map(expand, eval(syst, slv)) minus {0}; end do; end try; if slv <> NULL then Sol := eval(Sol, slv); end if; end if; if convert(Sol, 'set') = {0} then userinfo(2, LFS, "There is no rational solution: %1", time()-st); return NULL; end if; Sol := map(el-> dualnormal(u*el), Sol); userinfo(2, LFS, "The rational solution is found: %1", time()-st); else # denom(u) = 1 userinfo(2, LFS, "The polynomial solution is found: %1", time()-st); end if; elif u <> 1 then st1 := time(); Sol := map(el-> map(dualnormal, u*el), Sol); userinfo(2, LFS, sprintf("Time of dividing %a", P(p1['var'])/'U'(p1['var'])), time()-st1); userinfo(2, LFS, "The rational solution is found: %1", time()-st); else # u = 1 userinfo(2, LFS, "The polynomial solution is found: %1", time()-st); end if; return Sol; end proc: InducedRecurrenceForHomogeneous := proc(p1) local number_of_unknowns, number_of_equations, c, i, j, k, Syst, funcs, pp, rM, n, indts, deg, ncol, nlin, cns, vert_shft, MM, p, rr, B::nothing, v::nothing, y::nothing; if assigned(p1['induced_recurrence']) then return; end if; kernelopts(opaquemodules = false); if p1['input'] = 'matrices' then pp := table([ 'method' = 'ordinary', 'case' = p1['case'], 'q_par' = p1['q_par'], 'input' = 'matrices', 'var' = p1['var'], 'mat' = p1['source_system'], 'vec' = Vector['column'](p1['matsize']), 'matsize' = p1['matsize'], 'matdenom' = Vector(p1['matsize'], i -> primpart(lcm(seq(denom(p1['source_system'][i, k]), k = 1 .. p1['matsize'])), p1['var']))]); LinearFunctionalSystems:-Properties( pp['mat'], pp['var'], pp['case']) := pp; rM := LinearFunctionalSystems:-RecMat( pp['mat'], pp['var'], pp['case']); number_of_unknowns := pp['matsize']; else number_of_unknowns := LinearAlgebra:-ColumnDimension( p1['explicit_form'][0]); number_of_equations := LinearAlgebra:-RowDimension( p1['explicit_form'][0]); if p1['polycoeffs'] <> true then c := [seq(primpart(lcm( seq(seq(denom(p1['explicit_form'][j][i,k]), j=0..p1['order']), k=1..number_of_unknowns)), p1['var']), i = 1..number_of_equations)]; for i to number_of_equations do if c[i] <> 1 then for j from 0 to p1['order'] do p1['explicit_form'][j][i] := map(dualnormal, p1['explicit_form'][j][i]*c[i]); end do; end if; end do; p1['polycoeffs'] := true; end if; Syst := add(p1['explicit_form'][j]. map2(OreTools:-Apply, OrePoly(0$j, 1), p1['functions'], p1['Ring']), j = 0 .. p1['order']); Syst := convert(Syst, 'list'); funcs := convert(p1['functions'], 'list'); pp := table([ 'functions' = funcs, 'system0' = Syst, 'method' = 'ordinary', 'system' = Syst, 'case' = p1['case'], 'q_par' = p1['q_par'], 'shifts' = [({seq(i, i = 0..p1['order'])})$(nops(Syst))], 'qshifts' = [({seq(i, i = 0..p1['order'])})$(nops(Syst))], 'input' = 'lists', 'order' = p1['order'], 'f_names' = map2(op, 0, funcs), 'var' = p1['var']]); LinearFunctionalSystems:-Properties(pp['system0'], pp['functions']) := pp; rM := LinearFunctionalSystems:-RecMat(pp['system0'], pp['functions']); end if; indts := indets(rM[1], 'name'); indts := select(el -> convert(el, 'string') = "LinearFunctionalSystems/n", indts); if indts <> {} then LinearAlgebra:-Map(eval, rM[1], indts[1] = n); end if; for i from rM[3] to rM[2] while convert(rM[1][1..-1,(rM[2]-i)*number_of_unknowns+1 .. (rM[2]-i+1)*number_of_unknowns], 'set') = {0} do end do; rM[3] := i; userinfo(3, LFS, "The induced recurrence is %1+ ... + %2=0", B[rM[2]].v(n+rM[2]), B[rM[3]].v(n+rM[3])); p1['induced_recurrence'] := copy(rM[1][1..-1, 1 .. (rM[2]-rM[3]+1)*number_of_unknowns]), rM[2], rM[3], n, OreTools:-SetOreRing(n, 'shift'); return ; end proc: ###################################################################### # UniversalDenominatorForHomogeneousSystemInNormalForm # Output: the universal denominator for the homogeneous system ###################################################################### UniversalDenominatorForHomogeneousSystemInNormalForm := proc(pp) local rM, u, tt, k; if assigned(pp['part_of_denom']) then u := pp['part_of_denom']; elif pp['case'] in {'qdifference', 'qdifference'} then u := 1; end if; if pp['case'] = 'qdifference' then InducedRecurrenceForHomogeneous(pp); rM := copy(pp['induced_recurrence'][1]), pp['induced_recurrence'][2], pp['induced_recurrence'][3]; userinfo(3, LFS, "EG_lead for the induced recurrence"); EG_('lead', rM[1], rM[2]-rM[3], pp['matsize'], pp['matsize'], pp['induced_recurrence'][4], pp['induced_recurrence'][5]); tt := LinearAlgebra:-Determinant( rM[1][1 .. pp['matsize'], 1 .. pp['matsize']]); tt := qSolve(tt, pp['var'], pp['q_par'], pp['induced_recurrence'][4]); tt := min(0, min(tt) + rM[2]); u := u*pp['var']^tt; userinfo(3, LFS, "A lower bound of valuation of Laurent solutions of the system at %1 is %2", pp['var'] = 0, tt); elif pp['case'] = 'differential' then tt := table([ 'method' = 'ordinary', 'case' = pp['case'], 'input' = 'matrices', 'var' = pp['var'], 'mat' = pp['source_system'], 'vec' = Vector['column'](pp['matsize']), 'homogeneous' = true, 'matsize' = pp['matsize'], 'matdenom' = Vector(pp['matsize'], i -> primpart(lcm(seq(denom(pp['source_system'][i, k]), k = 1 .. pp['matsize'])), pp['var']))]); LinearFunctionalSystems:-Properties( tt['mat'], tt['var'], tt['case']) := tt; try u := LinearFunctionalSystems:-UniversalDenominator( pp['source_system'], pp['var'], pp['case']); catch: userinfo(3, LFS, StringTools[FormatMessage](lastexception[2..-1])); u := FAIL; end try; if u = FAIL then return FAIL; end if; end if; userinfo(3, 'LFS', "The universal denominator is", denom(u)); pp['universal_denominator'] := denom(u); return u; end proc: ###################################################################### # UniversalDenominatorForHomogeneous # Output: the universal denominator for the homogeneous system ###################################################################### UniversalDenominatorForHomogeneous := proc(pp) local S, x, OR, number_of_equations, number_of_unknowns, j, k, indts, M_l, leading_invers, V, M_t, trailing_invers, W, Res, q, y, rM, tt, d, i, a, c, p; if not assigned(pp['explicit_form']) then return UniversalDenominatorForHomogeneousSystemInNormalForm(pp); end if; x := pp['var']; OR := pp['Ring']; number_of_equations := LinearAlgebra:-RowDimension( pp['explicit_form'][0]); number_of_unknowns := LinearAlgebra:-ColumnDimension( pp['explicit_form'][0]); if pp['polycoeffs'] <> true then c := [seq(primpart(lcm(denom(pp['RHS_normal'][i]), seq(seq(denom(pp['explicit_form'][j][i,k]), j=0..pp['order']), k=1..number_of_unknowns)), x), i = 1..number_of_equations)]; for i to number_of_equations do if c[i] <> 1 then for j from 0 to pp['order'] do pp['explicit_form'][j][i] := map(dualnormal, pp['explicit_form'][j][i]*c[i]); end do; end if; end do; pp['polycoeffs'] := true; end if; pp['explicit_matrix'] := Matrix([seq(pp['explicit_form'][pp['order']-i], i=0..pp['order'])]); pp['EGl_explicit_matrix'] := copy(pp['explicit_matrix']); M_l := pp['EGl_explicit_matrix']; userinfo(3, LFS, "EG_lead for the given system"); EG_('lead', M_l, pp['order'], number_of_unknowns, number_of_equations, x, OR); if number_of_equations < number_of_unknowns or convert(M_l[number_of_unknowns], 'set') = {0} then userinfo(3, LFS, "The system is not full rank"); error "The system is not full rank"; end if; if pp['case'] in {'qdifference', 'difference'} then leading_invers := (M_l[1 .. number_of_unknowns, 1 .. number_of_unknowns])^(-1); V := lcm(convert(map(denom, leading_invers), 'set')[]); pp['EGt_explicit_matrix'] := copy(pp['explicit_matrix']); M_t := pp['EGt_explicit_matrix']; userinfo(3, LFS, "EG_trail for the given system"); EG_('trail', M_t, pp['order'], number_of_unknowns, number_of_equations, x, OR); trailing_invers := (M_t[1 .. number_of_unknowns, -number_of_unknowns .. -1])^(-1); W := lcm(convert(map(denom, trailing_invers), 'set')[]); end if; if pp['case'] = 'difference' then userinfo(3, LFS, "V = %1, W = %2", factor(V), factor(W)); Res := `LREtools/DENOM`(V, W, pp['order'], x); if has(Res, FAIL) then userinfo(3, 'LFS', "Wrong result of `LREtools/DENOM`: %1", Res); V := primpart(V, x); W := primpart(W, x); Res := 1/`LREtools/DENOM`(V, W, pp['order'], x); else Res := 1/Res; end if; elif pp['case'] = 'qdifference' then q := pp['q_par']; Res := 1/qDENOM(V, W, pp['order'], q, x); d := {x}; else # pp['case'] = 'differential' d := LinearAlgebra:-Determinant(M_l[1 .. number_of_unknowns, 1 .. number_of_unknowns]); d := map(el -> map2(op, 1, factors(el[1])[2])[], {sqrfree(d)[2][]}); userinfo(3, LFS, "A revealing polynomial is %1", mul(d)); if d = {} then userinfo(3, 'LFS', "The universal denominator is", 1); pp['universal_denominator'] := 1; return 1; end if; Res := 1; end if; if pp['case'] in {'qdifference', 'differential'} then for i to nops(d) do a := RootOf(d[i], x); userinfo(3, LFS, "For the system at %1", x = a); if a <> 0 then S := table(['case' = pp['case'], 'Ring' = pp['Ring'], 'explicit_form' = eval(copy(pp['explicit_form']), x = x + a), 'var' = x, 'order' = pp['order'], 'functions' = pp['functions'], 'homogeneous' = true, 'polycoeffs' = pp['polycoeffs']]); else S := pp; end if; InducedRecurrenceForHomogeneous(S); rM := copy(S['induced_recurrence'][1]), S['induced_recurrence'][2], S['induced_recurrence'][3]; userinfo(3, LFS, "EG_lead for the induced recurrence"); EG_('lead', rM[1], rM[2]-rM[3], number_of_unknowns, LinearAlgebra:-RowDimension(rM[1]), S['induced_recurrence'][4], S['induced_recurrence'][5]); p := LinearAlgebra:-Determinant(rM[1][1 .. number_of_unknowns, 1 .. number_of_unknowns]); if pp['case'] = 'qdifference' then tt := qSolve(p, x, q, S['induced_recurrence'][4]); else tt := Solve(p, S['induced_recurrence'][4]); end if; tt := min(0, min(tt) + rM[2]); if tt <> 0 then Res := Res*d[i]^tt; end if; userinfo(3, LFS, "A lower bound of valuation of Laurent solutions is %1", tt); end do; end if; userinfo(3, 'LFS', "The universal denominator is", denom(Res)); pp['universal_denominator'] := denom(Res); return Res; end proc: ###################################################################### qDENOM := proc(an, a0, n, q, x ) local u, A, B, c, k, i, h_s ; A := evala(Normal(subs( x=x*q^(-n), an ))) ; B := evala(Normal(a0)) ; kernelopts(opaquemodules = false); h_s := sort([op(LinearFunctionalSystems:-qTools:-qDispersion( B, A, q, x))]); if h_s = [FAIL] then return 1; end if; u := 1 ; for k from nops(h_s) by -1 to 1 do c := evala(Gcdex( A, subs(x=x*q^(h_s[k]),B),x)) ; u := u * mul(subs(x=x*q^(-i), c),i=0..h_s[k]) ; A := evala(Normal(A/c)) ; B := evala(Normal(B/subs(x=x*q^(-h_s[k]),c))); end do ; return u; end proc: # qDENOM ###################################################################### # EG_ ###################################################################### EG_ := proc(lt, ExplicitMatrix, n, m, rdim, x, Ring, Constraints := 0) # n - order # m - number of unknowns # rdim - number of equations local k, v, c, l_c, l_c1, R, den, rs, r, vv; c := 0; k := 1; while k <= rdim-c do if convert(ExplicitMatrix[k], 'set') = {0} then if k < rdim-c and convert(ExplicitMatrix[rdim-c], 'set') <> {0} then userinfo(3, LFS, "Swap the %-1 and the %-2 rows", k, rdim-c); ExplicitMatrix[k] := ExplicitMatrix[rdim-c]; ExplicitMatrix[rdim-c] := Vector['row']((n+1)*m); end if; c := c+1; else k := k+1; end if; end do; if c = rdim then return; elif c > 0 then procname(lt, ExplicitMatrix, n, m, rdim-c, x, Ring, Constraints); return; end if; if Constraints <> 0 then # for k to m do # Constraints[2][k] := Constraints[2][k] + # Shift_Column(lt, ExplicitMatrix, k, n, m, rdim, x, Ring); # end do; end if; for k to rdim do Shift_Row(lt, ExplicitMatrix, k, m, x, Ring); end do; if lt='lead' then v := LinearAlgebra:-NullSpace(ExplicitMatrix[1..rdim, 1..m]^%T); else v := LinearAlgebra:-NullSpace(ExplicitMatrix[1..rdim, -m..-1]^%T); end if; if v = {} then return; end if; v := v[1]; l_c := 0; for k to rdim do l_c1 := Row_Length(lt, ExplicitMatrix[k]); if lt = 'lead' then if v[k] <> 0 and l_c < l_c1 then c := k; l_c := l_c1; end if; else if v[k] <> 0 and l_c < l_c1 then c := k; l_c := l_c1; end if; end if; end do; userinfo(3, LFS, "The %1 reduction of the %-2 row", lt, c); den := lcm(map(denom,convert(v, 'set'))[]); if depends(den, x) then v := map(dualnormal, map(`*`, v, den)); end if; if Constraints <> 0 then rs := Solve(v[c], x); for r in rs do vv := convert(eval(ExplicitMatrix[c], x = r), 'list'); if convert(vv, 'set') = {0} then elif assigned(Constraints[1][r]) then Constraints[1][r] := [op(Constraints[1][r]),vv]; else Constraints[1][r] := [vv]; end if; end do; end if; R := map(dualnormal,(v^%T).ExplicitMatrix[1..rdim]); den := lcm(map(denom,convert(R, 'set'))[]); if depends(den, x) then R := map(dualnormal, den*R); end if; if convert(R, 'set') = {0} then if c < rdim then userinfo(3, LFS, "Swap the %-1 and the %-2 rows", c, rdim); ExplicitMatrix[c] := ExplicitMatrix[rdim]; end if; ExplicitMatrix[rdim] := R; procname(lt, ExplicitMatrix, n, m, rdim-1, x, Ring, Constraints); return; else ExplicitMatrix[c] := R; procname(lt, ExplicitMatrix, n, m, rdim, x, Ring, Constraints); return; end if; end proc: Shift_Column := proc(lt, C, k, n, m, rdim, x, Ring) # k -- unknown # n -- order # m - number of unknowns # rdim - number of equations local strt, stp, fnsh, i, j; if lt = 'lead' or OreTools:-Properties:-GetRingName(Ring) <> 'shift' then error "01"; end if; strt := `if`(lt = 'lead', 0, n); fnsh := `if`(lt = 'lead', n, 0); stp := `if`(lt = 'lead', 1, -1); for i from strt by stp to fnsh while convert(C[1..-1, m*i+k], 'set') = {0} do end do; if i = strt or i = fnsh + stp then return 0; end if; i := i; for j from strt by stp to strt-i do C[1..-1, m*j+k] := C[1..-1, m*(j-(strt-i))+k]; end do; for j from j by stp to fnsh do C[1..-1, m*j+k] := Vector(rdim); end do; i := (i-strt)*stp; if i = 1 then userinfo(3, LFS, "The %1 shift of the %-2 unknown", lt, k); else userinfo(3, LFS, "%1 %2 shifts of the %-3 unknown", i, lt, k); end if; return i; end proc; Shift_Row := proc(lt, C, k, m, x, Ring) local l, l1, l2, den, rng, dlt; try dlt := OreTools:-Properties:-Getdelta(Ring)(x, x); catch: end try; rng := `if`(lt = 'lead', 1 .. m, -m.. -1); l := Row_Length(lt, C[k]); while l <> 0 and convert(C[k][rng], 'set') = {0} do if lt = 'lead' and dlt <> 0 then den := C[k][l]; elif lt = 'trail' and dlt <> 0 then error "It's impossible to do trail shift for the case '%1'", OreTools:-Properties:-GetRingName(Ring); else den := 1; end if; if depends(den, x) then LinearAlgebra:-Map[(i,j) -> evalb(i=k)](dualnormal@`*`, C, 1/den); end if; if lt = 'lead' then userinfo(3, LFS, "The lead shift of the %-1 row", k); C[k] := [1] + map(OreTools:-Properties:-Getdelta(Ring), C[k], x); else l2 := 0; for l1 from LinearAlgebra:-ColumnDimension(C) to m by -m while convert(C[k][l1-m+1..l1], 'set') = {0} do l1-m+1..l1; l2 := l2 + 1; end do; if l2 = 1 then userinfo(3, LFS, "The trail shift of the %-1 row", k); else userinfo(3, LFS, "%1 trail shifts of the %-2 row", l2, k); end if; C[k,-l..-1] := OreTools:-Properties:-GetSigmaInverse(Ring)(el, x))@@l2) (C[k][-l..l1])>[1]; end if; LinearAlgebra:-Map[(i,j) -> evalb(i=k)](dualnormal, C); l := Row_Length(lt, C[k]); end do; den := lcm(map(denom,convert(C[k], 'set'))[]); if depends(den, x) then LinearAlgebra:-Map[(i,j) -> evalb(i=k)](dualnormal@`*`, C, den); end if; return; end proc: Row_Length := proc(lt, c) local k; for k from LinearAlgebra:-ColumnDimension(c) by -1 to 1 do if lt = 'lead' and c[k] <> 0 then return k; elif lt = 'trail' and c[-k] <> 0 then return k; end if; end do; return 0; end proc: ###################################################################### # Sdetermine # Output: a table with system's properties ###################################################################### Sdetermine := proc() local next_arg, m, A, b, x, y, i, err_str, q, nn, pp, opts, norm, S, dummy, term, indts, n, m2, k, shifts, dot_mul, cf, y_mul; option remember; next_arg := 1; if nargs < next_arg then error "The %-1 argument (a system) is expected", next_arg; end if; if args[next_arg]::'Matrix' then A := args[next_arg]; if not A::'Matrix'('square') then error "A square matrix is expected" elif not [entries(A, 'nolist')]::'list'('scalar') then error "A matrix with entries of type/scalar is expected" end if; m := LinearAlgebra:-RowDimension(A); if m = 0 then error "Non empty matrix is expected" end if; pp := table({ 'source_system' = A, 'normal_form' = Array(0..0, [A]), 'explicit_form' = Array(0..1, [-A, LinearAlgebra:-IdentityMatrix(m, compact=false)]), 'order' = 1, 'input' = 'matrices', 'matsize' = m}); userinfo(3, LFS, "The system is in the normal form"); next_arg := next_arg + 1; if nargs < next_arg then error "The %-1 argument of type Vector or type name is expected", next_arg; elif args[next_arg]::'Vector[column]' then b := args[next_arg]; if m <> LinearAlgebra:-RowDimension(b) then error "The %-1 argument is wrong, a Vector with %2 rows is expected", next_arg, m; elif not [entries(b, 'nolist')]::'list'('scalar') then error "The %-1 argument is wrong, a Vector with entries of type/scalar is expected", next_arg end if; pp['RHS'] := b; next_arg := next_arg+1; if nargs < next_arg or args[next_arg]::`=` and lhs(args[next_arg]) = 'Dummy' then error "The %-1 argument (an independent variable) is expected", next_arg; elif not args[next_arg]::{('name' = 'anything'^'name'), 'name'} then error "The %-1 argument of type name is expected", next_arg; end if; elif args[next_arg]::{('name' = 'anything'^'name'), 'name'} then pp['RHS'] := Vector(m); pp['homogeneous'] := true; userinfo(3, LFS, "The system is homogeneous"); else error "The %-1 argument of type Vector or type name is expected", next_arg; end if; if args[next_arg]::'name' then x := args[next_arg]; if not proper_name(x, {cases[], constants}) then error "The %-1 argument is wrong, an independent name cannot be '%2'", next_arg, x; end if; else # args[next_arg]::('name' = 'anything'^'name') x := lhs(args[next_arg]); q, nn := op(rhs(args[next_arg])); if not proper_name(x, {constants}) then error "The %-1 argument is wrong, an independent name cannot be '%2'", next_arg, x; elif not proper_name(q, {x, cases[]}) or not q::{'name', 'rational'} or q = -1 then error "The base for the 'qdifference' case is wrong, it cannot be %1", q; elif not proper_name(nn, {x, q, constants, cases[]}) then error "In the %-1 argument, the exponent %2 is wrong", next_arg, nn; elif nargs < next_arg+1 or args[next_arg+1]<>'qdifference' then error "The %-1 argument 'qdifference' is expected", next_arg+1; end if; end if; next_arg := next_arg + 1; pp['var'] := x; userinfo(3, LFS, "The independent variable is", x); if m = 1 then pp['functions'] := Vector['column']([y(x)]); else pp['functions'] := Vector['column']([seq(y[i](x), i = 1..m)]); end if; if nargs < next_arg or (not args[next_arg] in cases and not args[next_arg]::'qdifference'['anything']) then err_str := "The %-1 argument (one of "; for i in cases[1..-2] do err_str := cat(err_str, "'", i, "', "); end do: err_str := cat(err_str, "or '", cases[-1], "') is expected"); error err_str, next_arg; elif args[next_arg] = 'differential' then userinfo(3, LFS, "The case of the system is 'differential'"); pp['Ring'] := OreTools:-SetOreRing(x, 'differential'); pp['case'] := args[next_arg]; elif args[next_arg] = 'difference' then userinfo(3, LFS, "The case of the system is 'difference'"); pp['Ring'] := OreTools:-SetOreRing(x, 'shift'); pp['case'] := args[next_arg]; else if args[next_arg] = 'qdifference' and not assigned(q) then q := indets({A, pp['RHS']}, 'name'); q := select(proper_name, q, {x, cases[]}); if nops(q) <> 1 then error "The base for the 'qdifference' case is not recognized, use the %-1 argument in the form 'qdifference'[base]", next_arg; end if; q := q[]; elif args[next_arg]::'qdifference'['anything'] then q := op(args[next_arg]); if not q::{'name','rational'} or not proper_name(q, {x, cases[]}) or q in {1, -1} then error "The base for the 'qdifference' case is wrong, it cannot be %1", q; end if; end if; userinfo(3, LFS, "The case of the system is 'qdifference' with base", q); pp['Ring'] := SetQDifferenceRing(x, q); if assigned(nn) then pp['Ring'] := subsop(2 = ['qdifference', nn][], pp['Ring']); end if; pp['q_par'] := q; pp['case'] := 'qdifference'; end if; next_arg := next_arg + 1; m2 := m; elif args[next_arg]::'algebraic' or args[next_arg]::`=` and not lhs(args[next_arg]) in {'Dummy', 'rhs_normal'} then try S := `if`(args[next_arg]::`=`, (lhs-rhs)(args[next_arg]), args[next_arg]); catch: error "The form of the %-1 argument (a system) is not recognized", next_arg; end try; userinfo(3, LFS, "The system is in the equation form"); next_arg := next_arg + 1; if nargs < next_arg or args[next_arg]::`=` and lhs(args[next_arg]) in {'Dummy', 'rhs_normal'} then error "The %-1 argument (an unknown) is expected", next_arg; elif not args[next_arg]::'function'('name') then error "The %-1 argument of type function(name) is expected", next_arg end if; y := op(0, args[next_arg]); x := op(args[next_arg]); if not proper_name(x, {constants}) then error "The %-1 argument is wrong, an independent name cannot be '%2'", next_arg, x; elif not proper_name(y, {x, constants}) then error "The %-1 argument is wrong, %2 cannot be used for an unknown", next_arg, y(x); end if; pp := table({ 'source_system' = args[next_arg-1], 'var' = x}); userinfo(3, LFS, "The unknown function is", y(x)); next_arg := next_arg + 1; try shifts := check_unknown(S, y, x); catch: error; end try; S := convert(S + dummy, 'list'); pp['explicit_form'] := table('sparse'); n := -1; pp['RHS'] := 0; for term in S do if term = dummy then elif term::dependent(y) then indts := indets(term, `.`); indts := select(has, indts, y); if nops(indts) > 1 then error "The system's term %1 is wrong", term; elif nops(indts) = 1 then dot_mul := indts[1]; if nops(dot_mul) <> 2 then error "The `.`-element of two arguments is expected, but %1 has %2 arguments", dot_mul, nops(dot_mul); elif not has(op(2, dot_mul), y) or not term::'linear'(dot_mul) or not eval(term, dot_mul=dummy):: 'monomial'('anything',dummy) then error "The system's term %1 is wrong", term; end if; cf := coeff(term, dot_mul, 1) * op(1, dot_mul); if cf::'Matrix' then indts := {entries(cf, 'nolist')}; elif cf::'string' then error "The system has a wrong coefficient \"%1\"", cf; elif not cf::'scalar' then error "The system has a wrong coefficient %1", eval(cf); else indts := {cf}; end if; indts := remove(proper_expr, indts, {y}); if indts <> {} and indts[1]::'string' then error "The system has a wrong coefficient \"%1\"", indts[1]; elif indts <> {} then error "The system has a wrong coefficient %1", indts[1]; end if; y_mul := op(2, dot_mul); elif term::`*` then indts := convert(term, 'list'); indts, cf := selectremove(has, indts, y); cf := `*`(cf[]); if nops(indts) <> 1 then error "The system has a wrong term %1", term; elif not proper_expr(cf, {y}) then error "The system has a wrong coefficient %1", cf; end if; y_mul := indts[1]; else y_mul, cf := term, 1; end if; if y_mul = y(x) then k := 0; elif y_mul::{`^`, `+`} then error "The system's term %1 is wrong", term; elif not assigned(pp['Ring']) then try k, pp['Ring'] := what_is_case(y_mul, y, x, shifts, args[next_arg..-1]); catch: error; end try; if OreTools:-Properties:-GetRingName(pp['Ring']) in {'qdifference', 'qshift'} then pp['q_par'] := (OreTools:-Properties:-GetSigma( pp['Ring']))(x, x)/x; userinfo(3, LFS, "The case of the system is 'qdifference' with base %1", pp['q_par']); pp['case'] := 'qdifference'; elif OreTools:-Properties:-GetRingName(pp['Ring']) = 'shift' then userinfo(3, LFS, "The case of the system is 'difference'"); pp['case'] := 'difference'; else pp['case'] := 'differential'; userinfo(3, LFS, "The case of the system is '%1'", OreTools:-Properties:-GetRingName(pp['Ring'])); end if; else try k := OreTools:-Converters:-FromLinearEquationToOrePoly( y_mul, y, pp['Ring']); catch: userinfo(3, 'LFS', StringTools[FormatMessage](lastexception[2..-1])); k := OrePoly(0); end try; if nops(k) > 0 and op(-1, k) = 1 and {op(1..-2, k)} = {0} then k := nops(k)-1; else error "The system has a wrong element %1", y_mul; end if; end if; if not assigned(m) and cf::'Matrix' then m, m2 := LinearAlgebra:-Dimension(cf); if m = 0 or m2 = 0 then error "System's coefficients have to be non-empty matrices"; elif m <> m2 and {indices(pp['explicit_form'])} <> {} then error "System's coefficients have to be matrices with equal dimensions"; elif pp['RHS'] <> 0 and m <> 1 then error "The right-hand side is wrong, a Vector with %1 rows is expected", m; elif m = m2 then pp['explicit_form'] := map(el -> el*LinearAlgebra:-IdentityMatrix(m, compact=false), pp['explicit_form']); end if; pp['RHS'] := Vector(m, [pp['RHS']]); userinfo(3, LFS, "The number of equations is %1", m); userinfo(3, LFS, "The number of unknowns is %1", m2); elif assigned(m) and cf::'Matrix' and [m, m2] <> [LinearAlgebra:-Dimension(cf)] then error "System's coefficients have to be matrices with equal dimensions"; elif assigned(m) and not cf::'Matrix' and m <> m2 then error "System's coefficients have to be matrices with equal dimensions"; elif assigned(m) and not cf::'Matrix' then cf := cf*LinearAlgebra:-IdentityMatrix(m, compact=false); end if; pp['explicit_form'][k] := pp['explicit_form'][k] + cf; n := max(n, k); elif term::'Vector[column]' then # not dependent(term,y) indts := {entries(term, 'nolist')}; indts := remove(proper_expr, indts, {y}); if indts <> {} and (-indts[1])::'string' then error "The system has a wrong right-hand side \"%1\"", -indts[1]; elif indts <> {} then error "The system has a wrong right-hand side %1", -indts[1]; end if; if not assigned(m) then m := LinearAlgebra:-Dimension(term); if m = 0 then error "The right-hand side has to be non-empty vector"; end if; m2 := m; userinfo(3, LFS, "The number of equations is %1", m); userinfo(3, LFS, "The number of unknowns is %1", m2); pp['explicit_form'] := map(el -> el*LinearAlgebra:-IdentityMatrix(m, compact=false), pp['explicit_form']); pp['RHS'] := -term; elif LinearAlgebra:-Dimension(term) <> m then error "The right-hand side is wrong, a Vector with %1 rows is expected", m; else pp['RHS'] := pp['RHS'] - term; end if; # not term::'Vector[column]' elif assigned(m) and m <> 1 then error "The right-hand side is wrong, a Vector with %1 rows is expected", m; elif not proper_expr(term, {y}) then if (-term)::'string' then error "The system has a wrong right-hand side \"%1\"", -term; else error "The system has a wrong right-hand side %1", -term; end if; elif assigned(m) then pp['RHS'] := pp['RHS'] - ; else # not assigned(m) pp['RHS'] := pp['RHS'] - term; end if; end do; pp['order'] := n; userinfo(3, LFS, "The order of the system is %1", n); if not assigned(m) then m := 1; m2 := 1; userinfo(3, LFS, "The number of equations is %1", m); userinfo(3, LFS, "The number of unknowns is %1", m); pp['explicit_form'] := map(el -> Matrix(1, 1, 'fill' = el), pp['explicit_form']); pp['RHS'] := ; end if; if m2 = 1 then pp['functions'] := Vector['column']([y(x)]); else pp['functions'] := Vector['column']([seq(y[i](x), i = 1..m2)]); end if; pp['explicit_form'] := Array(0 .. n, {indices(pp['explicit_form'], 'pairs')}, fill=Matrix(m, m2)); elif args[next_arg]::`=` then error "The %-1 argument (a system) is expected", next_arg; else error "The form of the %-1 argument (a system) is not recognized", next_arg; end if; if not assigned(pp['homogeneous']) then if LinearAlgebra:-Equal(pp['RHS'], Vector(m)) then userinfo(3, LFS, "The system is homogeneous"); pp['homogeneous'] := true; end if; end if; opts := [args[next_arg..-1]]; if not opts::'list'(`=`) then error "optional arguments are expected to be of type equation, but received %0", remove(type, opts, `=`)[]; end if; if hasoption(opts, 'Dummy', 'i', 'opts') then S := args[1..next_arg-1], opts[]; else S := args[1..next_arg-1], 'Dummy'=true, opts[]; end if; if hasoption(opts, 'rhs_normal', 'norm', 'opts') then if pp['homogeneous'] = true then pp['RHS_normal'] := pp['RHS']; else try pp['RHS_normal'] := norm(copy(pp['RHS']), pp['Ring']); catch: error cat("The right-hand side normalizing: ", StringTools[FormatMessage](lastexception[2..-1])); end try; if pp['RHS_normal']::'Vector[column]' and LinearAlgebra:-RowDimension(pp['RHS_normal']) = m then pp['homogeneous'] := LinearAlgebra:-Equal( pp['RHS_normal'], Vector(m)); else pp['homogeneous'] := false; end if; userinfo(3, LFS, "The right-hand side is normalized"); if pp['homogeneous'] then userinfo(3, LFS, "The system is homogeneous"); else userinfo(3, LFS, "The system is inhomogeneous"); end if; end if; elif pp['homogeneous'] <> true then pp['homogeneous'] := false; userinfo(3, LFS, "The system is inhomogeneous"); end if; if pp['case'] = 'qdifference' and pp['input'] <> 'matrices' then hasoption(opts, x, 'indts', 'opts'); end if; if opts <> [] then error "Invalid optional arguments; first one is '%1'", opts[1]; end if; procname(S) := pp; return pp; end proc: what_is_case := proc(term, y, x, shift_indts) local k, Ring, indt, a, q, n, alpha, slv, opts; Ring := OreTools:-SetOreRing(x, 'differential'); try k := OreTools:-Converters:-FromLinearEquationToOrePoly(term, y, Ring); catch: k := OrePoly(0); end try; if nops(k) > 0 and op(-1, k) = 1 and {op(1..-2, k)} = {0} then return nops(k)-1, Ring; end if; Ring := OreTools:-SetOreRing(x, 'shift'); try k := OreTools:-Converters:-FromLinearEquationToOrePoly(term, y, Ring); catch: k := OrePoly(0); end try; if nops(k) > 0 and op(-1, k) = 1 and {op(1..-2, k)} = {0} then return nops(k)-1, Ring; end if; if op(0, term) <> y or has(op(term)/x, x) then error "The case of the system is not recognized by %1", term; end if; opts := select(type, [args[5..-1]], `=`); if hasoption(opts, x=('anything'^'name'), 'indt', 'opts') then q, n := op(indt); if not proper_name(q, {x, y}) or not q::{'name', 'rational'} or q = -1 then error "The base for the 'qdifference' case is wrong, it cannot be %1", q; end if; if not proper_name(n, {x, y, q, constants}) then error "The %1 option is wrong", x=indt; end if; Ring := SetQDifferenceRing(x, q); Ring := subsop(2 = ['qdifference', n][], Ring); OreTools:-Converters:-AddConversionRule('qdifference', eval(qdifferenceToOrePoly)); try k := OreTools:-Converters:-FromLinearEquationToOrePoly( term, y, Ring); catch: userinfo(3, 'LFS', StringTools[FormatMessage](lastexception[2..-1])); k := OrePoly(0); end try; if nops(k) > 0 and op(-1, k) = 1 and {op(1..-2, k)} = {0} then k := nops(k)-1; else error "The system has a wrong element %1 which is not (%2)-difference operator", term, q; end if; return k, Ring; end if; q := op(term)/x; if q::('anything'^'integer') then q, k := op(1, q), op(2, q); else k := 1; end if; if q = 0 or not proper_name(q, {x, y}) or not q::{'name', 'rational'} or q = -1 then error "The system has a wrong element %1", term; elif q::'name' then Ring := SetQDifferenceRing(x, q); OreTools:-Converters:-AddConversionRule('qdifference', eval(qdifferenceToOrePoly)); return k, Ring; end if; # q::'rational' for indt in shift_indts minus {term, y(x)} do a := op(indt)/x; if a = q then next; elif a::('anything'^'integer') then a := op(1, a); end if; if a::'name' then if has(a, {x, constants, y}) or has(y, a) then error "The system has a wrong element %1", indt; end if; error "The qdifference system with one base is expected but %1 and %2 have different bases", term, indt; elif not a::'rational' or a = -1 or a = 0 then error "The system has a wrong element %1", indt; elif a < 0 then error "Use an optional argument in the form %1=base^name to set base for the 'qdifference' case", x; end if; slv := solve(a = q^alpha, alpha); if slv = NULL then error "Use an optional argument in the form %1=base^name to set base for the 'qdifference' case", x; end if; slv := simplify(slv); if not slv::'rational' or slv < 0 then error "Use an optional argument in the form %1=base^name to set base for the 'qdifference' case", x; end if; if denom(slv) <> 1 then q := simplify(q^(1/denom(slv))); k := k*denom(slv); end if; end do; Ring := SetQDifferenceRing(x, q); OreTools:-Converters:-AddConversionRule('qdifference', eval(qdifferenceToOrePoly)); return k, Ring; end proc: qdifferenceToOrePoly := proc (expr, f, R) local x, q, a, k; x := OreTools:-Properties:-GetVariable(R); q := (OreTools:-Properties:-GetSigma(R))(x, x)/x; if expr = f(x) then return OrePoly(1); elif op(0, expr) <> f or nops(expr) <> 1 then error "%1 is not (%2)-difference operator", expr, q; end if; a := op(expr)/x; if a = q then return OrePoly(0, 1); elif a = 0 then error "%1 is not (%2)-difference operator", expr, q; elif q::'rational' and q > 0 and a::'rational' and a > 0 then k := simplify(log[q](a)); elif q::'rational' and q < 0 and a::'rational' and a < 0 then k := simplify(log[-q](-a)); if k::'posint' and k::'odd' then return OrePoly(0$k, 1); else error "%1 is not (%2)-difference operator", expr, q; end if; elif q::'rational' and q < 0 and a::'rational' and a > 0 then k := simplify(log[-q](a)); if k::'posint' and k::'even' then return OrePoly(0$k, 1); else error "%1 is not (%2)-difference operator", expr, q; end if; elif q::'name' and a::`^`(identical(q), 'posint') then k := op(2, a); else error "%1 is not (%2)-difference operator", expr, q; end if; if not k::'posint' then error "%1 is not (%2)-difference operator", expr, q; end if; return OrePoly(0$k, 1); end proc: check_unknown := proc(S, y, x) local indts, indts1; indts := indets(S, 'specfunc'(y)); if indts = {} then error "A linear system for %1 is expected", map(el->el(x), y); end if; indts, indts1 := selectremove(el -> nops(el) <> 1, indts); if nops(indts) = 1 then error "The system has a wrong element %1", indts[1]; elif nops(indts) > 1 then error "The system has wrong elements %1", indts; end if; return indts1; end proc: HGNormal := proc(b::Vector[column], OR) local case, m, opts, IsHypergeometricTerm, RNF, x, t1, q, n, define, b1, b2, i, j, term, InitPoint, Field, Nrml, indts; m := LinearAlgebra:-RowDimension(b); if not OR::'UnivariateOreRing' then return []; end if; x := OreTools:-Properties:-GetVariable(OR); opts := [args[3..-1]]; if not opts::'list'(`=`) then error "optional arguments are expected to be of type equation, but received %0", remove(type, opts, `=`)[]; end if; define := false; hasoption(opts, 'define_cf', 'define', 'opts'); hasoption(opts, 'constant_field', 'Field', 'opts'); hasoption(opts, 'Normal', 'Nrml', 'opts'); indts := {}; hasoption(opts, 'impossible_exponent', 'indts', 'opts'); if opts <> [] then error "Invalid optional arguments; first one is '%1'", opts[1]; end if; case := OreTools:-Properties:-GetRingName(OR); if case = 'differential' then if not assigned(Field) then Field := 'rational'; end if; if not assigned(Nrml) then Nrml := dualnormal; end if; IsHypergeometricTerm := proc(expr, x, R); DEtools['IsHyperexponential'](expr, x, R); end proc; RNF := proc(R, x) local F, V; F, V := DEtools['RationalCanonicalForm'][1](R, x); return V, F; end proc; InitPoint := proc(R, x) local i, Ri; if R = 0 then return {x = 0} end if; Ri := 0; for i from 0 while Ri = 0 do try Ri := eval(R, x = i); catch: end try; end do; return {x = i-1}; end proc; elif case = 'shift' then if not assigned(Field) then Field := 'rational'; end if; if not assigned(Nrml) then Nrml := proc(expr) dualnormal(value(expr)) end proc; end if; IsHypergeometricTerm := proc(expr, x, R); SumTools:-Hypergeometric:-IsHypergeometricTerm(expr, x, R); end proc; RNF := proc(R, x1) local z, r, s, u, v, F, V; z, r, s, u, v := SumTools:-Hypergeometric:-RationalCanonicalForm(R, x1); V := u/v; F := z*r/s; return V, F; end proc; InitPoint := proc(R, x) local Ri; Ri := LFS:-Solve(numer(R)*denom(R), x); if Ri = [] then return {x = 0} end if; Ri := {Ri[1]-1, map(`+`, Ri, 1)[]} minus {Ri[]}; return map2(`=`, x, Ri); end proc; elif case in {'qshift', 'qdifference'} then q := OreTools:-Properties:-GetSigma(OR)(x, x)/x; indts := indts minus {q}; if nops(OR) = 3 then n := op(3, OR); if indts <> {} and not proper_name(n, indts) then error "For %1, the exponent %2 is impossible, this name is used in the coefficients of the system", x=q^n, n; end if; else t1 := indets(b, '`^`'), indets(b, 'specfunc'({'Product','product'})); t1 := map2(op, 2, t1[1]), map2(op, 1, t1[2]), map2(op, [2, 2], t1[2]); t1 := indets([t1[1], t1[3]], 'name') minus indets(t1[2], 'name') minus {x, q, cases[], constants, 'shift', 'qshift'}; if nops(t1) > 1 or nops(t1) = 1 and indts <> {} and not proper_name(t1[1], indts) then error "The variable n (where %1) is not recognized", x=q^n; elif nops(t1) = 1 then n := t1[1]; end if; end if; IsHypergeometricTerm := SetQIsHypergeometricTerm(q, n); RNF := SetQRNF(q); InitPoint := SetQInitPoint(q, n); if not assigned(Field) then if not q::'name' then Field := 'rational'; else Field := 'ratpoly'('rational', q); end if; end if; if not assigned(Nrml) then Nrml := proc(expr) dualnormal(value(expr)) end proc; end if; elif not case in ['qdifference', 'shift', 'qshift'] then IsHypergeometricTerm := parse(sprintf( "proc() error \"The case '%a' is not determined\" end proc", case)); RNF := parse(sprintf( "proc() error \"The case '%a' is not determined\" end proc", case)); end if; for i to LinearAlgebra:-Dimension(b) do try b1[i], Field := HGNormal_step1(b[i], x, IsHypergeometricTerm, InitPoint, RNF, Nrml, Field, define); catch: error; end try; end do; t1 := map(op, convert(b1, 'list')); # step2 b2 := []; for i to m do for term in b1[i] do for j to nops(b2) do if b2[j][3] = term[3] then t1 := 1, OreTools:-Properties:-GetTheta1(OR); elif OreTools:-Properties:-Getdelta(OR)(x, x) = 0 then t1 := RNF(term[3]/b2[j][3], x); if t1[2] <> 1 then next; end if; else next; end if; try t1 := [t1[], map2(eval, Nrml(term[1]/b2[j][1]/t1[1]), b2[j][2])]; t1[3] := map(Nrml, t1[3]); catch: map2(Limit, term[1]/b2[j][1]/t1[1], b2[j][2]); next; end try; if nops(t1[3]) = 1 and t1[3][1]::Field and t1[3][1] <> 0 then b2[j][4][i] := b2[j][4][i] + t1[3][1]*term[2]*t1[1]; break; end if; end do; if j > nops(b2) then b2 := [b2[], [term[1], term[4], term[3], Vector(m, {i=term[2]})]]; end if; end do; end do; map(el -> LinearAlgebra:-Map(Nrml, el[4]), b2); b2 := remove(el -> LinearAlgebra:-Equal(el[4], Vector(m)), b2); if b2 = [] then b2 := Vector(m); else b2 := map(el-> [el[1],el[3],el[4]], b2), Field; end if; return b2; end proc: HGNormal_step1 := proc(term, x, IsHypergeometricTerm, InitPoint, RNF, Nrml, FIELD, define) local i, j, t1, F, V, R, dummy, indts, pnts, Field, res; Field := FIELD; if term::`+` then t1 := args[2..6]; res, Field := HGNormal_step1(op(1, term), t1, Field, define); t1, Field := HGNormal_step1(subsop(1=0, term), t1, Field, define); return [res[], t1[]], Field; elif term = 0 then return [], Field; elif term::'string' then error "\"%1\" is not hypergeometric", term; elif not term::'scalar' then error "%1 is not hypergeometric", term; end if; try t1 := IsHypergeometricTerm(term, x, 'R'); catch "The case": error; catch: t1 := false; end try; if not t1 then try t1 := expand(term); catch: t1 := term; end try; if t1 = term then error "%1 is not hypergeometric", term; elif t1 = 0 then return [], Field; end if; try t1 := procname(t1, args[2..-1]); catch: error "%1 is not hypergeometric", term; end try; return t1; end if; R := R; if define then indts := remove(type, indets(R, 'name') minus {x}, Field); if remove(proper_name, indts, {}) <> {} then error "A hypergeometric term in %1 over a field '%2' is expected but received %3", x, Field, term; elif indts <> {} then if Field = 'rational' then Field := 'ratpoly'('rational', indts); elif op(2, Field)::'set' then Field := 'ratpoly'('rational', op(2, Field) union indts); else Field := 'ratpoly'('rational', {op(2, Field)} union indts); end if; end if; end if; if not R::'ratpoly'(Field, x) then error "A hypergeometric term in %1 over a field '%2' is expected but received %3", x, Field, term; end if; pnts := [InitPoint(R, x)[]]; for i to nops(pnts) do try t1 := eval(term, pnts[i]); catch: try t1 := limit(term, pnts[i]); catch: t1 := FAIL; end try; end try; if t1::'SymbolicInfinity' or t1::'undefined' or t1 = exp(FAIL) or op(0, t1) = 'limit' and op(2, t1) = pnts[i] then t1 := FAIL; end if; try t1 := Nrml(t1); catch: t1 := FAIL; end try; if not t1::'scalar' then error "%1 is not hypergeometric", term; end if; pnts[i] := [pnts[i], t1]; end do; pnts := {pnts[]}; if map2(op, 2, pnts) = {FAIL} then error "%1 is not hypergeometric", term; end if; try V, F := RNF(R, x); catch: error; end try; t1 := term/V*dummy[1](x)*dummy[2](x); if define then t1 := selectremove(type, t1, 'ratpoly'('anything', x)); else t1 := selectremove(type, t1, 'ratpoly'(Field, x)); end if; try t1 := Nrml(t1[1])*t1[2]/dummy[1](x)/dummy[2](x); catch: error; end try; return [[t1, V, F, map2(op, 1, pnts)]], Field; end proc: ############################################################################### Solve := proc(P, x) local res, indts, cnt; if P::'polynom'('rational', x) then res := isolve(P); if res = {} or res = NULL then return []; else res := map2(op, [1,2], [res]); return sort(res); end if; elif P::'polynom'('anything', x) then indts := indets(P, {'name', 'RootOf'}); if P::'polynom'('rational', indts) then res := content(P, indts minus {x}); return procname(res, x); elif P::'ratpoly'('rational', indts) then res := content(numer(P), indts minus {x}); return procname(res, x); else res := primpart(P, x, 'cnt'); if cnt <> 1 then return procname(res, x); else WARNING("find all integer roots %1 for %2, some roots can be lost", x, P); end if; end if; end if; res := [{solve(P, x)}[]]; res := simplify(res, symbolic); if nops(res) = 1 and res[1]::RootOf then indts := indets(P, 'anything'^identical(x)); if indts <> {} then return qSolve(P, xx, op(1, indts[1]), x); end if; end if; return sort(select(type,res,'integer')); end proc: qSolve := proc(P, x, q, n) local e1, s, res, t, i; e1 := eval(P, x = q^n); res := [solve(e1)]; if res::list(set) then res := map(el -> `if`(lhs(el[1]) = n and rhs(el[1]) <> n, rhs(el[1]), `if`(lhs(el[2]) = n and rhs(el[2]) <> n, rhs(el[2]), NULL)), res); end if; res := simplify(res, symbolic); return sort(select(type,res,'integer')); end proc; ############################################################################### RNFQShift := proc(R1, x1, q) local z, r, s, u, v, V, F, alpha; z, r, s, u, v := QDifferenceEquations:-QRationalCanonicalForm[1](R1, q, x1); z := dualnormal(expand(z)); alpha := ldegree(numer(z), q) - ldegree(denom(z), q); z := dualnormal(z/q^alpha); V := x1^alpha*u/v; F := z*r/s; return V, F; end proc: InitPointQShift := proc(R, x, q, n) local Ri; Ri := LFS:-qSolve(numer(R)*denom(R), x, q, n); if Ri = [] then return {{x = 1, n = 0}}; end if; Ri := {Ri[1]-1, map(`+`, Ri, 1)[]} minus {Ri[]}; return map(el -> {x = q^el, n = el}, Ri); end proc; ############################################################################### # e.g., IsqHypergeometric(h,n,q^n=N,R) IsQHypergeometricTerm := proc(h,n,sub,RR) description "Check to see if the given expression is a q-hypergeometric term"; local R, dir, rat, i, tf; if not type(args[3],anything^name=name) then error "wrong type of arguments" elif hastype(h, tabular) then return false; end if; dir := 'up'; rat := false; for i from 4 to nargs do if member(args[i],{'up','down'}) then dir := args[i] elif type(args[i],'name') then rat := true end if; end do; (tf,R) := qratio(h,n,sub,dir); if tf then if rat then RR := R end if; return true else return false end if; end proc: # e.g., qratio(f,k,q^k=K,{'up','down'}); qratio := proc(f,k,sub,dir) local prod, rest, tf, sol, tf1, sol1, expr, R, i; kernelopts(opaquemodules = false); if has(f,{'Product','product'}) then if type(f,'specfunc'(anything,{'Product','product'})) then return qratio_prod(args) elif type(f,'specfunc'(anything,{'Product','product'})^integer) then (tf,sol) := qratio_prod(op(1,f),args[2..nargs]); if not tf then return (false, f) else return (true, sol^(op(2,f))) end if; elif type(f,`*`) then (prod,rest) := selectremove(type,f, {'specfunc'(anything,{'Product','product'}), ('specfunc'(anything,{'Product','product'}))^'integer'}); if type(prod,'specfunc'(anything,{'Product','product'})) then (tf,sol) := qratio_prod(prod,args[2..nargs]); if not tf then return (false, f) else (tf1,sol1) := procname(rest,args[2..nargs]); if not tf1 then return (false, f) else return (true, sol*sol1) end if; end if; elif prod <> 1 then for i to nops(prod) do (tf,sol[i]) := procname(op(i,prod),args[2..nargs]); if not tf then return (false, f) end if; end do; (tf1,sol1) := procname(rest,args[2..nargs]); if not tf1 then return (false, f) else return (true, `*`(seq(sol[i],i=1..nops(prod)))*sol1) end if; else for i to nops(rest) do (tf,sol[i]) := procname(op(i,rest),args[2..nargs]); if not tf then return (false, f) end if; end do; return (true, `*`(seq(sol[i],i=1..nops(prod)))) end if; else expr := QDifferenceEquations:-QSimplify(f); if expr = f then return false, f; end if; return procname(expr,args[2..nargs]) end if; end if; if dir = 'up' then expr := subs(k=k+1,f)/f else # dir = 'down' expr := f/subs(k=k-1,f) end if; R := SubsPower(sub,QDifferenceEquations:-QSimplify(QDifferenceEquations:-QSimpComb(expr))); if op([1,1], sub)::'name' and type(R,'ratpoly'(anything, {op(2,sub), op([1,1], sub)})) and not depends(R, k) then return true, R elif not op([1,1], sub)::'name' and type(R,'ratpoly'(anything, op(2,sub))) and not depends(R, k) then return true, R else return false, f end if; end proc: qratio_prod := proc(f,k,sub,dir) local func, ind, k0, ab, a, b, R; kernelopts(opaquemodules = false); func := op(1,f); ind := op([2,1],f); k0 := op([2,2,1],f); ab := op([2,2,2],f); if type(k0,'integer') and type(ab,'polynom'('integer',k)) and degree(ab,k)=1 then a := coeff(ab,k,1); b := coeff(ab,k,0); R := SubsPower(sub,eval(func,ind=k)); if op([1,1], sub)::name and not type(R,'ratpoly'(anything, {op(2,sub), op([1,1], sub)})) or not type(R,'ratpoly'(anything, op(2,sub))) or depends(R, k) then false, f elif not op([1,1], sub)::'name' and not type(R,'ratpoly'(anything, op(2,sub))) then false, f else if dir = 'up' then if a > 0 then if a = 1 and b = -1 then return true, R else return true, SubsPower(sub,product(func,ind=a*k+b+1..a*k+b+a)) end if; else return true, 1/SubsPower(sub,product(func,ind=a*k+b+a+1..a*k+b)) end if; else if a > 0 then if a = 1 and b = -1 then return true, SubsPower(sub,eval(func,ind=k-1)) else return true, SubsPower(sub,product(func,ind=a*k+b-a+1..a*k+b)) end if; else return true, 1/SubsPower(sub,product(func,ind=a*k+b+1..a*k+b-a)) end if; end if; end if; else return false, f end if; end proc: module() local recursive, `power/subs/sub`, `power/subs/subnames`; SubsPower := proc(subst, expr) local sub, subnames, f, z; if type(subst,equation) then sub := {subst} elif type(subst,{list(equation),set(equation)}) then sub := {op(subst)} else error "wrong type of arguments" end if; f:= subs(sub, expr); # Now we check if all powers were substituted... sub := select(type, sub, anything^name=anything); subnames := map(z -> op(2,lhs(z)), sub) intersect indets(f,name); if subnames <> {} then sub := select(unapply('has'('op'(1,z),subnames),z), sub); `power/subs/sub`:= [seq([op(op(1,z)),op(2,z)],z=sub)]; `power/subs/subnames`:= [seq(op(1,op(1,z)),z=sub)]; # Determine all power substitutions via recursive subnames := indets(f,anything^anything) minus\ indets(f,{anything^(-1),name^name,name^integer,\ (name^name)^integer,identical(I)}); f := eval(f,[seq(z=recursive(Normalizer(op(1,z)),\ expand(normal(op(2,z)))), z=subnames)]); end if; return f; end proc: recursive := proc(b,e) local z; if type(e,`+`) then return `*`(seq(procname(b,z),z=e)) elif type(b,`*`) then return map(procname,b,e); elif type(b,`^`) then return procname(op(1,b),expand(normal(op(2,b)*e))); elif member(b,`power/subs/subnames`,'z') then z := `power/subs/sub`[z]; if denom(numer(e)/z[2]) = 1 then return expand(z[3]^(e/z[2])); else return expand(b^e); end if; elif hastype(b,`^`) and has(b,`power/subs/subnames`) then return subs([seq(z=procname(op(z)),\ z=indets(b,`^`) minus indets(b,anything^(-1)))],b)^e; else return expand(b^e); end if; end proc: end module: SolveRec := module() local ModuleApply, DegreeBound, ConstructSolutionByTransformed, Polynomial, PolynomialCleanup, ResolveCeqs, SolveSystemV, SolveSystemW, _C::nothing; DegreeBound := proc (case, M0, sblk, ncol, nlin, t, x, deg0, q) local deg, p, rr, MM, TM, maxn, cns, vert_shft, k, i, j, mm; deg := -1; # Triangualize the trailing matrix and bound the degree of # the polynomial solution (by methods) cns := table(); vert_shft := Vector(sblk); EG_('trail', M0, ncol-1, sblk, nlin, x, OreTools:-SetOreRing(x, 'shift'), table({(1)=cns, (2)=vert_shft})); if convert(M0[sblk], 'set') = {0} then return FAIL; end if; MM := LinearAlgebra:-SubMatrix(M0, 1..sblk,sblk*(ncol-1)+1..sblk*ncol); p := dualnormal(LinearAlgebra:-LUDecomposition( MM, 'output'='determinant')); rr := Solve(p, x); if nops(rr)>0 then for i from nops(rr) by -1 to 1 do maxn := rr[i]; if assigned(cns[maxn]) then TM := Matrix([ [map(dualnormal,eval(MM,x=maxn))], seq([Vector['row'](sblk, [seq(cns[maxn][j][k], k=sblk*(ncol-1)+1..sblk*ncol)])], j=1..nops(cns[maxn])) ]); if LinearAlgebra:-Rank(TM)-infinity then deg := max(deg, maxn+max(seq(vert_shft[mm],mm=1..sblk))+t); end if; end if; return deg, [M0, cns, vert_shft, {}, {}], rr; end proc; Polynomial := proc(equations, unknownsV, unknownsW, cn) local eqns, vars, sols, unknowns, solV, solW, consts, k; unknowns := unknownsV union unknownsW; # Eliminate trivial equations simultaneously vars := unknowns; eqns := equations; sols := {}; while nops(vars)>0 do vars := eqns intersect unknownsV; sols := sols union map(`=`, vars, 0); eqns := subs(sols, equations); # Simplify equations as much as possible first eqns := map(PolynomialCleanup, eqns, unknowns); end do; solV, eqns := SolveSystemV(eqns minus {0}, unknownsV); # return FAIL if system W is not consistent solW := SolveSystemW(eqns minus {0}, unknownsW); sols := sols union solV; vars := unknownsV minus map(lhs, sols); consts := {seq(vars[k]=cn[k], k=1..nops(vars))}; return (eval(sols, consts) union consts,solW); end proc: SolveSystemV := proc(equations,unknowns) local eqns, k, eqn, i, sols, pivot, t, divtab, sol, var, vars; eqns := equations; for k while eqns <> {} do # Select the simplest equation to work with eqn := select(has, eqns, unknowns); eqn := ListTools:-FindMinimalElement([op(eqn)], (a,b)->length(a) {} do eqn := eqns[1]; for i from 2 to nops(eqns) while not type(eqn,'name') do # Select the simplest equation to work with if length(eqns[i]) < length(eqn) then eqn := eqns[i] end if end do; vars := indets(eqn) intersect unknowns; if vars = {} then return FAIL; else var := vars[1]; end if; eqns := eqns minus {eqn}; if not type(eqn,'polynom'('integer',vars)) then eqn := collect(eqn,vars); end if; pivot := coeff(eqn,var,1); for i from 2 to nops(vars) do # Select the simplest coefficient to eliminate t := coeff(eqn,vars[i],1); if type([pivot,t],'list(integer)') and abs(t) < abs(pivot) then pivot, var := t, vars[i]; elif type(t,`+`) and not type(pivot,`+`) then next; elif length(t) < length(pivot) then pivot, var := t, vars[i]; end if end do; eqn := subs(var=0,eqn); sol[k] := eqn, var, pivot; if pivot = 1 then eqns := map(dualnormal, subs(var=-eqn, eqns)) minus {0}; elif pivot = -1 then eqns := map(dualnormal, subs(var=eqn, eqns)) minus {0}; else divtab := table(); eqns := map( proc(eqn,var,vars,p,e,d) local r,s,t; r := coeff(eqn,var,1); if r = 0 then return eqn; end if; gcd(p, r, 's', 't'); r := subs(var=0, eqn); r := expand(s*r - t*e); if r = 0 then return NULL; elif assigned(d[eqn]) then divide(r,d[eqn],'r'); unassign('d[eqn]'); else r := primpart(r,vars) end if; if not type(s,'integer') then d[r] := s end if; return r; end proc, eqns, var, unknowns, pivot, eqn, 'divtab'); end if end do; sols := {}; for i from k-1 by -1 to 1 do userinfo(2,solve,`backsubstitution at:`,i); var := sol[i][2]; eqn := subs(sols,sol[i][1]); eqn := dualnormal(-1 / sol[i][3] * eqn); sols := sols union {var = eqn}; end do; return sols; end proc: PolynomialCleanup := proc(eqn,vars) local e; if type(eqn,'`^`') then e := 1; elif type(eqn,'`*`') then e := map(thisproc,args); elif type(eqn,'`+`') then e := expand(eqn); elif not type(eqn,'polynom(algnum)') and numer(eqn) <> eqn then e := thisproc(numer(eqn),vars); else e := eqn; end if; if e = 0 then 0 elif not has(e,vars) then 1 elif type(e,'name') then e elif type(e,'`+`') then primpart(e,vars) else thisproc(e,vars) end if end proc: ResolveCeqs := proc(ceqs::set, sol::Array, $) local sbs; if nops(ceqs)>0 then sbs := SolveTools:-Linear(ceqs, select(has, indets(ceqs), _C)); map['inplace'](z->eval(z,sbs), sol); end if; return NULL; end proc; ConstructSolutionByTransformed := proc(deg0, case, M, sblk, ncol, nlin, l, t, x, var, deg_min := 0) local deg, r, rr, maxVertShift, sol, ceqs, z::nothing, eq, m, m1, m2, i, vrs, rhsi, k, eqi, res, ceqs0, cns, mm, c, s, zeros, rec, j; if not case in ['differential','difference','qdifference'] then return NULL; end if; rec := M[1]; cns := eval(M[2]); maxVertShift := max(seq(M[3][mm],mm=1..sblk)); deg := deg0[1]; if type(deg0[2],'list') then rr := deg0[2]; end if; deg := max(deg-maxVertShift, -1); sol := Array(1..sblk, -maxVertShift+deg_min-l+t..deg); eq := Array(1..nlin); for k to nlin do eq[k] := dualnormal(add(rec[k,(ncol-1)*sblk+m]*z[m],m=1..sblk),'expanded'); end do; vrs := [seq(z[m],m=1..sblk)]; rec := LinearAlgebra:-SubMatrix(rec, 1..nlin, 1..sblk*(ncol-1)); rhsi := Vector['column'](nlin); zeros := 0; ceqs := {}; for i from deg-t by -1 to -l-maxVertShift+deg_min do for k to nlin do rhsi[k] := 0; # for m1 from max(1,i+l+1-deg) to min(ncol-1,ncol+i-t+maxVertShift-deg_min) do for m1 from max(1,i+l+1-deg) to ncol-1 do for m from 1 to sblk do if sol[m,ncol-m1+i+t]<>0 then rhsi[k] := rhsi[k]-eval(rec[k,(m1-1)*sblk+m],x=i) *sol[m,ncol-m1+i+t]; end if; end do; end do; rhsi[k] := dualnormal(rhsi[k],'expanded'); if case='qdifference' then rhsi[k] := map(simplify, rhsi[k], symbolic); end if; end do; eqi := [seq(dualnormal(eval(eq[m1],x=i)-rhsi[m1]),m1=1..nlin)]; if case='qdifference' then eqi := map(simplify, eqi, symbolic); end if; res, ceqs0 := Polynomial({op(eqi)},{op(vrs)}, select(has, indets(rhsi), {z, _C}), _C[i]); if assigned(rr) and ceqs0 = {} and map(rhs,res)={0} then zeros := zeros + 1; if zeros >= l-t+1 then k := -l-maxVertShift-1+deg_min; while nops(rr)>=1 do rr, r := rr[1..-2], rr[-1]; if r= -t-maxVertShift+deg_min then for m1 from 1 to nops(res) do m2 := op(lhs(res[m1])); sol[m2,i+t] := rhs(res[m1]); end do; else for m1 from 1 to nops(res) do if rhs(res[m1])<>0 then ceqs := ceqs union {rhs(res[m1])=0}; end if; end do; end if; ResolveCeqs(ceqs, sol); if case='qdifference' then sol := map(simplify, sol, symbolic); ceqs := map(simplify, ceqs, symbolic); end if; end do; ceqs := {}; for i in indices(cns,'nolist') do if (deg>=i+t) and (0<=i+l) then for m2 to nops(cns[i]) do rec := cns[i][m2]; c := 0; for m1 from max(1,ncol+i+t-deg) to min(ncol,ncol+i+t+maxVertShift-deg_min) do c := c + add(rec[(m1-1)*sblk+m]*sol[m,ncol-m1+i+t], m=1..sblk); end do; c := dualnormal(c); if c<>0 then ceqs := ceqs union {c=0}; end if; end do; end if; end do; ResolveCeqs(ceqs, sol); if case='qdifference' then sol := map(simplify, sol, symbolic); ceqs := map(simplify, ceqs, symbolic); end if; if deg_min = 0 then ceqs := {seq(seq(`if`(sol[k,m1]<>0,sol[k,m1]=0,NULL), m1 = -maxVertShift+deg_min..min(deg,-M[3][k]-1)), k = 1..sblk)}; ResolveCeqs(ceqs, sol); if case='qdifference' then sol := map(simplify, sol, symbolic); ceqs := map(simplify, ceqs, symbolic); end if; end if; # Form the solution s := Array(1..sblk); for i to sblk do if case='differential' or case='qdifference' then s[i] := add(var^(k+M[3][i])*sol[i,k], k=-l+t-M[3][i]+deg_min..deg); else # case='difference' s[i] := add(product(var-j,j=0..k-1+M[3][i])*sol[i,k], k=-l+t-M[3][i]+deg_min..deg); end if; end do; s := convert(s,'list'); vrs := select(has, indets(s), _C); return subs(seq(vrs[k]=_c[k],k=1..nops(vrs)), s); end proc; # ConstructSolutionByTransformed ####################################################################### ModuleApply := proc(case, rs, sblk, l, t, x, var, method, q, deg0) local M, nlin, ncol, deg, rr; # Construct the parameters: M := Matrix(rs); nlin := LinearAlgebra:-RowDimension(M); ncol := LinearAlgebra:-ColumnDimension(M) / sblk; # Compute the degree bound: deg := DegreeBound(case, M, sblk, ncol, nlin, t, x, deg0, q); if deg = FAIL then # The system has fewer equations then needed or incompatible if nargs=11 and type(args[11],'name') then assign(args[11], "Not full rank"); end if; return FAIL; end if; deg, M, rr := deg; # It can be changed due to removing dependent rows nlin := LinearAlgebra:-RowDimension(M[1]); # find the solution: return ConstructSolutionByTransformed([deg,rr], case, M, sblk, ncol, nlin, l, t, x, var); end proc; ####################################################################### end module; end module: