{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fa2add37-d5af-4239-886a-a1acd28e5960",
   "metadata": {},
   "outputs": [],
   "source": [
    "def intersection_number(F, G, point = (0,0)):\n",
    "    (x,y) = determine_variables(F, G)\n",
    "\n",
    "    # translate both curves to origin and calculate it there\n",
    "    F = F.subs({x:x + point[0], y:y + point[1]})\n",
    "    G = G.subs({x:x + point[0], y:y + point[1]})\n",
    "\n",
    "    # if $F(0,0)\\neq 0$ or $G(0,0)\\neq 0$ they don't intersect in the origin\n",
    "    if F.subs({x:0, y:0}) != 0 or G.subs({x:0, y:0}) != 0: return 0\n",
    "\n",
    "    # if $F$ or $G$ are zero they don't intersect properly\n",
    "    if F == 0 or G == 0: return Infinity\n",
    "\n",
    "    # we only look at factors of $x$\n",
    "    f = F.subs({y:0})\n",
    "    g = G.subs({y:0})\n",
    "\n",
    "    # $F$ contains a component $y=0$\n",
    "    if f == 0:\n",
    "        # $G$ contains a component $y=0$ too, no proper intersection\n",
    "        if g == 0: return infinity\n",
    "        # remove common $y^n$ in $F$, count degree of $x$ in $G$ and recurse\n",
    "        else:\n",
    "            f = F.quo_rem(y)[0]\n",
    "            return ldegree(g) + intersection_number(f, G)\n",
    "    # $F$ does not contain a component $y=0$\n",
    "    else:\n",
    "        # $G$ *does* contain a component $y=0$\n",
    "        if g == 0:\n",
    "            g = G.quo_rem(y)[0]\n",
    "            return ldegree(f) + intersection_number(F, g)\n",
    "        # we recurse as in condition (7) by removing factors of $x$\n",
    "        else:\n",
    "            a, b = f.lc(), g.lc()\n",
    "            r, s = f.degree(), g.degree()\n",
    "\n",
    "            # we drop the highest degree term\n",
    "            if r <= s:\n",
    "                return intersection_number(F, a*G - b*x^(s-r)*F)\n",
    "            else:\n",
    "                return intersection_number(b*F - a*x^(r-s)*G, G)\n",
    "\n",
    "def ldegree(F):\n",
    "    minimum = infinity\n",
    "\n",
    "    for (n, m) in F.dict():\n",
    "        minimum = min(minimum, n)\n",
    "\n",
    "    return minimum\n",
    "\n",
    "def determine_variables(F, G):\n",
    "    if len(F.variables()) == 2:\n",
    "        return F.variables()\n",
    "\n",
    "    if len(G.variables()) == 2:\n",
    "        return G.variables()\n",
    "\n",
    "    if len(F.variables()) == len(G.variables()) == 1:\n",
    "        if F.variables() == G.variables():\n",
    "            return (F.variable(0), 0)\n",
    "        else:\n",
    "            return (G.variable(0), F.variable(0))\n",
    "\n",
    "    return (0,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9cec71ba-3055-4771-8c8c-28e8ecc2d93b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 1): 1\n",
      "(0, 2): 3\n",
      "(0, 3): 2\n",
      "(0, 4): 4\n",
      "(0, 5): 6\n",
      "(1, 2): 2\n",
      "(1, 3): 2\n",
      "(1, 4): 3\n",
      "(1, 5): 6\n",
      "(2, 3): 4\n",
      "(2, 4): 7\n",
      "(2, 5): 10\n",
      "(3, 4): 6\n",
      "(3, 5): 8\n",
      "(4, 5): 14\n"
     ]
    }
   ],
   "source": [
    "P.<x,y> = PolynomialRing(QQ, 2)\n",
    "F = []\n",
    "F.append(y-x^2);\n",
    "F.append(y^2-x^3+x);\n",
    "F.append(y^2-x^3);\n",
    "F.append(y^2-x^3-x^2);\n",
    "F.append((x^2+y^2)^2+3*x^2*y-y^3);\n",
    "F.append((x^2+y^2)^3-4*x^2*y^2);\n",
    "\n",
    "for i in range(0,6):\n",
    "    for j in range(i + 1,6):\n",
    "        print('(%(i)d, %(j)d): %(m)d' % {'i': i, 'j': j, 'm': intersection_number(F[i], F[j])} )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ffbd7348-91a2-49fb-b7e2-eef7af8e926e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "intersection_number( x+y^2, x+y^2-x^3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "dcdf8d54-0638-4f1a-8c06-a8b27ce9d187",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.<x,y> = AffineSpace(CC, 2)\n",
    "X = A.subscheme([x+y^2])\n",
    "Y = A.subscheme([x+y^2-x^3])\n",
    "Q1 = Y([0,0])\n",
    "Q1.intersection_multiplicity(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "917c84cb-f6be-4202-9b4e-dee10c76c7e8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.<x,y,z> = AffineSpace(CC, 3)\n",
    "X = A.subscheme([y^2 - x^7*z])\n",
    "Q1 = X([1,1,1])\n",
    "Q1.multiplicity()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "c23af6ff-907e-4c79-ae21-d28b3d02ee4b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q2 = X([0,0,2])\n",
    "Q2.multiplicity() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "093fbda8-33a8-4f32-84c3-da86d53ed3af",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(x, y, z) |--> (-7*x^6*z, 2*y, -x^7)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[[x == 0, y == 0, z == r3]]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Note singularities all along z-axis: https://www.desmos.com/3d/tvj2vzsjoi\n",
    "f(x,y,z) = y^2 - x^7*z\n",
    "display( f.gradient() )\n",
    "solve( list(f.gradient()), (x,y,z) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e3b36460-c1ac-4758-a207-9da921e7fb08",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Gathmann example 2.13\n",
    "intersection_number( y^2-x^3, x^2-y^3 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f9841065-b209-4189-aacc-e79d463e4f64",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.<x,y> = AffineSpace(CC, 2)\n",
    "X = A.subscheme([ y^2-x^3 ])\n",
    "Y = A.subscheme([ x^2-y^3 ])\n",
    "Q1 = Y([0,0])\n",
    "Q1.intersection_multiplicity(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a70054e2-f126-4ac2-8725-3bd7cc302417",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q1 = Y([1,1])\n",
    "Q1.intersection_multiplicity(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "63ed262e-e748-4ef1-8977-b68b0f6801dd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "9a661325-24ac-425a-a782-c5a8747e0c30",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(x, y, z) |--> (-7*x^6*z, 2*y, -x^7)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd6020ba-8e1a-4741-b0cc-14888d699f11",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 10.4.beta5",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}