# Written by David Haws.

with(combinat);

###  testrec := proc(k)
###      local t;
###      t := k/2;
###      printf("t=%f\n",t);
###      if (t > 2) then
###          testrec(t);   
###      end if;
###  end proc;

### Recursive form of Katzman A using relation (1)
### Might be usefull for a better katzman precomputing routine.
katzmanrec := proc(n,d,i)
    local k, A;
###    printf ("n=%d, d=%d, i=%d\n",n,d,i);
    if (i < 0 or i > n*(d-1)) then
###       printf ("   Case i < 0 or i > n*(d-1)\n");
        return 0;
    end if;
    if (d = 1) then
###       printf ("   Case d=1\n");
        return 1;
    end if;
    if (n = 1) then
###       printf ("   Case n=1\n");
        return 1;
    end if;
    A := 0;
    for k from 1 to d do
        A := A + katzmanrec(n-1,d,i-d+k);
    end do;
    return A;
end proc;

### This is nearly the same as above, but it relies on precomputed data 
### if available.
katzmanrecP := proc(n,d,i)
    local k, A;
    global katzmanPrecompute, katzmanPrecomputeDim, katzmanCoeff;
    ### BEGIN BASE CASES ###
    if (i < 0 or i > n*(d-1)) then
        return 0;
    end if;
    if (d = 1) then
        return 1;
    end if;
    if (n = 1) then
        ### printf ("   Case n=1\n");
        return 1;
    end if;
    ### END   BASE CASES ###
   
    ### Now lets try to use precomputed data
    ### We precompute up to n+1 rank
    if (katzmanPrecompute = 1 and n+1 >= d) then
    ### OLD if (katzmanPrecompute = 1) then
        if (katzmanPrecomputeDim >= n-1) then
            ### printf ("Using precomputed data\n");
            A := 0;
            for k from 1 to d do
                    ### printf ("n=%d, d=%d, i-d+k=%d   \n",n,d,i-d+k);
                if (i-d+k >= 0 and i-d+k <= (n-1)*(d-1)) then
                    ### printf("    Value=%d\n",(katzmanCoeff[n-1,d])[i-d+k+1]);
                    A := A + (katzmanCoeff[n-1,d])[i-d+k+1];
                    ### printf (" ***\n");
                end if;
            end do;
            return A;
        end if;
    end if;
    ### printf ("Calculating Raw   n=%d, d=%d, \n",n,d); 
    A := katzman(n,d,i);
    return A;
    A := 0;
    for k from 1 to d do
        A := A + katzmanrecP(n-1,d,i-d+k);
    end do;
    return A;
end proc;

# This computes the Katzman A's. i.e. katzman(n,d,i) = A_{i}^{n,d}
katzman := proc(n,d,i)
#printf ("n=%d, d=%d, i=%d\n",n,d,i);
#if (i < 0) then
#    printf("katzman: Called with neg\n");
#end if;
#printf("i=%d\n",i);
if (katzmanPrecompute = 1) then
    if (katzmanPrecomputeDim >= n) then
        #printf ("Using precomputed data\n");
        if (i > n*(d-1)) then
            return 0;
        end if;
        return (katzmanCoeff[n,d])[i+1];
    else
        return coeff(expand(sum('t^k', 'k'=0..(d-1))^n),t,i)
    end if;
end if;
coeff(expand(sum('t^k', 'k'=0..(d-1))^n),t,i)
end proc;

# This computes the Katzman A's. i.e. katzman(n,d,i) = A_{i}^{n,d}
katzmanold := proc(n,d,i)
#if (i < 0) then
#    printf("katzman: Called with neg\n");
#end if;
#printf("i=%d\n",i);
coeff(expand(sum('t^k', 'k'=0..(d-1))^n),t,i)
end proc;

katzmanvecold := proc(n,d)
local i;
for i from 0 to n*(d-1) do
    printf ("%d, ",katzman(n,d,i));
end do;
printf("\n");
end proc;

katzmanvec := proc(n,d)
local i, TempArray;
TempArray := array(1..(n*(d-1)+1)):

for i from 0 to n*(d-1) do
    TempArray[i+1] := katzmanold(n,d,i):
end do;
return TempArray;
end proc;

### This uses the recursive and precompute katzman to find the vector
katzmanvecrec := proc(n,d)
    local i, TempArray;
    TempArray := array(1..(n*(d-1)+1)):

    for i from 0 to n*(d-1) do
        TempArray[i+1] := katzmanrecP(n,d,i):
    end do;
    return TempArray;
end proc;

calckatzman := proc(k)
    local A, i, l;
    A := array(1..k,1..k);
    for i from 1 to k do
        printf ("Calculating i=%d\n",i);
        for l from 1 to i do
            A[i,l] := katzmanvec(i,l);
        end do;
        return A;
    end do;
end proc;

