Computations of discrete orthogonal polynomials with Wolfram Mathematica
(Krisztina Kiss, Ilona Nagy)

Last modified: 13 July, 2016.

The file consists of the following chapters:

1. Norms of discrete orthogonal polynomials (short version but slow method)
2. Norms of discrete orthogonal polynomials (long version but fast method)
3. The (Subscript[p, m,j]) polynomials
4. The (Subscript[Ξ, 2]) matrix
5. Lambda matrix
6. (Subscript[Z, m]) matrices

In chapters 1 and 2 we provide two diferent methods to calculate the norm squares of the (Subscript[p, m,j]) polynomials for m=1,...,6 and m=1,...,5 and at the end of the chapters we give a list of the results in a product form. Chapters 1 and 2 correspond to section 7 Appendix. We used M instead of N in formula (2) because N is a protected symbol in Mathematica.
Chapter 3 corresponds to section 7 Appendix. Chapter 4 corresponds to Remark 1.
Chapter 5 corresponds to section 3 and 5.1. Chapter 6 corresponds to section 2.3 and section 7 Appendix.

1. Norms of discrete orthogonal polynomials (short version but slow method)

Definitions

Polynomials

n=100;
Maincoeff=Flatten[Table[a[i,j,j]->1,{i,1,n},{j,0,n}]];

alpha[m_,n_]:=(((2^n)Product[2i-1,{i,1,n}])/(n!))(1/(Product[M+i,{i,1,n}]))


p[m_,n_,x_]:=alpha[m,n]Sum[a[m,n,k](x^k),{k,0,n}]/.Maincoeff
p[m_,i_,0]:=alpha[m,0]a[m,i,0]/.Maincoeff

Scalar product for m=1,...,6

sum[1,i_,j_]:=sum[1,i,j]=Sum[p[1,i,x] p[1,j,x],{x,0,M-1}]/.Maincoeff
sum[2,i_,j_]:=sum[2,i,j]=Sum[p[2,i,x] p[2,j,x],{x2,0,M-1},{x,0,x2}]/.Maincoeff
sum[3,i_,j_]:=sum[3,i,j]=Sum[p[3,i,x] p[3,j,x],{x3,0,M-1},{x2,0,x3},{x,0,x2}]/.Maincoeff
sum[4,i_,j_]:=sum[4,i,j]=Sum[p[4,i,x] p[4,j,x],{x4,0,M-1},{x3,0,x4},{x2,0,x3},{x,0,x2}]/.Maincoeff
sum[5,i_,j_]:=sum[5,i,j]=Sum[p[5,i,x] p[5,j,x],{x5,0,M-1},{x4,0,x5},{x3,0,x4},{x2,0,x3},{x,0,x2}]/.Maincoeff
sum[6,i_,j_]:=sum[6,i,j]=Sum[p[6,i,x] p[6,j,x],{x6,0,M-1},{x5,0,x6},{x4,0,x5},{x3,0,x4},{x2,0,x3},{x,0,x2}]/.Maincoeff

List of the coefficients

coeff[m_,n_]:=Flatten[Table[a[m,j,i],{j,1,n},{i,0,j-1}],1]

List of scalar products for the orthogonal polynomials

eq[m_,n_]:=eq[m,n]=Flatten[Table[sum[m,i,j],{j,1,n},{i,0,j-1}],1]

Substitution rules for the coefficients

rules[m_,n_]:=rules[m,n]=Solve[eq[m,n]==0,coeff[m,n]];

Norm squares

normsquare[m_,n_]:=normsquare[m,n]=Table[sum[m,i,i],{i,0,n}]/.rules[m,n]//Factor

Calculations (up to degree 10)

The first output represents the total time to do the computation.

m=1

m1=AbsoluteTiming[normsquare[1,10]]

m=2

m2=AbsoluteTiming[normsquare[2,10]]

m=3

m3=AbsoluteTiming[normsquare[3,10]]

m=4

m4=AbsoluteTiming[normsquare[4,10]]

m=5

m5=AbsoluteTiming[normsquare[5,10]]

m=6

m6=AbsoluteTiming[normsquare[6,10]]

2. Norms of discrete orthogonal polynomials (long version but fast method)

Formulas

Quit;

ClearAll[S,list,listi,listj,listk,listl,listm,listn,listM,s,i,j,k,l,m,n,M, Maincoeff,p,alpha];
S[i_]:=Sum[(n^i),{n,0,s}]

list=Table[S[i],{i,0,16}];
listi=list/.s->i;
listj=list/.s->j;
listk=list/.s->k;
listl=list/.s->l;
listm=list/.s->m;
listn=list/.s->n;
listM=list/.s->M-1;

Maincoeff={a[1,0,0]->1,a[2,0,0]->1,a[3,0,0]->1,a[4,0,0]->1,a[5,0,0]->1,a[6,0,0]->1,a[7,0,0]->1,a[8,0,0]->1,
a[1,1,1]->1,a[1,2,2]->1,a[1,3,3]->1,a[1,4,4]->1,a[1,5,5]->1,a[1,6,6]->1,a[1,7,7]->1,a[1,8,8]->1,
a[2,1,1]->1,a[2,2,2]->1,a[2,3,3]->1,a[2,4,4]->1,a[2,5,5]->1,a[2,6,6]->1,a[2,7,7]->1,a[2,8,8]->1,
a[3,1,1]->1,a[3,2,2]->1,a[3,3,3]->1,a[3,4,4]->1,a[3,5,5]->1,a[3,6,6]->1,a[3,7,7]->1,a[3,8,8]->1,
a[4,1,1]->1,a[4,2,2]->1,a[4,3,3]->1,a[4,4,4]->1,a[4,5,5]->1,a[4,6,6]->1,a[4,7,7]->1,a[4,8,8]->1,
a[5,1,1]->1,a[5,2,2]->1,a[5,3,3]->1,a[5,4,4]->1,a[5,5,5]->1,a[5,6,6]->1,a[5,7,7]->1,a[5,8,8]->1,
a[6,1,1]->1,a[6,2,2]->1,a[6,3,3]->1,a[6,4,4]->1,a[6,5,5]->1,a[6,6,6]->1,a[6,7,7]->1,a[6,8,8]->1};

p[m_,n_,x_]:=alpha[m,n]Sum[a[m,n,k](x^k),{k,0,n}]
p[m_,i_,0]:=alpha[m,0]a[m,i,0]
alpha[m_,n_]:=(((2^n)Product[2i-1,{i,1,n}])/(n!))(1/(Product[M+i,{i,1,n}]))

m=1

Calculations

Program

sum1[y0_,z0_]:=Module[{y=y0,z=z0,poly0,poly1,scoeff,spowers,icoeff,ipowers},
poly0=y z;
scoeff=CoefficientList[poly0,s];
spowers=Take[listM,Length[scoeff]];

scoeff.spowers
]

Basis polynomials

ClearAll[poly10,poly20,poly21,poly30,poly31,poly32,poly40,poly41,poly42,poly43,poly50,poly51,poly52,poly53,poly54,poly60,poly61,poly62,poly63,poly64,poly65,poly70,poly71,poly72,poly73,poly74,poly75,poly76,poly80,poly81,poly82,poly83,poly84,poly85,poly86,poly87];
poly10=sum1[1,p[1,1,s]]; poly20=sum1[1,p[1,2,s]]; poly21=sum1[p[1,1,s],p[1,2,s]]; poly30=sum1[1,p[1,3,s]]; poly31=sum1[p[1,1,s],p[1,3,s]]; poly32=sum1[p[1,2,s],p[1,3,s]]; poly40=sum1[1,p[1,4,s]]; poly41=sum1[p[1,1,s],p[1,4,s]];
poly42=sum1[p[1,2,s],p[1,4,s]]; poly43=sum1[p[1,3,s],p[1,4,s]]; poly50=sum1[1,p[1,5,s]]; poly51=sum1[p[1,1,s],p[1,5,s]];
poly52=sum1[p[1,2,s],p[1,5,s]]; poly53=sum1[p[1,3,s],p[1,5,s]]; poly54=sum1[p[1,4,s],p[1,5,s]]; poly60=sum1[1,p[1,6,s]];
poly61=sum1[p[1,1,s],p[1,6,s]]; poly62=sum1[p[1,2,s],p[1,6,s]]; poly63=sum1[p[1,3,s],p[1,6,s]]; poly64=sum1[p[1,4,s],p[1,6,s]];
poly65=sum1[p[1,5,s],p[1,6,s]]; poly70=sum1[1,p[1,7,s]]; poly71=sum1[p[1,1,s],p[1,7,s]]; poly72=sum1[p[1,2,s],p[1,7,s]];
poly73=sum1[p[1,3,s],p[1,7,s]]; poly74=sum1[p[1,4,s],p[1,7,s]]; poly75=sum1[p[1,5,s],p[1,7,s]]; poly76=sum1[p[1,6,s],p[1,7,s]];
poly80=sum1[1,p[1,8,s]]; poly81=sum1[p[1,1,s],p[1,8,s]]; poly82=sum1[p[1,2,s],p[1,8,s]]; poly83=sum1[p[1,3,s],p[1,8,s]];
poly84=sum1[p[1,4,s],p[1,8,s]]; poly85=sum1[p[1,5,s],p[1,8,s]]; poly86=sum1[p[1,6,s],p[1,8,s]]; poly87=sum1[p[1,7,s],p[1,8,s]];

The coefficients as a function of M

rules1=Solve[{poly10,poly20,poly21,poly30,poly31,poly32,poly40,poly41,poly42,poly43,poly50,poly51,poly52,poly53,poly54,poly60,poly61,poly62,poly63,poly64,poly65,poly70,poly71,poly72,poly73,poly74,poly75,poly76,poly80,poly81,poly82,poly83,poly84,poly85,poly86,poly87}==0,{a[1,1,0],a[1,2,1],a[1,2,0],a[1,3,2],a[1,3,1],a[1,3,0],a[1,4,3],a[1,4,2],a[1,4,1],a[1,4,0],a[1,5,4],a[1,5,3],a[1,5,2],a[1,5,1],a[1,5,0],a[1,6,5],a[1,6,4],a[1,6,3],a[1,6,2],a[1,6,1],a[1,6,0],a[1,7,6],a[1,7,5],a[1,7,4],a[1,7,3],a[1,7,2],a[1,7,1],a[1,7,0],a[1,8,7],a[1,8,6],a[1,8,5],a[1,8,4],a[1,8,3],a[1,8,2],a[1,8,1],a[1,8,0]}][[1]]

Norm squares (for m=1 up to degree 8)

m1normsquare={sum1[1,1],
sum1[p[1,1,s],p[1,1,s]],
sum1[p[1,2,s],p[1,2,s]],
sum1[p[1,3,s],p[1,3,s]],
sum1[p[1,4,s],p[1,4,s]],
sum1[p[1,5,s],p[1,5,s]],
sum1[p[1,6,s],p[1,6,s]],
sum1[p[1,7,s],p[1,7,s]],
sum1[p[1,8,s],p[1,8,s]]}/.rules1/.Maincoeff//Factor

m=2

Calculations

Program

sum2[y0_,z0_]:=Module[{y=y0,z=z0,poly0,poly1,scoeff,spowers,icoeff,ipowers},
poly0=y z;
scoeff=CoefficientList[poly0,s];
spowers=Take[listi,Length[scoeff]];

poly1=scoeff.spowers;(*sum 1*)
icoeff=CoefficientList[poly1,i];
ipowers=Take[listM,Length[icoeff]]//Expand;

icoeff.ipowers//Expand(*sum 2*)
]

Basis polynomials

ClearAll[poly10,poly20,poly21,poly30,poly31,poly32,poly40,poly41,poly42,poly43,poly50,poly51,poly52,poly53,poly54,poly60,poly61,poly62,poly63,poly64,poly65,poly70,poly71,poly72,poly73,poly74,poly75,poly76,poly80,poly81,poly82,poly83,poly84,poly85,poly86,poly87];
poly10=sum2[1,p[2,1,s]]; poly20=sum2[1,p[2,2,s]]; poly21=sum2[p[2,1,s],p[2,2,s]]; poly30=sum2[1,p[2,3,s]]; poly31=sum2[p[2,1,s],p[2,3,s]]; poly32=sum2[p[2,2,s],p[2,3,s]]; poly40=sum2[1,p[2,4,s]]; poly41=sum2[p[2,1,s],p[2,4,s]]; poly42=sum2[p[2,2,s],p[2,4,s]]; poly43=sum2[p[2,3,s],p[2,4,s]]; poly50=sum2[1,p[2,5,s]]; poly51=sum2[p[2,1,s],p[2,5,s]]; poly52=sum2[p[2,2,s],p[2,5,s]]; poly53=sum2[p[2,3,s],p[2,5,s]]; poly54=sum2[p[2,4,s],p[2,5,s]]; poly60=sum2[1,p[2,6,s]]; poly61=sum2[p[2,1,s],p[2,6,s]]; poly62=sum2[p[2,2,s],p[2,6,s]]; poly63=sum2[p[2,3,s],p[2,6,s]]; poly64=sum2[p[2,4,s],p[2,6,s]]; poly65=sum2[p[2,5,s],p[2,6,s]]; poly70=sum2[1,p[2,7,s]]; poly71=sum2[p[2,1,s],p[2,7,s]]; poly72=sum2[p[2,2,s],p[2,7,s]]; poly73=sum2[p[2,3,s],p[2,7,s]]; poly74=sum2[p[2,4,s],p[2,7,s]]; poly75=sum2[p[2,5,s],p[2,7,s]]; poly76=sum2[p[2,6,s],p[2,7,s]]; poly80=sum2[1,p[2,8,s]]; poly81=sum2[p[2,1,s],p[2,8,s]]; poly82=sum2[p[2,2,s],p[2,8,s]]; poly83=sum2[p[2,3,s],p[2,8,s]]; poly84=sum2[p[2,4,s],p[2,8,s]]; poly85=sum2[p[2,5,s],p[2,8,s]]; poly86=sum2[p[2,6,s],p[2,8,s]]; poly87=sum2[p[2,7,s],p[2,8,s]];

The coefficients as a function of M

rules2=Solve[{poly10,poly20,poly21,poly30,poly31,poly32,poly40,poly41,poly42,poly43,poly50,poly51,poly52,poly53,poly54,poly60,poly61,poly62,poly63,poly64,poly65,poly70,poly71,poly72,poly73,poly74,poly75,poly76,poly80,poly81,poly82,poly83,poly84,poly85,poly86,poly87}==0,{a[2,1,0],a[2,2,1],a[2,2,0],a[2,3,2],a[2,3,1],a[2,3,0],a[2,4,3],a[2,4,2],a[2,4,1],a[2,4,0],a[2,5,4],a[2,5,3],a[2,5,2],a[2,5,1],a[2,5,0],a[2,6,5],a[2,6,4],a[2,6,3],a[2,6,2],a[2,6,1],a[2,6,0],a[2,7,6],a[2,7,5],a[2,7,4],a[2,7,3],a[2,7,2],a[2,7,1],a[2,7,0],a[2,8,7],a[2,8,6],a[2,8,5],a[2,8,4],a[2,8,3],a[2,8,2],a[2,8,1],a[2,8,0]}][[1]]

Norm squares (for m=2 up to degree 7)

