////////////////////////////////////////////////////////////////////////////////
version="1.0";

category="Algebraic Geometry";
info="
LIBRARY:    schubert.lib   Proceduces for Intersection Theory

PROCEDURES:
    projectiveSpace(n)     Chow ring of projective space of dimension n
    grassmannian(k,n)      Chow ring of Grassmannian of k-planes on an n-space
    projectiveBundle(f,d)  Chow ring of projective bundles
    grassmannBundle(k,f,d) Chow ring of Grassmann bundles
    integral(f)            compute intersection number
    dual(f)                dual of a vector bundle
    tensorBundle(f,g,d)    tensor of two vector bundles
    symmPower(n,f,d)       n-th symmetric power of a vector bundle
    chern(f,d)             total Chern class of a vector bundle
    segre(f,d)             total Segre class of a vector bundle
    rankBundle(f)          rank of a vector bundle

KEYWORDS:        Intersection Theory, Enumerative Geometry, Schubert Calculus
";

////////////////////////////////////////////////////////////////////////////////

LIB "general.lib";

////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
/// Procedures concerned with abstract varieties ///////////////////////////////
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////

proc projectiveSpace(int n)
"USAGE: projectiveSpace(n); n int
RETURN:     ring
PURPOSE:    create a polynomial ring in one variable h, where h is the hyperplane
            class of this projective space. The monomial ordering of this ring is
            'wp' (weighted reverse lexicographical ordering) with vector (1). 
            Moreover, we export an ideal (say 'rels') generated by h^{n+1} and 
            the quotient ring modulo this ideal should be its Chow ring.
EXAMPLE: example projectiveSpace; shows examples
"
{
    ring CR = 0, (h), wp(1);
    int d = n;
    poly u = 1;
    poly v = 1 + h;
    list l;
    int i;
    for (i=1;i<=n;i++)
    {
        l=insert(l,1,size(l));
        u=u+(-1)^i*schur(l,v);
    }
    ideal rels = h^(n+1);
    rels = std(rels);
    poly s = reduce(logg(u,d)+n,rels);
    poly q = reduce(logg(v,d)+1,rels);
    export (s,q,rels);
    return (CR);
}
example
{
    "EXAMPLE:"; echo=2;
    def P3 = projectiveSpace(3);
    setring P3;
    s; q; rels;
}

proc grassmannian(int k, int n)
"USAGE: grassmannian(k,n); k int, n int
RETURN:     ring
PURPOSE:    create a polynomial ring in n variables, in which the first k 
            variables s(1),...,s(k) will be the Chern classes of the 
            tautological subbundle on Grassmannian G(k,n) and the next n-k 
            variables q(1),...,q(n-k) will be the Chern classes of the 
            tautological quotient bundle on G(k,n). The monomial ordering 
            of this ring is 'wp' (weighted reverse lexicographical ordering) 
            with vector (1..k,1..n-k).
            Moreover, we export the Chern characters of the tautological 
            subbundle and quotient bundle on G(k,n) (say 's' and 'q'). They 
            are the polynomials. One more thing, we export an ideal 
            (say 'rels') generated by the relations between the Chern 
            variables and the quotient ring modulo this ideal should be Chow 
            ring of G(k,n).
EXAMPLE: example grassmannian; shows examples
"
{
    ring CR = 0, (q(1..n-k)), wp(1..n-k);
    int d = k*(n-k);
    int i,j,h;
    poly v = 1;
    poly u = 1;
    for (j=1;j<=n-k;j++) {v=v+q(j);}
    list l;
    for (i=1;i<=k;i++)
    {
        l=insert(l,1,size(l));
        u=u+(-1)^i*schur(l,v);
    }
    ideal rels;
    for (h=k+1;h<=n;h++)
    {
        l=insert(l,1,size(l));
        rels=rels,schur(l,v);
    }
    rels = std(rels);
    poly s = reduce(logg(u,d)+k,rels);
    poly q = reduce(logg(v,d)+n-k,rels);
    export (s,q,rels);
    return (CR);
}
example
{
    "EXAMPLE:"; echo=2;
    def G24 = grassmannian(2,4);
    setring G24;
    s; q; rels;
}