### This assumes that katzmanCoeff has been allocated as array(1..k,1..k);
calckatzmanrec := proc(k)
    local A, i, l;
    global katzmanPrecompute, katzmanPrecomputeDim, katzmanCoeff;
    katzmanPrecompute := 1;
    katzmanPrecomputeDim := 0;
    ### i = dimension
    for i from 1 to k do
        printf ("Calculating i=%d\n",i);
        ### l is up to k
        for l from 1 to i+2 do
            ###  if (l mod k/10 = 0) then
            ###      printf (" %d\n",l);
            ###  end if;
            ### printf("    Calculating i=%d <---> l=%d\n",i,l);
            ### printf ("i=%d, l=%d\n",i,l);
            ### print (katzmanvecrec(i,l));
            if (l <= i+2) then
                katzmanCoeff[i,l] := katzmanvecrec(i,l);
            else
                katzmanCoeff[i,l] := katzmanvec(i,l);
            end if;
        end do;
        ### printf ("\n");
        katzmanPrecomputeDim := i;
    end do;
end proc;

# This computes the pth coefficient of the h-vector(d,n)
hcoeff := proc(d,n,p)
local tempcoeff,s,j,k;
tempcoeff := 0:

for s from 0 to d-1 do
    for j from 0 to s do
        for k from 0 to j do
#            printf("s=%d, j=%d, k=%d\n",s,j,k);
#            printf("    n-j=%d, d-s=%d, (p-k)(d-s)=%d   (+-1)=%d\n",n-j,d-s,(p-k)*(d-s),(-1)^(s+j+k));
            tempcoeff := tempcoeff + (-1)^(s+j+k)*numbcomb(n,s)*numbcomb(s,j)*numbcomb(j,k)*katzman(n-j,d-s,(p-k)*(d-s)):
        end do;
    end do;
end do;
tempcoeff;
end proc;

# This computes the h-vector using the expression exactly as it appears
# in katzman
hsumunmod := proc (d,n)
local s,j,l,stemp,jtemp,ltemp,kcount, hvecunmod;
hvecunmod := 0:
stemp := 0:
jtemp := 0:
ltemp := 0:
kcount := 0:
for s from 0 to d-1 do
    for j from 0 to s do
        for l from 0 to n-j do
            ltemp := ltemp + katzman(n-j,d-s,l*(d-s))*t^l:
            kcount := kcount + 1:
        end do;
        jtemp := jtemp + (-1)^j*numbcomb(s,j)*(1-t)^j*ltemp:
        ltemp := 0;
    end do;
    stemp := stemp + (-1)^s*numbcomb(n,s)*jtemp:
    jtemp := 0:
end do;
    #printf ("katzman called %d times\n",kcount);
    stemp;
end proc;

# This computes the h-vector polynomial from the rewritten epxression
# that appeared in katzman paper
hsum := proc (d,n)
local hvec, kcount,s,j,k,l;
hvec := 0:
kcount := 0:
for s from 0 to d-1 do
    for j from 0 to s do
        for k from 0 to j do
            for l from 0 to n-j do
            hvec := hvec + ((-1)^(s+j+k)*numbcomb(n,s)*numbcomb(s,j)*numbcomb(j,k))*katzman(n-j,d-s,l*(d-s))*t^(l+k):
            kcount := kcount + 1:
#            printf("hsum: pow=%d\n",(l+k));
            end do;
        end do;
    end do;
end do;
#printf ("katzman called %d times\n",kcount);
hvec;
end proc;

sympoly := proc (p)
local deg, halfpt, i;
deg := degree(p):

halfpt := floor((deg+1)/2):
for i from 0 to halfpt do
    if (coeff(p,t,i) <> coeff(p,t,deg-i)) then
        printf ("Not symmetric at i=%d\n",i);
        i := halfpt:
    end if;
end do;


end proc;

# This procedure checks if the polynomial p is unimodal
unimodal := proc (p)
local inc, oldcoeff, i;
# Test to see if the polynomial p is unimodal
# Returns 1 if unimodal, 0 else

inc := 1:   # Signifies that it is still increasing
oldcoeff := 0:
for i from 0 to degree(p) do
    if (abs(coeff(p,t,i)) <= abs(oldcoeff)) then
        # We are decreasing
        # Not increasing any more
        inc := 0;
        # This must be the peak
    else
        # we are increasing then, problem if inc == 0
        if (inc = 0) then
            return (0);
        end if;
    end if;
    oldcoeff := abs(coeff(p,t,i));    
end do;
return (1);
end proc;

# This will compute hsum(d,n) with d=0 to a, n=0 to b
htest := proc (a,b)
local n,dlim,d;
for n from 1 to b do
    dlim := min(a,n-1);
    for d from 1 to dlim do
        printf("d=%d, n=%d\n",d,n);
        print(hsum(d,n));
    end do;
    printf("<****************> n=%d <****************>\n",n);
end do;
end proc;

# This will compute hsum(d,n) with d=0 to a, n=0 to b
htestunmod := proc (a,b)
local n,dlim,d;
for n from 1 to b do
    dlim := min(a,n-1);
    for d from 1 to dlim do
        printf("d=%d, n=%d\n",d,n);
        print(simplify(hsumunmod(d,n)));
    end do;
    printf("<****************> n=%d <****************>\n",n);
end do;
end proc;

