As the title of this blog suggests math, it's time to put this in practice. For one of my bachelor projects I had to do exercises from Fulton's Algebraic Curves. People who know me in real life are familiar with my total lack of perseverance for actual calculations. And for one of the exercise I had to do calculations for 15 intersection numbers (fair enough, I didn't have to do them all, but I like to do things properly). So instead of doing the calculations by hand, I decided to implement this in Sage. No other implementation of this is known to me, except for Magma's which is non-free in all meanings of the word.

What is it about? I refer you to the Wikipedia article and Fulton's excellent book. But it concerns how plane curves intersect and how their intersections behave in the local rings to these points. They lead to Max Noether's theorem, also known as the AF+BG theorem and Bézout's theorem which is a global statement relating the local information to the entire curve. Interesting stuff, by all means :).

So, having written an implementation for myself, I decided to share this. 

def intersection_number(F, G, point = (0,0)):
    (x,y) = determine_variables(F, G)

    # translate both curves to origin and calculate it there
    F = F.subs({x:x + point[0], y:y + point[1]})
    G = G.subs({x:x + point[0], y:y + point[1]})

    # if $F(0,0)\neq 0$ or $G(0,0)\neq 0$ they don't intersect in the origin
    if F.subs({x:0, y:0}) != 0 or G.subs({x:0, y:0}) != 0: return 0

    # if $F$ or $G$ are zero they don't intersect properly
    if F == 0 or G == 0: return Infinity

    # we only look at factors of $x$
    f = F.subs({y:0})
    g = G.subs({y:0})

    # $F$ contains a component $y=0$
    if f == 0:
        # $G$ contains a component $y=0$ too, no proper intersection
        if g == 0: return infinity
        # remove common $y^n$ in $F$, count degree of $x$ in $G$ and recurse
        else:
            f = F.quo_rem(y)[0]
            return ldegree(g) + intersection_number(f, G)
    # $F$ does not contain a component $y=0$
    else:
        # $G$ *does* contain a component $y=0$
        if g == 0:
            g = G.quo_rem(y)[0]
            return ldegree(f) + intersection_number(F, g)
        # we recurse as in condition (7) by removing factors of $x$
        else:
            a, b = f.lc(), g.lc()
            r, s = f.degree(), g.degree()

            # we drop the highest degree term
            if r <= s:
                return intersection_number(F, a*G - b*x^(s-r)*F)
            else:
                return intersection_number(b*F - a*x^(r-s)*G, G)

def ldegree(F):
    minimum = infinity

    for (n, m) in F.dict():
        minimum = min(minimum, n)

    return minimum

def determine_variables(F, G):
    if len(F.variables()) == 2:
        return F.variables()

    if len(G.variables()) == 2:
        return G.variables()

    if len(F.variables()) == len(G.variables()) == 1:
        if F.variables() == G.variables():
            return (F.variable(0), 0)
        else:
            return (G.variable(0), F.variable(0))

    return (0,0)

The algorithm is pretty straightforward: based on Fulton's proof of the unicity of the intersection number all the operations I perform are related to reductions in the local rings. The only non-obvious parts might be ldegree and determine_variables. The first determines the lowest degree of the first, assuming the polynomial is expressed in 2 variables while the latter determines which variables the algorithm can use. It has no a priori means of determining which variables are used (unless by supplying them as arguments, which would be ugly). I'm sure this can look better, but for the moment it works. I've got a big red niceify and contact Sage developers note on my list by the way. For the moment it's implemented for algebraic schemes, but I'm not yet sure whether it's as it should be.

In case you wish to use it, the following code solves all exercises from Fulton's Algebraic Curves: 

import intersection_number

P. = PolynomialRing(QQ, 2)
F = []
F.append(y-x^2);
F.append(y^2-x^3+x);
F.append(y^2-x^3);
F.append(y^2-x^3-x^2);
F.append((x^2+y^2)^2+3*x^2*y-y^3);
F.append((x^2+y^2)^3-4*x^2*y^2);

for i in range(0,6):
    for j in range(i + 1,6):
        print '(%(i)d, %(j)d): %(m)d' % {'i': i, 'j': j, \
              'm': intersection_number.intersection_number(F[i], F[j])}

To use this: perform sage intersection_number.sage first and then you can use it by calling sage example.sage.