////////////////////////////////////////////////////////////////////////////////

proc projectiveBundle(poly f, int d)
"USAGE: projectiveBundle(f,d); f poly, d int
RETURN:     ring
PURPOSE:    create a polynomial ring which we work on.
            we export an ideal (say 'rels') generated by the relations between 
            the Chern variables and the quotient ring modulo this ideal should 
            be its Chow ring.
EXAMPLE: example projectiveBundle; shows examples
"
{
    def R = basering;
    setring R;
    int r = rankBundle(f);
    poly c = chern(f,d);
    ring P = 0, (z), wp(1);
    def CR = P + R;
    setring CR;
    int dd = r-1 + d;
    poly c = imap(R,c);
    ideal rels = imap(R,rels);
    poly g = z^r;
    int i;
    for (i=1;i<=r;i++) {g=g+z^(r-i)*part(c,i);}
    rels=std(rels,g);
    poly u = 1 + z;
    poly Q = reduce(logg(u,dd)+1,rels);
    export (rels,Q);
    return (CR);
}
example
{
    "EXAMPLE:"; echo=2;
    def G35 = grassmannian(3,5);
    setring G35;
    int d = 6;
    poly f = symmPower(2,dual(s),6);
    def PB = projectiveBundle(f,6);
    setring PB;
    int d = 14;
    Q; rels;
}
////////////////////////////////////////////////////////////////////////////////

proc grassmannBundle(int k, poly f, int d)
"USAGE: grassmannBundle(k,f,d); k int, f poly, d int
RETURN:     ring
PURPOSE:    create a polynomial ring which we work on.
            we export the Chern characters of the tautological subbundle and 
            quotient bundle on this Grassmann bundle (say 'S' and 'Q'). 
            They are the polynomials. One more thing, we also export an ideal 
            (say 'rels') generated by the relations between the Chern variables 
            and the quotient ring modulo this ideal should be its Chow ring.
EXAMPLE: example grassmannBundle; shows examples
"
{
    def R = basering;
    setring R;
    int r = rankBundle(f);
    ring G = 0, (S(1..k),Q(1..r-k)), wp(1..k,1..r-k);
    def CR = G + R;
    setring CR;
    int dd = k*(r-k)+d;
    poly f = imap(R,f);
    ideal rels = imap(R,rels);
    poly u = 1;
    poly v = 1;
    int i,j,h;
    for (i=1;i<=k;i++) {u=u+S(i);}  
    for (j=1;j<=r-k;j++) {v=v+Q(j);}
    poly g = u*v;
    ideal I = part(g,1)-part(expp(f,r),1);
    for (h=2;h<=r;h++) {I = I,part(g,h)-part(expp(f,r),h);}
    def J = I + rels;
    ideal rels = std(J);
    poly S = reduce(logg(u,dd)+k,rels);
    poly Q = reduce(logg(v,dd)+r-k,rels);
    export (rels,S,Q);
    return (CR);
}
example
{
    "EXAMPLE:"; echo=2;
    def G25 = grassmannian(2,5);
    setring G25;
    int d = 6;
    poly f = dual(symmPower(2,q,6));
    def GB = grassmannBundle(5,f,6);
    setring GB;
    int d = 14;
    S; Q; rels;
}

////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
/// Auxilary Static Procedures in this Library /////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
/// - part
/// - parts
/// - logg
/// - expp
/// - adams
/// - wedges
/// - integral
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////

static proc part(poly f, int n)
{
    int i;
    poly p;
    for (i=1;i<=size(f);i++)
    {
        if (deg(f[i])==n)
        {
            p=p+f[i];
        }
    }
    return (p);
}

