# 250715 RK(8,11) without rk5.txt with(LinearAlgebra): # RK_U123_(n) with 1-2-3simplifying assumption # B_s(r1,r2); coordinate-wise multiplication # delta_, weights global variables # calcPsi_(A, p); calculate vectors psi for coefficions of matrix A: psi[1], ... # OrderConditions(A, p); Equations, order contions # RK_8_11(c5,c6,c7,c9,c10) main function # --------------------------------------------------------------------- # --------------------------------------------------------------------- # RK_U123_(n) with 1st & 2nd & 3rd (dim M3'=2) simplifying assumption, b2=0, a[i,3]=..., a[i,2]=..., a[n,i]=... RK_U123_ := proc(n::integer)::Matrix; local A, i,j, eq; global a,b,c; # A := RK_setABC(n); A := Matrix(n+1,n+1); for i from 2 to n do for j from 1 to i-1 do A[i,j]:=a[i,j]; end do; end do; for j from 1 to n do A[n+1,j]:=b[j]; end do; c[1]:=0; c[n+1]:=1; for i from 2 to n+1 do A[i,1]:=c[i] - add(A[i,j], j=2..i-1); a[i,1]:=A[i,1]; end do; b[1]:=A[n+1,1]; a[n,1]:=b[1]; # dim M3'=2 ---------- for i from 4 to n+1 do A[i,3] := (c[i]^2*(c[i]-c[3]) - add(A[i,j]*c[j]*(3*c[j]-2*c[3]), j=4..i-1) )/c[3]^2; a[i,3] := A[i,3]; end do; b[3] := A[n+1,3]; # UPR(2) -------------- for i from 3 to n do A[i,2] := (c[i]^2/2 - add(A[i,j]*c[j], j=3..i-1) )/c[2]; a[i,2] := A[i,2]; end do; # special calculation for the last row (b[4]:=...) # b[2]=0 eq := add(A[i,j]*c[j], j=3..n) - 1/2; b[4] := solve(eq,b[4]); A[n+1,4] := b[4]; # UPR(1) -------------- for j from 4 to n-1 do A[n,j] := (b[j]*(1-c[j]) - add(b[i]*A[i,j], i=j+1..n-1))/b[n]; a[n,j] := A[n,j]; end do; # a[n,1] -------------- for i from 2 to n+1 do A[i,1] :=c[i] - add(A[i,j], j=2..i-1); a[i,1] := A[i,1]; end do: b[1]:=A[n+1,1]; a[n,1]:=A[n+1,1]; simplify(A); end proc: # --------------------------------------------------------------------- delta_ := Vector(200, [1,1,2,1,6,2,1,3,24, 8, 4,12, 6, 3, 4, 2,1, # 1..17 120,40,20,60,30,15,20,10, 5,24, 8, 4,12,12, 6, 6, 3, 4, 2, 1, # 18..37 720,240, 120, 360, 180, 90, 120, 60, 30, 144, 48, 24, 72, 72, 36, 36, # 38..53 18, 24, 12, 6, 120, 40, 20, 60, 30, 15, 20, 10, 5, 48, 16, 8, 24, 36, # 54..71 18, 9, 24, 8, 4, 12, 12, 6, 8, 6, 3, 4, 2, 1]); # 72..85 # delta[86..200] = 0! fill after! weights_:= Vector(200, [0,1,2,2,3,3,3,3, 4, 4, 4, 4, 4, 4, 4, 4, 4, # 1..17 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, # 18..37 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, # 38..53 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, # 54..71 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, # 72..85 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, # 86..200 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7 ]); # ---- coordinate-wise multiplication ------------------------------ B_s := proc(r1::Vector, r2::Vector)::Vector; simplify(Vector[column](Array(r1)*Array(r2))); end proc; # --------------------------------------------------------------------- # calculate vectors psi for coefficions of matrix A: psi[1], ... # for each i, psi[i] - vector! calcPsi_ := proc( A::Matrix, p::integer)::Vector; local n1, neq, e, psi, i; global weights, delta_, HasA2e, OneLeg; n1:=RowDimension(A); #print( 777, n1); e:=Vector(n1, [1$n1]); if p<1 then neq := 1; end if; if p=2 then neq := 2; end if; if p=3 then neq := 4; end if; if p=4 then neq := 8; end if; if p=5 then neq := 17; end if; if p=6 then neq := 37; end if; if p=7 then neq := 85; end if; if p=8 then neq :=200; end if; if p>8 then ERROR(`in calcPsi: p>8`); fi; psi:=Vector(neq); psi[ 1]:= e: # e if p<=1 then return psi; end if; psi[ 2]:= simplify(A . e): # Ae if p<=2 then return psi; end if; psi[ 3]:= simplify(A . psi[2]): # A2e psi[ 4]:= B_s( psi[2], psi[2] ): # Ae*Ae if p<=3 then return psi; end if; # * # | # * * * * # | | \ / # * * * * # t1 t2 t3 t4 psi[ 5]:= simplify(A . psi[3]): # A3e psi[ 6]:= B_s( psi[3], psi[2] ): # A2e*Ae psi[ 7]:= B_s( psi[4], psi[2] ): # Ae*Ae*Ae psi[ 8]:= simplify(A . psi[4]): # A(Ae*Ae) if p<=4 then return psi; end if; # * # | # * * * * # | | \ / # * * * * * * * # | \ / \|/ | # * * * * # t5 t6 t7 t8 psi[ 9]:= simplify(A . psi[5]): # A4e psi[10]:= simplify(A . psi[6]): # A(A2e*Ae) psi[11]:= simplify(A . psi[7]): # A(Ae*Ae*Ae) psi[12]:= simplify(A . psi[8]): # A2(Ae*Ae) psi[13]:= B_s( psi[5], psi[2] ): # A3e*Ae psi[14]:= B_s( psi[8], psi[2] ): # A(Ae*Ae)*Ae psi[15]:= B_s( psi[3], psi[3] ): # A2e*A2e psi[16]:= B_s( psi[3], psi[4] ): # A2e*Ae*Ae psi[17]:= B_s( psi[4], psi[4] ): # Ae*Ae*Ae*Ae # * # | # * * * * * # | | \ / | # * * * * * * * * * * * * * # | \ / \|/ | | \ / | | | # * * * * * * * * * * * * * * ** * # | | | | \ / |/ \ / \|/ \||/ # * * * * * * * * * # t9 t10 t11 t12 t13 t14 t15 t16 t17 # вес = 5 -------------------------------------------- if p<=5 then return psi; end if; # one leg psi[18]:= simplify(A . psi[ 9]): # A5e psi[19]:= simplify(A . psi[10]): # A2(A2e*Ae) psi[20]:= simplify(A . psi[11]): # A2(Ae*Ae*Ae) psi[21]:= simplify(A . psi[12]): # A3(Ae*Ae) psi[22]:= simplify(A . psi[13]): # A(A3e*Ae) psi[23]:= simplify(A . psi[14]): # A(A(Ae*Ae)*Ae) psi[24]:= simplify(A . psi[15]): # A(A2e*A2e) psi[25]:= simplify(A . psi[16]): # A(A2e*Ae*Ae) psi[26]:= simplify(A . psi[17]): # A(Ae*Ae*Ae*Ae) # * # | # * * * * * # | | \ / | # * * * * * * * * * * * * * # | \ / \|/ | | \ / | | | # * * * * * * * * * * * * * * ** * # | | | | \ / |/ \ / \|/ \||/ # * * * * * * * * * # | | | | | | | | | # * * * * * * * * * # t18 t19 t20 t21 t22 t23 t24 t25 t26 # two legs psi[27]:= B_s( psi[ 9], psi[2 ] ): # A4e*Ae psi[28]:= B_s( psi[10], psi[2 ] ): # A(A2e*Ae)*Ae psi[29]:= B_s( psi[11], psi[2 ] ): # A(Ae*Ae*Ae)*Ae psi[30]:= B_s( psi[12], psi[2 ] ): # A2(Ae*Ae)*Ae psi[31]:= B_s( psi[ 5], psi[3 ] ): # A3e*A2e psi[32]:= B_s( psi[ 8], psi[3 ] ): # A(Ae*Ae)*A2e # three legs psi[33]:= B_s( psi[13], psi[2 ] ): # A3e*Ae*Ae psi[34]:= B_s( psi[14], psi[2 ] ): # A(Ae*Ae)*Ae*Ae psi[35]:= B_s( psi[15], psi[2 ] ): # A2e*A2e*Ae # * # | # * * * * * # | \ / | | # * * * * * * * * * # \|/ \|/ \|/ # * * * # t33 t34 t35 # four legs psi[36]:= B_s( psi[16], psi[2 ] ): # A2e*Ae*Ae*Ae # * # | # * ** * # \||/ # * # t36 # five legs psi[37]:= B_s( psi[17], psi[2 ] ): # Ae*Ae*Ae*Ae*Ae # * *** * # \|||/ # * # t37 if p<=6 then return psi; end if; # вес = 6 -------------------------------------------- # Дерево Ae*Ae*Ae*Ae*Ae*Ae номер 181! # one leg psi[38]:= simplify(A . psi[18]): # A6e psi[39]:= simplify(A . psi[19]): # A3(A2e*Ae) psi[40]:= simplify(A . psi[20]): # A3(Ae*Ae*Ae) psi[41]:= simplify(A . psi[21]): # A4(Ae*Ae) psi[42]:= simplify(A . psi[22]): # A2(A3e*Ae) psi[43]:= simplify(A . psi[23]): # A2(A(Ae*Ae)*Ae) psi[44]:= simplify(A . psi[24]): # A2(A2e*A2e) psi[45]:= simplify(A . psi[25]): # A2(A2e*Ae*Ae) psi[46]:= simplify(A . psi[26]): # A2(Ae*Ae*Ae*Ae) psi[47]:= simplify(A . psi[27]): # A(A4e*Ae) psi[48]:= simplify(A . psi[28]): # A((A2e*Ae)*Ae) psi[49]:= simplify(A . psi[29]): # A((Ae*Ae*Ae)*Ae) psi[50]:= simplify(A . psi[30]): # A(A2(Ae*Ae)*Ae) psi[51]:= simplify(A . psi[31]): # A(A3e*A2e) psi[52]:= simplify(A . psi[32]): # A(A(Ae*Ae)*A2e) psi[53]:= simplify(A . psi[33]): # A(A3e*Ae*Ae) psi[54]:= simplify(A . psi[34]): # A(A(Ae*Ae)*Ae*Ae) psi[55]:= simplify(A . psi[35]): # A(A2e*A2e*Ae) psi[56]:= simplify(A . psi[36]): # A(A2*Ae*Ae*Ae) psi[57]:= simplify(A . psi[37]): # A(Ae*Ae*Ae*Ae*Ae) # two legs psi[58]:= B_s( psi[18], psi[2 ]): # A5e*Ae psi[59]:= B_s( psi[19], psi[2 ]): # A2(A2e*Ae)*Ae psi[60]:= B_s( psi[20], psi[2 ]): # A2(Ae*Ae*Ae)*Ae psi[61]:= B_s( psi[21], psi[2 ]): # A3(Ae*Ae)*Ae psi[62]:= B_s( psi[22], psi[2 ]): # A(A3e*Ae)*Ae psi[63]:= B_s( psi[23], psi[2 ]): # A(A(Ae*Ae)*Ae)*Ae psi[64]:= B_s( psi[24], psi[2 ]): # A(A2e*A2e)*Ae psi[65]:= B_s( psi[25], psi[2 ]): # A(A2e*Ae*Ae)*Ae psi[66]:= B_s( psi[26], psi[2 ]): # A(Ae*Ae*Ae*Ae)*Ae psi[67]:= B_s( psi[ 9], psi[3 ]): # A4e*A2e psi[68]:= B_s( psi[10], psi[3 ]): # A(A2e*Ae)*A2e psi[69]:= B_s( psi[11], psi[3 ]): # A(Ae*Ae*Ae)*A2e psi[70]:= B_s( psi[12], psi[3 ]): # A2(Ae*Ae)*A2e psi[71]:= B_s( psi[5], psi[5 ]): # A3e*A3e psi[72]:= B_s( psi[8], psi[5 ]): # A(Ae*Ae)*A3e psi[73]:= B_s( psi[8], psi[8 ]): # A(Ae*Ae)*A(Ae*Ae) # three legs psi[74]:= B_s( psi[27], psi[2 ]): # A4e*Ae*Ae psi[75]:= B_s( psi[28], psi[2 ]): # A(A2e*Ae)*Ae*Ae psi[76]:= B_s( psi[29], psi[2 ]): # A(Ae*Ae*Ae)*Ae*Ae psi[77]:= B_s( psi[30], psi[2 ]): # A2(Ae*Ae)*Ae*Ae psi[78]:= B_s( psi[31], psi[2 ]): # A3e*A2e*Ae psi[79]:= B_s( psi[32], psi[2 ]): # A(Ae*Ae)*A2e*Ae psi[80]:= B_s( psi[15], psi[3 ]): # A2e*A2e*A2e # four legs psi[81]:= B_s( psi[33], psi[2 ]): # A3e*Ae*Ae*Ae psi[82]:= B_s( psi[34], psi[2 ]): # A(Ae*Ae)*Ae*Ae*Ae psi[83]:= B_s( psi[35], psi[2 ]): # A2e*A2e*Ae*Ae # five legs psi[84]:= B_s( psi[36], psi[2 ] ): # A2e*Ae*Ae*Ae*Ae # six legs psi[85]:= B_s( psi[37], psi[2 ] ): # Ae*Ae*Ae*Ae*Ae*Ae if p<=7 then return psi; end if; # now L7: # L7 = A(L6) + (48 pcs) # Ae*L6 + (48 pcs) # A2e*L5 (where no Ae*.., indexes 18..26, 31, 32 total 11 pcs) # +A3e*L4 (where no Ae*.., A2e* indexes 9..12 total 4 pcs) # +A(Ae*Ae)*L4 (where no Ae*.., A2e* indexes 9..12 total 4 pcs) # total 48+48+11+4+4=115 pcs for i from 0 to 47 do: # A(L6) 86..133 psi[86+i] := simplify(A . psi[38+i]): delta_[86+i] := delta_[38+i]*7: end: for i from 0 to 47 do: # Ae*L6 134..181 psi[134+i] := B_s(psi[2], psi[38+i]): delta_[134+i] := delta_[38+i]: end: for i from 0 to 8 do: # A2e*L5 182..190 psi[182+i] := B_s(psi[3], psi[18+i]): delta_[182+i] := 2*delta_[18+i]: end: psi[191] := B_s(psi[3], psi[31]): delta_[191] := 2*delta_[31]: psi[192] := B_s(psi[3], psi[32]): delta_[192] := 2*delta_[32]: for i from 0 to 3 do: # A3e*L4 193..196 psi[193+i] := B_s(psi[5], psi[9+i]): delta_[193+i] := delta_[5]*delta_[9+i]: end: for i from 0 to 3 do: # A(Ae**Ae)*L4 197..200 psi[197+i] := B_s(psi[8], psi[9+i]): delta_[197+i] := delta_[8]*delta_[9+i]: end: return psi; end proc: # --------------------------------------------------------------------- OrderConditions := proc(A::Matrix, p::integer)::Vector; local n1, psi, neq, eqs, i,j; global delta_, weigths_; n1:=RowDimension(A); psi := calcPsi_(A, p); neq := Dimension(psi); eqs := Vector(neq); for i from 1 to neq do eqs[i]:=add(A[n1,j]*psi[i][j], j=1..n1-1)-1/(1+weights_[i])/delta_[i]; end do; for i from 1 to neq do eqs[i]:=factor(eqs[i]); end do; eqs; end proc; # ----------------------------------------------------------------------------- # simpify A,b,c simp := proc(a,b,c) local i,j,n; #global a,b,c,A,eqs,n, neq; n:=RowDimension(A)-1; for i from 1 to n do b[i]:=factor(b[i]); end; for i from 2 to n do c[i]:=factor(c[i]); end; for i from 2 to n do for j from 1 to i-1 do a[i,j] := factor(a[i,j]); end; end; end: # ---------------------------------------------- RK_8_11 := proc(c5,c6,c7,c9,c10) global a,b,c, A, eqs, neq, n,p; local i,sol, u,car, f0,f1,den, co74,co75,co1,co2,co3,simm4, j, a_75, a74,t74,de741,de742; a:='a'; b:='b'; c:='c'; c[5]:=c5: c[6]:=c6: c[7]:=c7: c[9]:=c9: c[10]:=c10: c[3]:=2/3: c[2]:=2/3*c[3]: b[2]:=0: # Fixed! c[11]:=1: c[12]:=1: n:=11: p:=8: A:=RK_U123_(n): eqs := OrderConditions(A,p): neq := Dimension(eqs); for i from 1 to 200 do if not i in [13,28,63,33,72,27,62,74,75,147,58, 59,134,154,135] then eqs[i]:=0; end: end: c[4]:= c[5]*(56*c[6]^2*c[5]^3-63*c[5]^2*c[6]^2-9*c[5]^2*c[6]+c[5]^2+27*c[6]^2*c[5]+c[6]*c[5]-3*c[6]^2)/ (168*c[6]^2*c[5]^4-168*c[6]*c[5]^4-357*c[6]^2*c[5]^3+393*c[5]^3*c[6]+3*c[5]^3+263*c[5]^2*c[6]^2- 303*c[5]^2*c[6]-6*c[5]^2-63*c[6]^2*c[5]+78*c[6]*c[5]+4*c[5]+6*c[6]^2-8*c[6]); c[4]:=factor(c[4]); c[8] := 1/2*(3*c[7]+14*c[5]*c[7]^2-28*c[5]*c[7]+3*c[5]+14*c[7]*c[5]^2+14*c[6]*c[7]^2- 28*c[6]*c[7]+3*c[6]-168*c[6]*c[5]*c[7]^2+231*c[6]*c[5]*c[7]-28*c[6]*c[5]+ 98*c[6]*c[5]^2*c[7]^2-168*c[6]*c[7]*c[5]^2+14*c[6]*c[5]^2+14*c[6]^2*c[7]+ 98*c[6]^2*c[5]*c[7]^2-168*c[6]^2*c[5]*c[7]+14*c[6]^2*c[5]+98*c[6]^2*c[7]*c[5]^2) / (6+21*c[5]^2+21*c[6]^2-84*c[5]*c[7]^2+77*c[5]*c[7]+98*c[5]^2*c[7]^2-84*c[6]*c[7]^2+ 77*c[6]*c[7]+77*c[6]*c[5]-84*c[6]*c[5]^2-84*c[6]^2*c[7]-84*c[6]^2*c[5]+ 98*c[6]^2*c[7]^2+98*c[6]^2*c[5]^2-84*c[7]*c[5]^2+364*c[6]*c[5]*c[7]^2- 294*c[6]*c[5]*c[7]-490*c[6]*c[5]^2*c[7]^2+364*c[6]*c[7]*c[5]^2- 490*c[6]^2*c[5]*c[7]^2+364*c[6]^2*c[5]*c[7]-490*c[6]^2*c[7]*c[5]^2+ 735*c[6]^2*c[5]^2*c[7]^2+21*c[7]^2-21*c[7]-21*c[5]-21*c[6]); c[8]:=factor(c[8]); #c8n := 3*c[7]+14*c[5]*c[7]^2-28*c[5]*c[7]+3*c[5]+14*c[7]*c[5]^2+14*c[6]*c[7]^2- # 28*c[6]*c[7]+3*c[6]-168*c[6]*c[5]*c[7]^2+231*c[6]*c[5]*c[7]-28*c[6]*c[5]+ # 98*c[6]*c[5]^2*c[7]^2-168*c[6]*c[7]*c[5]^2+14*c[6]*c[5]^2+14*c[6]^2*c[7]+ # 98*c[6]^2*c[5]*c[7]^2-168*c[6]^2*c[5]*c[7]+14*c[6]^2*c[5]+98*c[6]^2*c[7]*c[5]^2; #c8n := factor(c8n); #c8d := 6+21*c[5]^2+21*c[6]^2-84*c[5]*c[7]^2+77*c[5]*c[7]+98*c[5]^2*c[7]^2-84*c[6]*c[7]^2+ # 77*c[6]*c[7]+77*c[6]*c[5]-84*c[6]*c[5]^2-84*c[6]^2*c[7]-84*c[6]^2*c[5]+ # 98*c[6]^2*c[7]^2+98*c[6]^2*c[5]^2-84*c[7]*c[5]^2+364*c[6]*c[5]*c[7]^2- # 294*c[6]*c[5]*c[7]-490*c[6]*c[5]^2*c[7]^2+364*c[6]*c[7]*c[5]^2- # 490*c[6]^2*c[5]*c[7]^2+364*c[6]^2*c[5]*c[7]-490*c[6]^2*c[7]*c[5]^2+ # 735*c[6]^2*c[5]^2*c[7]^2+21*c[7]^2-21*c[7]-21*c[5]-21*c[6]; #c8d := factor(c8d); #print(777, factor(c[8]), c8n, c8d, factor(c8n/c8d/2)); a[5,4] := c[5]^2*(3*c[3]-2*c[5])/(6*c[4]*(-c[4]+c[3])); f0 := -1/375*(-120*c[3]+64)/(c[3]-c[5])/c[5]; f1 := -1/375*(375*c[3]*c[4]-375*c[4]^2)/(c[3]-c[5])/c[5]; a[6,5]:= 125/16*f0*c[6]^2*(1-c[6]) + f1*a[6,4]; a_75:=(27*c[6]*c[5]^5-30*c[6]*c[5]^4-8*c[6]*c[5]^3+84*c[6]^2*c[5]^6-217*c[6]^2*c[5]^5+126*c[6]^2*c[5]^4+ 48*c[6]^2*c[5]^3+2*c[6]^2*c[5]^2-84*c[6]^3*c[5]^5+231*c[6]^3*c[5]^4-207*c[6]^3*c[5]^3- 9*c[6]^3*c[5]^2+2*c[6]^3*c[5]-14*c[6]^4*c[5]^3+54*c[6]^4*c[5]^2-6*c[6]^4*c[5]-3*c[5]^5+ 4*c[5]^4+12*c[7]*c[6]^4+27*c[7]*c[5]^5-39*c[7]*c[5]^4+4*c[7]*c[5]^3-16*c[7]*c[6]^3+ 12*c[7]^2*c[6]^2-7*c[7]^2*c[5]^5+27*c[7]^2*c[5]^4-3*c[7]^2*c[5]^3-9*c[7]^2*c[6]^3- 327*c[7]*c[6]*c[5]^5+423*c[7]*c[6]*c[5]^4-12*c[7]*c[6]*c[5]^3-756*c[7]*c[6]^2*c[5]^6+ 2877*c[7]*c[6]^2*c[5]^5-2695*c[7]*c[6]^2*c[5]^4+360*c[7]*c[6]^2*c[5]^3- 84*c[7]*c[6]^2*c[5]^2-4704*c[7]*c[6]^3*c[5]^6+12348*c[7]*c[6]^3*c[5]^5-14175*c[7]*c[6]^3*c[5]^4+ 8874*c[7]*c[6]^3*c[5]^3-2232*c[7]*c[6]^3*c[5]^2+291*c[7]*c[6]^3*c[5]+4704*c[7]*c[6]^4*c[5]^6- 12348*c[7]*c[6]^4*c[5]^5+13230*c[7]*c[6]^4*c[5]^4-7266*c[7]*c[6]^4*c[5]^3+ 1714*c[7]*c[6]^4*c[5]^2-216*c[7]*c[6]^4*c[5]-2*c[7]*c[6]*c[5]^2+8*c[7]*c[6]^2*c[5]+ 63*c[7]^2*c[6]*c[5]^5-145*c[7]^2*c[6]*c[5]^4-153*c[7]^2*c[6]*c[5]^3+4900*c[7]^2*c[6]^2*c[5]^6- 13545*c[7]^2*c[6]^2*c[5]^5+14175*c[7]^2*c[6]^2*c[5]^4-6251*c[7]^2*c[6]^2*c[5]^3+ 1602*c[7]^2*c[6]^2*c[5]^2-4704*c[7]^2*c[6]^3*c[5]^6+12348*c[7]^2*c[6]^3*c[5]^5- 12033*c[7]^2*c[6]^3*c[5]^4+4410*c[7]^2*c[6]^3*c[5]^3-1163*c[7]^2*c[6]^3*c[5]^2+ 162*c[7]^2*c[6]^3*c[5]-196*c[7]^2*c[6]^4*c[5]^4+756*c[7]^2*c[6]^4*c[5]^3- 84*c[7]^2*c[6]^4*c[5]^2+60*c[7]^2*c[6]*c[5]^2- 216*c[7]^2*c[6]^2*c[5]-6*c[7]^2*c[6]*c[5]) / (-c[5]^2+9*c[6]*c[5]^2+c[6]*c[5]+28*c[6]^2*c[5]^3-63*c[6]^2*c[5]^2+ 9*c[6]^2*c[5]-c[6]^2+28*c[6]^3*c[5]^2)/24; a[7,5]:=a_75*(c[7]-c[5])*c[7]/(c[6]-c[5])/(c[5]-1)/(c[5]-1/2)/c[5]^2/(7*c[5]^2-7*c[5]+1); a74:= -168*c6*c5^7+501*c6*c5^6-513*c6*c5^5+174*c6*c5^4-4*c6*c5^3+504*c6^2*c5^7-1449*c6^2*c5^6+ 1433*c6^2*c5^5-441*c6^2*c5^4-18*c6^2*c5^3+4*c6^2*c5^2+42*c6^3*c5^5-190*c6^3*c5^4+ 126*c6^3*c5^3-12*c6^3*c5^2+81*c7*c5^6-162*c7*c5^5+72*c7*c5^4+32*c7^2*c6^2-24*c7^2*c6^3+ 168*c7^2*c5^7-777*c7^2*c5^6+1214*c7^2*c5^5-732*c7^2*c5^4+168*c7^2*c5^3-16*c7^2*c5^2+ 12*c7^3*c5-24*c7^3*c6+18*c7^3*c6^2+189*c7^3*c5^5-515*c7^3*c5^4+423*c7^3*c5^3-126*c7^3*c5^2- 216*c7*c6*c5^3+97083*c7*c6^2*c5^6+33189*c7*c6^2*c5^4-79857*c7*c6^2*c5^5-58716*c7*c6^2*c5^7+ 3363*c7*c6*c5^5-720*c7*c6*c5^4-3627*c7*c6*c5^6+1260*c7*c6*c5^7-6876*c7*c6^2*c5^3+ 864*c7*c6^2*c5^2+58338*c7*c6^3*c5^5-23076*c7*c6^3*c5^4+4602*c7*c6^3*c5^3-594*c7*c6^3*c5^2+ 24*c7*c6*c5^2-48*c7*c6^2*c5+14112*c7*c6^2*c5^8+36*c7*c6^3*c5-76146*c7*c6^3*c5^6+ 51156*c7*c6^3*c5^7-14112*c7*c6^3*c5^8+33279*c7^2*c6*c5^5-14990*c7^2*c6*c5^4+4044*c7^2*c6*c5^3- 32928*c7^2*c6^2*c5^7+56455*c7^2*c6^2*c5^6-65079*c7^2*c6^2*c5^5+51903*c7^2*c6^2*c5^4- 23682*c7^2*c6^2*c5^3+5244*c7^2*c6^2*c5^2+51646*c7^2*c6^3*c5^5-43134*c7^2*c6^3*c5^4+ 19170*c7^2*c6^3*c5^3-4076*c7^2*c6^3*c5^2-432*c7^2*c6*c5^2-624*c7^2*c6^2*c5+9408*c7^2*c6^2*c5^8+ 468*c7^2*c6^3*c5-34104*c7^2*c6^3*c5^6+9408*c7^2*c6^3*c5^7+16*c7^2*c6*c5-9408*c7^2*c6*c5^8- 34104*c7^3*c6*c5^6+49203*c7^3*c6*c5^5-35916*c7^3*c6*c5^4+14666*c7^3*c6*c5^3-14112*c7^3*c6^2*c5^7+ 45864*c7^3*c6^2*c5^6-58135*c7^3*c6^2*c5^5+35532*c7^3*c6^2*c5^4-12141*c7^3*c6^2*c5^3+ 2812*c7^3*c6^2*c5^2-588*c7^3*c6^3*c5^5+2660*c7^3*c6^3*c5^4-1764*c7^3*c6^3*c5^3+ 168*c7^3*c6^3*c5^2-3519*c7^3*c6*c5^2-351*c7^3*c6^2*c5+450*c7^3*c6*c5-45339*c7^2*c6*c5^6- 9*c5^6+18*c5^5-8*c5^4+9408*c7^3*c6*c5^7+32592*c7^2*c6*c5^7; t74:=56*c[6]^2*c[5]^3-63*c[5]^2*c[6]^2-9*c[5]^2*c[6]+c[5]^2+27*c[6]^2*c[5]+c[6]*c[5]-3*c[6]^2; de741:= (-c[5]^2+9*c[6]*c[5]^2+c[6]*c[5]+28*c[6]^2*c[5]^3-63*c[6]^2*c[5]^2+ 9*c[6]^2*c[5]-c[6]^2+28*c[6]^3*c[5]^2); de742:=3*c[5]^3-12*c[5]^2+8*c[5]-16*c[6]+813*c[6]*c[5]^3-336*c[6]*c[5]^4-609*c[6]*c[5]^2+156*c[6]*c[5]- 525*c[6]^2*c[5]^3+168*c[6]^2*c[5]^4+445*c[6]^2*c[5]^2-117*c[6]^2*c[5]+12*c[6]^2; a[7,4]:=a74/(c[5]-1)/(2*c[5]-1)/(7*c[5]^2-7*c[5]+1)/c[4]^2*t74*c[7]/de741/de742/12; # now b[:] sol:=solve({b[3],b[4]}, {b[5],b[6]}); assign(sol); b[2]; b[3]; b[4]; for j from 1 to 3 do u[j]:=factor( add(b[i]*A[i,j], i=j+1..n) - (1-c[j])*b[j] ); end; b[7]:=solve(u[1],b[7]); u[2]:=factor(u[2]); u[3]:=factor(u[3]); for i from 1 to n do b[i]:=factor(b[i]); end; # ---- carcas ------------------------------------------- car[1]:=add(b[i]*c[i] , i=2..n) - 1/2; car[2]:=add(b[i]*c[i]^2, i=2..n) - 1/3; car[3]:=add(b[i]*c[i]^3, i=2..n) - 1/4; car[4]:=add(b[i]*c[i]^4, i=2..n) - 1/5; car[5]:=add(b[i]*c[i]^5, i=2..n) - 1/6; car[6]:=add(b[i]*c[i]^6, i=2..n) - 1/7; car[7]:=add(b[i]*c[i]^7, i=2..n) - 1/8; sol:=solve({car[4],car[5],car[6],car[7]},{b[11],b[10],b[9],b[8]}); assign(sol); #for i from 1 to 11 do printf("b[%d]=%a\n", i, b[i]); end; den := (-c[6]+c[5])*(91*c[5]*c[6]-52*c[5]-52*c[6]+27)*(-c[6]+c[3])*c[6]; co74 := -c[4]*(91*c[5]*c[6]-52*c[5]-52*c[6]+27)*(-c[6]+c[5])*(-c[4]+c[3]); co75 := -c[5]*(91*c[5]*c[6]-52*c[5]-52*c[6]+27)*(-c[6]+c[5])*(-c[5]+c[3]); co1:= c[5]*c[6] * (5*c[6]-4) * (52*c[5]-27) * (-2+3*c[3])/6; co2:=-26*c[5]^2*c[3]+27/2*c[3]*c[5]+187/2*c[6]^2*c[3]-135/2*c[3]*c[6]+455/3*c[5]^2*c[6]^2- 364/3*c[5]^2*c[6]+36*c[6]-455/2*c[3]*c[5]^2*c[6]^2+455/2*c[5]^2*c[3]*c[6]-91/2*c[6]^2*c[3]*c[5]-45*c[6]^2; co3:=-130*c[6]^2*c[3]+104*c[3]*c[6]-91/3*c[5]^2*c[6]-364/3*c[5]*c[6]^2+364/3*c[5]*c[6]- 182*c[5]*c[3]*c[6]-9*c[5]-181/3*c[6]+455/2*c[6]^2*c[3]*c[5]+52/3*c[5]^2+208/3*c[6]^2; a[7,6] := (co1*c[7] + co2*c[7]^2 + co3*c[7]^3 + co74*a[7,4] + co75*a[7,5])/den; a[6,4]:=-(168*c[5]^5*c[6]^2-168*c[5]^5*c[6]-441*c[5]^4*c[6]^2+21*c[5]^3*c[6]^3+501*c[5]^4*c[6]+ 374*c[5]^3*c[6]^2-95*c[5]^2*c[6]^3-9*c[5]^4-522*c[5]^3*c[6]-12*c[5]^2*c[6]^2+ 63*c[5]*c[6]^3+18*c[5]^3+192*c[5]^2*c[6]-72*c[5]*c[6]^2-6*c[6]^3-8*c[5]^2- 12*c[5]*c[6]+8*c[6]^2)*(168*c[5]^4*c[6]^2-168*c[5]^4*c[6]-357*c[5]^3*c[6]^2+ 393*c[5]^3*c[6]+263*c[5]^2*c[6]^2+3*c[5]^3-303*c[5]^2*c[6]-63*c[5]*c[6]^2- 6*c[5]^2+78*c[5]*c[6]+6*c[6]^2+4*c[5]-8*c[6])^2 * c[6] / ((12*(-1+c[5]))*(2*c[5]-1)*(7*c[5]^2-7*c[5]+1) * (168*c[5]^4*c[6]^2-336*c[5]^4*c[6]-525*c[5]^3*c[6]^2+813*c[5]^3*c[6]+445*c[5]^2*c[6]^2+ 3*c[5]^3-609*c[5]^2*c[6]-117*c[5]*c[6]^2-12*c[5]^2+156*c[5]*c[6]+12*c[6]^2+8*c[5]-16*c[6])* (56*c[5]^3*c[6]^2-63*c[5]^2*c[6]^2-9*c[5]^2*c[6]+27*c[5]*c[6]^2+ c[5]^2+c[5]*c[6]-3*c[6]^2)*c[5]^2); simp(a,b,c): a[10,9]:=solve(eqs[ 13], a[10,9]): if length(a[10,9])=1 then return -1; end; simp(a,b,c): a[10,8]:=solve(eqs[ 28], a[10,8]): if length(a[10,8])=1 then return -1; end; simp(a,b,c): a[10,7]:=solve(eqs[ 63], a[10,7]): if length(a[10,7])=1 then return -1; end; simp(a,b,c): a[ 9,8]:=solve(eqs[ 33], a[ 9,8]): if length(a[ 9,8])=1 then return -1; end; simp(a,b,c): a[ 8,7]:=solve(eqs[ 72], a[ 8,7]): if length(a[ 8,7])=1 then return -1; end; simp(a,b,c): a[10,6]:=solve(eqs[ 27], a[10,6]): if length(a[10,6])=1 then return -1; end; simp(a,b,c): a[10,4]:=solve(eqs[ 62], a[10,4]): if length(a[10,4])=1 then return -1; end; simp(a,b,c): a[ 9,7]:=solve(eqs[ 74], a[ 9,7]): if length(a[ 9,7])=1 then return -1; end; simp(a,b,c): a[ 9,6]:=solve(eqs[ 75], a[ 9,6]): if length(a[ 9,6])=1 then return -1; end; simp(a,b,c): a[10,5]:=solve(eqs[147], a[10,5]): if length(a[10,5])=1 then return -1; end; simp(a,b,c): a[ 8,6]:=solve(eqs[ 58], a[ 8,6]): if length(a[ 8,6])=1 then return -1; end; simp(a,b,c): a[ 8,5]:=solve(eqs[ 59], a[ 8,5]): if length(a[ 8,5])=1 then return -1; end; simp(a,b,c): a[ 9,5]:=solve(eqs[134], a[ 9,5]): if length(a[ 9,5])=1 then return -1; end; simp(a,b,c): a[ 8,4]:=solve(eqs[154], a[ 8,4]): if length(a[ 8,4])=1 then return -1; end; simp(a,b,c): a[ 9,4]:=solve(eqs[135], a[ 9,4]): if length(a[ 9,4])=1 then return -1; end; simp(a,b,c): for i from 2 to n+1 do for j from 1 to i-1 do A[i,j]:=factor(A[i,j]); end: end: return A; end proc: