#######################################################
## make-win-package tcore
## Mon Jun 12 09:41:50 EDT 2023
#######################################################
#######################################################

printf("BEGIN tcore package\n");
printf("THIS VERSION DATED Mon Jun 12 09:41:50 EDT 2023 \n");
tcore:=table();
tcore[PHI1]:=proc(ptn,p)  #ptn ->[pcore, pquotient]
	local tcr,tqu,bgv:
	tcr:=tcore[tcoreofptn](ptn,p);
	tqu:=tcore[tquot](ptn,p);
	bgv:=[tcr,tqu];
	RETURN(bgv):
end:


#HERE
##testptn:=[1, 1, 2, 4, 4, 5, 6, 6, 6, 7, 7, 7, 7, 7, 8, 9, 9, 10, 10, 13, 16, 17];
##tcr:=fcoreofptn(testptn);
##tqu:=tquot(testptn,5);
##bgv:=[tcr,tqu];
##
################################################################################
###11/03/03:
##
##ptns4:=partition(4):
##tqu1:=[[2],[2],[2],[2],[2]]:
##tqu2:=[[1+1],[2],[2],[2],[2]]:
##tqu3:=[[1+1],[1+1],[2],[2],[2]]:
##tqu4:=[[1+1],[2],[1+1],[2],[2]]:
##
##for j from 1 to 5 do
##  bgv:=[ptns4[j],tqu1]:
##  lprint(bgv);
##  ptn:=invphi1(bgv,5):
##  lprint(ptn);
##  lprint(tcrank(ptn,5),srank(ptn));
##od:
##
###++++++++++++++++++++++++++++++++++++++++++++++++++++
tcore[addrimcell]:=proc(L) 
      local n,i,j,pp,diagarray,np,Li,Ni,Nj,newdiagarray:
## Input = L = [diagarray,[i,j]]
## [i,j] is position of current cell
## We put * in next cell on the rim
      diagarray:=L[1]:
      i:=L[2][1]:
      j:=L[2][2]:
      np:=nops([indices(diagarray)]):
      ## Find position of next cell
      Li:=nops(diagarray[i]):
      if j < Li then
         Ni:=i:
         Nj:=j+1:
      else
         Ni:=i-1:
         Nj:=j:
      fi:
      newdiagarray:=diagarray:
      newdiagarray[Ni][Nj]:="*":
      RETURN([newdiagarray,[Ni,Nj]]):
end:        


tcore[addrimthook]:=proc(ptn,j,L,t)
local len,nmp,cj,nj,prenewptn,newptn,k:
    # add rim t-hook starting from jth part
    len:=0:
    nmp:=tcore[np](ptn):
    if j>nmp then
       #add dummy zeros
       prenewptn:=[op(ptn),seq(0,k=1..(j-nmp))]:
    else
       prenewptn:=ptn:
    fi:
    ################################################################
    #11/11/03: Convert newptn to array. Apparently maple complains
    #          if list length > 100.
    #Eg:> b:=[seq(0,i=1..101)];
    #   > b[101]:=b[101]+1;    
    #Error, assigning to a long list, please use arrays
    ################################################################
    newptn:=convert(prenewptn,array):	
    cj:=j:  #current part
    while len<L*t do
       if len=0 then
          newptn[cj]:=newptn[cj]+1:
       else
          if cj=1 then
             newptn[cj]:=newptn[cj]+1:
          else
             nj:=cj-1:
             if newptn[nj]<newptn[cj] then
                newptn[nj]:=newptn[nj]+1:
                cj:=nj:
             else
                 newptn[cj]:=newptn[cj]+1:
             fi:
           fi:
        fi:
        len:=len+1:
     od;
 #convert back to list
 RETURN(convert(newptn,list)):
end:
##ptn:=[11,7,4,4,2,2,2,0,0,0,0,0,0,0,0,0];findhookinpos(ptn,5,4,2);
##  ptn2:=addrimthook(ptn,16,2,5);findhookinpos(ptn2,5,4,2);
##  ptn3:=addrimthook(ptn2,12,2,5);findhookinpos(ptn3,5,4,1);
##  ptn4:=addrimthook(ptn3,5,1,5);findhookinpos(ptn4,5,4,1);
##  ptn5:=addrimthook(ptn4,3,1,5);findhookinpos(ptn5,5,4,1);
##  ptn6:=addrimthook(ptn5,2,1,5);findhookinpos(ptn6,5,3,3);
##  ptn7:=addrimthook(ptn6,22,3,5);findhookinpos(ptn7,5,3,3);
##  ptn8:=addrimthook(ptn7,18,3,5);findhookinpos(ptn8,5,3,2);
##  ptn9:=addrimthook(ptn8,11,2,5);findhookinpos(ptn9,5,2,3);
##  ptn10:=addrimthook(ptn9,16,3,5);findhookinpos(ptn10,5,1,5);
##  ptn8:=addrimthook(ptn7,18,3,5);findhookinpos(ptn8,5,3,2);
##  ptn11:=addrimthook(ptn10,20,5,5);findhookinpos(ptn11,5,1,2);
##  ptn12:=addrimthook(ptn11,7,2,5);findhookinpos(ptn12,5,1,1);
##  ptn13:=addrimthook(ptn12,2,1,5);
tcore[avec2nvec]:=proc(avec)
 local t,N4,nvec:
 t:=nops(avec):
 N4 := [[1, -1, 0, 0, 0], [0, 1, -1, 0, 0], [0, 0, 1, -1, 0], 
        [0, 0, 0, 1, -1], [1, 1, 0, -1, -1]]: 
 if t=5 then
    nvec:=veccombo(avec,N4):
    RETURN(nvec):
else
    ERROR("Only t=5 implemented");
fi:
end:
 tcore[aveccyc]:=proc(av)
    local t,newav,j:
    # cyclically permute components of avec
    t:=nops(av):
    newav:=[av[t],seq(av[j],j=1..t-1)]:
 RETURN(newav):
 end:
  tcore[darray2ptn]:=proc(L)
        local n,i,j,pp,np,jp,ptn:
  ## Input = L = diagarray
        np:=nops([indices(L)]):
        ## Find position of next cell
        ptn:=[]:
        j:=np:
        while j >= 1 do
          ptn:=[op(ptn),nops(L[j])]:
          j := j - 1:
        od:
        RETURN(ptn):
  end:
  tcore[findcell]:=proc(darray,t,r,k)
  local np,f,i,darray0,j,parti,L;
  np:=nops([indices(darray)]):
  f:=1:
  i:=1:
  darray0:=array(1..np+t):
  ##darray0 is darray plus some extre empty parts
  for j from 1 to np do
  darray0[j]:=darray[j]:
  od:
  for j from np+1 to np+t do
  darray0[j]:=[]:
  od:
  while f=1 do
    parti:=nops(darray0[i]):
    j:=t*(r-1)+k+i:
    if j<parti then
      L:=[i,j]:
    else
      f:=0:
    fi:
    i:=i+1:
  od:
  RETURN(L);
  end:
tcore[findhookinpos]:=proc(ptn,t,w,p)
   local nmp,mnp,find,j,fj,prenewptn,newptn,k;
    # find part number where can insert rim pt-hook from word W_w
    nmp:=tcore[np](ptn): 

    #add p*t dummy zeros just in case
 
 prenewptn:=[op(ptn),seq(0,k=1..(p*t))]:
 newptn:=convert(prenewptn,array): #Convert to array
     
 ## FIND pth consecutive E in W_w
 mnp:=nops(newptn):
 find := 0: j:=1: fj:=0: 
 while find<p do
   if modp(newptn[j]-j,t)=w then
      fj:=j:
      find:=find+1:
   fi:
   j:=j+1:
 od:
 RETURN(fj);
end:
tcore[freqtab]:=proc(L)
	local fset,i,M,j:
	fset:={}:
	for i from 1 to nops(L) do
	   fset:=fset union {L[i]}:
	   if type(M[L[i]],integer) then
	      M[L[i]]:=M[L[i]]+1:
	   else
	      M[L[i]]:=1:
	   fi:
	od:
	for j in fset do
	lprint(j,M[j]);
	od;
end:
tcore[invphi1]:=proc(bigvec,t)   
    #Inverse of GSK phi1 bijection
