Revision history [back]

It seems to work, thanks. Code like

m=matrix.identity(nn+1)
nn1_half=(nn+1)/2
for j in range(nn1_half) :
  m=m*elementary_matrix(nn+1,row1=j,row2=nn-j)
vecordinv_matrix=m # a (nn+1)-sized matrix with 1 on the anti-diagonal and 0 as other entries
#
A=matrix(CC,nn+1)
A[0,0:4]=vector((N(m11_tot(aa,nn)),N(m12_tot(aa,nn)),N(m13_tot(aa,nn)),N(m14_tot(aa,nn))))
A[1,0:5]=vector((N(m12_tot(aa,nn)),N(m22_tot(aa,nn)),N(m23_tot(aa,nn)),N(m24_tot(aa,nn)), \
                 N(m25_tot(aa,nn))))
A[2,0:6]=vector((N(m13_tot(aa,nn)),N(m23_tot(aa,nn)),N(m33_tot(aa,nn)),N(m34_tot(aa,nn)), \
                 N(m35_tot(aa,nn)),N(m36_tot(aa,nn))))
A[3,0:7]=vector((N(m14_tot(aa,nn)),N(m24_tot(aa,nn)),N(m34_tot(aa,nn)),N(m33_tot(aa,nn)), \
                 N(m34_tot(aa,nn)),N(m35_tot(aa,nn)),N(m36_tot(aa,nn))))
A[4,1:8]=vector((N(m25_tot(aa,nn)),N(m35_tot(aa,nn)),N(m34_tot(aa,nn)),N(m33_tot(aa,nn)), \
                 N(m34_tot(aa,nn)),N(m35_tot(aa,nn)),N(m36_tot(aa,nn))))
for j in [5..nn-5] :
  A[j,j-3:j+4]=vector((N(m36_tot(aa,nn)),N(m35_tot(aa,nn)),N(m34_tot(aa,nn)),N(m33_tot(aa,nn)), \
                       N(m34_tot(aa,nn)),N(m35_tot(aa,nn)),N(m36_tot(aa,nn))))
A_aux=matrix(CC,nn+1,5)
for j in range(5) :
  A_aux[:,j]=vecordinv_matrix*A[j,:].transpose()
  A[nn-j,:]=A_aux[:,j].transpose()
M=nn*A

correctly outputs the desired results. Differently from Matlab/Octave, a:b means from a to b-1 and not from a to b (the right extremum is excluded). Also, matrix has to be declared with the right type of entries (complex, CC, in this case), otherwise they are assumed to be integers, and in the case they are not, an error message is displayed (something that has the meaning of attempted coercition of non integral number in an integer).