//////////////////////////////////////////////////////////////////////////
version="version schubert.lib 4.0.0.0 Jun_2013 "; // $Id$

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

AUTHOR:     Hiep Dang,  email: hiep@mathematik.uni-kl.de

OVERVIEW:

    We implement two new classes (variety and sheaf) and methods for computing
    with them. Here a variety is represented by the dimension and the Chow ring.
    A sheaf on a variety is represented by the Chern character.

    In particular, we implement the concrete varieties such as projective spaces
    , Grassmannians, and projective bundles. Finally, the most important thing
    is a method for computing the intersection numbers (degrees of 0-cycles).

PROCEDURES:
    makeVariety(int,ideal)            create a variety
    productVariety(variety,variety)   product of two varieties
    ChowRing(variety)                 create the Chow ring of a variety
    dimension(variety)                dimension of a variety
    relations(variety)                relations of a variety
    Grassmannian(int,int)             create a Grassmannian as a variety
    projectiveSpace(int)              create a projective space as a variety
    projectiveBundle(sheaf)           create a projective bundle as a variety
    integral(variety,poly)            degree of a 0-cycle on a variety
    makeSheaf(variety,poly)           create a sheaf
    ChernCharacter(sheaf)             the Chern character of a sheaf
    totalChernClass(sheaf)            the total Chern class of a sheaf
    ChernClass(sheaf,int)             the k-th Chern class of a sheaf
    topChernClass(sheaf)              the top Chern class of a sheaf
    totalSegreClass(sheaf)            the total Segre class of a sheaf
    dualSheaf(sheaf)                  the dual of a sheaf
    tensorSheaf(sheaf,sheaf)          the tensor of two sheaves
    symmetricPowerSheaf(sheaf,int)    the k-th symmetric power of a sheaf
    minusSheaf(sheaf,sheaf)           the quotient of two sheaves
    plusSheaf(sheaf,sheaf)            the direct sum of two sheaves

REFERENCES:

    Hiep Dang,  Intersection theory with applications to the computation of
                Gromov-Witten invariants, Ph.D thesis, TU Kaiserslautern, 2013.

    Sheldon Katz and Stein A. Stromme, Schubert-A Maple package for intersection
                theory and enumerative geometry, 1992.

    Daniel R. Grayson, Michael E. Stillman, Stein A. Stromme, David Eisenbud and
    Charley Crissman, Schubert2-A Macaulay2 package for computation in
                intersection theory, 2010.

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

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

LIB "general.lib";
LIB "homolog.lib";

////////////////////////////////////////////////////////////////////////////////
/// create new classes: varieties and sheaves //////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

proc mod_init()
{
    newstruct("variety","int dimension, ring baseRing, ideal relations");
    newstruct("sheaf","variety currentVariety, poly ChernCharacter");

    system("install","variety","print",variety_print,1);
    system("install","variety","*",productVariety,2);
    system("install","sheaf","print",sheaf_print,1);
    system("install","sheaf","*",tensorSheaf,2);
    system("install","sheaf","+",plusSheaf,2);
    system("install","sheaf","-",minusSheaf,2);
}

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

proc part(poly f, int n)
"USAGE:     part(f,n); f poly, n int
RETURN:     poly
PURPOSE:    return the homogeneous component of degree n of the polynomial f.
EXAMPLE:    example part; shows examples
"
{
    int i;
    poly p;
    for (i=1;i<=size(f);i++)
    {
        if (deg(f[i])==n) {p=p+f[i];}
    }
    return (p);
}
example
{
    "EXAMPLE:"; echo=2;
    ring r = 0,(x,y,z),wp(1,2,3);
    poly f = 1+x+x2+x3+x4+y+y2+y3+z+z2+xy+xz+yz+xyz;
    part(f,0);
    part(f,1);
    part(f,2);
    part(f,3);
    part(f,4);
    part(f,5);
    part(f,6);
}

proc parts(poly f, int i, int j)
"USAGE:     parts(f,i,j); f poly, i int, j int
RETURN:     poly
PURPOSE:    return a polynomial which is the sum of the homogeneous components
            of degree from i to j of the polynomial f.
EXAMPLE:    example parts; shows examples
"
{
    int k;
    poly p;
    for (k=i;k<=j;k++)
    {
        p=p+part(f,k);
    }
    return (p);
}
example
{
    "EXAMPLE:"; echo=2;
    ring r = 0,(x,y,z),wp(1,2,3);
    poly f = 1+x+x2+x3+x4+y+y2+y3+z+z2+xy+xz+yz+xyz;
    parts(f,2,4);
}
proc logg(poly f, int n)
"USAGE:     logg(f,n); f poly, n int
RETURN:     poly
PURPOSE:    computing the Chern character from the total Chern class.
EXAMPLE:    example logg; shows examples
"
{
    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);
}
example
{
    "EXAMPLE:"; echo=2;
    ring r = 0,(x,y),wp(1,2);
    poly f = 1+x+y;
    logg(f,4);
}