local tc,ntc,ptn,tq,j,tqj,r,ntqj,k,p,posi,newptn,nn,optn;
    tc:=bigvec[1]:
    ntc:=nops(tc):
    ptn:=[seq(tc[ntc-j+1],j=1..ntc)]:
    tq:=bigvec[2]:
    for j from t by -1 to 1 do
        tqj:=tq[j];
        r:=j-1:  #residue
	ntqj:=nops(tqj):
	if ntqj>0 then
	  for k from ntqj by -1 to 1 do
 	    p:=tqj[k]:
	    posi:=tcore[findhookinpos](ptn,t,r,p); 
	    ##print("j=",j,"p=",p,"posi=",posi);
            newptn:=tcore[addrimthook](ptn,posi,p,t);
	    ##print("newptn=",newptn);
	    ptn:=newptn:
	  od:
	fi:
    od:
    nn:=tcore[np](ptn):
    optn:=[seq(ptn[nn-j+1],j=1..nn)]:
    RETURN(optn):
end:
tcore[ispos]:=proc(n)
    # Returns true if n > 0
        if n>0 then
           true;
        else
           false;
        fi:
end:

tcore[istcore]:=proc(ptn,p) ##determines if ptn is a p-core
       local x,i,R,i1:
        x:=0:
        for i from 0 to p-1 do
           R[i]:=tcore[rvec](ptn,p,i):
        od:
        for i from 0 to p-1 do
           i1:=modp(i+1,p):
           x:=x + R[i]^2 - R[i]*R[i1]:
        od:
        if x=R[0] then
           RETURN(true);
        else
           RETURN(false);
        fi:
end:
tcore[lp] :=proc(ptn) #largest part
     ptn[np(ptn)]:
end:
tcore[makebiw]:=proc(ptn,t,mj)
local np,reg,i,k,ip,j,a,b,r,hh,pp;
## Make the bi-infinite words W[0],W[1],...W[t-1]
## with j from -mj to mj
## EXAMPLE:
##  > ptn := [1, 1, 2, 4, 4, 5, 6, 6, 6];
##  > tcore[makebiw](ptn,5,3);
##     -3-2-1 0 1 2 3
##  W0  E E E E N E N
##  W1  E E N N E N N
##  W2  E E E N N N N
##  W3  E E E E E N N
##  W4  E E N E E N N
##  
        np:=nops(ptn):
        reg:=array(-100..100,0..t-1):
        for i from -100 to 100 do
          for k from 0 to t-1 do
             reg[i,k]:=" N":
          od:
	od:
	ip:=[seq(ptn[np-k],k=0..(np-1)),seq(0,k=0..(mj+1)*t)]:
	for j from 1 to nops(ip) do
           a:=ip[j]:
           b:=modp(a-j,t):
	   r:=floor( (a-j)/t ) + 1:
           reg[r,b]:=" E":
	od:
        hh:=cat("   ",seq(r,r=-mj..-1),seq(cat(" ",r),r=0..mj)):
        printf("%s\n",hh);
	for i from 0 to t-1 do
           pp:=cat("W",i," ",seq(reg[r,i],r=-mj ..mj)):
           printf("%s\n",pp);
	od:
        RETURN():
end:
  tcore[markrimhookV2]:=proc(L,N)      
        local n,i,j,pp,diagarray,k,newL;
  ## Input = L = [diagarray,[i,j]]
  ## Mark a rim-hook of length N starting at position [i,j]
  ## [i,j] is position of current cell
        diagarray:=L[1]: i:=L[2][1]: j:=L[2][2]:
        for k from 1 to N do
        ##lprint("k=",k);
             ##lprint("i=",i,"j=",j);
             newL:=tcore[addrimcell]([diagarray,[i,j]]):
             diagarray:=newL[1]:
             i:=newL[2][1]:
             j:=newL[2][2]:
        od:
        RETURN(diagarray):
  end:
tcore[nep]:=proc(ptn)
    local x,i,p:
    #number of even parts
    x:=0:
    for i from 1 to nops(ptn) do
      p:=ptn[i]:
      if modp(p,2)=0 then
         x:=x+1:
      fi:
    od:
RETURN(x):
end:
tcore[nepo]:=proc(ptn)
    local y:
    y:=nep(ptn):
    if modp(y,2)=1 then
       RETURN(true):
    else
       RETURN(false):
    fi:
end:
tcore[np] :=proc(ptn) #number of parts
    nops(ptn):
end:
tcore[numnepo]:=proc(n)
    local ptns,x,j;
    ptns:=combinat[partition](n):
    x:=0:
    for j from 1 to nops(ptns) do
       if nepo(ptns[j]) then
          x:=x+1:
       fi:
    od:
    RETURN(x):
