# ************************************************************************* # # adaptation of EP code (version 05.11.97) to difference case # # global super-reduction algorithm for systems with rational function # # coefficients (p-adic version). # # # # # # ************************************************************************* # # g_super_reduce(A, x, p, lambda [,T ,invT ] [,k2 ]) # INPUT: A - a square matrix with rational function entries # x - a name # p - an irreducible polynomial or infinity # lambda - a name # T - (optional) a name # invT - (optional) a name # k2 - (optional) a positive integer # OUTPUT: The list of a super-irreducible form of A at the point p and the # rho-polygon of A at p. # If T is passed as argument, T and invT are assigned so that # T[A] := T^(-1) (A &* T - T') is the computed super (k2-) # irreducible form. # g_super_reduce := proc(A, x, p, lambda) local t0, i, p1, dp, s, k, k2, n, mk, T, invT, ntab, valtab, M, G, theta, theta0, rk, P, rhopoly, min_val; global Ti; t0 := time(); if (nargs < 4) or (nargs > 7) or not type(x, 'name') or not type(lambda, 'name') or not type(A, 'matrix'(ratpoly('anything', x), 'square')) then error "wrong number or type of arguments" fi; Ti := 0; k2 := infinity; n := linalg['rowdim'](A); G := linalg['matrix'](n, n, 0); T := Id(n); invT := copy(T); # If p = 1/x (which represents infinity), do change of variable x -> 1/x if p = 1/x then M := normalm(-1/x^2*subs(x = 1/x, eval(A))); p1 := x; for i from 1 to n do valtab[i] := valtab[i]-2 od else error "3rd argument has not correct type" fi; k := 1; userinfo(1, 'super_reduce', "enter super_reduce, reduction step = 1"); dp := 1; rhopoly := NULL; s := infinity; mk := 0; while (k <= min(s-1, k2)) and (mk < n) do userinfo(3, 'super_reduce', `computing G-matrix`); g_sort_cols(M, x, p1, n, k, lambda, ntab, valtab, min_val, G, false, T, invT); userinfo(1, 'super_reduce', `valuations:`, sort(convert(valtab, 'list'))); userinfo(4, 'super_reduce', `ntab:`, convert(ntab, 'list')); s := -valtab[1]; if s = 0 then break fi; mk := add(ntab[i], i=1..k); theta := evala(Rem(linalg['det'](G), p1, x)); userinfo(5, 'super_reduce', `reduction step:`, k, `theta:`, theta); if type(min_val, 'name') then theta0 := rem(linalg['charpoly'](linalg['submatrix'](G, 1..mk, 1..mk), lambda), p1, x); theta0 := evala(Normal( theta0/lambda^ldegree(theta0, lambda))); if degree(theta0, lambda) > 0 then min_val := -s fi fi; while (theta = 0) and (s > 1) do rk := prank(linalg['submatrix'](G, 1..n, 1..mk), x, p1); if rk < mk then pgauss(linalg['submatrix'](G, 1..n, 1..mk), x, p1, 'P'); userinfo(4, 'super_reduce', `columns are linearly dependent`); userinfo(5, 'super_reduce', `kernel dim = `, mk-rk); g_column_rank(M, x, p1, n, k, ntab, valtab, mk, G, P, T, invT); userinfo(3, 'super_reduce', `recomputing G-matrix`); g_sort_cols(M, x, p1, n, k, lambda, ntab, valtab, min_val, G, false, T, invT) else userinfo(4, 'super_reduce', `columns are linearly independent`); g_qtcd(M, x, p1, n, k, min_val, ntab, valtab, mk, lambda, G, T, invT); g_sort_cols(M, x, p1, n, k, lambda, ntab, valtab, min_val, G, false, T, invT) fi; userinfo(1, 'super_reduce', `valuations:`, sort(convert(valtab, 'list'))); userinfo(4, 'super_reduce', `ntab:`, convert(ntab, 'list')); userinfo(4, 'super_reduce', `mk:`, mk); theta := expand(evala(Rem(linalg['det'](G), p1, x))); userinfo(5, 'super_reduce', `reduction step:`, k, `theta:`, theta); s := -valtab[1]; mk := add(ntab[i], i=1..k); if type(min_val, 'name') then theta0 := rem(linalg['charpoly'](linalg['submatrix'](G, 1..mk, 1..mk), lambda), p1, x); theta0 := normal( theta0/lambda^ldegree(theta0, lambda)); if degree(theta0, lambda) > 0 then min_val := -s fi fi; od; theta := primpart(theta, lambda); if k < s-1 then theta := normal(theta/lambda^ldegree(theta, lambda)) fi; if k = 1 then # Moser irreducible if s = 1 then ntab := [n]; userinfo(1, 'super_reduce', `matrix is regular singular at`, p1); rhopoly := [0, primpart(rem(subs(lambda = lambda*dp, theta0), p1, x), lambda)*lambda^(n-degree(theta0, lambda))]; userinfo(1, 'super_reduce', `slope:`, 0, `indicial equation:`, rhopoly[2]); break else # invariant of biggest slope (Katz approximation) theta := primpart(rem(subs(lambda = -dp*lambda, theta), p1, x), lambda); theta0 := primpart(rem( theta0 , p1, x), lambda); rhopoly := [s-2, collect(theta , lambda)], [s-1, collect(theta0, lambda)]; userinfo(1, 'super_reduce', `Katz approximation:`, s-1, `invariant:`, op(2, [rhopoly])[2]); userinfo(1, 'super_reduce', `slope `, s-2, `invariant:`, op(1, [rhopoly])[2]) fi; else theta := primpart(rem(subs(lambda = -dp*lambda, theta), p1, x), lambda); rhopoly := [s-k-1, collect(theta, lambda)], rhopoly; userinfo(1, 'super_reduce', `slope `, s-k-1, `invariant:`, op(1, [rhopoly])[2]) fi; if s > 1 then k:= k+1; min_val := valtab[1]; userinfo(1, 'super_reduce', `increase reduction step to: `, k) fi od; if mk = n then for i from k+1 to s-1 do userinfo(1, 'super_reduce', `increase reduction step to: `, i); rhopoly := [s-i-1, 1], rhopoly; userinfo(1, 'super_reduce', `slope:`, s-i-1, `invariant:`, op(1, [rhopoly])[2]); od fi; userinfo(1, 'super_reduce', `time: `,time()-t0); if s <= 0 then error "Did not know this point is reachable"; fi; userinfo(5, 'super_reduce', `rho-polygon:`, [rhopoly]); [ [seq(max(seq(ldegree(denom(T[i,s]),x)-ldegree(numer(T[i,s]),x), s=1..n)),i=1..n)], [] , [rhopoly]] end: # g_sort_cols(M, x, p, n, k, lambda, ntab, valtab, s, G, flag, # [T, invT]) # INPUT: M - a rational function square matrix # x - a name # p - an irreducible polynomial # n - a positive integer (the dimension of M) # k - a positive integer (the reduction step) # lambda - a name # ntab - a name or an array of integers # valtab - a name or an array of integers # s - a name or an integer # G - a name or a square matrix of dim. n # flag - a boolean # T - (optional) a name or a square matrix of dim. n # invT - (optional) a name or a square matrix of dim. n # OUTPUT: M with the n columns sorted in increasing p-valuations. T and invT # are assigned to the corresponding permutation matrices, G is # assigned to the matrix consisting of the mk leading vectors of # p-valuation s,s-1,..,s-(k-1) and the left n-mk leading vectors of # p-valuation >= s-k. valtab and ntab are modified to the actual # p-valuations of the columns and numbers of columns having same # p-valuations resp. If s is an integer, it is assumed to be the # valuation of the matrix M at p. If flag = true, the first n_0+..+ # n_(k-1) columns of G will not be recomputed. # g_sort_cols := proc(M, x, p, n, k, lambda, ntab, valtab, s, G, flag, T, invT) local d, cc, col, v_mat, v_col, v, i, j, h, m, tmp, M_col, G_col, perm, a, U, V; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if type(s, 'name') then v_mat := infinity else v_mat := s fi; for i from 0 to k do col[i] := NULL od; # t0:=time(); if (k = 1) or ((k > 1) and not flag) then if type(s, 'name') then ntab := array([0, 0]) ##seq(0, i=1..n)]) else ntab := array(1..-s, [seq(0, i=1..-s)]) fi; valtab := array(1..n); G := linalg['matrix'](n, n, 0); for j from 1 to n do v_col := infinity; cc := NULL; for i from 1 to n do v := pval((M[i, j]), x, p); if v < v_col then v_col := v; cc := i elif v = v_col then cc := cc, i fi od; # cc is the seq of coefficients of the leading vector valtab[j] := v_col; if v_mat = infinity then v_mat := v_col; d := k else d := simplify(v_mat-v_col) fi; if d > 0 then for h from k-1 to d by -1 do col[h] := col[h-d] od; for h from d-1 to 1 by -1 do col[h] := NULL od; col[0] := [j, [cc]]; v_mat := v_col elif -d <= k-1 then col[v_col-v_mat] := col[v_col-v_mat], [j, [cc]]; fi; od; m := 0; G_col := 1; perm := NULL; for h from 0 to k-1 do tmp := [col[h]]; for j from 1 to nops(tmp) do M_col := tmp[j][1]; if M_col <> G_col then perm := perm, [M_col, G_col] fi; if p = x then for i in tmp[j][2] do G[i, G_col] := subs(x = 0, normal(x^(-v_mat-h)*M[i, M_col])) od else for i in tmp[j][2] do a := normal(p^(-v_mat-h)*M[i, M_col]); gcdex(denom(a), p, x, 'U', 'V'); G[i, G_col] := rem(numer(a)*U, p, x) od; fi; G_col := G_col+1 od; ntab[h+1] := nops(tmp); m := m+ntab[h+1]; od; ntab[k+1] := n-m; perm := [perm]; for tmp in perm do d := valtab[tmp[1]]; valtab[tmp[1]] := valtab[tmp[2]]; valtab[tmp[2]] := d; d := op(tmp); M := linalg['swapcol'](M, d); M := linalg['swaprow'](M, d); G := linalg['swaprow'](G, d); if nargs = 13 then T := linalg['swapcol'](T, d); invT := linalg['swaprow'](invT, d) fi od else m := add(ntab[i+1], i=0..k-1); fi; # Rest of columns have valuations v_mat+k if p = x then for j from m+1 to n do for i from 1 to n do G[i, j] := subs(x = 0, normal(p^(-v_mat-k)*M[i, j])) od od else for j from m+1 to n do for i from 1 to n do a := normal(p^(-v_mat-k)*M[i, j]); gcdex(denom(a), p, x, 'U', 'V'); G[i, j] := rem(numer(a)*U, p, x) od od fi; for i from m+1 to n do G[i, i] := G[i, i] + lambda od; #AVOIR ntab := convert(ntab, 'array'); valtab := convert(valtab, 'array'); NULL; end: # g_column_rank(M, x, p, n, k, ntab, valtab, mk, G, P, T, invT) # INPUT: M - a rational function matrix # x - a name # p - an irreducible polynomial or infinity # n - a positive integer (dimension of M) # k - a positive integer (the reduction step) # ntab - a list of integers # vals - a list of integers # mk - an integer # P - a matrix # T - a matrix # invT - a matrix # OUTPUT: A matrix M' with mu_(k,p)(M') < mu_(k,p)(M) # g_column_rank := proc(M, x, p, n, k, ntab, valtab, mk, G, P, T, invT) local TT, invTT, i, j, mk1, c; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; userinfo(2, 'super_reduce', `reducing column rank`); P := linalg['transpose'](pcolechelon(linalg['transpose'](P), x, p)); mk1 := mk-ntab[k]; c := linalg['coldim'](P); for j from c to 1 by -1 do # eventual permutation on P, M, T and invT if P[mk-c+j, j] = 0 then i := mk-c+j; while P[i, j] = 0 do i := i-1 od; userinfo(5, 'super_reduce', `permute rows`,i,mk-c+j); P := linalg['swaprow'](P, i, mk-c+j); M := linalg['swapcol'](M, i, mk-c+j); M := linalg['swaprow'](M, i, mk-c+j); G := linalg['swaprow'](G, i, mk-c+j); if nargs = 12 then T := linalg['swapcol'](T, i, mk-c+j); invT := linalg['swaprow'](invT, i, mk-c+j) fi fi od; # generate constant transformation which reduces row rank for i from 1 to mk1 do P := linalg['mulrow'](P, i, p^(k+valtab[1]-valtab[i]-1)) od; TT := linalg['copyinto'](P, Id(n), 1, mk-c+1); invTT := linalg['inverse'](TT); if nargs = 12 then T := normalm(T &* TT); invT := normalm(invTT &* invT) fi; M := normalm(subs(x=x/(x+1),eval(invTT)) &* (M &* TT + 1/x^2*subs(x=x/(x+1),eval(TT))-1/x^2*TT)); end: # g_qtcd(M, x, p, n, k, s, ntab, valtab, mk, lambda, G, T, invT); # INPUT: M - a rational function matrix # x - a name # p - an irreducible polynomial or infinity # n - a positive integer (dimension of M) # k - a positive integer (the reduction step) # s - a negative integer or a name # ntab - an array of integers # valtab - an array of integers # mk - an integer # lambda - a name # G - a matrix (the G-matrix of M) # T - a matrix polynomial # invT - a matrix rational function # OUTPUT: A matrix M' whose G-matrix G' has reduced row rank # g_qtcd := proc(M, x, p, n, k, s, ntab, valtab, mk, lambda, G, T, invT) local GG, q, mk1, Gq, r, Y, i, P, invP, D, invD, U, V; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; userinfo(2, 'super_reduce', `enter qtcd-algorithm`); q := 0; mk1 := copy(mk); do GG := subs(lambda = 0, normalm(G)); userinfo(5, 'super_reduce', `q =`, q); Gq := linalg['submatrix'](GG, 1..n-q, 1..n-q); r := prank(linalg['submatrix'](Gq, 1..mk1, 1..n-q), x, p); if r < mk1 then userinfo(5, 'super_reduce', `rows are linearly dependent`); break fi; pgauss(linalg['transpose'](Gq), x, p, 'P'); Y := linalg['subvector'](P, 1..n-q, 1); if Y[n-q] = 0 then userinfo(5, 'super_reduce', `permute matrix and solution vector Y`); i := n-q-1; while Y[i] = 0 do i := i-1 od; Y[n-q] := Y[i]; Y[i] := 0; M := linalg['swaprow'](M, n-q, i); M := linalg['swapcol'](M, n-q, i); if nargs = 13 then invT := linalg['swaprow'](invT, n-q, i); T := linalg['swapcol'](T, n-q, i) fi; G := linalg['swaprow'](G, n-q, i); G := linalg['swapcol'](G, n-q, i); fi; # Normalizing Y[n-q] -> yields unimodular transformation matrix if degree(p, x) > 1 then gcdex(Y[n-q], p, x, 'U', 'V'); Y := map(rem, normalm(U * Y), p, x) else Y := normalm(1/Y[n-q]*Y); fi; P := Id(n); for i from 1 to n-q do P[n-q, i] := Y[i] od; invP := linalg['inverse'](P); M := normalm(subs(x=x/(x+1), P) &* (M &* invP + 1/x^2*subs(x=x/(x+1),eval(invP))-1/x^2*invP)); if nargs = 13 then T := normalm(T &* invP); invT := normalm(P &* invT) fi; if q = n-mk1 then break fi; userinfo(3, 'super_reduce', `recomputing G-matrix`); if nargs = 13 then g_sort_cols(M, x, p, n, k, lambda, ntab, valtab, s, G, false, T, invT) else g_sort_cols(M, x, p, n, k, lambda, ntab, valtab, s, G, false) fi; userinfo(1, 'super_reduce', `valuations:`, sort(convert(valtab, 'list'))); userinfo(4, 'super_reduce', `ntab:`, convert(ntab, 'list')); q := q+1 od; GG := subs(lambda = 0, normalm(G)); userinfo(5, 'super_reduce', `exit qtcd with q = `, q); userinfo(3, 'super_reduce', `applying diagonal transformation:`, mk, n-mk-q, q); D := linalg['diag'](seq(p, i=1..mk), seq(1, i=1..n-mk-q), seq(p, i=1..q)); invD := linalg['diag'](seq(1/p, i=1..mk), seq(1, i=1..n-mk-q), seq(1/p, i=1..q)); M := normalm(subs(x=x/(x+1),eval(invD)) &* (M &* D + 1/x^2*subs(x=x/(x+1),eval(D))-1/x^2*D)); if nargs = 13 then T := normalm(T &* D); invT := normalm(invD &* invT) fi; end: # pgauss(A, x, p, P) # INPUT: A - a m x n - matrix polynomial of degree < degree(p) # x - a name # p - an irreducible polynomial # P - (optional) a name # invP - (optional) a name # OUTPUT: the algebraic rank r of A as a matrix of K[x]/(p). I implemented # this since Maple cannot do this on matrices containing RootOfs. # Let N be the nullspace of A as matrix over K[x]/(p). Then P will # be assigned to a basis of N # pgauss := proc(A, x, p, P) local A1, T, i, j, k, l, pivot, n, m, r, t; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; m := linalg['rowdim'](A); n := linalg['coldim'](A); if ((p = x) or (p = infinity)) and not has(convert(A, 'set'), 'RootOf') then A1 := linalg['gausselim'](linalg['augment'](linalg['transpose'](A), linalg['diag'](seq(1, i=1..n))), 'r'); k := 1; i := 1; r := 0; while (k <= m) and (i <= n) do if A1[i, k] = 0 then k := k+1 else i := i+1; r := r+1 fi od; if (nargs = 4) and (r < n) then P := linalg['transpose'](linalg['submatrix'](A1, r+1..n, m+1..m+n)); fi; RETURN(r); fi; A1 := linalg['transpose'](A); t := n; T := linalg['diag'](seq(1, i=1..n)); i := 1; j := 1; r := 0; while (j <= m) and (i <= n) do pivot := 0; for k from i to n do if A1[k, j] <> 0 then pivot := A1[k, j]; A1 := linalg['swaprow'](A1, i, k); T := linalg['swaprow'](T, i, k); r := r+1; break; fi; od; if pivot = 0 then j := j+1; else if p = x then # p = x => there must be algebraic extensions for k from i+1 to n do if A1[k, j] <> 0 then for l from j+1 to m do A1[k, l] := evala(Normal(A1[k, l]-A1[i, l]*A1[k, j]/pivot)) od; for l from 1 to t do T[k, l] := evala(Normal(T[k, l] - T[i, l]*A1[k, j]/pivot)); od; A1[k, j] := 0; fi od else error "This part is only for the differential case"; # gcdex(pivot, p, x, 'U', 'V'); # for k from i+1 to n do # if A1[k, j] <> 0 then # for l from j+1 to m do # A1[k, l] := evala(Rem(A1[k, l]-U*A1[i, l]*A1[k, j], p, x)); # od; # for l from 1 to t do # T[k, l] := evala(Expand(T[k, l] - U*T[i, l]*A1[k, j])); # od; # A1[k, j] := 0; # fi # od fi; j := j+1; i := i+1: fi; od; if nargs = 4 then if r < n then P := linalg['transpose'](linalg['submatrix'](T, r+1..n, 1..n)); P := map(rem, normalm(P), p, x); fi; fi; r end: # prank(A, x, p) # INPUT: A - a m x n - matrix polynomial of degree < degree(p) # x - a name # p - an irreducible polynomial # OUTPUT: the algebraic rank r of A as a matrix of K[x]/(p). Same comment # as for pgauss. # prank := proc(A, x, p) local A1, U, V, i, j, k, l, pivot, n, m, r; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; n := linalg['rowdim'](A); m := linalg['coldim'](A); A1 := copy(A); # if (RootOfs = {}) and ((p = x) or (p = infinity)) then return linalg['rank'](A) fi; i := 1; j := 1; r := 0; while (j <= m) and (i <= n) do pivot := 0; for k from i to n do if A1[k, j] <> 0 then pivot := A1[k, j]; A1 := linalg['swaprow'](A1, i, k); r := r+1; break; fi; od; if pivot = 0 then j := j+1; else gcdex(pivot, p, x, 'U', 'V'); A1[i, j] := 1; for k from j+1 to m do A1[i, k] := rem(A1[i, k]*U, p, x) od; for k from i+1 to n do if A1[k, j] <> 0 then for l from j+1 to m do A1[k, l] := evala( Rem(A1[k, l] - A1[i, l]*A1[k, j], p, x)) od; A1[k, j] := 0; fi od; j := j+1; i := i+1 fi; od; r end: # pcolechelon(A, x, p) # INPUT: A - a m x n - matrix polynomial of degree < degree(p) # x - a name # p - an irreducible polynomial # OUTPUT: A as a matrix of K[x]/(p) transformed to column echelon form # pcolechelon := proc(A, x, p) local i, j, k, l, p1, pivot, n, m, A1, U, V; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; n := linalg['rowdim'](A); m := linalg['coldim'](A); A1 := copy(A); if (p = infinity) or (p = 1/x) then p1 := x else p1 := p fi; i := n; j := m; while (j >= 1) and (i >= 1) do pivot := 0; for k from i to 1 by -1 do if A1[k, j] <> 0 then pivot := A1[k, j]; A1 := linalg['swaprow'](A1, i, k); break; fi; od; if pivot = 0 then j := j-1; else gcdex(pivot, p1, x, 'U', 'V'); A1[i, j] := 1; for k from j-1 to 1 by -1 do A1[i, k] := rem(A1[i, k]*U, p1, x) od; for k from i-1 to 1 by -1 do if A1[k, j] <> 0 then for l from j-1 to 1 by -1 do A1[k, l] := evala( Rem(A1[k, l] - A1[i, l]*A1[k, j], p1, x)); od; A1[k, j] := 0; fi od; j := j-1; i := i-1 fi; od; normalm(A1) end: # pval(f, x, p [, lc]) # INPUT: f - a a rational function # x - a name # p - a polynomial or 1/x # lc - (optional) a name # OUTPUT: the order of f at p. If lc is passed as argument, it is # assigned to the leading coefficient of f at p. # pval := proc(f, x, p) option remember; poly_pval(expand(numer(f)), x, p)- poly_pval(expand(denom(f)), x, p) end: # poly_pval(a, x, p [,lc ]) # INPUT: a - a polynomial # x - a name # p - a polynomial or 1/x # lc - (optional) a name # OUTPUT: the order of a at p # poly_pval := proc(a, x, p, lc) local a1, ord, temp,q; option remember; a1 := collect(a, x); if not type(a, 'polynom') then error "polynomial expected" fi; if a = 0 then if nargs = 4 then lc := 0 fi; return infinity fi; if (p = 1/x) or (p = infinity) then if nargs = 4 then lc := lcoeff(a1, x) fi; return -degree(a, x) fi; if p = x then if nargs = 4 then lc := tcoeff(a1, x) fi; return ldegree(a, x) fi; ord := 0; temp := a; if indets(a, 'RootOf') union indets(p, 'RootOf') = {} then while divide(temp, p, 'q') do temp := diff(temp, x); ord := ord+1 od else while evala(Divide(temp, p, 'q')) do temp := q; ord := ord+1 od fi; if nargs = 4 then if ord = 0 then lc := evala(Rem(a, p, x)) else lc := evala(Rem(q, p, x)) fi fi; ord end: # normalm(M) # INPUT: M - a matrix # OUTPUT: same as evalm, but entries of M are normalized # normalm := proc(M) map(normal, evalm(M)) end: # Id(n) # INPUT: n - an integer # OUTPUT: The n x n identity matrix # Id := proc(n) local i; linalg['diag'](seq(1, i = 1..n)) end: