> restart; > with(LinearAlgebra): > with(Involutive): > pol:=proc(a,n) > t^n+add((-1)^i*a[i]*t^(n-i),i=1..n); > end proc: > pol(a,2); 2 t - a[1] t + a[2] > REL:=proc(a,n,b,m,c) > local A,B,C; > A:=CompanionMatrix(pol(a,n),t); > B:=CompanionMatrix(pol(b,m),t); > C:=Matrix(convert(Matrix(n,n,(i,j)->A[i,j]*B),listlist)); > map(i->coeff(CharacteristicPolynomial(C,t)-pol(c,n*m),t,n*m-i),[$1..n* > m]); > end proc: > REL1:=proc(a,n,b,m,c) > local A,B,C; > A:=CompanionMatrix(pol(a,n-1)*(t-1),t); > B:=CompanionMatrix(pol(b,m),t); > C:=Matrix(convert(Matrix(n,n,(i,j)->A[i,j]*B),listlist)); > map(i->coeff(CharacteristicPolynomial(C,t)-pol(c,n*m),t,n*m-i),[$1..n* > m]); > end proc: > # > L:=REL(a,2,b,3,c); L := [-a[1] b[1] + c[1], 2 2 -2 a[2] b[2] + b[1] a[2] + a[1] b[2] - c[2], 3 3 a[1] a[2] b[3] - b[1] a[2] a[1] b[2] - a[1] b[3] + c[3], 2 2 2 2 -2 b[3] a[2] b[1] + a[2] b[1] a[1] b[3] + a[2] b[2] 2 3 2 - c[4], -a[2] a[1] b[3] b[2] + c[5], a[2] b[3] - c[6]] > sL:=solve({op(L)},map(i->c[i],{$1..6})); 3 2 sL := {c[1] = a[1] b[1], c[6] = a[2] b[3] , 2 c[5] = a[2] a[1] b[3] b[2], 2 2 c[2] = -2 a[2] b[2] + b[1] a[2] + a[1] b[2], 3 c[3] = -3 a[1] a[2] b[3] + b[1] a[2] a[1] b[2] + a[1] b[3], c[4] = 2 2 2 2 -2 b[3] a[2] b[1] + a[2] b[1] a[1] b[3] + a[2] b[2] } > > L1:=subs([c[6]=1,b[3]=1,a[2]=1],L)[1..-2]; 2 2 L1 := [-a[1] b[1] + c[1], -2 b[2] + b[1] + a[1] b[2] - c[2], 3 3 a[1] - b[1] a[1] b[2] - a[1] + c[3], 2 2 -2 b[1] + b[1] a[1] + b[2] - c[4], -a[1] b[2] + c[5]] > v := [ c[5] = 10, c[4] = 8, c[3] = 6, c[2] = 4, c[1] = 2, > b[2]=2,b[1]=1,a[1]=1]; v := [c[5] = 10, c[4] = 8, c[3] = 6, c[2] = 4, c[1] = 2, b[2] = 2, b[1] = 1, a[1] = 1] > > Verf:=proc(v,gg) > local an1,an2,aa; > global J; > aa:=lhs(v[-1]); > J:=InvolutiveBasisGINV(J,[op(v[1..-2]),aa=gg]): > an1:=nops(map(r->if not has(r,aa) then r end if,J)); > an2:=nops(map(r->if not > has(LeadingMonomial(r,[op(v[1..-2]),aa=gg]),aa) then r end if,J)); > return(an1,an2,gg,nops(J)); > end proc: > J:=copy(L1); 2 2 J := [-a[1] b[1] + c[1], -2 b[2] + b[1] + a[1] b[2] - c[2], 3 3 a[1] - b[1] a[1] b[2] - a[1] + c[3], 2 2 -2 b[1] + b[1] a[1] + b[2] - c[4], -a[1] b[2] + c[5]] > v := [c[5] = 10, c[4] = 8, c[3] = 6, c[2] = 4, c[1] = 2, b[2] = 2, > b[1] = 1,a[1]=1]; v := [c[5] = 10, c[4] = 8, c[3] = 6, c[2] = 4, c[1] = 2, b[2] = 2, b[1] = 1, a[1] = 1] # sequence of degrees for a[1]: 1,5,10 > Verf(v,1);Verf(v,5);Verf(v,10); 0, 5, 1, 5 35, 37, 5, 95 48, 48, 10, 105 > J:=select(r->not has(r,a[1]),J): > nops(J); 48 > v := [ c[5] = 10, c[4] = 8, c[3] = 6, c[2] = 4, c[1] = 2, b[2] = 2, > b[1] = 1]; v := [c[5] = 10, c[4] = 8, c[3] = 6, c[2] = 4, c[1] = 2, b[2] = 2, b[1] = 1] # sequence of degrees for b[1]: 1,4,8,16,24 > Verf(v,1);Verf(v,4);Verf(v,8); 0, 17, 1, 31 0, 13, 4, 25 0, 16, 8, 37 > Verf(v,16);Verf(v,20); 16, 23, 16, 77 25, 27, 20, 83 > Verf(v,24); 28, 28, 24, 83 > J:=select(r->not has(r,b[1]),J): > nops(J);indets(J); 28 {b[2], c[1], c[2], c[3], c[4], c[5]} > v := [ c[5] = 10, c[4] = 8, c[3] = 6, c[2] = 4, c[1] = 2, b[2] = 2]; v := [c[5] = 10, c[4] = 8, c[3] = 6, c[2] = 4, c[1] = 2, b[2] = 2] # sequence of degrees for b[2]: 4,10,14,20,25,30,35,40 > Verf(v,4);Verf(v,10);Verf(v,14); 0, 15, 4, 25 0, 12, 10, 91 0, 16, 14, 84 > Verf(v,20);Verf(v,30); Warning, resulting involutive basis is big; reading it may take a while... 6, 19, 20, 148 Warning, resulting involutive basis is big; reading it may take a while... 21, 21, 30, 164 > > J:=select(r->not has(r,b[2]),J): > v:=map(i->c[7-i]=7-i,[$2..6]); v := [c[5] = 5, c[4] = 4, c[3] = 3, c[2] = 2, c[1] = 1] Warning, resulting involutive basis is big; reading it may take a while... > DEG:=proc(p)degree(subs(map(i->c[i]=t^i,[$1..8]),LeadingMonomial(p,v)) > ); > end proc: > > map(DEG,J); [19, 20, 21, 21, 22, 22, 23, 23, 23, 24, 24, 25, 25, 25, 26, 26, 27, 27, 28, 30, 31] > ErzDet1:=copy(J): > save v, > ErzDet1,"/home/plesken/CharlesMatrix.d/paper.d/Sechs.d/ErzDet1.m"; > read("/home/plesken/CharlesMatrix.d/paper.d/Sechs.d/ErzDet1.m"): > J:=InvolutiveBasisGINV(J,v): > AssertInvBasis(J,v): > PolHilbertSeries(t); 2 3 4 5 6 7 8 1 + 5 t + 15 t + 35 t + 68 t + 115 t + 173 t + 244 t + 327 t 9 10 / 420 93 11 \ + 422 t + t |----- + -------- + --------| |1 - t 2 3| \ (1 - t) (1 - t) / > DEG:=proc(p) > degree(subs(map(i->c[i]=t^i,[$1..8]),LeadingMonomial(p,v))); > end proc: > map(DEG,J); [19, 20, 21, 21, 22, 22, 23, 23, 23, 24, 24, 25, 25, 25, 26, 26, 27, 27, 28, 30, 31] > nops(select(irreduc,J)); 21 > nops(J); 21 > J1:=subs(map(i->c[i]=c[i]/t^i,[$1..7]),J): > J1:=map(numer,J1): > J1:=algsubs(t^6=c[6],J1): > v:=map(i->c[7-i]=7-i,[$1..6]); v := [c[6] = 6, c[5] = 5, c[4] = 4, c[3] = 3, c[2] = 2, c[1] = 1] > map(DEG, J1); [19, 20, 21, 21, 22, 22, 23, 23, 23, 24, 24, 25, 25, 25, 26, 26, 27, 27, 28, 30, 31] > indets(J1); {c[1], c[2], c[3], c[4], c[5], c[6]} > Jf:=InvolutiveBasisGINV(J1,v): Warning, resulting involutive basis is big; reading it may take a while... > map(DEG, Jf); [19, 20, 21, 21, 22, 22, 23, 23, 23, 24, 24, 25, 25, 25, 26, 26, 27, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31, 32, 32, 33, 33, 34, 34, 35, 35, 36, 36, 37, 37, 38, 39, 41] > > > # Minimal generating set by inspection > > AssertInvBasis(Jf,v): > Hf:=evala(subs(map(i->c[i]=t^i,[$1..6]),SubmoduleBasis(v))); 19 11 9 10 15 7 8 13 Hf := t (1 + t + t - 4 t - 4 t + 3 t - 3 t - 4 t + 3 t 14 18 2 3 4 19 21 22 / + 2 t - t + 2 t + 2 t + 3 t - t - t + t ) / ( / 6 5 4 3 2 (-1 + t ) (-1 + t ) (-1 + t ) (-1 + t ) (-1 + t ) (t - 1)) > series(Hf,t=0,45); 19 20 21 22 23 24 25 26 t + 2 t + 5 t + 9 t + 17 t + 25 t + 40 t + 55 t + 27 28 29 30 31 32 78 t + 103 t + 137 t + 175 t + 224 t + 277 t + 33 34 35 36 37 38 344 t + 418 t + 507 t + 605 t + 721 t + 848 t + 39 40 41 42 43 996 t + 1158 t + 1343 t + 1546 t + 1775 t + 2024 44 45 t + O(t ) > interface(rtablesize=27); 24 > Transpose(<<[$1..25]>|DEG(Jf[i]),[$1..25])>>); [1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 , 11 , 12 , 13 , 14 , 15 , 16 , 17 , 18 , 19 , 20 , 21 , 22 , 23 , 24 , 25] [19 , 20 , 21 , 21 , 22 , 22 , 23 , 23 , 23 , 24 , 24 , 25 , 25 , 25 , 26 , 26 , 27 , 27 , 27 , 28 , 28 , 29 , 29 , 30 , 30] > J20:=AssertInvBasis(Jf,v, "extract"=[1,2]): > PolInvReduce(J1[3]); > J24:=AssertInvBasis(Jf,v, "extract"=1..11): > LeadingMonomial(PolInvReduce(J1[12]),v); 2 3 c[4] c[3] c[2] c[6] > LeadingMonomial(PolInvReduce(J1[13]),v); 2 3 c[4] c[3] c[2] c[6] > LeadingMonomial(PolInvReduce(J1[14]),v); 2 2 c[5] c[4] c[3] c[6] > coeff(coeff(PolInvReduce(J1[12]),c[4],2),c[3],3); 3 2 4 6 c[5] c[1] + 61/8 c[2] c[1] - 17/8 c[2] c[1] - 2 c[5] c[2] c[1] 3 2 - 3/4 c[2] c[1] + c[2] c[6] > coeff(coeff(PolInvReduce(J1[13]),c[4],2),c[3],3); 3 2 4 6 -10 c[5] c[1] - 305/4 c[2] c[1] + 85/4 c[2] c[1] 3 2 + 20 c[5] c[2] c[1] + 15/2 c[2] c[1] - 10 c[2] c[6] > LeadingMonomial(PolInvReduce(J1[13]+10*J1[12]),v); 1 # Hence, 2 generators 25 instead of 3. > DEG(Jf[15]); 26 > J25:=AssertInvBasis(Jf,v, "extract"=1..14): > LeadingMonomial(PolInvReduce(J1[15]),v); 2 2 c[5] c[4] c[2] c[6] > LeadingMonomial(PolInvReduce(J1[16]),v); 1 # Hence, one generator of degree 26 instead of 2. # # > J26:=AssertInvBasis(Jf,v, "extract"=1..16): > LeadingMonomial(PolInvReduce(J1[17]),v); 2 2 c[5] c[4] c[3] c[6] > LeadingMonomial(PolInvReduce(J1[18]),v); 2 2 c[5] c[4] c[3] c[6] > coeff(coeff(PolInvReduce(J1[17]),c[6],2),c[5],1); 2 2 2 2 2 2 4 -4 c[4] c[3] + 4 c[3] c[2] + 2 c[4] c[1] - 5/8 c[3] c[1] 3 4 6 2 6 + 23/4 c[2] c[1] - 11/4 c[4] c[1] - 41/2 c[2] c[1] 49 7 8 + -- c[3] c[1] + 7/4 c[2] c[1] + 2 c[4] c[3] c[2] c[1] 16 2 2 2 2 - 3/2 c[3] c[2] c[1] - 2 c[4] c[2] c[1] 3 2 3 - 11/4 c[4] c[3] c[1] + 13/4 c[3] c[2] c[1] 4 5 10 + 2 c[4] c[2] c[1] - 31/4 c[3] c[2] c[1] - 161/8 c[1] > coeff(coeff(PolInvReduce(J1[18]),c[6],2),c[5],1); 2 2 2 2 2 4 2 68 c[4] c[3] - 68 c[3] c[2] - 37 c[4] c[1] - 3 c[2] c[1] 2 4 3 4 6 - 27/8 c[3] c[1] - 717/8 c[2] c[1] + 265/4 c[4] c[1] 4377 2 6 3351 7 567 8 + ---- c[2] c[1] - ---- c[3] c[1] + --- c[2] c[1] 16 32 16 3 - 43 c[4] c[3] c[2] c[1] + 9 c[3] c[2] c[1] 2 2 2 2 + 21 c[3] c[2] c[1] + 145/4 c[4] c[2] c[1] 3 2 3 + 199/4 c[4] c[3] c[1] - 487/8 c[3] c[2] c[1] 4 2625 5 6647 10 - 93/2 c[4] c[2] c[1] + ---- c[3] c[2] c[1] + ---- c[1] 16 16 > LeadingMonomial(PolInvReduce(17*J1[17]+J1[18]),v); 3 3 c[4] c[3] c[6] # Hence, 2 generators of degree 27, one less than in Janet basis. > J27:=AssertInvBasis(Jf,v, "extract"=1..19): > LeadingMonomial(PolInvReduce(J1[18]),v); 1 # Hence, no generators in degree 28 > J29:=AssertInvBasis(Jf,v, "extract"=1..23): > LeadingMonomial(PolInvReduce(J1[-2]),v); 2 2 2 c[5] c[4] c[6] # Hence, one generator of degree 30, one less than in Janet basis > J30:=AssertInvBasis(Jf,v, "extract"=1..25): > LeadingMonomial(PolInvReduce(J1[-1]),v); 2 3 2 c[4] c[3] c[2] c[6] > map(DEG,Jf[1..26]); [19, 20, 21, 21, 22, 22, 23, 23, 23, 24, 24, 25, 25, 25, 26, 26, 27, 27, 27, 28, 28, 29, 29, 30, 30, 31] > J30n:=AssertInvBasis(Jf[1..25],v): > LeadingMonomial(PolInvReduce(J1[-1],J30n,v),v); 2 3 2 c[4] c[3] c[2] c[6] > > > > > > > # Minimal generating set > > ssystem("date"); [0, "Fri May 2 23:59:00 PDT 2008 "] > > > ithprime(1599); 13487 > InvolutiveOptions("char", %); 0 > > > AssertInvBasis(Jfin, v): > > Jfin_multvar := map(i->i[5], P_T_List): > save_P_T_List := P_T_List: > save_PolTab_Var := PolTab_Var: > > var := map(lhs, v); > nvar := nops(var); > n_Jfin := nops(Jfin); var := [c[6], c[5], c[4], c[3], c[2], c[1]] nvar := 6 n_Jfin := 42 > > > Jfin_degs := map(DEG, Jfin); Jfin_degs := [19, 20, 21, 21, 22, 22, 23, 23, 23, 24, 24, 25, 25, 25, 26, 26, 27, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31, 32, 32, 33, 33, 34, 34, 35, 35, 36, 36, 37, 37, 38, 39, 41] > > degs := [Jfin_degs[1]]: > first_deg := [1]: > for i from 2 to nops(Jfin_degs) do > if Jfin_degs[i] > degs[-1] then > degs := [op(degs), Jfin_degs[i]]; > first_deg := [op(first_deg), i]; > elif Jfin_degs[i] < degs[-1] then > ERROR(`not sorted!`); > fi; > od: > first_deg := [op(first_deg), nops(Jfin_degs)+1]: > first_deg; [1, 2, 3, 5, 7, 10, 12, 15, 17, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 41, 42, 43] > nops(degs); 22 > degs; [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 41] > > Jfin_lm := LeadingMonomial(Jfin, v): > > JJ := [[1]]: > for i from 2 to nops(degs) do > print("degree", degs[i]); > > # in the present context, the following 3 lines are equivalent to > # AssertInvBasis(Jfin, v, "extract"=1..(first_deg[i]-1)): > # but are much faster > P_T_List := save_P_T_List[1..(first_deg[i]-1)]: > PolTab_Var := save_PolTab_Var[1..(first_deg[i]-1)]: > pnumber_t := first_deg[i]-1: > > ind := [seq(select(j->Jfin_degs[j] = degs[i]-k, > [$1..(first_deg[i]-1)]), k=1..6)]; > LL := []: > for k from 1 to 6 do > for l from 1 to nops(ind[k]) do > if Jfin_multvar[ind[k][l]][6-k+1] = 0 then > LL := [op(LL), var[6-k+1]*Jfin[ind[k][l]]]; > fi: > od: > od: > > if LL = [] then > JJ := [op(JJ), [$first_deg[i]..(first_deg[i+1]-1)]]: > else > RR := []; > RL := []; > RC := []; > for j from 1 to nops(LL) do > gaussred := true; > print ("j", j); > p := LL[j]; > while p <> 0 and gaussred do > gaussred := false; > p := PolInvReduce(p): > > lm := LeadingMonomial(p, v, "C"): > lc := coeffs(lm, var, 'mon'); > lm := mon; > > k := 1; > while k <= nops(RL) and lm <> RL[k] do > k := k+1; > od; > if k <= nops(RL) then > print("Gauss reduction"); > p := normal(p - lc/RC[k]*RR[k]); > gaussred := true; > fi; > od; > if p <> 0 then > RR := [op(RR), p]: > RL := [op(RL), lm]: > RC := [op(RC), lc]: > fi; > od; > JJ := [op(JJ), select(j->not member(Jfin_lm[j], RL), > [$first_deg[i]..(first_deg[i+1]-1)])]: > fi; > od: "degree", 20 "degree", 21 "degree", 22 "degree", 23 "degree", 24 "j", 1 "j", 2 "degree", 25 "j", 1 "j", 2 "j", 3 "degree", 26 "j", 1 "j", 2 "j", 3 "j", 4 "j", 5 "degree", 27 "j", 1 "j", 2 "j", 3 "j", 4 "j", 5 "Gauss reduction" "j", 6 "j", 7 "Gauss reduction" "degree", 28 "j", 1 "j", 2 "j", 3 "j", 4 "j", 5 "j", 6 "j", 7 "degree", 29 "j", 1 "j", 2 "j", 3 "j", 4 "Gauss reduction" "j", 5 "Gauss reduction" "j", 6 "Gauss reduction" "j", 7 "degree", 30 "j", 1 "j", 2 "j", 3 "j", 4 "j", 5 "j", 6 "degree", 31 "j", 1 "j", 2 "Gauss reduction" "j", 3 "j", 4 "Gauss reduction" "j", 5 "j", 6 "degree", 32 "j", 1 "j", 2 "j", 3 "j", 4 "j", 5 "j", 6 "degree", 33 "j", 1 "j", 2 "j", 3 "j", 4 "Gauss reduction" "j", 5 "Gauss reduction" "j", 6 "j", 7 "Gauss reduction" "degree", 34 "j", 1 "j", 2 "j", 3 "j", 4 "j", 5 "degree", 35 "j", 1 "j", 2 "j", 3 "Gauss reduction" "j", 4 "j", 5 "Gauss reduction" "j", 6 "degree", 36 "j", 1 "j", 2 "j", 3 "j", 4 "j", 5 "degree", 37 "j", 1 "j", 2 "Gauss reduction" "j", 3 "j", 4 "Gauss reduction" "j", 5 "degree", 38 "j", 1 "j", 2 "j", 3 "j", 4 "degree", 39 "j", 1 "j", 2 "j", 3 "j", 4 "Gauss reduction" "j", 5 "Gauss reduction" "degree", 41 "j", 1 "j", 2 "j", 3 "j", 4 "Gauss reduction" > > > ssystem("date"); [0, "Sat May 3 01:09:29 PDT 2008 "] > > > > nops(JJ); 22 > map(nops, JJ); [1, 1, 2, 2, 3, 2, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] > JJ; [[1], [2], [3, 4], [5, 6], [7, 8, 9], [10, 11], [12], [15], [17], [20], [], [25], [], [], [], [], [], [], [], [], [], []] > > indi := map(op, JJ); indi := [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 15, 17, 20, 25] > map(i->Jfin_degs[i], indi); [19, 20, 21, 21, 22, 22, 23, 23, 23, 24, 24, 25, 26, 27, 28, 30] > nops(indi); 16 > > # verify in characteristic 0: > > InvolutiveOptions("char", 0); 13487 > Jminerz := map(i->Jfin[i], indi): > > Jneu := InvolutiveBasisGINV(Jminerz, v): > AssertInvBasis(Jneu, v): > Fneu := FactorModuleBasis(v): > nops(Jneu); 42 > simplify(F - Fneu); 0 > > > > > > #