end:
tcore[nvec2alphavec]:=proc(nvec)  #only for t=5
local n0,n1,n2,n3,n4,n5,n6,a3,a0,a2,a1,a4,a5,a6,t;
    t:=nops(nvec):
    if t=5 then
	n0:=nvec[1]:
	n1:=nvec[2]:
	n2:=nvec[3]:
	n3:=nvec[4]:
	n4:=nvec[5]:
	a3:=1/5*n0+1/5+3/5*n2+2/5*n1+4/5*n3; 
        a0:=-1/5*n3+1/5*n0+1/5-2/5*n2-3/5*n1; 
        a2:=-1/5*n1-2/5*n3-3/5*n0+2/5+1/5*n2; 
        a1:=-4/5*n2-1/5*n1-2/5*n3-3/5*n0+2/5;
	a4:=1-a0-a1-a2-a3:
	RETURN([a0,a1,a2,a3,a4]):
    fi:
    if t=7 then
	n0:=nvec[1]:
	n1:=nvec[2]:
	n2:=nvec[3]:
	n3:=nvec[4]:
	n4:=nvec[5]:
	n5:=nvec[6]:
	n6:=nvec[7]:
	a3 := 6/7*n1+5/7*n3+1/7*n4+3/7*n0-1/7+2/7*n2+4/7*n5; 
        a4 := -1/7*n4-3/7*n0+1/7-2/7*n2-4/7*n5+2/7*n3+1/7*n1; 
        a0 := -1/7*n5+1/7*n0-2/7*n4+2/7-4/7*n2-3/7*n3-5/7*n1; 
        a2 := 3/7*n2+6/7*n5+1/7*n0+5/7*n4+2/7+4/7*n3+2/7*n1; 
	a1 := -2/7*n0+3/7+1/7*n2+2/7*n5-3/7*n4-1/7*n3-4/7*n1; 
        a5 := -5/7*n3-1/7*n4-3/7*n0+1/7-2/7*n2-4/7*n5+1/7*n1;
	a6:=1-a0-a1-a2-a3-a4-a5:
	RETURN([a0,a1,a2,a3,a4,a5,a6]):
    fi:
end:
tcore[nvec2ptn]:=proc(nvec)
     # Returns partition (t-core) with given n-vector
	local t,X,j,i,pparts,k,ppartsA,cnvec,ppartsB,partsA,pp2a,partsB;
	t:=nops(nvec):
	##print("t=",t);
	X:=0:
	for j from 1 to t do
	   i:=j-1:
	   if nvec[j]>0 then
	      X:=X+nvec[j]:
	   fi:
	od:
	## X = size of Durfee square
	##print("size of Durfee square = ",X);
	ppartsA:=[]:
	for j from 1 to t do
	   i:=j-1:
	   if nvec[j]>0 then
	      for k from 1 to nvec[j] do
	          ppartsA:=[op(ppartsA),t*(k-1)+i]:
	      od:
	   fi:
	od:
	ppartsA:=sort(ppartsA):
##	
## now look at conjugate ptn
##	
	for j from 1 to t do
	   cnvec[j]:=-nvec[t-j+1]:
	od:
	ppartsB:=[]:
	for j from 1 to t do
	   i:=j-1:
	   if cnvec[j]>0 then
	      for k from 1 to cnvec[j] do
	          ppartsB:=[op(ppartsB),t*(k-1)+i]:
	      od:
	   fi:
	od:
	ppartsB:=sort(ppartsB):
	##print("FROB PTN:");
	##print(ppartsA);
	##print(ppartsB);
	##print("");
	##print([seq(ppartsA[j]+X-j+1,j=1..X)]);
        partsA:=[seq(ppartsA[j]+X-j+1,j=1..X)];
	pp2a:=[seq(ppartsB[j]-j+1,j=1..X)];
	## remove zero parts
	pp2a:=select(ispos,pp2a);
	##if pp2a[1]=0 then
	## pp2a:=[seq(ppartsB[j]-j+1,j=2..X)];
        ##fi:
        partsB:=combinat[conjpart](pp2a);                          
	##print([seq(ppartsA[j]+X-j+1,j=1..X)]);
	##RETURN([seq(pparts2[j]+X-j+1,j=1..X)]);
	RETURN([op(partsB),op(partsA)]);
end:
tcore[printdarray]:=proc(L) 
      local n,i,j,pp,np,jp:
## Input = L = diagarray
      np:=nops([indices(L)]):
      ## Find position of next cell
      for i from 1 to np do
        jp:=nops(L[i]):
        pp:=cat(seq(cat(L[i][j]," ") ,j=1..jp)):
        printf("%s\n",pp);
      od:
      RETURN():
