Using Singular
Singular is a free and open-source computer algebra system for polynomial computation. It can be either downloaded or run in a browser from its official website: https://www.singular.uni-kl.de/
The Singular script generates the relations between the parameters which describe the target molecule in the Henderson-Hasselbalch model and the coefficients of its binding polynomial. Additionally, on this website, we present some sample experiments that can be run on its output.
Generating the binding polynomial equations
The script consists of a single main function, generateRingAndIdeal. Here is an example of how it is called:
bash-3.2$ Singular ligands.sing
Singular ligands.sing
SINGULAR / Development
A Computer Algebra System for Polynomial Computations / version 4.1.0
0<
by: W. Decker, G.-M. Greuel, G. Pfister, H. Schoenemann \ Nov 2016
FB Mathematik der Universitaet, D-67653 Kaiserslautern \
> int n1 = 4;
> int n2 = 3;
> def r = generateRingAndIdeal(n1,n2);
> r;
// coefficients: QQ
// number of vars : 31
// block 1 : ordering dp
// : names a(1)(1) a(1)(2) a(1)(3) a(2)(1) a(2)(2) a(2)(3) a(3)(1) a(3)(2) a(3)(3) a(4)(1) a(4)(2) a(4)(3) g(1) g(2) g(3) g(4) h(1) h(2) h(3) w(1)(1) w(1)(2) w(1)(3) w(2)(1) w(2)(2) w(2)(3) w(3)(1) w(3)(2) w(3)(3) w(4)(1) w(4)(2) w(4)(3)
// block 2 : ordering C
> setring r;
> I;
I[1]=g(1)*h(1)*w(1)(1)+g(1)*h(2)*w(1)(2)+g(1)*h(3)*w(1)(3)+g(2)*h(1)*w(2)(1)+g(2)*h(2)*w(2)(2)+g(2)*h(3)*w(2)(3)+g(3)*h(1)*w(3)(1)+g(3)*h(2)*w(3)(2)+g(3)*h(3)*w(3)(3)+g(4)*h(1)*w(4)(1)+g(4)*h(2)*w(4)(2)+g(4)*h(3)*w(4)(3)-a(1)(1)
I[2]=g(1)*h(1)*h(2)*w(1)(1)*w(1)(2)+g(1)*h(1)*h(3)*w(1)(1)*w(1)(3)+g(1)*h(2)*h(3)*w(1)(2)*w(1)(3)+g(2)*h(1)*h(2)*w(2)(1)*w(2)(2)+g(2)*h(1)*h(3)*w(2)(1)*w(2)(3)+g(2)*h(2)*h(3)*w(2)(2)*w(2)(3)+g(3)*h(1)*h(2)*w(3)(1)*w(3)(2)+g(3)*h(1)*h(3)*w(3)(1)*w(3)(3)+g(3)*h(2)*h(3)*w(3)(2)*w(3)(3)+g(4)*h(1)*h(2)*w(4)(1)*w(4)(2)+g(4)*h(1)*h(3)*w(4)(1)*w(4)(3)+g(4)*h(2)*h(3)*w(4)(2)*w(4)(3)-a(1)(2)
I[3]=g(1)*h(1)*h(2)*h(3)*w(1)(1)*w(1)(2)*w(1)(3)+g(2)*h(1)*h(2)*h(3)*w(2)(1)*w(2)(2)*w(2)(3)+g(3)*h(1)*h(2)*h(3)*w(3)(1)*w(3)(2)*w(3)(3)+g(4)*h(1)*h(2)*h(3)*w(4)(1)*w(4)(2)*w(4)(3)-a(1)(3)
I[4]=g(1)*g(2)*h(1)*w(1)(1)*w(2)(1)+g(1)*g(2)*h(2)*w(1)(2)*w(2)(2)+g(1)*g(2)*h(3)*w(1)(3)*w(2)(3)+g(1)*g(3)*h(1)*w(1)(1)*w(3)(1)+g(2)*g(3)*h(1)*w(2)(1)*w(3)(1)+g(1)*g(3)*h(2)*w(1)(2)*w(3)(2)+g(2)*g(3)*h(2)*w(2)(2)*w(3)(2)+g(1)*g(3)*h(3)*w(1)(3)*w(3)(3)+g(2)*g(3)*h(3)*w(2)(3)*w(3)(3)+g(1)*g(4)*h(1)*w(1)(1)*w(4)(1)+g(2)*g(4)*h(1)*w(2)(1)*w(4)(1)+g(3)*g(4)*h(1)*w(3)(1)*w(4)(1)+g(1)*g(4)*h(2)*w(1)(2)*w(4)(2)+g(2)*g(4)*h(2)*w(2)(2)*w(4)(2)+g(3)*g(4)*h(2)*w(3)(2)*w(4)(2)+g(1)*g(4)*h(3)*w(1)(3)*w(4)(3)+g(2)*g(4)*h(3)*w(2)(3)*w(4)(3)+g(3)*g(4)*h(3)*w(3)(3)*w(4)(3)-a(2)(1)
I[5]=g(1)*g(2)*h(1)*h(2)*w(1)(1)*w(1)(2)*w(2)(1)*w(2)(2)+g(1)*g(2)*h(1)*h(3)*w(1)(1)*w(1)(3)*w(2)(1)*w(2)(3)+g(1)*g(2)*h(2)*h(3)*w(1)(2)*w(1)(3)*w(2)(2)*w(2)(3)+g(1)*g(3)*h(1)*h(2)*w(1)(1)*w(1)(2)*w(3)(1)*w(3)(2)+g(2)*g(3)*h(1)*h(2)*w(2)(1)*w(2)(2)*w(3)(1)*w(3)(2)+g(1)*g(3)*h(1)*h(3)*w(1)(1)*w(1)(3)*w(3)(1)*w(3)(3)+g(2)*g(3)*h(1)*h(3)*w(2)(1)*w(2)(3)*w(3)(1)*w(3)(3)+g(1)*g(3)*h(2)*h(3)*w(1)(2)*w(1)(3)*w(3)(2)*w(3)(3)+g(2)*g(3)*h(2)*h(3)*w(2)(2)*w(2)(3)*w(3)(2)*w(3)(3)+g(1)*g(4)*h(1)*h(2)*w(1)(1)*w(1)(2)*w(4)(1)*w(4)(2)+g(2)*g(4)*h(1)*h(2)*w(2)(1)*w(2)(2)*w(4)(1)*w(4)(2)+g(3)*g(4)*h(1)*h(2)*w(3)(1)*w(3)(2)*w(4)(1)*w(4)(2)+g(1)*g(4)*h(1)*h(3)*w(1)(1)*w(1)(3)*w(4)(1)*w(4)(3)+g(2)*g(4)*h(1)*h(3)*w(2)(1)*w(2)(3)*w(4)(1)*w(4)(3)+g(3)*g(4)*h(1)*h(3)*w(3)(1)*w(3)(3)*w(4)(1)*w(4)(3)+g(1)*g(4)*h(2)*h(3)*w(1)(2)*w(1)(3)*w(4)(2)*w(4)(3)+g(2)*g(4)*h(2)*h(3)*w(2)(2)*w(2)(3)*w(4)(2)*w(4)(3)+g(3)*g(4)*h(2)*h(3)*w(3)(2)*w(3)(3)*w(4)(2)*w(4)(3)-a(2)(2)
I[6]=g(1)*g(2)*h(1)*h(2)*h(3)*w(1)(1)*w(1)(2)*w(1)(3)*w(2)(1)*w(2)(2)*w(2)(3)+g(1)*g(3)*h(1)*h(2)*h(3)*w(1)(1)*w(1)(2)*w(1)(3)*w(3)(1)*w(3)(2)*w(3)(3)+g(2)*g(3)*h(1)*h(2)*h(3)*w(2)(1)*w(2)(2)*w(2)(3)*w(3)(1)*w(3)(2)*w(3)(3)+g(1)*g(4)*h(1)*h(2)*h(3)*w(1)(1)*w(1)(2)*w(1)(3)*w(4)(1)*w(4)(2)*w(4)(3)+g(2)*g(4)*h(1)*h(2)*h(3)*w(2)(1)*w(2)(2)*w(2)(3)*w(4)(1)*w(4)(2)*w(4)(3)+g(3)*g(4)*h(1)*h(2)*h(3)*w(3)(1)*w(3)(2)*w(3)(3)*w(4)(1)*w(4)(2)*w(4)(3)-a(2)(3)
I[7]=g(1)*g(2)*g(3)*h(1)*w(1)(1)*w(2)(1)*w(3)(1)+g(1)*g(2)*g(3)*h(2)*w(1)(2)*w(2)(2)*w(3)(2)+g(1)*g(2)*g(3)*h(3)*w(1)(3)*w(2)(3)*w(3)(3)+g(1)*g(2)*g(4)*h(1)*w(1)(1)*w(2)(1)*w(4)(1)+g(1)*g(3)*g(4)*h(1)*w(1)(1)*w(3)(1)*w(4)(1)+g(2)*g(3)*g(4)*h(1)*w(2)(1)*w(3)(1)*w(4)(1)+g(1)*g(2)*g(4)*h(2)*w(1)(2)*w(2)(2)*w(4)(2)+g(1)*g(3)*g(4)*h(2)*w(1)(2)*w(3)(2)*w(4)(2)+g(2)*g(3)*g(4)*h(2)*w(2)(2)*w(3)(2)*w(4)(2)+g(1)*g(2)*g(4)*h(3)*w(1)(3)*w(2)(3)*w(4)(3)+g(1)*g(3)*g(4)*h(3)*w(1)(3)*w(3)(3)*w(4)(3)+g(2)*g(3)*g(4)*h(3)*w(2)(3)*w(3)(3)*w(4)(3)-a(3)(1)
I[8]=g(1)*g(2)*g(3)*h(1)*h(2)*w(1)(1)*w(1)(2)*w(2)(1)*w(2)(2)*w(3)(1)*w(3)(2)+g(1)*g(2)*g(3)*h(1)*h(3)*w(1)(1)*w(1)(3)*w(2)(1)*w(2)(3)*w(3)(1)*w(3)(3)+g(1)*g(2)*g(3)*h(2)*h(3)*w(1)(2)*w(1)(3)*w(2)(2)*w(2)(3)*w(3)(2)*w(3)(3)+g(1)*g(2)*g(4)*h(1)*h(2)*w(1)(1)*w(1)(2)*w(2)(1)*w(2)(2)*w(4)(1)*w(4)(2)+g(1)*g(3)*g(4)*h(1)*h(2)*w(1)(1)*w(1)(2)*w(3)(1)*w(3)(2)*w(4)(1)*w(4)(2)+g(2)*g(3)*g(4)*h(1)*h(2)*w(2)(1)*w(2)(2)*w(3)(1)*w(3)(2)*w(4)(1)*w(4)(2)+g(1)*g(2)*g(4)*h(1)*h(3)*w(1)(1)*w(1)(3)*w(2)(1)*w(2)(3)*w(4)(1)*w(4)(3)+g(1)*g(3)*g(4)*h(1)*h(3)*w(1)(1)*w(1)(3)*w(3)(1)*w(3)(3)*w(4)(1)*w(4)(3)+g(2)*g(3)*g(4)*h(1)*h(3)*w(2)(1)*w(2)(3)*w(3)(1)*w(3)(3)*w(4)(1)*w(4)(3)+g(1)*g(2)*g(4)*h(2)*h(3)*w(1)(2)*w(1)(3)*w(2)(2)*w(2)(3)*w(4)(2)*w(4)(3)+g(1)*g(3)*g(4)*h(2)*h(3)*w(1)(2)*w(1)(3)*w(3)(2)*w(3)(3)*w(4)(2)*w(4)(3)+g(2)*g(3)*g(4)*h(2)*h(3)*w(2)(2)*w(2)(3)*w(3)(2)*w(3)(3)*w(4)(2)*w(4)(3)-a(3)(2)
I[9]=g(1)*g(2)*g(3)*h(1)*h(2)*h(3)*w(1)(1)*w(1)(2)*w(1)(3)*w(2)(1)*w(2)(2)*w(2)(3)*w(3)(1)*w(3)(2)*w(3)(3)+g(1)*g(2)*g(4)*h(1)*h(2)*h(3)*w(1)(1)*w(1)(2)*w(1)(3)*w(2)(1)*w(2)(2)*w(2)(3)*w(4)(1)*w(4)(2)*w(4)(3)+g(1)*g(3)*g(4)*h(1)*h(2)*h(3)*w(1)(1)*w(1)(2)*w(1)(3)*w(3)(1)*w(3)(2)*w(3)(3)*w(4)(1)*w(4)(2)*w(4)(3)+g(2)*g(3)*g(4)*h(1)*h(2)*h(3)*w(2)(1)*w(2)(2)*w(2)(3)*w(3)(1)*w(3)(2)*w(3)(3)*w(4)(1)*w(4)(2)*w(4)(3)-a(3)(3)
I[10]=g(1)*g(2)*g(3)*g(4)*h(1)*w(1)(1)*w(2)(1)*w(3)(1)*w(4)(1)+g(1)*g(2)*g(3)*g(4)*h(2)*w(1)(2)*w(2)(2)*w(3)(2)*w(4)(2)+g(1)*g(2)*g(3)*g(4)*h(3)*w(1)(3)*w(2)(3)*w(3)(3)*w(4)(3)-a(4)(1)
I[11]=g(1)*g(2)*g(3)*g(4)*h(1)*h(2)*w(1)(1)*w(1)(2)*w(2)(1)*w(2)(2)*w(3)(1)*w(3)(2)*w(4)(1)*w(4)(2)+g(1)*g(2)*g(3)*g(4)*h(1)*h(3)*w(1)(1)*w(1)(3)*w(2)(1)*w(2)(3)*w(3)(1)*w(3)(3)*w(4)(1)*w(4)(3)+g(1)*g(2)*g(3)*g(4)*h(2)*h(3)*w(1)(2)*w(1)(3)*w(2)(2)*w(2)(3)*w(3)(2)*w(3)(3)*w(4)(2)*w(4)(3)-a(4)(2)
I[12]=g(1)*g(2)*g(3)*g(4)*h(1)*h(2)*h(3)*w(1)(1)*w(1)(2)*w(1)(3)*w(2)(1)*w(2)(2)*w(2)(3)*w(3)(1)*w(3)(2)*w(3)(3)*w(4)(1)*w(4)(2)*w(4)(3)-a(4)(3)
The command generateRingAndIdeals takes two integers \(n_1\), \(n_2\) as input and returns a ring in the following variables:
\(a_{i,j}\) represent the coefficient of the bivariate binding polynomial in front of the monomial with exponent vector \((i,j)\)
\(g_i\) represent the site energy of the \(i\)-th site for the first ligand type
\(h_j\) represent the site energy of the \(j\)-th site for the second ligand type
\(w_{i,j}\) represent the interaction energy between the \(i\)-th site for the first ligand type and the \(j\)-th site for the second ligand type
The ring contains an ideal \(I\), which is generated by the relations between the variables.
Moreover, it is possible to provide two additional integers:
bash-3.2$ Singular ligands.sing
SINGULAR / Development
A Computer Algebra System for Polynomial Computations / version 4.1.0
0<
by: W. Decker, G.-M. Greuel, G. Pfister, H. Schoenemann \ Nov 2016
FB Mathematik der Universitaet, D-67653 Kaiserslautern \
> int n1 = 4;
> int n2 = 3;
> int b = 99;
> int c = 101;
> def r = generateRingAndIdeal(n1,n2,b,c);
> r;
// coefficients: ZZ/101
// number of vars : 12
// block 1 : ordering dp
// : names w(1) w(2) w(3) w(4) w(5) w(6) w(7) w(8) w(9) w(10) w(11) w(12)
// block 2 : ordering C
> setring r;
> I;
I[1]=2*w(1)+35*w(2)+w(3)-42*w(4)-28*w(5)-21*w(6)+4*w(7)-31*w(8)+2*w(9)-28*w(10)+15*w(11)-14*w(12)+30
I[2]=w(1)*w(2)+26*w(1)*w(3)-50*w(2)*w(3)-21*w(4)*w(5)-41*w(4)*w(6)+40*w(5)*w(6)+2*w(7)*w(8)-49*w(7)*w(9)+w(8)*w(9)-14*w(10)*w(11)+40*w(10)*w(12)-7*w(11)*w(12)+29
I[3]=13*w(1)*w(2)*w(3)+30*w(4)*w(5)*w(6)+26*w(7)*w(8)*w(9)+20*w(10)*w(11)*w(12)+39
I[4]=-11*w(1)*w(4)-41*w(2)*w(5)+45*w(3)*w(6)-23*w(1)*w(7)-22*w(4)*w(7)-49*w(2)*w(8)+19*w(5)*w(8)+39*w(3)*w(9)-11*w(6)*w(9)-41*w(1)*w(10)-48*w(4)*w(10)+19*w(7)*w(10)+40*w(2)*w(11)-32*w(5)*w(11)-21*w(8)*w(11)+30*w(3)*w(12)-24*w(6)*w(12)-41*w(9)*w(12)+44
I[5]=45*w(1)*w(2)*w(4)*w(5)-42*w(1)*w(3)*w(4)*w(6)-28*w(2)*w(3)*w(5)*w(6)+39*w(1)*w(2)*w(7)*w(8)-11*w(4)*w(5)*w(7)*w(8)+4*w(1)*w(3)*w(7)*w(9)+17*w(4)*w(6)*w(7)*w(9)-31*w(2)*w(3)*w(8)*w(9)+45*w(5)*w(6)*w(8)*w(9)+30*w(1)*w(2)*w(10)*w(11)-24*w(4)*w(5)*w(10)*w(11)-41*w(7)*w(8)*w(10)*w(11)-28*w(1)*w(3)*w(10)*w(12)-18*w(4)*w(6)*w(10)*w(12)+45*w(7)*w(9)*w(10)*w(12)+15*w(2)*w(3)*w(11)*w(12)-12*w(5)*w(6)*w(11)*w(12)+30*w(8)*w(9)*w(11)*w(12)-10
I[6]=-21*w(1)*w(2)*w(3)*w(4)*w(5)*w(6)+2*w(1)*w(2)*w(3)*w(7)*w(8)*w(9)-42*w(4)*w(5)*w(6)*w(7)*w(8)*w(9)-14*w(1)*w(2)*w(3)*w(10)*w(11)*w(12)-9*w(4)*w(5)*w(6)*w(10)*w(11)*w(12)-28*w(7)*w(8)*w(9)*w(10)*w(11)*w(12)+19
I[7]=-25*w(1)*w(4)*w(7)+17*w(2)*w(5)*w(8)+38*w(3)*w(6)*w(9)-27*w(1)*w(4)*w(10)+17*w(1)*w(7)*w(10)+47*w(4)*w(7)*w(10)-18*w(2)*w(5)*w(11)+45*w(2)*w(8)*w(11)-36*w(5)*w(8)*w(11)+37*w(3)*w(6)*w(12)-42*w(3)*w(9)*w(12)-27*w(6)*w(9)*w(12)+43
I[8]=38*w(1)*w(2)*w(4)*w(5)*w(7)*w(8)-22*w(1)*w(3)*w(4)*w(6)*w(7)*w(9)+19*w(2)*w(3)*w(5)*w(6)*w(8)*w(9)+37*w(1)*w(2)*w(4)*w(5)*w(10)*w(11)-42*w(1)*w(2)*w(7)*w(8)*w(10)*w(11)-27*w(4)*w(5)*w(7)*w(8)*w(10)*w(11)-48*w(1)*w(3)*w(4)*w(6)*w(10)*w(12)+19*w(1)*w(3)*w(7)*w(9)*w(10)*w(12)+5*w(4)*w(6)*w(7)*w(9)*w(10)*w(12)-32*w(2)*w(3)*w(5)*w(6)*w(11)*w(12)-21*w(2)*w(3)*w(8)*w(9)*w(11)*w(12)+37*w(5)*w(6)*w(8)*w(9)*w(11)*w(12)-18
I[9]=-11*w(1)*w(2)*w(3)*w(4)*w(5)*w(6)*w(7)*w(8)*w(9)-24*w(1)*w(2)*w(3)*w(4)*w(5)*w(6)*w(10)*w(11)*w(12)-41*w(1)*w(2)*w(3)*w(7)*w(8)*w(9)*w(10)*w(11)*w(12)-48*w(4)*w(5)*w(6)*w(7)*w(8)*w(9)*w(10)*w(11)*w(12)+16
I[10]=-43*w(1)*w(4)*w(7)*w(10)+5*w(2)*w(5)*w(8)*w(11)+29*w(3)*w(6)*w(9)*w(12)-48
I[11]=29*w(1)*w(2)*w(4)*w(5)*w(7)*w(8)*w(10)*w(11)+47*w(1)*w(3)*w(4)*w(6)*w(7)*w(9)*w(10)*w(12)-36*w(2)*w(3)*w(5)*w(6)*w(8)*w(9)*w(11)*w(12)-28
I[12]=-27*w(1)*w(2)*w(3)*w(4)*w(5)*w(6)*w(7)*w(8)*w(9)*w(10)*w(11)*w(12)+46
If \(b>0\), the procedure will insert random integers in \([1,b]\) into \(a_{i,j}\) and \(g_i\), \(h_i\). This makes the \(n_1+n_2\) equations involving \(a_{i,0}\) and \(g_i\) resp. \(a_{0,j}\) and \(h_j\) obsolete, which are hence omitted.
If \(c>0\), the characteristic of the ground field is set to be the next bigger prime than \(c\). This speeds up computations but can also falsify results.
To write the generators of the ideal \(I\) into a file binding.sing, type (note that “:w” creates a new file or overwrites exisiting files while “:a” creates a new file or appends to existing files)
> write(":w binding.sing",I);
Computing mixed volumes in gfan
The mixed volumes of the generators of the ideal I can be computed as follows:
list P;
for (int i=1; i<=size(I); i++)
{
P[i] = newtonPolytope(I[i]);
vertices(P[i]);
}
mixedVolume(P);
Computing number of roots
The number of roots of the ideal I, counted with multiplicity, equals the dimension of the quotient ring \(\mathbb C[x]/I\) as a \(\mathbb C\) vector space. It can be computed as follows (hint: add “option(prot);” before the Groebner basis computation to receive some output during the computation):
ideal gbI = groebner(I);
vdim(gbI);
In-depth description of the script
First, we construct some helper functions to generate the polynomials that determine a(i)(j) in terms of the remaining variables.
WW(I,J) takes integers vectors I, J; it returns the product of all w(i)(j) with i in I and j in J
GG(I) takes an integer vector I; it returns the product of all g(i) with i in I
HH(J) takes an integer vector J; it returns the product of all h(j) with j in J
proc WW(intvec I, intvec J)
{
poly WIJ = 1;
int i,j;
for (i=1; i<=size(I); i++)
{
for (j=1; j<=size(J); j++)
{
WIJ = WIJ*w(I[i])(J[j]);
}
}
return (WIJ);
}
proc GG(intvec I)
{
poly GI = 1;
int i;
for (i=1; i<=size(I); i++)
{
GI = GI*g(I[i]);
}
return (GI);
}
proc HH(intvec J)
{
poly HJ = 1;
int j;
for (j=1; j<=size(J); j++)
{
HJ = HJ*h(J[j]);
}
return (HJ);
}
Moreover, we define a helper function which returns the binary expansion of an integer, e.g. input 5 yields output (1,0,1). This allows us to enumerate over all subsets of {1,2,…,n} by iterating over all integers 0,…,2^n-1.
proc intToBinary(int intN)
{
intvec binN;
for (int i=30; i>=0; i--)
{
if (intN >= 2^i)
{
intN = intN - 2^i;
binN[i+1] = 1;
}
}
intvec b;
for (i=1; i<=size(binN); i++)
{
if (binN[i]==1)
{
b[size(b)+1] = i;
}
}
return (intvec(b[2..size(b)]));
}
First, it begins with declaring a ring with the necessary variales. Here
int n1=5;
int n2=2;
ring r = 0,(a(1..n1)(1..n2),g(1..n1),h(1..n2),w(1..n1)(1..n2)),dp;
Finally, we construct some helper lists.
L is a list with n1 entries and its k-th entry contains all subsets of \(\{1,\ldots,n_1\}\) of cardinality k
M is a list with n2 entries and its l-th entry contains all subsets of \(\{1,\ldots,n_2\}\) of cardinality l
list L;
for (int i=1; i<=n1; i++)
{
L[i] = list();
}
intvec Li;
for (i=1; i<=2^n1-1; i++)
{
Li = intToBinary(i);
L[size(Li)] = L[size(Li)] + list(Li);
}
list M;
for (int j=1; j<=n2; j++)
{
M[j] = list();
}
intvec Mj;
for (j=1; j<=2^n2-1; j++)
{
Mj = intToBinary(j);
M[size(Mj)] = M[size(Mj)] + list(Mj);
}
Using this we can construct the polynomials.
int l,m;
ideal I;
poly gij;
for (i=1; i<=n1; i++)
{
for (j=1; j<=n2; j++)
{
gij = -a(i)(j);
for (l=1; l<=size(L[i]); l++)
{
for (m=1; m<=size(M[j]); m++)
{
gij = gij + poly(GG(L[i][l])*HH(M[j][m])*WW(L[i][l],M[j][m]));
}
}
I[size(I)+1] = gij;
}
}