m2normsquare={sum2[1,1],
sum2[p[2,1,s],p[2,1,s]],
sum2[p[2,2,s],p[2,2,s]],
sum2[p[2,3,s],p[2,3,s]],
sum2[p[2,4,s],p[2,4,s]],
sum2[p[2,5,s],p[2,5,s]],
sum2[p[2,6,s],p[2,6,s]],
sum2[p[2,7,s],p[2,7,s]]}/.rules2/.Maincoeff//Factor

m=3

Calculations

Program

sum3[y0_,z0_]:=Module[{y=y0,z=z0,poly0,poly1,poly2,scoeff,spowers,icoeff,ipowers,jcoeff,jpowers},
poly0=y z;
scoeff=CoefficientList[poly0,s];
spowers=Take[listi,Length[scoeff]];

poly1=scoeff.spowers;(*sum 1*)
icoeff=CoefficientList[poly1,i];
ipowers=Take[listj,Length[icoeff]]//Expand;

poly2=icoeff.ipowers;(*sum 2*)
jcoeff=CoefficientList[poly2,j];
jpowers=Take[listM,Length[jcoeff]]//Expand;

jcoeff.jpowers//Expand(*sum 3*)
]

Basis polynomials

ClearAll[poly10,poly20,poly21,poly30,poly31,poly32,poly40,poly41,poly42,poly43,poly50,poly51,poly52,poly53,poly54,
poly60,poly61,poly62,poly63,poly64,poly65,poly70,poly71,poly72,poly73,poly74,poly75,poly76,poly80,poly81,poly82,poly83,poly84,poly85,poly86,poly87];
poly10=sum3[1,p[3,1,s]]; poly20=sum3[1,p[3,2,s]]; poly21=sum3[p[3,1,s],p[3,2,s]]; poly30=sum3[1,p[3,3,s]];
poly31=sum3[p[3,1,s],p[3,3,s]]; poly32=sum3[p[3,2,s],p[3,3,s]]; poly40=sum3[1,p[3,4,s]]; poly41=sum3[p[3,1,s],p[3,4,s]];
poly42=sum3[p[3,2,s],p[3,4,s]]; poly43=sum3[p[3,3,s],p[3,4,s]]; poly50=sum3[1,p[3,5,s]]; poly51=sum3[p[3,1,s],p[3,5,s]];
poly52=sum3[p[2,3,s],p[3,5,s]]; poly53=sum3[p[3,3,s],p[3,5,s]]; poly54=sum3[p[3,4,s],p[3,5,s]]; poly60=sum3[1,p[3,6,s]];
poly61=sum3[p[3,1,s],p[3,6,s]]; poly62=sum3[p[3,2,s],p[3,6,s]]; poly63=sum3[p[3,3,s],p[3,6,s]]; poly64=sum3[p[3,4,s],p[3,6,s]];
poly65=sum3[p[3,5,s],p[3,6,s]]; poly70=sum3[1,p[3,7,s]]; poly71=sum3[p[3,1,s],p[3,7,s]]; poly72=sum3[p[3,2,s],p[3,7,s]];
poly73=sum3[p[3,3,s],p[3,7,s]]; poly74=sum3[p[3,4,s],p[3,7,s]]; poly75=sum3[p[3,5,s],p[3,7,s]]; poly76=sum3[p[3,6,s],p[3,7,s]];

The coefficients as a function of M

rules3=Solve[{poly10,poly20,poly21,poly30,poly31,poly32,poly40,poly41,poly42,poly43,poly50,poly51,poly52,poly53,poly54,poly60,poly61,poly62,poly63,poly64,poly65,poly70,poly71,poly72,poly73,poly74,poly75,poly76}==0,{a[3,1,0],a[3,2,1],a[3,2,0],a[3,3,2],a[3,3,1],a[3,3,0],a[3,4,3],a[3,4,2],a[3,4,1],a[3,4,0],a[3,5,4],a[3,5,3],a[3,5,2],a[3,5,1],a[3,5,0],a[3,6,5],a[3,6,4],a[3,6,3],a[3,6,2],a[3,6,1],a[3,6,0],a[3,7,6],a[3,7,5],a[3,7,4],a[3,7,3],a[3,7,2],a[3,7,1],a[3,7,0]}][[1]]

Norm squares (for m=3 up to degree 7)

m3normsquare={sum3[1,1],
sum3[p[3,1,s],p[3,1,s]],
sum3[p[3,2,s],p[3,2,s]],
sum3[p[3,3,s],p[3,3,s]],
sum3[p[3,4,s],p[3,4,s]],
sum3[p[3,5,s],p[3,5,s]],
sum3[p[3,6,s],p[3,6,s]],
sum3[p[3,7,s],p[3,7, s]]}/.rules3/.Maincoeff//Factor

m=4

Calculations

Program

sum4[y0_,z0_]:=Module[{y=y0,z=z0,poly0,poly1,poly2,poly3,scoeff,spowers,icoeff,ipowers,jcoeff,jpowers,kcoeff,kpowers},
poly0=y z;
scoeff=CoefficientList[poly0,s];
spowers=Take[listi,Length[scoeff]];

poly1=scoeff.spowers;(*sum 1*)
icoeff=CoefficientList[poly1,i];
ipowers=Take[listj,Length[icoeff]]//Expand;

poly2=icoeff.ipowers;(*sum 2*)
jcoeff=CoefficientList[poly2,j];
jpowers=Take[listk,Length[jcoeff]]//Expand;

poly3=jcoeff.jpowers;(*sum 3*)
kcoeff=CoefficientList[poly3,k];
kpowers=Take[listM,Length[kcoeff]]//Expand;

kcoeff.kpowers//Expand(*sum 4*)
]

Basis polynomials

ClearAll[poly10,poly20,poly21,poly30,poly31,poly32,poly40,poly41,poly42,poly43,poly50,poly51,poly52,poly53,poly54,poly60,poly61,poly62,poly63,poly64,poly65,poly70,poly71,poly72,poly73,poly74,poly75,poly76,poly80,poly81,poly82,poly83,poly84,poly85,poly86,poly87];
poly10=sum4[1,p[4,1,s]]; poly20=sum4[1,p[4,2,s]]; poly21=sum4[p[4,1,s],p[4,2,s]]; poly30=sum4[1,p[4,3,s]]; poly31=sum4[p[4,1,s],p[4,3,s]]; poly32=sum4[p[4,2,s],p[4,3,s]]; poly40=sum4[1,p[4,4,s]]; poly41=sum4[p[4,1,s],p[4,4,s]]; poly42=sum4[p[4,2,s],p[4,4,s]]; poly43=sum4[p[4,3,s],p[4,4,s]]; poly50=sum4[1,p[4,5,s]]; poly51=sum4[p[4,1,s],p[4,5,s]]; poly52=sum4[p[4,2,s],p[4,5,s]]; poly53=sum4[p[4,3,s],p[4,5,s]]; poly54=sum4[p[4,4,s],p[4,5,s]];
poly60=sum4[1,p[4,6,s]]; poly61=sum4[p[4,1,s],p[4,6,s]]; poly62=sum4[p[4,2,s],p[4,6,s]]; poly63=sum4[p[4,3,s],p[4,6,s]];
poly64=sum4[p[4,4,s],p[4,6,s]]; poly65=sum4[p[4,5,s],p[4,6,s]]; poly70=sum4[1,p[4,7,s]]; poly71=sum4[p[4,1,s],p[4,7,s]];
poly72=sum4[p[4,2,s],p[4,7,s]]; poly73=sum4[p[4,3,s],p[4,7,s]]; poly74=sum4[p[4,4,s],p[4,7,s]]; poly75=sum4[p[4,5,s],p[4,7,s]];poly76=sum4[p[4,6,s],p[4,7,s]];

The coefficients as a function of M

rules4=Solve[{poly10,poly20,poly21,poly30,poly31,poly32,poly40,poly41,poly42,poly43,poly50,poly51,poly52,poly53,poly54,poly60,poly61,poly62,poly63,poly64,poly65,poly70,poly71,poly72,poly73,poly74,poly75,poly76}==0,{a[4,1,0],a[4,2,1],a[4,2,0],a[4,3,2],a[4,3,1],a[4,3,0],a[4,4,3],a[4,4,2],a[4,4,1],a[4,4,0],a[4,5,4],a[4,5,3],a[4,5,2],a[4,5,1],a[4,5,0],a[4,6,5],a[4,6,4],a[4,6,3],a[4,6,2],a[4,6,1],a[4,6,0],a[4,7,6],a[4,7,5],a[4,7,4],a[4,7,3],a[4,7,2],a[4,7,1],a[4,7,0]}][[1]]

Norm squares (for m=4 up to degree 7)

m4normsquare={sum4[1,1],
sum4[p[4,1,s],p[4,1,s]],
sum4[p[4,2,s],p[4,2,s]],
sum4[p[4,3,s],p[4,3,s]],
sum4[p[4,4,s],p[4,4,s]],
sum4[p[4,5,s],p[4,5,s]],
sum4[p[4,6,s],p[4,6,s]]}/.rules4/.Maincoeff//Factor

m=5

Calculations

Program

sum5[y0_,z0_]:=Module[{y=y0,z=z0,poly0,poly1,poly2,poly3,poly4,scoeff,spowers,icoeff,ipowers,jcoeff,jpowers,kcoeff,kpowers,lcoeff,lpowers},
poly0=y z;
scoeff=CoefficientList[poly0,s];
spowers=Take[listi,Length[scoeff]];

poly1=scoeff.spowers;(*sum 1*)
icoeff=CoefficientList[poly1,i];
ipowers=Take[listj,Length[icoeff]]//Expand;

poly2=icoeff.ipowers;(*sum 2*)
jcoeff=CoefficientList[poly2,j];
jpowers=Take[listk,Length[jcoeff]]//Expand;

poly3=jcoeff.jpowers;(*sum 3*)
kcoeff=CoefficientList[poly3,k];
kpowers=Take[listl,Length[kcoeff]]//Expand;

poly4=kcoeff.kpowers;(*sum 4*)
lcoeff=CoefficientList[poly4,l];
lpowers=Take[listM,Length[lcoeff]]//Expand;

lcoeff.lpowers//Expand
]

Basis polynomials

