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[addrimthook]:=proc(ptn,j,L,t)
local len,nmp,cj,nj,prenewptn,newptn:
    # 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[findhookinpos]:=proc(ptn,t,w,p)
   local nmp,mnp,find,j,fj,prenewptn,newptn;
    # 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[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[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[nep]:=proc(ptn)
    #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:
nepo:=proc(ptn)
    y:=nep(ptn):
    if modp(y,2)=1 then
       RETURN(true):
    else
       RETURN(false):
    fi:
end:

numnepo:=proc(n)
    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[np] :=proc(ptn) #number of parts
    nops(ptn):
end:

freqtab:=proc(L)
     #Nice for making tables
	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:

ispos:=proc(n)
    # This function is needed by nvec2ptn
    # Returns true if n > 0
        if n>0 then
           true;
        else
           false;
        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[ptn2nvec]:=proc(ptn,p)
    # 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)
    [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;
	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[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[tcoreofptn]:=proc(ptn,p) ##compute t-core of ptn
##	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[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:
save( tcore, 
freqtab
, 
ispos
, 
"c:\\Program Files\\Maple 13\\mylib\\tcore.m");
lprint("EXIT RESTART AND TEST tcore PACKAGE");
#++++++++++++++++++++++++++++++++++++++++++++++++++++