////////////////////////////////////////////////////////////////////////////////

static proc parts(poly f, int i, int j)
{
    int k;
    poly p;
    for (k=i;k<=j;k++)
    {
        p=p+part(f,k);
    }
    return (p);
}

////////////////////////////////////////////////////////////////////////////////

static proc logg(poly f, int n)
{
    poly p;
    int i,j,k,m;
    if (n==0) {p=0;}
    if (n==1) {p=part(f,1);}
    else
    {
        list l=-part(f,1);
        for (j=2;j<=n;j++)
        {
            poly q;
            for (k=1;k<j;k++)
            {
                q=q+part(f,k)*l[j-k];
            }
            q=-j*part(f,j)-q;
            l=insert(l,q,size(l));
            kill q;
        }
        for (m=1;m<=n;m++)
        {
            p=p+1/factorial(m)*(-1)^m*l[m];
        }
    }
    return (p);
}

////////////////////////////////////////////////////////////////////////////////

static proc expp(poly f, int n)
{
    poly p;
    int i,j,k;
    if (deg(f)==0) {p=1;}
    else
    {
        list l=1;
        for (i=1;i<=n;i++)
        {
            poly q;
            for (j=1;j<=i;j++)
            {
                q=q+factorial(j)*(-1)^(j-1)*l[i-j+1]*part(f,j)/i;
            }
            l=insert(l,q,size(l));
            kill q;
        }
        for (k=1;k<=size(l);k++)
        {
            p=p+l[k];
        }
    }
    return (p);
}

////////////////////////////////////////////////////////////////////////////////

static proc adams(poly f, int n)
{
    poly p;
    int i;
    for (i=0;i<=deg(f);i++)
    {
        p=p+n^i*part(f,i);
    }
    return (p);
}

////////////////////////////////////////////////////////////////////////////////

static proc wedges(int n, poly f, int d)
{
    int i,j;
    list l;
    if (n==0) {l=1;}
    if (n==1) {l=1,f;}
    else
    {
        l=1,f;
        for (i=2;i<=n;i++)
        {
            poly q;
            for (j=1;j<=i;j++)
            {
                q=q+((-1)^(i-j))*parts(l[j]*adams(f,i-j+1),0,d)/i;
            }
            l=insert(l,q,size(l));
            kill q;
        }
    }
    return (l);
}

////////////////////////////////////////////////////////////////////////////////

static proc schur(list p, poly f)
{
    int i,j;
    int n = size(p);
    matrix M[n][n];
    for (i=1;i<=n;i++)
    {
        for (j=1;j<=n;j++)
        {
            M[i,j] = part(f,p[i]+j-i);
        }
    }
    return (det(M));
}

////////////////////////////////////////////////////////////////////////////////

proc integral(poly f)
"USAGE: integral(f); f poly
RETURN:     int
PURPOSE:    compute the intersection number. 
            That is the degree of top Chern class.
EXAMPLE: example integral; shows examples
"
{
    return (leadcoef(reduce(f,rels)));
}
example
{
    "EXAMPLE:"; echo=2;
    def G25 = grassmannian(2,5);
    setring G25;
    integral(q(2)^3);
}

////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
/// Procedures concerned with vector bundles ///////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////

proc symmPower(int n, poly f, int d)
"USAGE: symmPower(n,f,d); n int, f poly, d int
RETURN:     poly
INPUT:      n -- integer
            f -- Chern character of a given vector bundle
            d -- dimension of current variety.
OUTPUT:     a polynomial which is the Chern character of the n-th symmetric 
            power of the vector bundle defined by 'f'
