#####################################################################
# This Maple V script computes the Hodge numbers of a smooth complete 
# intersection and Betti numbers of the moduli space of vector bundles
# It's based on a formulas  given inof Hirzebruch given in SGA7 exp
# XI, and
# [Z] D. Zagier, "Elementary aspects of the Verlinde formula
# and of the Harder-Narasimhan-Atiyah-Bott  formula." Proc.
# of Hirzebruch's 65 birthday conference.
# See the file "formulas.pdf" for more background.
#
# To compute the middle Hodge numbers (i.e h^pq, p+q = dim X) of
# a smooth c.i. of degree (d_1,..d_k) in P^n, type:
#     hodge([d_1,... d_k],n);
# to obtain:
#     h^dimX,0, h^(dimX-1),1,...
#
# (Recall that by weak Lefschetz, the nonmiddle hodge numbers 
# can be computed trivially as  h^pq = delta_pq)
#
# To calculate k terms of the Poincare polynomial for the  moduli space M
# of rank n, and fixed determinant degree d vector bundles over a
# genus g curve when (n,d)=1,
# use poincare_vb_trunc(n,d,g,k). Use  poincare_vb(n,d,g) to calculate
# the full polynomial which has degree k = 2dim M = 2(n^2-1)(g-1)
#
# Donu Arapura dvb@math.purdue.edu
#####################################################################



 # Auxillary functions for hodge_ci
 # In the interests of efficiency, the algorithm  discards terms of high degree;
 # consequently the code is a bit convoluted.
 # truncat(p,level) yields p if it's a monomial otherwise it truncates terms
 # with deg > level. 


        truncat := proc(p,level) local total, m, seq,q;
               total := 0;
               q := expand(p);
               if nops([coeffs(q)]) = 1 then
                  total := q
                else
                  seq := op(q);
                  for m in seq do
                     if( degree(m,x)+degree(m,y) <= level) then
                       total := total+m
                     fi
                   od
                fi;
                total;
                end;

 # hinv(p,n) computes the truncation of the sum of the
 # 1st n-terms of the geom. series 1/(1+p)
      hinv := proc(p,n) local i,total,prod;
       total :=1;
       prod := 1;
       for i to n do 
        prod := truncat( -p*prod,n);
        total := total + prod;
       od;
       total;
       end;

  # hnum computes numerator of a rational function needed below.
       hnum := proc(d,level) 
         local total,i,j, top; 
          total := 0;
          top := max(d,level); 
          for i from 0 to top-1 do for j from 0 to top-2-i do 
              total := total+binomial(d-1,i+j+1)*x**i*y**j od od; 
          total 
        end;

   # 1 + hdenminus1 computes denominator.
       hdenminus1 := proc(d,level)
        local total,i,j, top; 
         total := 0;
          top := max(d, level); 
         for i to top do for j to top-i do 
             total := total-binomial(d,i+j)*x**i*y**j od od; 
         total 
        end;

   # hquot = hnum/(1+hdenminus1) expanded as a truncated series (upto level)
      hquot := proc(a,level) hnum(a,level)*hinv(hdenminus1(a,level),level) end;

   # hpoly gives (truncated) generating function for primitive
   # hodge numbers [SGA7 part II,p54]
       hpoly := proc(degs,level)
        local p,i;
          p := 1;
          for i in degs do 
             p := truncat(p*(1+(1+x)*(1+y)*hquot(i,level)), level)
           od;
          expand((p-1)*hinv((x+y +x*y), level));
        end;

   
   # hodge1 computes sequence h^pq with p+q = level
       hodge1 := proc(degs,level) 
         local htemp, hseq,i; 
          hseq := NULL;
          htemp := hpoly(degs,level);
          for i from 0 to level do
            if (2*i = level) then
              hseq := hseq, (coeff(coeff(htemp,x,i),y,level-i)+1);
            else
               hseq := hseq, (coeff(coeff(htemp,x,i),y,level-i)); 

             fi;
           od;
         end; 

   
# Main procedure:
      hodge_ci := proc(degs, n); hodge1(degs, n - nops(degs)) end;

      fast_hodge_ci = hodge_ci;

###########################################
# This batch of functions will be used to
# construction a list of ordered paritions of an
# integer.

successor :=
proc(lis) local i,s, len, newlis;
s := {};
len := nops(lis);
for i from 1 to len do
newlis := lis;
newlis[i] := newlis[i]+1;
s := s union {newlis};
od;
s;
end;

Successor :=
proc(st) local i,s, len;
len := nops(st);
s := {};
for i from 1 to len do
s := successor(st[i]) union s;
od;
s;
end;

opart1 :=
proc(k,n) local i,s, lis;
lis := [seq(1,i=1..k)];
s := {lis};
for i from 1 to n-k do
s := Successor(s);
od;
s;
end;

opartitions :=
proc(n) local s,i;
s := {};
for i from 1 to n do
s := s union opart1(i,n);
od;
convert(s,list);
end;

####################
# Below I implemented Zagier's formulas. I've
# tried to follow his notation as far as posible
# There is a small typo in [Z, p 458], the definition given
# for M_{g-1} should be for M_g

P :=
proc(n,g) local i, prd;
prd := 1;
for i from 1 to n do
 prd := prd*(1+t^(2*i-1))^(2*g)/(1-t^(2*i))^2;
od;
prd*(1-t^(2*n));
end;

# bracket(x) := <x> in [Z]
bracket := x -> 1 + floor(x) -x;

# Sum elements of a list
sumlist :=
proc(lis) local i, len;
len := nops(lis);
sum(lis[i], i =1..len);
end;

M :=
proc(lis,lambda) local i, len, sm,Sm;
len := nops(lis);
sm := 0;
Sm :=0;
for i from 1 to len-1 do
 sm := lis[i]+sm;
 Sm := (lis[i]+lis[i+1])*bracket(sm*lambda)+Sm;
od;
Sm;
end;

Mg :=
proc(lis,lambda,g) local i,j, sm, len;
len := nops(lis);
sm := 0;
for i from 1 to len do
for j from i+1 to len do
 sm := lis[i]*lis[j] + sm;
od; od;
M(lis, lambda) + (g-1)*sm;
end;

# Summand in [Z, Cor. p458]
# lis = [n_1,...,n_k]

Summand :=
proc(lis,n, d,g) local prd,len,i;
len := nops(lis);
prd := -(-1)^len*t^(2*Mg(lis, d/n,g));
for i from 1 to len-1 do
 prd := prd/(1-t^(2*lis[i]+2*lis[i+1]));
od;
for i from 1 to len do
 prd := prd*P(lis[i], g);
od;
prd;
end;

# First approx to [Z, Cor. p458]
# lis = list of partitions of n

Poincare1 :=
proc(lis, n,d,g) local i,sm, len, Lis;
len := nops(lis);
Lis := map(L-> Summand(L,n,d,g), lis);
sm := sum(Lis[i], i=1..len);
sm/P(1,g);
end;


poincare_vb_trunc :=
proc(n,d,g, numterms);
convert(series(Poincare1(opartitions(n), n,d,g),t, numterms+1), polynom);
end;

poincare_vb :=
(n,d,g) -> poincare_vb_trunc(n,d,g, 2*(n^2-1)*(g-1));