polysem.lib

version="1.0";
category="Miscellaneous";
info="
LIBRARY:   polysem.lib  Polynomial semantics of linear logic
AUTHOR:    Daniel Murfet

PROCEDURES:
";

// NOTE: We include a version of matrix.lib which suppresses some
// unnecessary output from the procedure "rowred"
LIB "linalg.lib";
LIB "matrix.lib";
LIB "ring.lib";

////////////////////////////////////////////////////////////////////
// USAGE GUIDE
//
	
////////////////////////////////////////////////////////////////////
// sem_binaryseq
//
// Given a binary sequence S encoded as an intvec with entires 0,1
// returns the product of the given matrices A, B with A standing
// for 0 and B standing for 1, e.g. given 001 returns A * A * B

proc sem_binaryseq(matrix A, matrix B, intvec S)
{
	if( ncols(A) != nrows(A) || ncols(B) != nrows(B) || ncols(A) != ncols(B) )
	{
		print("[sem_binaryseq] Bad matrices, exiting.");
		return();
	}	
	
	int k = ncols(A);
	matrix tS[k][k] = unitmat(k);

	int i;
	for(i=size(S);i>=1;i--)
	{
		if( S[i] == 0 )
		{
			tS = A * tS;
		}
		else
		{
			tS = B * tS;
		}
	}
	
	return(tS);
}

////////////////////////////////////////////////////////////////////
// matrix_to_vect
//
// Given a matrix of polynomials, returns it as a list of coefficients.
// We must be given the total weight of monomials L to truncate at.
//
// The ordering in which we write out the coefficients is as follows.
// The outer loop is the set of monomials of weight <= L, then the row,
// then the column. The monomials are ordered by grevlex.

proc matrix_to_vect(matrix A, int L)
{
	int numVars = nvars(basering);
	list monoms = monomialdict(L, numVars);
	
	list vectorCoeffs = matrix_to_vect_givenmonoms(A,L,monoms);
	
	return(vectorCoeffs);
}

proc matrix_to_vect_givenmonoms(matrix A, int L, list monoms)
{
	int numVars = nvars(basering);
	int Q = size(monoms);

	poly xprod = 1;
	int i,j;
	for(i=1; i<=numVars; i++)
	{
    	xprod = xprod * var(i);
	}    
    
	list vectorCoeffs;
	matrix koffer;
	int a, b;
	for(a=1;a<=Q;a++)
	{
    	for(i=1;i<=k;i++)
	    {
        	for(j=1;j<=k;j++)
	        {
    	    	koffer = coef(A[i,j], xprod);
        		int found = 0;
        	
	        	for(b=1; b<=ncols(koffer); b++)
    	        {   
					if( leadexp(koffer[1,b]) == monoms[a] )
					{
						found = 1;
						vectorCoeffs = vectorCoeffs + list( number(koffer[2,b]) );   
					}
				}
            
				if( found == 0 )
				{
					vectorCoeffs = vectorCoeffs + list( 0 );
				}
        	}
    	}
	}

	return(vectorCoeffs);
}

////////////////////////////////////////////////////////////////////
// monomialdict
//
// Given integers k >= 1 and nv >= 1 returns all intvecs of length
// nv with sum <= k.
//
// The order is by
//    - sum, from zero up to k
//    - for a fixed sum, in lexicographic order
//
// That is, the monomials are ordered according to grevlex

proc monomialdict(int k, int nv)
{
 	int i;
 	list monoms;
 	for(i=0; i<=k; i++)
 	{
		monoms = monoms + partitions(i, nv);
	}
	
	return(monoms);
}

////////////////////////////////////////////////////////////////////
// partitions
//
// Given integers N and M with N >= 0 and M > 0 returns all elements
// in NN^M which sum to N (here NN means non-negative integers) as a list
// of intvecs.

proc partitions(int N, int M)
{
	if( N < 0 || M <= 0 )
	{
		print("[partitions] Bad integers, exiting.");
		return();
	}
	
	if( M == 1 )
	{
		return(list(intvec(N)));
	}
	
	list P;
    int i, j;
    
    for( i=0; i<=N; i++)
    {
    	list miniP = partitions(N-i, M-1);
    
    	for( j = 1; j<=size(miniP); j++ )
    	{
    		P = P + list( concat_intvec( intvec(i), miniP[j] ) );
    	}
    }
    
    return(P);
}

////////////////////////////////////////////////////////////////////
// displayIndices

proc displayIndices(int k)
{
	int numVars = nvars(basering);
	list monoms = monomialdict(k, numVars);
 	
 	print("Monomials of total degree <= " + string(k));
	print("");
 	
 	int i;
 	for(i=1;i<=size(monoms);i++)
 	{
 		print(string(i) + "    " + string(monomial(monoms[i])) + "    [" + string(monoms[i]) + "]");
 	}
}

////////////////////////////////////////////////////////////////////
// seqsfrominterval
//
// Given an integer N >=1 and M >= 1 return the set of all sequences
// of integers of length M with entries in 1,...,N

proc seqsfrominterval(int N, int M)
{
	if( M == 1 )
	{
		list L;
		int i;
		for(i=1;i<=N;i++)
		{
			L = L + list(intvec(i));
		}
		return(L);
	}
	
	list L_small = seqsfrominterval(N, M-1);
	
	list L;
	int i,j,k;
	intvec v;
	
	for(i=1;i<=N;i++)
	{
		for(j=1;j<=size(L_small);j++)
		{
			v[1] = i;
			
			for(k=1;k<=M-1;k++)
			{
				v[k+1] = (L_small[j])[k];
			}
			
			L = L + list(v);
		}
	}
	
	return(L);
}

proc concat_intvec( intvec L, intvec R )
{
	int n = size(L);
	int m = size(R);
	intvec ret;
	
	int i;
	for(i=1;i<=n;i++)
	{
		ret[i] = L[i];
	}
	
	for(i=1;i<=m;i++ )
	{
		ret[n+i] = R[i];
	}
	
	return(ret);
}