proc expp(poly f, int n)
"USAGE:     expp(f,n); f poly, n int
RETURN:     poly
PURPOSE:    computing the total Chern class from the Chern character.
EXAMPLE:    example expp; shows examples
"
{
    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);
}
example
{
    "EXAMPLE:"; echo=2;
    ring r = 0,(x),dp;
    poly f = 3+x;
    expp(f,3);
}

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));
}

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

proc variety_print(variety V)
{
    "A variety of dimension", V.dimension;
}

proc makeVariety(int d, ideal i)
{
    def R = basering;
    variety V;
    V.dimension = d;
    V.baseRing = R;
    V.relations = i;
    return(V);
}

proc dimension(variety V)
{
    return (V.dimension);
}

proc ChowRing(variety V)
{
    def R = V.baseRing;
    setring R;
    ideal rels = V.relations;
    qring CR = std(rels);
    return (CR);
}

proc relations(variety V)
{
    def R = V.baseRing;
    setring R;
    ideal i = V.relations;
    return (i);
}

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

proc Grassmannian(int k, int n, list #)
"USAGE:     Grassmannian(k,n); k int, n int
RETURN:     variety
PURPOSE:    create a variety as a Grassmannian G(k,n). This variety has
            diemnsion k(n-k) and its Chow ring is the quotient ring of a
            polynomial ring in n-k variables q(1),...,q(n-k) which are the
            Chern classes of the tautological quotient bundle on the
            Grassmannian G(k,n), modulo the ideal generated by n-k polynomials
            which come from the Giambelli formula. The monomial ordering of this
            ring is 'wp' 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 'subBundle' and 'quotientBundle').
