###  WPROGS.txt  used to check paper "Weighted Partition Identities and
###  Divisor Sums"
###  05.16.16


with(combinat):
with(numtheory):

DC:=proc(ptn)
   # partitions into distinct parts
   local sptn:
   sptn:=convert(ptn,set):
   if nops(ptn)=nops(sptn) then
   RETURN(true):
   else
   RETURN(false):
   fi:
end proc:


smp:=proc(ptn)
   # smallest part
   RETURN(ptn[1]):
end:
numsmp:=proc(ptn)
local c,s,j,nm:
   c:=0: s:=ptn[1]: nm:=nops(ptn):
   for j from 1 to nm do if ptn[j]=s then c:=c+1: fi: od:
   RETURN(c):
end:

# ffwstat:=proc(n)
#  # Fokkink-Fokkink-Wang partition stat
#  local ptns,dptns,ptn,x:
#   ptns:=partition(n):
#   dptns:=select(DC,ptns):
#   x:=add( (-1)^(nops(ptn)+1)*smp(ptn),ptn in dptns):
#   RETURN(x):
# end:

# 12.23.16
largp:=ptn->ptn[nops(ptn)]:
# zffwstat:=proc(n)
#  # Fokkink-Fokkink-Wang partition stat
#  local ptns,dptns,ptn,x,LP:
#   ptns:=partition(n):
#   dptns:=select(DC,ptns):
#   x:=add( (-1)^(nops(ptn)+1)*z^(largp(ptn)-smp(ptn))*add(z^j,j=1..smp(ptn)),ptn in dptns):
#   ##x:=add( (-1)^(nops(ptn)+1)*smp(ptn),ptn in dptns):
#   RETURN(x):
# end:
ffwstatz1:=proc(n)
 # Fokkink-Fokkink-Wang partition stat
 local ptns,dptns,ptn,x,LP:
  ptns:=partition(n):
  dptns:=select(DC,ptns):
  x:=map(ptn->(-1)^(nops(ptn))*smp(ptn)*(smp(ptn)-2*largp(ptn)-1)/2, dptns):
  x2:=convert(x,`+`):
  ##x:=add( (-1)^(nops(ptn)+1)*smp(ptn),ptn in dptns):
  RETURN(x2):
end:
# ffwstat2:=proc(n)
#  # Fokkink-Fokkink-Wang partition stat
#  local ptns,dptns,ptn,x:
#   ptns:=partition(n):
#   dptns:=select(DC,ptns):
#   x:=add( (-1)^(smp(ptn)+1)*numsmp(ptn),ptn in dptns):
#   RETURN(x):
# end:
fftable:=proc(n)
 # Fokkink-Fokkink-Wang partition stat
 local ptns,dptns,ptn,x:
  ptns:=partition(n):
  dptns:=select(DC,ptns):
  for ptn in dptns do
  lprint(ptn,(-1)^(nops(ptn)+1)*smp(ptn));
  od:
  ##x:=add( (-1)^(nops(ptn)+1)*smp(ptn),ptn in dptns):
  RETURN():
end:
ffw1:=ptn->(-1)^nops(ptn)*smp(ptn)*(smp(ptn)-2*largp(ptn)-1)/2:
fftable1:=proc(n)
 # Fokkink-Fokkink-Wang partition stat
 local ptns,dptns,ptn,x:
  ptns:=partition(n):
  dptns:=select(DC,ptns):
  for ptn in dptns do
  lprint(ptn,nops(ptn),smp(ptn),largp(ptn), ffw1(ptn));
  od:
  ##x:=add( (-1)^(nops(ptn)+1)*smp(ptn),ptn in dptns):
  RETURN():
end:

numdivs := n -> nops(divisors(n)):
 zdivs:=n->add(z^k, k in divisors(n));

###############################################################################
## omega1a and omega1b

nup:=ptn->nops(ptn):  # number of  parts
lap:=ptn->ptn[nup(ptn)]: # larggest part
oset:=N->{seq(2*k-1,k=1..trunc((N+1)/2))}:

PC1:=proc(ptn)
## Test ptn in P[o]
local j,ptn2,con2,oset1,ptns,ptns2;
ptns:=convert(ptn,set):
if modp(lap(ptn),2)=0 then
  ptns2:=ptns minus {lap(ptn)}:
else
  ptns2:=ptns:
fi: ##print("ptns2=",ptns2);
con2:=map(x->modp(x,2),ptns2):
if con2<>{1} then RETURN(false):
else
  oset1:=oset(lap(ptn)):
  ptns:=convert(ptn,set):
  if `subset`(oset1,ptns) then RETURN(true): else RETURN(false): fi: 
fi:  
end:

lop:=proc(ptn) ##largest odd partition
local x,p,j:
   x:=1:
   for j from 1 to nops(ptn) do
     p:=ptn[j]: if modp(p,2)=1 and p>x then x:=p:
   fi: od:
RETURN(x):
end:

###############################################################################
## omega2a and omega2b
eset:=N->{seq(2*k,k=1..trunc((N)/2))}:

PC2:=proc(ptn)
## Test ptn in P[e]
local j,ptn2,con2,eset1,ptns,ptns2;
ptns:=convert(ptn,set):
  if modp(lap(ptn),2)=1 then
    ptns2:=ptns minus {lap(ptn)}:
  else
    ptns2:=ptns:
  fi: ##print("ptns2=",ptns2);
  con2:=map(x->modp(x,2),ptns2):
  if con2<>{0} and con2<>{} then RETURN(false):
  else
    eset1:=eset(lap(ptn)):
    ptns:=convert(ptn,set):
    if `subset`(eset1,ptns) then RETURN(true): else RETURN(false): fi: 
  fi:
end:

lep:=proc(ptn) ##largest even part
local x,p,j:
x:=0:
for j from 1 to nops(ptn) do
  p:=ptn[j]: if modp(p,2)=0 and p>x then x:=p:
fi: od:
RETURN(x):
end:

###############################################################################
ptn2freqs0:=proc(ptn)
local n,l,freq,i,j,pt:
        n:=nops(ptn):
if ptn=[] then
   RETURN([]):
else
        l:=ptn[n]:
        freq:=array(1..l):
        for i from 1 to l do
           freq[i]:=0:
        od:
        for j from 1 to n do
           pt:=ptn[j]:
           freq[pt]:=freq[pt]+1:
        od:
        RETURN(convert(freq,list)):
fi:
end:

ptn2freqs:=proc(ptn)
  # convert to frequence notation
  local flist,fout,j:
  flist:=ptn2freqs0(ptn):
  fout:=[]:
  for j from 1 to nops(flist) do
     if flist[j]<>0 then
        fout:=[op(fout),cat(j)^flist[j]]:
     fi:
  od:
  RETURN(fout):
end: