{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Reduziert den Ausdruck <term> im Quotienten-Ring modulo der Gleichung der elliptischen Kurve.\n",
    "# @param E .. elliptische Kurve; Parameter a1,a2,a3,a4 und a6 dürfen c enthalten\n",
    "# @param term .. Element aus Q(x,y,c) = Quotient aus zwei Polynomen mit\n",
    "#                Variablen x und y und Parameter c\n",
    "# @returns term modulo Q(x,y,c)/I, wobei Q(x,y,c) = Polynom/Polynom ist \n",
    "#          und I das Ideal der Gleichung der elliptischen Kurve E\n",
    "def ReduceInQuotientRing(term, E):\n",
    "    a1, a2, a3, a4, a6 = E.a_invariants();\n",
    "    \n",
    "    Q.<x,y,c> = QQ['x','y','c'];\n",
    "    J = Q.ideal(y^2 + a1*x*y + a3*y - x^3 - a2*x^2 - a4*x - a6);\n",
    "    R.<X,Y,C> = QuotientRing(Q, J);\n",
    "    \n",
    "    # In <term> werden die Variablen ersetzt gemäß x -> X, y -> Y, c -> C.\n",
    "    termXYZ = sage_eval(str(term), locals={'x':X, 'y':Y, 'c':C});\n",
    "    \n",
    "    # Reduzieren im Ideal modulo der Gleichung der elliptischen Kurve.\n",
    "    reduc = R(termXYZ).reduce(J.gens());\n",
    "    \n",
    "    # Rückbenennung X -> x, Y -> y, C -> c\n",
    "    return sage_eval(str(reduc), locals={'X':x, 'Y':y, 'C':c});"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Bestimmt die Gleichung der Gerade durch die Punkte Q und R, die auf der elliptischen Kurve E liegen.\n",
    "# Vorlage: _line_ in https://github.com/sagemath/sage/blob/1961f9430ee403bb5632853af3e18d2d2c858187/src/sage/\n",
    "#                             schemes/elliptic_curves/ell_point.py\n",
    "# @param E   .. elliptische Kurve; die Koeffizienten a1,a2,a3,a4 und a6 dürfen \n",
    "#               Parameter c enthalten\n",
    "# @param Q,R .. Zwei Punkte der elliptischen Kurve\n",
    "def EllipticLine(E, Q, R):\n",
    "    if (Q == E(0)) or (R == E(0)):\n",
    "        if Q == R:\n",
    "            return E.base_field().one()\n",
    "        if Q == E(0):\n",
    "            return x - R[0]\n",
    "        if R == E(0):\n",
    "            return x - Q[0]\n",
    "    elif Q != R:\n",
    "        if Q[0] == R[0]:\n",
    "            return x - Q[0]\n",
    "        else:\n",
    "            m = (R[1] - Q[1])  / (R[0] - Q[0])\n",
    "            return y - Q[1] - m * (x - Q[0])\n",
    "    else:\n",
    "        a1, a2, a3, a4, a6 = E.a_invariants()\n",
    "        numerator = (3*Q[0]**2 + 2*a2*Q[0]+a4-a1*Q[1])\n",
    "        denominator = (2*Q[1] + a1*Q[0] + a3)\n",
    "        if denominator == 0:\n",
    "            return x - Q[0]\n",
    "        else:\n",
    "            m = numerator / denominator\n",
    "            return y - Q[1] - m * (x - Q[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Vorlage: _miller_ in https://github.com/sagemath/sage/blob/1961f9430ee403bb5632853af3e18d2d2c858187/src/sage/\n",
    "#                             schemes/elliptic_curves/ell_point.py\n",
    "# @param E .. elliptische Kurve; die Koeffizienten a1,a2,a3,a4 und a6 dürfen Parameter c enthalten\n",
    "# @param P .. Punkt der elliptischen Kurve (Anwendungsfall: P ist n-Torsionspunkt von E)\n",
    "# @param n .. natürliche Zahl (Anwendungsfall: n ist die Ordnung von P)\n",
    "# @return Gleichung für rationale Funktion f:E -> CC in den Variablen x, y mit Parameter c\n",
    "#         mit Divisor div(f) = n[P] - [nP] - (n-1)[OO]. Falls P ein n-Torsionspunkt von E\n",
    "#         ist, hat f also den Divisor div(f) = n[P] - n[OO].\n",
    "def WeilFunction(E, P, n):\n",
    "    \n",
    "    if not P.curve() is E:\n",
    "            raise ValueError(\"point must be on the elliptic curve\")\n",
    "    \n",
    "    if P == E(0):\n",
    "        raise ValueError(\"P cannot be the basepoint of the elliptic curve\")\n",
    "        \n",
    "    if n.is_zero():\n",
    "        raise ValueError(\"n must be nonzero.\")\n",
    "    \n",
    "    n_is_negative = False\n",
    "    if n < 0:\n",
    "        n = n.abs()\n",
    "        n_is_negative = True\n",
    "     \n",
    "    one = E.base_field().one()\n",
    "    t = one\n",
    "    V = P\n",
    "    S = 2*V\n",
    "    nbin = n.bits()\n",
    "    i = n.nbits() - 2\n",
    "    while i > -1:\n",
    "        S = 2*V\n",
    "        ell = EllipticLine(E, V, V)\n",
    "        vee = EllipticLine(E, S, -S)\n",
    "        t = (t**2)*(ell/vee)\n",
    "        V = S\n",
    "        if nbin[i] == 1:\n",
    "            S = V + P\n",
    "            ell = EllipticLine(E, V, P)\n",
    "            vee = EllipticLine(E, S, -S)\n",
    "            t = t*(ell/vee)\n",
    "            V = S\n",
    "        i = i-1\n",
    "    \n",
    "    if n_is_negative:\n",
    "        vee = EllipticLine(E, V, -V)\n",
    "        t = 1/(t*vee)\n",
    "        \n",
    "    # Im Quotientenring reduzieren, damit eine nennerfreie Darstellung entsteht\n",
    "    return ReduceInQuotientRing(t, E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Drückt x^n für n >= 3 durch kleinere x-Potenzen aus.\n",
    "# @param E   .. elliptische Kurve; die Koeffizienten a1,a2,a3,a4 und a6 dürfen \n",
    "#               Parameter c enthalten\n",
    "# @param n   .. natürliche Zahl\n",
    "def XPower(n, E):\n",
    "    a1, a2, a3, a4, a6 = E.a_invariants();\n",
    "    \n",
    "    if n <= 2:\n",
    "        return x^n\n",
    "    if n >= 3:\n",
    "        return expand(x^(n-3) * (y^2 + a1*x*y + a3*y - a2*x^2 - a4*x - a6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Ersetzt in <term> alle Vorkommnisse von <search> durch (<repl>).\n",
    "# Unterstützt werden die Variablen x, y und c.\n",
    "def ReplaceSubTerm(term, search, repl):\n",
    "    term_str = str(term);\n",
    "    term_repl = term_str.replace(str(search), \"(\" + str(repl) + \")\");\n",
    "    return sage_eval(term_repl, locals={'x':x, 'y':y, 'c':c})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Ersetzt in <term> alle Vorkommnisse von <x^n> durch Potenzen x^2, x, 1.\n",
    "# Parameter term: Polynom in x, y unc c.\n",
    "# @param E    .. elliptische Kurve; die Koeffizienten a1,a2,a3,a4 und a6 dürfen \n",
    "#                Parameter c enthalten\n",
    "# @param term .. Polynom in x, y unc c.\n",
    "def ReplaceHighXMonoms(term, E):\n",
    "    # Polynomring als \"Monom-Wörterbuch\"\n",
    "    R.<x,y,c> = sage.rings.polynomial.multi_polynomial_ring.MPolynomialRing_polydict(QQ, 3, order='lex');\n",
    "    \n",
    "    # term als Polynom in x, y, c interpretieren\n",
    "    poly = R(term)\n",
    "    \n",
    "    # Maximalgrad von x im Polynom <term> bestimmen.\n",
    "    max_x_degree = poly.degree(x)\n",
    "        \n",
    "    # Iteration von n bis 3\n",
    "    for e in range(max_x_degree, 2, -1):         \n",
    "        # Ersetze x^e durch kleinere Potenzen.\n",
    "        term = ReplaceSubTerm(term, x^e, XPower(e, E));   \n",
    "        # Ausmultiplizieren, damit im nächsten Schritt x^(e-1) ersetzt werden kann.\n",
    "        term = expand(term); \n",
    "    return term;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Interpretiert <term> als Polynom in x und y und faktorisiert die Koeffizienten dieses Polynoms.\n",
    "# Rückgabe erfolgt als String.\n",
    "# @param term .. Polynom in x, y und c\n",
    "def FactorCoefficientsXYPoly(term):\n",
    "    # Polynomring als \"Monom-Wörterbuch\"\n",
    "    R.<x,y,c> = sage.rings.polynomial.multi_polynomial_ring.MPolynomialRing_polydict(QQ, 3, order='lex');\n",
    "    \n",
    "    # term als Polynom in x, y, c interpretieren\n",
    "    poly = R(term)\n",
    "    \n",
    "    # Höchsten Grad von x und y bestimmen\n",
    "    degree_x = poly.degree(x)\n",
    "    degree_y = poly.degree(y)\n",
    "    \n",
    "    # Leeres Polynom und leeren Polynom-String anlegen.\n",
    "    new_poly = 0\n",
    "    new_poly_str = \"\"\n",
    "    \n",
    "    # Über alle Monome x^a y^b iterieren\n",
    "    for b in range(degree_y, -1, -1):\n",
    "        for a in range(degree_x, -1, -1):\n",
    "            # Koeffizient von x^a y^b lesen\n",
    "            coeff = poly.coefficient({x:a, y:b})\n",
    "            \n",
    "            # Koeffizient faktorisieren, wenn nicht 0.\n",
    "            if (coeff != 0):\n",
    "                coeff = factor(coeff);\n",
    "                \n",
    "                # Neues Polynom um Monom x^a y^b mit faktorisiertem Koeffizienten erweitern. \n",
    "                new_poly = new_poly + coeff * x^a * y^b;\n",
    "                \n",
    "                # Polynomstring um Monom x^a y^b mit faktorisiertem Koeffizienten erweitern\n",
    "                new_poly_str = new_poly_str + \"+ (\" + str(coeff) + \") * \" + str(x^a * y^b)\n",
    "                \n",
    "                # Übersichtliche formatierte Konsolenausgabe erzeugen\n",
    "                print \"+ (\", coeff, \") * \", x^a * y^b\n",
    "                \n",
    "    # Überprüfen, ob die beiden Polynome übereinstimmen.\n",
    "    diff = term - sage_eval(new_poly_str, locals={'x':x, 'y':y, 'c':c})\n",
    "    if diff != 0:\n",
    "        print \"Error: Wrong polynomial calculated. Diff = \", diff\n",
    "            \n",
    "    return new_poly_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Berechnet die Weil-Funktion für X_1(4),\n",
    "# reduziert die Potenzen x^3 und höher bzgl. der Gleichung der elliptischen Kurve\n",
    "# und faktorisiert die Koeffizienten, um eine kompakte Darstellung der Weil-Funktion zu erhalten.\n",
    "def CalcWeilFunction4():\n",
    "    # Elliptische Kurve X_1(4) und 4-Torsionspunkt P definieren\n",
    "    x, y, c = var('x y c');\n",
    "    E_4 = EllipticCurve([1, c, c, 0, 0]);\n",
    "    P_4 = E_4(0, 0);\n",
    "    \n",
    "    # Weil-Funktion mittels Miller-Algorithmus berechnen\n",
    "    # (Funktion wird gleich im Quotientenring reduziert, \n",
    "    #  damit eine nennerfreie Darstellung entsteht)\n",
    "    f = WeilFunction(E_4, P_4, 4)\n",
    "    \n",
    "    # Potenzen x^n für n >= 3 ersetzen durch x^2, x, 1\n",
    "    f = ReplaceHighXMonoms(f, E_4)\n",
    "    \n",
    "    # f als Polynom in x,y sehen und Koeffizienten faktorisieren\n",
    "    f_str = FactorCoefficientsXYPoly(f)\n",
    "    \n",
    "    return f_str;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+ ( -1 ) *  y\n",
      "+ ( 1 ) *  x^2\n",
      "CPU times: user 34 ms, sys: 66 ms, total: 100 ms\n",
      "Wall time: 99.8 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'+ (-1) * y+ (1) * x^2'"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time Weil4_str = CalcWeilFunction4(); \n",
    "Weil4_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Berechnet die Weil-Funktion für X_1(5),\n",
    "# reduziert die Potenzen x^3 und höher bzgl. der Gleichung der elliptischen Kurve\n",
    "# und faktorisiert die Koeffizienten, um eine kompakte Darstellung der Weil-Funktion zu erhalten.\n",
    "def CalcWeilFunction5():\n",
    "    # Elliptische Kurve X_1(5) und 5-Torsionspunkt P definieren\n",
    "    x, y, c = var('x y c');\n",
    "    E_5 = EllipticCurve([c+1, c, c, 0, 0]);\n",
    "    P_5 = E_5(0, 0);\n",
    "    \n",
    "    # Weil-Funktion mittels Miller-Algorithmus berechnen\n",
    "    # (Funktion wird gleich im Quotientenring reduziert, \n",
    "    #  damit eine nennerfreie Darstellung entsteht)\n",
    "    f = WeilFunction(E_5, P_5, 5)\n",
    "    \n",
    "    # Potenzen x^n für n >= 3 ersetzen durch x^2, x, 1\n",
    "    f = ReplaceHighXMonoms(f, E_5)\n",
    "    \n",
    "    # f als Polynom in x,y sehen und Koeffizienten faktorisieren\n",
    "    f_str = FactorCoefficientsXYPoly(f)\n",
    "    \n",
    "    return f_str;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+ ( 1 ) *  x*y\n",
      "+ ( 1 ) *  y\n",
      "+ ( -1 ) *  x^2\n",
      "CPU times: user 1.4 s, sys: 584 ms, total: 1.98 s\n",
      "Wall time: 1.99 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'+ (1) * x*y+ (1) * y+ (-1) * x^2'"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time Weil5_str = CalcWeilFunction5(); \n",
    "Weil5_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Berechnet die Weil-Funktion für X_1(6),\n",
    "# reduziert die Potenzen x^3 und höher bzgl. der Gleichung der elliptischen Kurve\n",
    "# und faktorisiert die Koeffizienten, um eine kompakte Darstellung der Weil-Funktion zu erhalten.\n",
    "def CalcWeilFunction6():\n",
    "    # Elliptische Kurve X_1(6) und 6-Torsionspunkt P definieren\n",
    "    x, y, c = var('x y c');\n",
    "    E_6 = EllipticCurve([1-c, \n",
    "                         -c*(c+1), \n",
    "                         -c*(c+1), \n",
    "                         0, 0]);\n",
    "    P_6 = E_6(0, 0);\n",
    "    \n",
    "    # Weil-Funktion mittels Miller-Algorithmus berechnen\n",
    "    # (Funktion wird gleich im Quotientenring reduziert, \n",
    "    #  damit eine nennerfreie Darstellung entsteht)\n",
    "    f = WeilFunction(E_6, P_6, 6)\n",
    "    \n",
    "    # Potenzen x^n für n >= 3 ersetzen durch x^2, x, 1\n",
    "    f = ReplaceHighXMonoms(f, E_6)\n",
    "    \n",
    "    # f als Polynom in x,y sehen und Koeffizienten faktorisieren\n",
    "    f_str = FactorCoefficientsXYPoly(f)\n",
    "    \n",
    "    return f_str;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+ ( 1 ) *  y^2\n",
      "+ ( (-1) * (c + 1) ) *  x*y\n",
      "+ ( (-1) * (c + 1)^2 ) *  y\n",
      "+ ( (c + 1)^2 ) *  x^2\n",
      "CPU times: user 405 ms, sys: 11 ms, total: 416 ms\n",
      "Wall time: 417 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'+ (1) * y^2+ ((-1) * (c + 1)) * x*y+ ((-1) * (c + 1)^2) * y+ ((c + 1)^2) * x^2'"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time Weil6_str = CalcWeilFunction6(); \n",
    "Weil6_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Berechnet die Weil-Funktion für X_1(7),\n",
    "# reduziert die Potenzen x^3 und höher bzgl. der Gleichung der elliptischen Kurve\n",
    "# und faktorisiert die Koeffizienten, um eine kompakte Darstellung der Weil-Funktion zu erhalten.\n",
    "def CalcWeilFunction7():\n",
    "    # Elliptische Kurve X_1(7) und 7-Torsionspunkt P definieren\n",
    "    x, y, c = var('x y c');\n",
    "    E_7 = EllipticCurve([1-c-c^2, \n",
    "                         c^2*(c+1), \n",
    "                         c^2*(c+1), \n",
    "                         0, 0]);\n",
    "    P_7 = E_7(0, 0);\n",
    "    \n",
    "    # Weil-Funktion mittels Miller-Algorithmus berechnen\n",
    "    # (Funktion wird gleich im Quotientenring reduziert, \n",
    "    #  damit eine nennerfreie Darstellung entsteht)\n",
    "    f = WeilFunction(E_7, P_7, 7)\n",
    "    \n",
    "    # Potenzen x^n für n >= 3 ersetzen durch x^2, x, 1\n",
    "    f = ReplaceHighXMonoms(f, E_7)\n",
    "    \n",
    "    # f als Polynom in x,y sehen und Koeffizienten faktorisieren\n",
    "    f_str = FactorCoefficientsXYPoly(f)\n",
    "    \n",
    "    return f_str;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+ ( c - 1 ) *  y^2\n",
      "+ ( 1 ) *  x^2*y\n",
      "+ ( (-1) * c^3 ) *  x*y\n",
      "+ ( c^4 ) *  y\n",
      "+ ( (-1) * c^4 ) *  x^2\n",
      "CPU times: user 8.52 s, sys: 214 ms, total: 8.73 s\n",
      "Wall time: 8.75 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'+ (c - 1) * y^2+ (1) * x^2*y+ ((-1) * c^3) * x*y+ (c^4) * y+ ((-1) * c^4) * x^2'"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time Weil7_str = CalcWeilFunction7(); \n",
    "Weil7_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Berechnet die Weil-Funktion für X_1(8),\n",
    "# reduziert die Potenzen x^3 und höher bzgl. der Gleichung der elliptischen Kurve\n",
    "# und faktorisiert die Koeffizienten, um eine kompakte Darstellung der Weil-Funktion zu erhalten.\n",
    "def CalcWeilFunction8():\n",
    "    # Elliptische Kurve X_1(8) und 8-Torsionspunkt P definieren\n",
    "    x, y, c = var('x y c');\n",
    "    E_8 = EllipticCurve([1-2*c^2, \n",
    "                         -c*(2*c+1)*(c+1)^2, \n",
    "                         -c*(2*c+1)*(c+1)^3, \n",
    "                         0, 0]);\n",
    "    P_8 = E_8(0, 0);\n",
    "    \n",
    "    # Weil-Funktion mittels Miller-Algorithmus berechnen\n",
    "    # (Funktion wird gleich im Quotientenring reduziert, \n",
    "    #  damit eine nennerfreie Darstellung entsteht)\n",
    "    f = WeilFunction(E_8, P_8, 8)\n",
    "    \n",
    "    # Potenzen x^n für n >= 3 ersetzen durch x^2, x, 1\n",
    "    f = ReplaceHighXMonoms(f, E_8)\n",
    "    \n",
    "    # f als Polynom in x,y sehen und Koeffizienten faktorisieren\n",
    "    f_str = FactorCoefficientsXYPoly(f)\n",
    "    \n",
    "    return f_str;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+ ( 1 ) *  x*y^2\n",
      "+ ( (2) * (c + 3/2) * (c + 1)^3 ) *  y^2\n",
      "+ ( (-2) * (c + 1)^2 ) *  x^2*y\n",
      "+ ( (-4) * (c + 1/2)^2 * (c + 1)^4 ) *  x*y\n",
      "+ ( (-4) * (c + 1/2)^2 * (c + 1)^7 ) *  y\n",
      "+ ( (4) * (c + 1/2)^2 * (c + 1)^6 ) *  x^2\n",
      "CPU times: user 1.16 s, sys: 12 ms, total: 1.17 s\n",
      "Wall time: 1.17 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'+ (1) * x*y^2+ ((2) * (c + 3/2) * (c + 1)^3) * y^2+ ((-2) * (c + 1)^2) * x^2*y+ ((-4) * (c + 1/2)^2 * (c + 1)^4) * x*y+ ((-4) * (c + 1/2)^2 * (c + 1)^7) * y+ ((4) * (c + 1/2)^2 * (c + 1)^6) * x^2'"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time Weil8_str = CalcWeilFunction8(); \n",
    "Weil8_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Berechnet die Weil-Funktion für X_1(9),\n",
    "# reduziert die Potenzen x^3 und höher bzgl. der Gleichung der elliptischen Kurve\n",
    "# und faktorisiert die Koeffizienten, um eine kompakte Darstellung der Weil-Funktion zu erhalten.\n",
    "def CalcWeilFunction9():\n",
    "    # Elliptische Kurve X_1(9) und 9-Torsionspunkt P definieren\n",
    "    x, y, c = var('x y c');\n",
    "    E_9 = EllipticCurve([c^3 + c^2 + 1, \n",
    "                         c^2 * (c+1) * (c^2+c+1), \n",
    "                         c^2 * (c+1) * (c^2+c+1), \n",
    "                         0, 0]);\n",
    "    P_9 = E_9(0, 0);\n",
    "    \n",
    "    # Weil-Funktion mittels Miller-Algorithmus berechnen\n",
    "    # (Funktion wird gleich im Quotientenring reduziert, \n",
    "    #  damit eine nennerfreie Darstellung entsteht)\n",
    "    f = WeilFunction(E_9, P_9, 9)\n",
    "    \n",
    "    # Potenzen x^n für n >= 3 ersetzen durch x^2, x, 1\n",
    "    f = ReplaceHighXMonoms(f, E_9)\n",
    "    \n",
    "    # f als Polynom in x,y sehen und Koeffizienten faktorisieren\n",
    "    f_str = FactorCoefficientsXYPoly(f)\n",
    "    \n",
    "    return f_str;\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+ ( 1 ) *  y^3\n",
      "+ ( (c - 1) * (c^2 + c + 1) ) *  x*y^2\n",
      "+ ( (c^2 + c + 1)^2 * (c^3 + 2*c - 1) ) *  y^2\n",
      "+ ( (-2) * (c - 1/2) * (c^2 + c + 1)^2 ) *  x^2*y\n",
      "+ ( c^4 * (c^2 + c + 1)^3 ) *  x*y\n",
      "+ ( c^4 * (c^2 + c + 1)^4 ) *  y\n",
      "+ ( (-1) * c^4 * (c^2 + c + 1)^4 ) *  x^2\n",
      "CPU times: user 26.4 s, sys: 52 ms, total: 26.4 s\n",
      "Wall time: 26.7 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'+ (1) * y^3+ ((c - 1) * (c^2 + c + 1)) * x*y^2+ ((c^2 + c + 1)^2 * (c^3 + 2*c - 1)) * y^2+ ((-2) * (c - 1/2) * (c^2 + c + 1)^2) * x^2*y+ (c^4 * (c^2 + c + 1)^3) * x*y+ (c^4 * (c^2 + c + 1)^4) * y+ ((-1) * c^4 * (c^2 + c + 1)^4) * x^2'"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time Weil9_str = CalcWeilFunction9(); \n",
    "Weil9_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Berechnet die Weil-Funktion für X_1(10),\n",
    "# reduziert die Potenzen x^3 und höher bzgl. der Gleichung der elliptischen Kurve\n",
    "# und faktorisiert die Koeffizienten, um eine kompakte Darstellung der Weil-Funktion zu erhalten.\n",
    "def CalcWeilFunction10():\n",
    "    # Elliptische Kurve X_1(10) und 10-Torsionspunkt P definieren\n",
    "    x, y, c = var('x y c');\n",
    "    E_10 = EllipticCurve([-c^3 - 2*c^2 + 4*c + 4, \n",
    "                          (c+1) * (c+2) * c^3, \n",
    "                          (c+1) * (c+2) * (c^2+6*c+4) * c^3, \n",
    "                          0, 0]);\n",
    "    P_10 = E_10(0, 0);\n",
    "    \n",
    "    # Weil-Funktion mittels Miller-Algorithmus berechnen\n",
    "    # (Funktion wird gleich im Quotientenring reduziert, \n",
    "    #  damit eine nennerfreie Darstellung entsteht)\n",
    "    f = WeilFunction(E_10, P_10, 10)\n",
    "    \n",
    "    # Potenzen x^n für n >= 3 ersetzen durch x^2, x, 1\n",
    "    f = ReplaceHighXMonoms(f, E_10)\n",
    "    \n",
    "    # f als Polynom in x,y sehen und Koeffizienten faktorisieren\n",
    "    f_str = FactorCoefficientsXYPoly(f)\n",
    "    \n",
    "    return f_str;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+ ( (2) * (c^2 - 2*c - 2) ) *  y^3\n",
      "+ ( 1 ) *  x^2*y^2\n",
      "+ ( (-2) * (c + 1/2) * c^4 ) *  x*y^2\n",
      "+ ( c^6 * (c^3 + 16*c^2 + 22*c + 8) ) *  y^2\n",
      "+ ( (-3) * (c + 2/3) * c^6 ) *  x^2*y\n",
      "+ ( (c + 1)^2 * c^10 ) *  x*y\n",
      "+ ( (-1) * (c + 1)^2 * c^12 * (c^2 + 6*c + 4) ) *  y\n",
      "+ ( (c + 1)^2 * c^12 ) *  x^2\n",
      "CPU times: user 12.8 s, sys: 89 ms, total: 12.9 s\n",
      "Wall time: 12.9 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'+ ((2) * (c^2 - 2*c - 2)) * y^3+ (1) * x^2*y^2+ ((-2) * (c + 1/2) * c^4) * x*y^2+ (c^6 * (c^3 + 16*c^2 + 22*c + 8)) * y^2+ ((-3) * (c + 2/3) * c^6) * x^2*y+ ((c + 1)^2 * c^10) * x*y+ ((-1) * (c + 1)^2 * c^12 * (c^2 + 6*c + 4)) * y+ ((c + 1)^2 * c^12) * x^2'"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time Weil10_str = CalcWeilFunction10(); \n",
    "Weil10_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Berechnet die Weil-Funktion für X_1(12),\n",
    "# reduziert die Potenzen x^3 und höher bzgl. der Gleichung der elliptischen Kurve\n",
    "# und faktorisiert die Koeffizienten, um eine kompakte Darstellung der Weil-Funktion zu erhalten.\n",
    "def CalcWeilFunction12():\n",
    "    # Elliptische Kurve X_1(12) und 12-Torsionspunkt P definieren\n",
    "    x, y, c = var('x y c');\n",
    "    E_12 = EllipticCurve([6*c^4 - 8*c^3 + 2*c^2 + 2*c - 1,\n",
    "                          -c * (c-1)^2 * (2*c-1) * (2*c^2-2*c+1) * (3*c^2-3*c+1),\n",
    "                          -c * (c-1)^5 * (2*c-1) * (2*c^2-2*c+1) * (3*c^2-3*c+1),\n",
    "                          0, 0]);\n",
    "    P = E_12(0, 0);\n",
    "    \n",
    "    # Miller-Algorithmus durchführen\n",
    "    f1 = 1;\n",
    "    f2 = f1^2 * EllipticLine(E_12, P, P) / EllipticLine(E_12, 2*P, -2*P); \n",
    "    f2 = factor(f2);\n",
    "    f3 = f2 * EllipticLine(E_12, 2*P, P) / EllipticLine(E_12, 3*P, -3*P);\n",
    "    f3 = factor(f3);\n",
    "    f6 = f3^2 * EllipticLine(E_12, 3*P, 3*P) / EllipticLine(E_12, 6*P, -6*P);\n",
    "    f6 = factor(f6);\n",
    "    f12a = f6^2 * EllipticLine(E_12, 6*P, 6*P);\n",
    "    f12a = factor(f12a);\n",
    "    f12 = f12a / EllipticLine(E_12, 12*P, -12*P);\n",
    "    f12 = factor(f12);\n",
    "    \n",
    "    # Im Quotientenring reduzieren, damit eine nennerfreie Darstellung \n",
    "    # entsteht\n",
    "    f = ReduceInQuotientRing(f12, E_12);\n",
    "    \n",
    "    # Potenzen x^n für n >= 3 ersetzen durch x^2, x, 1\n",
    "    f = ReplaceHighXMonoms(f, E_12)\n",
    "    \n",
    "    # f als Polynom in x,y sehen und Koeffizienten faktorisieren\n",
    "    f_str = FactorCoefficientsXYPoly(f)\n",
    "    \n",
    "    return f_str;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+ ( 1 ) *  y^4\n",
      "+ ( (12) * (c^2 - 4/3*c + 1/2) * (c^2 - c + 1/2) ) *  x*y^3\n",
      "+ ( (144) * (c - 1)^4 * (c^2 - c + 1/2)^2 * (c^4 - 5/2*c^3 + 8/3*c^2 - 49/36*c + 5/18) ) *  y^3\n",
      "+ ( (60) * (c - 1)^2 * (c^2 - 6/5*c + 2/5) * (c^2 - c + 1/2)^2 ) *  x^2*y^2\n",
      "+ ( (1008) * (c - 1)^4 * (c^2 - 8/7*c + 5/14) * (c^2 - c + 1/3)^2 * (c^2 - c + 1/2)^3 ) *  x*y^2\n",
      "+ ( (1728) * (c - 1)^8 * (c^2 - c + 1/3)^2 * (c^2 - c + 1/2)^4 * (c^4 - 7/2*c^3 + 19/4*c^2 - 11/4*c + 7/12) ) *  y^2\n",
      "+ ( (2592) * (c - 1)^6 * (c^2 - 10/9*c + 1/3) * (c^2 - c + 1/3)^2 * (c^2 - c + 1/2)^4 ) *  x^2*y\n",
      "+ ( (10368) * (c - 1/2)^2 * (c - 1)^8 * (c^2 - c + 1/3)^4 * (c^2 - c + 1/2)^5 ) *  x*y\n",
      "+ ( (-20736) * (c - 1/2)^2 * (c - 1)^13 * (c^2 - c + 1/3)^4 * (c^2 - c + 1/2)^6 ) *  y\n",
      "+ ( (20736) * (c - 1/2)^2 * (c - 1)^10 * (c^2 - c + 1/3)^4 * (c^2 - c + 1/2)^6 ) *  x^2\n",
      "CPU times: user 13min 35s, sys: 14.3 s, total: 13min 49s\n",
      "Wall time: 13min 53s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'+ (1) * y^4+ ((12) * (c^2 - 4/3*c + 1/2) * (c^2 - c + 1/2)) * x*y^3+ ((144) * (c - 1)^4 * (c^2 - c + 1/2)^2 * (c^4 - 5/2*c^3 + 8/3*c^2 - 49/36*c + 5/18)) * y^3+ ((60) * (c - 1)^2 * (c^2 - 6/5*c + 2/5) * (c^2 - c + 1/2)^2) * x^2*y^2+ ((1008) * (c - 1)^4 * (c^2 - 8/7*c + 5/14) * (c^2 - c + 1/3)^2 * (c^2 - c + 1/2)^3) * x*y^2+ ((1728) * (c - 1)^8 * (c^2 - c + 1/3)^2 * (c^2 - c + 1/2)^4 * (c^4 - 7/2*c^3 + 19/4*c^2 - 11/4*c + 7/12)) * y^2+ ((2592) * (c - 1)^6 * (c^2 - 10/9*c + 1/3) * (c^2 - c + 1/3)^2 * (c^2 - c + 1/2)^4) * x^2*y+ ((10368) * (c - 1/2)^2 * (c - 1)^8 * (c^2 - c + 1/3)^4 * (c^2 - c + 1/2)^5) * x*y+ ((-20736) * (c - 1/2)^2 * (c - 1)^13 * (c^2 - c + 1/3)^4 * (c^2 - c + 1/2)^6) * y+ ((20736) * (c - 1/2)^2 * (c - 1)^10 * (c^2 - c + 1/3)^4 * (c^2 - c + 1/2)^6) * x^2'"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time Weil12_str = CalcWeilFunction12(); \n",
    "Weil12_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 7.0",
   "language": "",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}