end:        
tcore[ptn2nvec]:=proc(ptn,p)
    local k:
    # Computes n-vector of given p-core
	[seq(tcore[rvec](ptn,p,k)-tcore[rvec](ptn,p,k+1),k=0..p-1)];
end:
tcore[ptn2rvec]:=proc(ptn,p) #compute r-vector of partition (mod p)
    local k:
    [seq(tcore[rvec](ptn,p,k),k=0..p-1)]:
end:
tcore[ptnnorm]:=proc(ptn)
   convert(ptn,`+`); #sum of parts
end:
tcore[randpcore]:=proc(p,num)
    # Generates a random p-core
local x,i,a,L,j;
	x:=0:
	for i from 1 to p-1 do
	  a:=modp(rand(),num)-trunc(num/2);
	  L[i]:=a:
	  x:=a+x:
        od:
	L[p]:=-x:
	L:=[seq(L[j],j=1..p)];
	tcore[nvec2ptn](L);   
end:
  tcore[removerimhook]:=proc(darraymarked)
     local np,ndarray,i,L,k,j,f,m;
     # remove hook corresponding to marks in darray
     np:=nops([indices(darraymarked)]):
     ndarray:=array(1..np):
     for i from 1 to np do
        L:=darraymarked[i]:
        # find first mark in row
        k:=infty:
        j:=1:
        f:=0:
        while f=0 and j<=nops(L) do
            if L[j]="*" then
               k:=j-1:
               f:=1:
            else
               j:=j+1:
            fi:
        od:
        if f=0 then
           k:=nops(L):
        fi:
        if k>0 then
           ndarray[i]:=[seq(L[m],m=1..k)]:
        else
           ndarray[i]:=[]:
        fi:
     od:
     RETURN(ndarray):
  end:
tcore[rvec]:=proc(ptn,p,k)  ##number of nodes in p-residue diagram of ptn
                     ##colored k.
        local x,m,nn,j:
        x:=0:
        nn:=tcore[np](ptn):
        for m from 1 to nn do
           j:=nn-m+1:
           x := x + trunc( (ptn[j] + modp(p-m-k,p))/p):
        od:
        RETURN(x):
end:
tcore[tcorechanges]:=proc()
printf("**************************************************************\n");
printf("*\n");
printf("*\n");
printf("*   partitions package version 0.1 -  Novewmber, 2008\n");
printf("*   partitions package version 0.2 -  Fri Jun  9 17:38:57 EDT 2023\n");
printf("*  This version tested on MAPLE 2022\n");
printf("*\n");
printf("*\n");
printf("*  Changes since previous version 0.1\n");
printf("*\n");
printf("*     * NEW FUNCTIONS:\n");
printf("*    addrimcell, avec2nvec, aveccyc, darray2ptn, findcell, freqtab\n");
printf("*    ispos, makebiw, markrimhookV2, nepo, numnepo, nvec2alphavec\n");
printf("*    printdarray, removerimhook, tquot, tresdiag, tresdiag2array\n");
printf("*    veccombo, tcorepversion, tcorechanges\n");
printf("*\n");
printf("*  Wed Jan 4 16:14:30 EST 2012: public version on qseries.org\n");
printf("*\n");
printf("*     * Fixed some bugs\n");
printf("        Made maple txt help files\n");
printf("*\n");
printf("*\n");
printf("*\n");
printf("**************************************************************\n");
RETURN():
end:
tcore[tcoreofptn]:=proc(ptn,p) ##compute t-core of ptn
        local km,nvec,k;
##	local r0,r1,r2,r3,r4,nvec;
##	r0:=rvec(ptn,5,0):
##        r1:=rvec(ptn,5,1):
##        r2:=rvec(ptn,5,2):
##        r3:=rvec(ptn,5,3):
##        r4:=rvec(ptn,5,4):
##        nvec:=[rvec(ptn,t,
##	nvec:=[r0-r1,r1-r2,r2-r3,r3-r4,r4-r0]:
        km:=k->modp(k,p):
        nvec:=[seq(tcore[rvec](ptn,p,km(k))-tcore[rvec](ptn,p,km(k+1)),k=0..(p-1))];
	tcore[nvec2ptn](nvec);
