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