ClearAll[poly10,poly20,poly21,poly30,poly31,poly32,poly40,poly41,poly42,poly43,poly50,poly51,poly52,poly53,poly54,poly60,poly61,poly62,poly63,poly64,poly65,poly70,poly71,poly72,poly73,poly74,poly75,poly76,poly80,poly81,poly82,poly83,poly84,poly85,poly86,poly87];
poly10=sum5[1,p[5,1,s]]; poly20=sum5[1,p[5,2,s]]; poly21=sum5[p[5,1,s],p[5,2,s]]; poly30=sum5[1,p[5,3,s]];
poly31=sum5[p[5,1,s],p[5,3,s]]; poly32=sum5[p[5,2,s],p[5,3,s]]; poly40=sum5[1,p[5,4,s]]; poly41=sum5[p[5,1,s],p[5,4,s]]; poly42=sum5[p[5,2,s],p[5,4,s]]; poly43=sum5[p[5,3,s],p[5,4,s]]; poly50=sum5[1,p[5,5,s]]; poly51=sum5[p[5,1,s],p[5,5,s]]; poly52=sum5[p[5,2,s],p[5,5,s]]; poly53=sum5[p[5,3,s],p[5,5,s]]; poly54=sum5[p[5,4,s],p[5,5,s]];
poly60=sum5[1,p[5,6,s]]; poly61=sum5[p[5,1,s],p[5,6,s]]; poly62=sum5[p[5,2,s],p[5,6,s]]; poly63=sum5[p[5,3,s],p[5,6,s]];
poly64=sum5[p[5,4,s],p[5,6,s]];poly65=sum5[p[5,5,s],p[5,6,s]];

The coefficients as a function of M

rules5=Solve[{poly10,poly20,poly21,poly30,poly31,poly32,poly40,poly41,poly42,poly43,
poly50,poly51,poly52,poly53,poly54,poly60,poly61,poly62,poly63,poly64,poly65}==0,{a[5,1,0],a[5,2,1],a[5,2,0],a[5,3,2],a[5,3,1],a[5,3,0],a[5,4,3],a[5,4,2],a[5,4,1],a[5,4,0],a[5,5,4],a[5,5,3],a[5,5,2],a[5,5,1],a[5,5,0],a[5,6,5],a[5,6,4],a[5,6,3],a[5,6,2],a[5,6,1],a[5,6,0]}][[1]]

Norm squares (for m=5 up to degree 7)

m5normsquare={sum5[1,1],
sum5[p[5,1,s],p[5,1,s]],
sum5[p[5,2,s],p[5,2,s]],
sum5[p[5,3,s],p[5,3,s]],
sum5[p[5,4,s],p[5,4,s]],
sum5[p[5,5,s],p[5,5,s]],
sum5[p[5,6,s],p[5,6,s]]}/.rules5/.Maincoeff//Factor

Norm squares

m1normsquare

m2normsquare

m3normsquare

m4normsquare

m5normsquare

3. The (Subscript[p, m,j]) polynomials

m=1, j=0,...,5

Table[p[1,i,x]/.rules1/.Maincoeff,{i,0,5}]

m=2, j=0,...,5

Table[p[2,i,x]/.rules2/.Maincoeff,{i,0,5}]

m=3, j=0,...,5

Table[p[3,i,x]/.rules3/.Maincoeff,{i,0,5}]

4. The (Subscript[Ξ, 2]) matrix

G1 and G2 respectively denote the matrices G(1,(Subscript[ν, 1])) and G(2,(Subscript[ν, 2])) in Remark 1. The (Subscript[Ξ, 2]) matrix is given in list form and matrix form.

ClearAll[r,q]
r[M_,m_,i_]:=m!Binomial[M-1+m-i,m];
q[M_,m_,j_,i_]:=r[M,m-1,i]p[m,j,i]/.rules1/.rules2/.Maincoeff

G1={CoefficientList[p[1,0,x]/.rules1/.Maincoeff,x]~Join~{0,0,0,0,0,0},
CoefficientList[p[1,1,x]/.rules1/.Maincoeff,x]~Join~{0,0,0,0,0},
CoefficientList[p[1,2,x]/.rules1/.Maincoeff,x]~Join~{0,0,0,0},
CoefficientList[p[1,3,x]/.rules1/.Maincoeff,x]~Join~{0,0,0},
CoefficientList[p[1,4,x]/.rules1/.Maincoeff,x]~Join~{0,0},
CoefficientList[p[1,5,x]/.rules1/.Maincoeff,x]~Join~{0},
CoefficientList[p[1,6,x]/.rules1/.Maincoeff,x]}

G2={CoefficientList[q[M,2,0,x]/.rules2/.Maincoeff,x]~Join~{0,0,0,0,0},
CoefficientList[q[M,2,1,x]/.rules2/.Maincoeff,x]~Join~{0,0,0,0},
CoefficientList[q[M,2,2,x]/.rules2/.Maincoeff,x]~Join~{0,0,0},
CoefficientList[q[M,2,3,x]/.rules2/.Maincoeff,x]~Join~{0,0},
CoefficientList[q[M,2,4,x]/.rules2/.Maincoeff,x]~Join~{0},
CoefficientList[q[M,2,5,x]/.rules2/.Maincoeff,x]}//Factor