hunitest := proc (a,b)
local uplim,d,n;
# This will test for unimodality from d=0 to a and n=0 to b
# Note that we need d <= n

for n from 1 to b do
    uplim := min(a,n-1);
    printf("    d from 1 to %d\n",uplim);
    for d from 1 to uplim do
        if (unimodal(hsum(d,n)) = 0) then
            printf("Found non-unimodal");
            print(n,d);
        end if;
    end do;
    printf("<****************> n=%d <****************>\n",n);
end do;
end proc;

hunitestunmod := proc (a,b)
local uplim,d,n;
# This will test for unimodality from d=0 to a and n=0 to b
# Note that we need d <= n

for n from 1 to b do
    uplim := min(a,n-1);
    printf("    d from 1 to %d\n",uplim);
    for d from 1 to uplim do
        if (unimodal(simplify(hsumunmod(d,n))) = 0) then
            printf("Found non-unimodal");
            print(n,d);
        end if;
    end do;
    printf("<****************> n=%d <****************>\n",n);
end do;
end proc;

hrandunitest := proc (b)
local n,d;
# This will test for unimodality from d=0 to a and n=0 to b
# Note that we need d <= n

for n from 1 to b do
    d := (rand() mod n-1) + 1;
    printf ("   d=%d\n",d);
    if (unimodal(hsum(d,n)) = 0) then
        printf("Found non-unimodal");
        print(n,d);
    end if;
    printf("<****************> n=%d <****************>\n",n);
end do;
end proc;

clist := proc(d,n)
local s,j,k,p;
#sum := 0:
#vec :=
for s from 0 to d-1 do
    for j from 0 to s do
        for k from 0 to j do
            printf ("s=%d, j=%d, k=%d, (-1)^(s+j+k)=%d coef=%d\n",s,j,k,(-1)^(s+j+k),numbcomb(n,s)*numbcomb(s,j)*numbcomb(j,k));
            printf ("                                          (");
            for p from 0 to n do
                printf ("%d, ",katzman(n-j,d-s,(p-k)*(d-s)));
            end do;
            printf (")\n\n");
#            sum := sum + katzman(n-j,d-s,(p-k)*(d-s))*(-1)^(s+j+k):
        end do;
    end do;
end do;
#printf ("Sum=%d\n",sum);
end proc;

termlist := proc(d)
local s,j,k;
#sum := 0:
#vec :=
for s from 0 to d-1 do
    for j from 0 to s do
        for k from 0 to j do
            printf ("s=%d, j=%d, k=%d, (-1)^(s+j+k)=%d (n %d) (%d %d) (%d %d)\n",s,j,k,(-1)^(s+j+k),s,s,j,j,k);
            printf ("                                          (");
            printf ("A^{n-%d,%d}_{(p-%d)(%d)}, ",j,d-s,k,d-s);
            printf (")\n\n");
#            sum := sum + katzman(n-j,d-s,(p-k)*(d-s))*(-1)^(s+j+k):
        end do;
    end do;
end do;
#printf ("Sum=%d\n",sum);
end proc;

d3 := proc(n)
local p,d;
for p from 0 to n do
    if (p <> 2) then
        for d from 0 to 2 do
            printf("%20d - %d*%20d = %20d\n",katzman(n,3,3*p-2+d),n,numbcomb(n,2*p-1-2+d), katzman(n,3,3*p-2+d) - n*numbcomb(n,2*p-1-2+d));
        end do;
    else
        for d from 0 to 2 do
            printf("%20d - %d*%20d + %d = %20d\n",katzman(n,3,3*p-2+d),n,numbcomb(n,2*p-1-2+d),numbcomb(n,2), katzman(n,3,3*p-2+d) - n*numbcomb(n,2*p-1-2+d) + numbcomb(n,2));
        end do;
    end if;
end do;

end proc;


katsum := proc(n,d)
local k,c,p,l;
for k from 0 to n do
    c := numbcomb(n,k):
    printf("%2d * ",c);
    for p from 1 to k do
        printf("   0, "):
    end do;
    for l from 0 to (k)*(d-2) do
        printf("%4d, ", katzman(k,d-1,l)):
    end do;
    printf ("\n");
end do;
end proc;

katsum2 := proc(n,d)
local k,c,p,l;
for k from 0 to n do
    c := numbcomb(n,k):
    for p from 1 to k do
        printf("  0,"):
    end do;
    for l from 0 to (k)*(d-2) do
        printf("%3d,", c*katzman(k,d-1,l)):
    end do;
    printf ("\n");
end do;


end proc;

termcoeff := proc(l,n)
    local q, sum1;
    sum1 := 0:
    for q from 0 to l do
#        sum1 := sum1 + numbcomb(n,q)*katzman(q,2,l-q);
        sum1 := sum1 + numbcomb(n,q)*numbcomb(q,l-q);;
    end do;
    sum1:
end proc;

termcoeff2 := proc(l,n)
    local q, sum1;
    sum1 := 0:
    for q from 0 to l do
        printf("%d * %d + ",numbcomb(n,q),numbcomb(q,l-q));
    end do;
    printf("\n");
end proc;