EXAMPLE: example symmPower(n,f,d); shows examples
"
{
    int i,j,m;
    int r = rankBundle(f);
    poly result;
    list s,w;
    if (n==0) {result=1;}
    if (n==1) {result=f;}
    else
    {
        s = 1,f;
        w = wedges(n,f,d);
        for (i=2;i<=n;i++)
        {
            if (i<=r) {m=i;}
            else {m=r;}
            poly q;
            for (j=1;j<=m;j++)
            {
                q = q + ((-1)^(j+1))*parts(w[j+1]*s[i-j+1],0,d);
            }
            s = insert(s,q,size(s));
            kill q;
        }
        result = s[n+1];
    }
    return (result);
}
example
{
    "EXAMPLE:"; echo=2;
    def G24 = grassmannian(2,4);
    setring G24;
    int d = 4;
    symmPower(3,q,d);
}
////////////////////////////////////////////////////////////////////////////////

proc dual(poly f)
"USAGE: dual(f); f poly
RETURN:     poly
INPUT:      f -- Chern character of a vector bundle
OUTPUT:     a polynomial which is the Chern character of the dual of the vector
            bundle defined by 'f'.
EXAMPLE: example dual(f); shows examples
"
{
    return (adams(f,-1));
}
example
{
    "EXAMPLE:"; echo=2;
    def G24 = grassmannian();
    setring G24;
    dual(s);
}

////////////////////////////////////////////////////////////////////////////////

proc tensorBundle(poly f, poly g, int d)
"USAGE: tensorBundle(f,g,d); f poly, g poly, d int
RETURN:     poly
INPUT:      f, g -- Chern characters of two vector bundles
            d -- dimension of current variety. 
OUTPUT:     a polynomial which is the Chern character of the tensor of two 
            vector bundles defined by 'f' and 'g' respectively.
EXAMPLE: example tensorBundle(f,g,d); shows examples
"
{
    return (parts(f*g,0,d));
}
example
{
    "EXAMPLE:"; echo=2;
    def G34 = grassmannian(3,4);
    setring G34;
    int d = 3;
    tensorBundle(s,q,d);
}

////////////////////////////////////////////////////////////////////////////////

proc chern(poly f, int d)
"USAGE: chern(f,d); f poly, d int
RETURN:     poly
INPUT:      f -- Chern character of a given vector bundle
            d -- dimension of current variety
OUTPUT:     a polynomial which is the total Chern class of the vector bundle 
            defined by 'f'.
EXAMPLE: example chern(f,d); shows examples
"
{
    poly c = reduce(expp(f,d),rels);
    return (c);
}
example
{
    "EXAMPLE:"; echo=2;
    def G24 = grassmannian(2,4);
    setring G24;
    int d = 4;
    chern(q,d);
}

////////////////////////////////////////////////////////////////////////////////

proc segre(poly f, int d)
"USAGE: segre(f,d); f poly, d int
RETURN:     poly
INPUT:      f -- Chern character of a given vector bundle
            d -- dimension of current variety
OUTPUT:     a polynomial which is the total Segre class of the vector bundle 
            defined by 'f'.
EXAMPLE: example segre(f,d); shows examples
"
{
    poly g = dual(chern(f,d));
    if (g==1) {return (1);}
    else
    {
        poly t,h;
        int i,j;
        for (i=0;i<=d;i++) {t = t + (1-g)^i;}
        for (j=0;j<=d;j++) {h = h + part(t,j);}
        return (reduce(h,rels));
    }
}
example
{
    "EXAMPLE:"; echo=2;
    def G24 = grassmannian(2,4);
    setring G24;
    int d = 4;
    segre(q,d);
}

////////////////////////////////////////////////////////////////////////////////

proc rankBundle(poly f)
"USAGE: rankBundle(f); f poly
RETURN:     poly
INPUT:      f -- Chern character of a given vector bundle
OUTPUT:     an integer which is the rank of the vector bundle defined by 'f'.
EXAMPLE: example rankBundle(f); shows examples
"
{
    return (int(part(f,0)));
}
example
{
    "EXAMPLE:"; echo=2;
    def G24 = grassmannian(2,4);
    setring G24;
    int d = 4;
    rankBundle(q);
}