G2.Inverse[G1]//Factor

G2.Inverse[G1]//Factor//MatrixForm

5. Lambda matrix

Case m=1

Creating the lambda matrix

c[1,l_,1]:=p[1,l,M-1]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff
c[1,l_,0]:=-p[1,l,-1]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff
szumma[l_]:=Sum[L[l,s]p[1,s,j],{s,0,l}]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff

l=0

CoefficientList[szumma[0],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff

lambda0={1,-1,1,0,0,0,0,0,0,0}

l=1

CoefficientList[szumma[1],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff

righthandside=CoefficientList[szumma[1],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff
lefthandside=CoefficientList[p[1,1,j-1],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff
sol1=Solve[righthandside==lefthandside,{L[1,0],L[1,1]}]

lambda1={c[1,1,1],c[1,1,0],L[1,0],L[1,1],0,0,0,0,0,0}/.sol1[[1]]//Factor

l=2

s1=CoefficientList[szumma[2],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff
s2=CoefficientList[p[1,2,j-1],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff
sol2=Solve[s1==s2,{L[2,0],L[2,1],L[2,2]}]

lambda2={c[1,2,1],c[1,2,0],L[2,0],L[2,1],L[2,2],0,0,0,0,0}/.sol2[[1]]//Factor

l=3

s1=CoefficientList[szumma[3],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff
s2=CoefficientList[p[1,3,j-1],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff
sol3=Solve[s1==s2,{L[3,0],L[3,1],L[3,2],L[3,3]}]

lambda3={c[1,3,1],c[1,3,0],L[3,0],L[3,1],L[3,2],L[3,3],0,0,0,0}/.sol3[[1]]//Factor

l=4

s1=CoefficientList[szumma[4],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff;
s2=CoefficientList[p[1,4,j-1],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff;
sol4=Solve[s1==s2,{L[4,0],L[4,1],L[4,2],L[4,3],L[4,4]}]

lambda4={c[1,4,1],c[1,4,0],L[4,0],L[4,1],L[4,2],L[4,3],L[4,4],0,0,0}/.sol4[[1]]//Factor

l=5

s1=CoefficientList[szumma[5],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff;
s2=CoefficientList[p[1,5,j-1],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff;
sol5=Solve[s1==s2,{L[5,0],L[5,1],L[5,2],L[5,3],L[5,4],L[5,5]}]

lambda5={c[1,5,1],c[1,5,0],L[5,0],L[5,1],L[5,2],L[5,3],L[5,4],L[5,5],0,0}/.sol5[[1]]//Factor

l=6

s1=CoefficientList[szumma[6],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff;
s2=CoefficientList[p[1,6,j-1],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff;
sol6=Solve[s1==s2,{L[6,0],L[6,1],L[6,2],L[6,3],L[6,4],L[6,5],L[6,6]}]

lambda6={c[1,6,1],c[1,6,0],L[6,0],L[6,1],L[6,2],L[6,3],L[6,4],L[6,5],L[6,6],0}/.sol6[[1]]//Factor

l=7

s1=CoefficientList[szumma[7],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff;
s2=CoefficientList[p[1,7,j-1],j]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff;
sol7=Solve[s1==s2,{L[7,0],L[7,1],L[7,2],L[7,3],L[7,4],L[7,5],L[7,6],L[7,7]}]

lambda7={c[1,7,1],c[1,7,0],L[7,0],L[7,1],L[7,2],L[7,3],L[7,4],L[7,5],L[7,6],L[7,7]}/.sol7[[1]]//Factor

The lambda matrix in list form:

{lambda0,lambda1,lambda2,lambda3,lambda4,lambda5,lambda6,lambda7}

6. (Subscript[Z, m]) matrices

Notations: In formula (3) (Subscript[r, N,m-1])(i)  is denoted by r[M_,m_,i_]. In formula (28) (Subscript[(OverTilde[q]), N,m+j-2])(i) is denoted by qq[M_,m_,j_,i_]. In formula (29) (Subscript[ζ, j,l]^m) is denoted by xi[m,j,l] and (Subscript[(UnderBar[ζ]), j]^m) is denoted by XI[m,j].

m≥1

ClearAll[r,qq,M,c]
r[M_,m_,i_]:=m!Binomial[M-1+m-i,m];
qq[M_,m_,j_,i_]:=r[M,m-1,i-1]p[m,j,i-1]-r[M,m-1,i]p[m,j,i]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff

ClearAll[szumma2];
szumma2[m_,j_,i_]:=Sum[xi[m,j,l]p[1,l,i],{l,0,m+j-2}]/.rules1/.rules2/.rules3/.rules4/.rules5/.Maincoeff;

ClearAll[c,M];
c[m_,j_,1]:=r[M,m-1,M-1]p[m,j,M-1]/.rules1/.rules2/.rules3/.rules4/.rules5;
c[m_,j_,0]:=-r[M,m-1,-1]p[m,j,-1]/.rules1/.rules2/.rules3/.rules4/.rules5;

Substitution rules for the coefficients (Subscript[ζ, j,l]^m) for m=1,...,5.

ClearAll[l,m];
Zrule1[m_,l_]:=Solve[CoefficientList[qq[M,m,l,j],j]==CoefficientList[szumma2[m,l,j],j],Table[xi[m,l,i],{i,0,l-1}]]

ClearAll[l,m];
Zrule2[m_,l_]:=Solve[CoefficientList[qq[M,m,l,j],j]==CoefficientList[szumma2[m,l,j],j],Table[xi[m,l,i],{i,0,l}]]

ClearAll[l,m];
Zrule3[m_,l_]:=Solve[CoefficientList[qq[M,m,l,j],j]==CoefficientList[szumma2[m,l,j],j],Table[xi[m,l,i],{i,0,l+1}]]

ClearAll[l,m];
Zrule4[m_,l_]:=Solve[CoefficientList[qq[M,m,l,j],j]==CoefficientList[szumma2[m,l,j],j],Table[xi[m,l,i],{i,0,l+2}]]

ClearAll[l,m];
Zrule5[m_,l_]:=Solve[CoefficientList[qq[M,m,l,j],j]==CoefficientList[szumma2[m,l,j],j],Table[xi[m,l,i],{i,0,l+3}]]

m=1

Calculation of the row vector (Subscript[(UnderBar[ζ]), j]^1):

XI[1,0]={c[1,0,1],c[1,0,0],0,0,0,0}/.Zrule1[1,0]/.Maincoeff//Factor

XI[1,1]={c[1,1,1],c[1,1,0],xi[1,1,0],0,0,0}/.Zrule1[1,1]/.Maincoeff//Factor

XI[1,2]={c[1,2,1],c[1,2,0],xi[1,2,0],xi[1,2,1],0,0}/.Zrule1[1,2]/.Maincoeff//Factor

XI[1,3]={c[1,3,1],c[1,3,0],xi[1,3,0],xi[1,3,1],xi[1,3,2],0}/.Zrule1[1,3]/.Maincoeff//Factor

XI[1,4]={c[1,4,1],c[1,4,0],xi[1,4,0],xi[1,4,1],xi[1,4,2],xi[1,4,3]}/.Zrule1[1,4]/.Maincoeff//Factor

The (Subscript[Z, 1]) matrix in list form and matrix form:

Table[XI[1,i][[1]],{i,0,4}]//Factor

Table[XI[1,i][[1]],{i,0,4}]//Factor//MatrixForm

m=2

Calculation of the row vector (Subscript[(UnderBar[ζ]), j]^2):

XI[2,0]={c[2,0,1],c[2,0,0],xi[2,0,0],0,0,0}/.Zrule2[2,0]/.Maincoeff//Factor

XI[2,1]={c[2,1,1],c[2,1,0],xi[2,1,0],xi[2,1,1],0,0}/.Zrule2[2,1]/.Maincoeff//Factor

XI[2,2]={c[2,2,1],c[2,2,0],xi[2,2,0],xi[2,2,1],xi[2,2,2],0}/.Zrule2[2,2]/.Maincoeff//Factor

XI[2,3]={c[2,3,1],c[2,3,0],xi[2,3,0],xi[2,3,1],xi[2,3,2],xi[2,3,3]}/.Zrule2[2,3]/.Maincoeff//Factor

The (Subscript[Z, 2]) matrix in list form and matrix form:

Table[XI[2,i][[1]],{i,0,3}]//Factor

Table[XI[2,i][[1]],{i,0,3}]//Factor//MatrixForm

m=3

Calculation of the row vector (Subscript[(UnderBar[ζ]), j]^3):

XI[3,0]={c[3,0,1],c[3,0,0],xi[3,0,0],xi[3,0,1],0,0}/.Zrule3[3,0]/.Maincoeff//Factor

XI[3,1]={c[3,1,1],c[3,1,0],xi[3,1,0],xi[3,1,1],xi[3,1,2],0}/.Zrule3[3,1]/.Maincoeff//Factor

XI[3,2]={c[3,2,1],c[3,2,0],xi[3,2,0],xi[3,2,1],xi[3,2,2],xi[3,2,3]}/.Zrule3[3,2]/.Maincoeff//Factor

The (Subscript[Z, 3]) matrix in list form and matrix form:

Table[XI[3,i][[1]],{i,0,2}]//Factor

Table[XI[3,i][[1]],{i,0,2}]//Factor//MatrixForm

m=4

Calculation of the row vector (Subscript[(UnderBar[ζ]), j]^4):

XI[4,0]={c[4,0,1],c[4,0,0],xi[4,0,0],xi[4,0,1],xi[4,0,2],0}/.Zrule4[4,0]/.Maincoeff//Factor

XI[4,1]={c[4,1,1],c[4,1,0],xi[4,1,0],xi[4,1,1],xi[4,1,2],xi[4,1,3]}/.Zrule4[4,1]/.Maincoeff//Factor

The (Subscript[Z, 4]) matrix in list form:

Table[XI[4,i][[1]],{i,0,1}]//Factor

m=5

The (Subscript[Z, 5]) matrix in list form:

XI[5,0]={c[5,0,1],c[5,0,0],xi[5,0,0],xi[5,0,1],xi[5,0,2],xi[5,0,3]}/.Zrule5[5,0]/.Maincoeff//Factor