EXAMPLE:    example Grassmannian; shows examples
"
{
    string q;
    if (size(#)==0) {q = "q";}
    else
    {
        if (typeof(#[1]) == "string") {q = #[1];}
        else {Error("invalid optional argument");}
    }
    variety G;
    G.dimension = k*(n-k);
    execute("ring r = 0,("+q+"(1..n-k)),wp(1..n-k);");
    setring r;
    G.baseRing = r;
    int i,j;
    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);
    }
    l=insert(l,1,size(l));
    ideal rels = schur(l,v);
    int h = k+2;
    while (h<=n)
    {
        l=insert(l,1,size(l));
        rels = rels,schur(l,v);
        h++;
    }
    G.relations = rels;
    int d = k*(n-k);
    poly subBundle = reduce(logg(u,d)+k,std(rels));
    poly quotientBundle = reduce(logg(v,d)+n-k,std(rels));
    export (subBundle,quotientBundle);
    kill u,v,d,l,rels;
    return (G);
}
example
{
    "EXAMPLE:"; echo=2;
    def G24 = Grassmannian(2,4);
    G24;
    def r = G24.baseRing;
    subBundle;
    quotientBundle;
    G24.dimension;
    G24.relations;
    ChowRing(G24);
}

proc projectiveSpace(int n, list #)
"USAGE:     projectiveSpace(n); n int
RETURN:     variety
PURPOSE:    create a variety as a projective space of dimension n. The Chow ring
            is the quotient ring in one variable h modulo the ideal generated by
            h^(n+1).
EXAMPLE:    example projectiveSpace; shows examples
"
{
    string h;
    if (size(#)==0) {h = "h";}
    else
    {
        if (typeof(#[1]) == "string") {h = #[1];}
        else {Error("invalid optional argument");}
    }
    variety P;
    P.dimension = n;
    execute("ring r = 0, ("+h+"), wp(1);");
    setring r;
    P.baseRing = r;
    ideal rels = var(1)^(n+1);
    P.relations = rels;
    poly u = 1;
    poly v = 1 + var(1);
    list l;
    int i;
    for (i=1;i<=n;i++)
    {
        l=insert(l,1,size(l));
        u=u+(-1)^i*schur(l,v);
    }
    poly subBundle = reduce(logg(u,n)+n,std(rels));
    poly quotientBundle = reduce(logg(v,n)+1,std(rels));
    export(subBundle,quotientBundle);
    kill rels,u,v,l;
    return (P);
}
example
{
    "EXAMPLE:"; echo=2;
    def P = projectiveSpace(3);
    P;
    P.dimension;
    def r = P.baseRing;
    setring r;
    P.relations;
    ChowRing(P);
}

proc projectiveBundle(sheaf S, list #)
"USAGE:     projectiveBundle(S); S sheaf
RETURN:     variety
PURPOSE:    create a variety which we work on.
EXAMPLE:    example projectiveBundle; shows examples
"
{
    string z;
    if (size(#)==0) {z = "z";}
    else
    {
        if (typeof(#[1]) == "string") {z = #[1];}
        else {Error("invalid optional argument");}
    }
    variety A;
    def B = S.currentVariety;
    def R = B.baseRing;
    setring R;
    ideal rels = B.relations;
    int r = rankSheaf(S);
    A.dimension = r - 1 + B.dimension;
    poly c = totalChernClass(S);
    execute("ring P = 0, ("+z+"), wp(1);");
    def CR = P + R;
    setring CR;
    A.baseRing = CR;
    poly c = imap(R,c);
    ideal rels = imap(R,rels);
    poly g = var(1)^r;
    int i;
    for (i=1;i<=r;i++) {g=g+var(1)^(r-i)*part(c,i);}
    A.relations = rels,g;
    poly u = 1 + var(1);
    poly f = logg(u,A.dimension)+1;
    poly QuotientBundle = reduce(f,std(A.relations));
    export (QuotientBundle);
    kill f,rels;
    return (A);
}
example
{
    "EXAMPLE:"; echo=2;
    def G35 = Grassmannian(3,5);
    def R = G35.baseRing;
    setring R;
    def S = makeSheaf(G35, subBundle);
    def B = symmetricPowerSheaf(dualSheaf(S),2);
    def PB = projectiveBundle(B);
    PB;
    def P = PB.baseRing;
    setring P;
    QuotientBundle;
    ChowRing(PB);
}

proc productVariety(variety U, variety V)
{
    //def br = basering;
    def ur = U.baseRing; setring ur;
    ideal ii1 = U.relations;
    def vr = V.baseRing; setring vr;
    ideal ii2 = V.relations;
    variety W;
    W.dimension = U.dimension + V.dimension;
    def temp = ringtensor(ur,vr);
    setring temp;
    W.baseRing = temp;
    ideal i1 = imap(ur,ii1);
    ideal i2 = imap(vr,ii2);
    W.relations = i1 + i2;
    //setring br;
    return (W);
}

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

proc integral(variety V, poly f)
"USAGE:     integral(V,f); V variety, f poly
RETURN:     int
PURPOSE:    compute the intersection numbers.
EXAMPLE:    example integral; shows examples
"
{
    def R = V.baseRing;
    setring R;
    ideal rels = V.relations;
    return (leadcoef(reduce(f,std(rels))));
}
example
{
    "EXAMPLE:"; echo=2;
    def G24 = Grassmannian(2,4);
    def R = G24.baseRing;
    setring R;
    integral(G24,q(1)^4);
}

proc SchubertClass(list p)
"USAGE:     SchubertClass(p); p list
RETURN:     poly
PURPOSE:    compute the Schubert class on a Grassmannian.
EXAMPLE:    example SchubertClass; shows examples
"
{
    def R = basering;
    setring R;
    poly f = 1;
    if (size(p) == 0) {return (f);}
    int i;
    for (i=1;i<=nvars(R);i++)
    {
        f = f + var(i);
    }
    return (schur(p,f));
}
example
{
    "EXAMPLE:"; echo=2;
    def G = Grassmannian(2,4);
    def r = G.baseRing;
    setring r;
    list p = 1,1;
    SchubertClass(p);
}

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

proc dualPartition(int k, int n, list p)
{
    while (size(p) < k)
    {
        p = insert(p,0,size(p));
    }
    int i;
    list l;
    for (i=1;i<=size(p);i++)
    {
        l[i] = n-k-p[size(p)-i+1];
    }
    return (l);
}

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

////////////////////////////////////////////////////////////////////////////////
////////// Procedures concerned with sheaves ///////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

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

proc sheaf_print(sheaf S)
{
    def V = S.currentVariety;
    def R = V.baseRing;
    setring R;
    poly f = S.ChernCharacter;
    "A sheaf of rank ", int(part(f,0));
}

proc makeSheaf(variety V, poly ch)
{
    def R = basering;
    sheaf S;
    S.currentVariety = V;
    S.ChernCharacter = ch;
    return(S);
}

proc currentVariety(sheaf S)
{
    return (S.currentVariety);
}

proc ChernCharacter(sheaf S)
{
    return (S.ChernCharacter)
}

proc rankSheaf(sheaf S)
"USAGE:     rankSheaf(S); S sheaf
RETURN:     int
INPUT:      S -- a sheaf
OUTPUT:     an integer which is the rank of the sheaf S.
EXAMPLE:    example rankSheaf(S); shows examples
"
{
    def V = S.currentVariety;
    def R = V.baseRing;
    setring R;
    poly f = S.ChernCharacter;
    return (int(part(f,0)));
}
example
{
    "EXAMPLE:"; echo=2;
    def G24 = Grassmannian(2,4);
    def R = G24.baseRing;
    setring R;
    def S = makeSheaf(G24,subBundle);
    rankSheaf(S);
}

proc totalChernClass(sheaf S)
"USAGE:     totalChernClass(S); f sheaf
RETURN:     poly
INPUT:      S -- a sheaf
OUTPUT:     a polynomial which is the total Chern class of the sheaf S
EXAMPLE:    example totalChernClass(S); shows examples
"
{
    def V = S.currentVariety;
    int d = V.dimension;
    def R = V.baseRing;
    setring R;
    poly ch = S.ChernCharacter;
    poly f = expp(ch,d);
    ideal rels = std(V.relations);
    return (reduce(f,rels));
}

proc ChernClass(sheaf S, int i)
{
    return (part(totalChernClass(S),i));
}

proc topChernClass(sheaf S)
"USAGE:     topChernClass(S); S sheaf
RETURN:     poly
INPUT:      S -- a sheaf
OUTPUT:     the top Chern class of the sheaf S
EXAMPLE:    example topChernClass(S); shows examples
"
{
    return (ChernClass(S,rankSheaf(S)));
}
example
{
    "EXAMPLE:"; echo=2;
    def G24 = Grassmannian(2,4);
    def R = G24.baseRing;
    setring R;
    def S = makeSheaf(G24,quotientBundle);
    def B = symmetricPowerSheaf(S,3);
    topChernClass(B);
}

proc totalSegreClass(sheaf S)
"USAGE:     totalSegreClass(S); S sheaf
RETURN:     poly
INPUT:      S -- a sheaf
OUTPUT:     a polynomial which is the total Segre class of the sheaf S.
EXAMPLE:    example totalSegreClass(S); shows examples
"
{
    //def D = dualSheaf(S);
    def V = S.currentVariety;
    def R = V.baseRing;
    setring R;
    poly f = totalChernClass(S);
    poly g;
    int d = V.dimension;
    ideal rels = std(V.relations);
    if (f == 1) {return (1);}
    else
    {
        poly t,h;
        int i,j;
        for (i=0;i<=d;i++) {t = t + (1-f)^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);
    def R = G24.baseRing;
    setring R;
    def S = makeSheaf(G24,subBundle);
    totalSegreClass(S);
}

proc dualSheaf(sheaf S)
"USAGE:     dualSheaf(S); S sheaf
RETURN:     sheaf
INPUT:      S -- a sheaf
OUTPUT:     the dual of the sheaf S
EXAMPLE:    example dualSheaf(S); shows examples
"
{
    def V = S.currentVariety;
    int d = V.dimension;
    def R = V.baseRing;
    setring R;
    poly ch = S.ChernCharacter;
    poly f = adams(ch,-1);
    sheaf D;
    D.currentVariety = V;
    D.ChernCharacter = f;
    return (D);
}
example
{
    "EXAMPLE:"; echo=2;
    def G24 = Grassmannian(2,4);
    def R = G24.baseRing;
    setring R;
    def S = makeSheaf(G24,subBundle);
    dualSheaf(S);
}

proc tensorSheaf(sheaf A, sheaf B)
"USAGE:     tensorSheaf(A,B); A sheaf, B sheaf
RETURN:     sheaf
INPUT:      A, B -- two sheaves
OUTPUT:     the tensor of two sheaves A and B
EXAMPLE:    example tensorSheaf(A,B); shows examples
"
{
    sheaf S;
    def V1 = A.currentVariety;
    def V2 = B.currentVariety;
    def R1 = V1.baseRing;
    setring R1;
    poly c1 = A.ChernCharacter;
    def R2 = V2.baseRing;
    setring R2;
    poly c2 = B.ChernCharacter;
    if (nvars(R1) < nvars(R2))
    {
        S.currentVariety = V2;
        poly c = imap(R1,c1);
        poly f = parts(c*c2,0,V2.dimension);
        S.ChernCharacter = f;
        return (S);
    }
    else
    {
        setring R1;
        S.currentVariety = V1;
        poly c = imap(R2,c2);
        poly f = parts(c1*c,0,V1.dimension);
        S.ChernCharacter = f;
        return (S);
    }
}
example
{
    "EXAMPLE:"; echo=2;
    def G34 = Grassmannian(3,4);
    def R = G34.baseRing;
    setring R;
    def S = makeSheaf(G34,subBundle);
    def Q = makeSheaf(G34, quotientBundle);
    S*Q;
}

proc symmetricPowerSheaf(sheaf S, int n)
"USAGE:     symmetricPowerSheaf(S,n); S sheaf, n int
RETURN:     sheaf
INPUT:      S -- a sheaf
            n -- an integer
OUTPUT:     the n-th symmetric power of the sheaf S
EXAMPLE:    example symmetricPowerSheaf(S,n); shows examples
"
{
    def V = S.currentVariety;
    def R = V.baseRing;
    setring R;
    int r = rankSheaf(S);
    int d = V.dimension;
    int i,j,m;
    poly f = S.ChernCharacter;
    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];
    }
    sheaf A;
    A.currentVariety = V;
    A.ChernCharacter = result;
    return (A);
}
example
{
    "EXAMPLE:"; echo=2;
    def G24 = Grassmannian(2,4);
    def R = G24.baseRing;
    setring R;
    def S = makeSheaf(G24,quotientBundle);
    def B = symmetricPowerSheaf(S,3);
    B;
}

proc minusSheaf(sheaf A, sheaf B)
{
    sheaf S;
    def V1 = A.currentVariety;
    def V2 = B.currentVariety;
    def R1 = V1.baseRing;
    setring R1;
    poly c1 = A.ChernCharacter;
    def R2 = V2.baseRing;
    setring R2;
    poly c2 = B.ChernCharacter;
    if (nvars(R1) < nvars(R2))
    {
        S.currentVariety = V2;
        poly c = imap(R1,c1);
        S.ChernCharacter = c - c2;
        return (S);
    }
    else
    {
        setring R1;
        S.currentVariety = V1;
        poly c = imap(R2,c2);
        S.ChernCharacter = c1 - c;
        return (S);
    }
}

proc plusSheaf(sheaf A, sheaf B)
{
    sheaf S;
    def V1 = A.currentVariety;
    def V2 = B.currentVariety;
    def R1 = V1.baseRing;
    setring R1;
    poly c1 = A.ChernCharacter;
    def R2 = V2.baseRing;
    setring R2;
    poly c2 = B.ChernCharacter;
    if (nvars(R1) < nvars(R2))
    {
        S.currentVariety = V2;
        poly c = imap(R1,c1);
        S.ChernCharacter = c + c2;
        return (S);
    }
    else
    {
        setring R1;
        S.currentVariety = V1;
        poly c = imap(R2,c2);
        S.ChernCharacter = c1 + c;
        return (S);
    }
}

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

proc geometricMultiplicity(ideal I, ideal J)
"USAGE:     geometricMultiplicity(I,J); I,J = ideals
RETURN:     int
INPUT:      I, J -- two ideals
OTPUT:      an integer which is the intersection multiplicity of two subvarieties
            defined by the ideals I, J at the origin.
EXAMPLE:    example geometricMultiplicity(I,J); shows examples
"
{
    def R = basering;
    setring R;
    ideal K = I + J;
    int v = vdim(std(K));
    return (v);
}
example
{
    "EXAMPLE:"; echo=2;
    ring r = 0, (x,y,z,w), ds;
    ideal I = xz,xw,yz,yw;
    ideal J = x-z,y-w;
    geometricMultiplicity(I,J);
}

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

proc serreMultiplicity(ideal I, ideal J)
"USAGE:     serreMultiplicity(I,J); I,J = ideals
RETURN:     int
INPUT:      I, J -- two ideals
OUTPUT:     the intersection multiplicity (defined by J. P. Serre)
            of two subvarieties defined by the ideals I, J at the origin.
EXAMPLE:    example serreMultiplicity(I,J); shows examples
"
{
    def R = basering;
    setring R;
    int i = 0;
    int s = 0;
    module m = std(Tor(i,I,J));
    int t = sum(hilb(m,2));
    kill m;
    while (t != 0)
    {
        s = s + ((-1)^i)*t;
        i++;
        module m = std(Tor(i,I,J));
        t = sum(hilb(m,2));
        kill m;
    }
    return (s);
}
example
{
    "EXAMPLE:"; echo=2;
    ring r = 0, (x,y,z,w), ds;
    ideal I = xz,xw,yz,yw;
    ideal J = x-z,y-w;
    serreMultiplicity(I,J);
}