end:
tcore[tcorepversion]:=proc()
printf("****************************************************\n");
printf("*\n");
printf("*  tcore package version 0.2  \n");
printf("*  Fri Jun  9 17:38:57 EDT 2023\n");
printf("*  This version tested on MAPLE 2022\n"); 
printf("*\n");
printf("*  Please report any problems to fgarvan@ufl.edu\n");
printf("*  Previous versions:\n");
printf("      NONE           \n");
printf("*\n");
printf("*  Please report any problems to fgarvan@ufl.edu\n");
printf("*  Previous versions:\n");
printf("        0.1 - Nov 2008 \n");
printf("****************************************************\n");
RETURN():
end:
tcore[tcores]:=proc(p,n) #p-cores of n
        local ptns,i,L,nmp,j,ptn,rr,pcore,pcores:
        ptns:=combinat[partition](n):
        pcore:=ptn->tcore[istcore](ptn,p):
        select(pcore,ptns):
end:
tcore[tcrank]:=proc(ptn,p)  ##This gives tcore crank of a partition
                     ##Takes values 0,1,...,p-1
        local h,x,m,i,j,lam:
        h:=t->(t - (p-1)/2)^(p-3):
        x:=0:
        m:=tcore[np](ptn):
        for i from 1 to m do
           j:=m-i+1:
           lam:=ptn[j]:
           x := x + modp(  h(lam-i)-h(i-1), p):
        od:
        RETURN(modp(x,p)):
end:
tcore[tquot]:=proc(ptn,t) ##compute t-quotient of ptn
local nps,ptnz,lambdabar,lamhat,i,muibar,lamhati,j,ki,k,nik;
	nps:=tcore[np](ptn):
	ptnz:=[seq(ptn[nps-j+1],j=1..nps),seq(0,i=1..t)]:
	##print("ptnz=",ptnz);
	lambdabar:=[seq(ptnz[j]-j,j=1..(nps+t))]:
	##print("lambdabar=",lambdabar);
	lamhat:=[]:
	for i from 0 to (t-1) do
	    muibar:=[]:                      
	    lamhati:=[]:
	    for j from 1 to nps+t do
	      if modp(lambdabar[j],t)=i then
		muibar:=[op(muibar),lambdabar[j]]:
	      fi:
	    od:
	    ki:=nops(muibar):
	    for k from 1 to (ki-1) do
		nik:=floor(muibar[k]/t) - floor(muibar[k+1]/t) - 1:
	        if nik>0 then
	          lamhati:=[op(lamhati),seq(k,j=1..nik)]:
	        fi:
	    od:
	    lamhat:=[op(lamhat),lamhati]:
	od:
	RETURN(lamhat):
end:
tcore[tresdiag]:=proc(ptn,t) #t-residue diagram of partition
      local n,i,j,pp:
      n:=nops(ptn):
      i:=1:
      for j from n by -1 to 1 do
        pp:=cat(seq(cat(modp(j-i,t)," ") ,j=1..ptn[j])):
        printf("%s\n",pp);
        i := i + 1:
      od:
      RETURN():
end:        


tcore[tresdiag2array]:=proc(ptn,t) #t-residue diagram of partition
      local n,i,j,pp,diagarray:
      diagarray:=array(1..nops(ptn)):
      n:=nops(ptn):
      i:=1:
      for j from n by -1 to 1 do
        diagarray[i]:=[seq(modp(j-i,t),j=1..ptn[j])]:
        i := i + 1:
      od:
      RETURN(diagarray):
end:        


tcore[veccombo]:=proc(C,VL)
   local DV,newV,j,c,k:
   # If C=(c1,c2,..cm) and VL is a list vectors [v1,v2, ...  vm]
   # RETURNS  c1*v1 + c2*v2 + ... + cm*vm
   if nops(C)<>nops(VL) then
     ERROR(nops(C),"<>",nops(VL));
   else
     DV:=nops(VL[1]):
     newV:=[]:
     for j from 1 to DV do
        c:=add(C[k]*VL[k][j],k=1..nops(C)):
        newV:=[op(newV),c]:
     od:
     RETURN(newV):
   fi:
 end:
printf("TABLE TYPE tcore = %a\n",type(tcore,table));
printf("WARNING: tcore package will be save to location given by\n");
printf("Homelib\n");                                                    
savelib( tcore   
,
FileTools:-JoinPath([Homelib,"tcore.mla"]));
printf("END tcore package\n");
## mylib  above must be changed to the name of the directory
## in which you want stuff stored