class WeilObject(object):
# Konstruktor: Gibt elliptische Kurve mit vorberechneter Weil-Funktion zurück.
# @param n .. Natürliche Zahl aus dem Bereich {4,...,10,12}
# self._E = Elliptische Kurve X_1(n) in den Variablen x, y und Parameter c
# self._P = n-Torsionspunkt (0,0) auf E
# self._n = natürliche Zahl n
# self._f = Weil-Funktion E -> C mit Divisor div(f)=n{P} - n[OO]
# self._s = Wohlformatierter String der Weil-Funktion f
#
def __init__(self, n):
self._n = n
x, y, c = var('x y c');
if n == 4:
self._E = EllipticCurve([1, c, c, 0, 0]); # X_1(4)
self._P = self._E(0,0); # 4-Torsionspunkt (0,0)
self._f = x^2 - y; # TODO: http://trac.sagemath.org/5930
self._s = "x^2 - y";
return;
if n == 5:
self._E = EllipticCurve([c+1, c, c, 0, 0]);
self._P = self._E(0,0);
self._f = -x^2 + x*y + y;
self._s = "-x^2 + x*y + y";
return;
if n == 6:
self._E = EllipticCurve([1-c, -c*(c+1), -c*(c+1), 0, 0]);
self._P = self._E(0,0);
self._f = y^2 - (c + 1)*x*y - (c + 1)^2*y + (c + 1)^2*x^2;
self._s = "y^2 - (c + 1)*x*y - (c + 1)^2*y + (c + 1)^2*x^2";
return;
if n == 7:
self._E = EllipticCurve([1-c-c^2, c^2*(c+1), c^2*(c+1), 0, 0]);
self._P = self._E(0,0);
self._f = (c - 1)*y^2 + x^2*y - c^3*x*y + c^4*y - c^4*x^2;
self._s = "(c - 1)*y^2 + x^2*y - c^3*x*y + c^4*y - c^4*x^2";
return;
if n == 8:
self._E = EllipticCurve([1-2*c^2, -c*(2*c+1)*(c+1)^2, -c*(2*c+1)*(c+1)^3, 0, 0]);
self._P = self._E(0,0);
self._f = (x*y^2 + (2*c + 3)*(c + 1)^3*y^2 - 2*(c + 1)^2*x^2*y - (2*c + 1)^2*(c + 1)^4*x*y
- (2*c + 1)^2*(c + 1)^7*y + (2*c + 1)^2*(c + 1)^6*x^2);
self._s = ("x*y^2 \n+ (2*c + 3)*(c + 1)^3*y^2 \n- 2*(c + 1)^2*x^2*y \n- (2*c + 1)^2*(c + 1)^4*x*y \n" +
"- (2*c + 1)^2*(c + 1)^7*y \n+ (2*c + 1)^2*(c + 1)^6*x^2");
return;
if n == 9:
self._E = EllipticCurve([c^3 + c^2 + 1, c^2 * (c+1) * (c^2+c+1),
c^2 * (c+1) * (c^2+c+1), 0, 0]);
self._P = self._E(0,0);
self._f = (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)*(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 - c^4*(c^2 + c + 1)^4*x^2);
self._s = ("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)*(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- c^4*(c^2 + c + 1)^4*x^2");
return;
if n == 10:
self._E = 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]);
self._P = self._E(0,0);
self._f = (2*(c^2 - 2*c - 2)*y^3 + x^2*y^2 - (2*c + 1)*c^4*x*y^2
+ (c^3 + 16*c^2 + 22*c + 8)*c^6*y^2 - (3*c + 2)*c^6*x^2*y
+ (c + 1)^2*c^10*x*y - (c + 1)^2*(c^2 + 6*c + 4)*c^12*y
+ (c + 1)^2*c^12*x^2);
self._s = ("2*(c^2 - 2*c - 2)*y^3 \n+ x^2*y^2 \n- (2*c + 1)*c^4*x*y^2 \n" +
"+ (c^3 + 16*c^2 \n+ 22*c + 8)*c^6*y^2 \n- (3*c + 2)*c^6*x^2*y \n" +
"+ (c + 1)^2*c^10*x*y \n- (c + 1)^2*(c^2 + 6*c + 4)*c^12*y \n" +
"+ (c + 1)^2*c^12*x^2");
return;
if n == 12:
self._E = 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]);
self._P = self._E(0,0);
self._f = (y^4
+ (6*c^2 - 8*c + 3)*(2*c^2 - 2*c + 1)*x*y^3
+ (c - 1)^4*(2*c^2 - 2*c + 1)^2*(36*c^4 - 90*c^3 + 96*c^2 - 49*c + 10)*y^3
+ 3*(c - 1)^2*(5*c^2 - 6*c + 2)*(2*c^2 - 2*c + 1)^2*x^2*y^2
+ (c - 1)^4*(14*c^2 - 16*c + 5)*(3*c^2 - 3*c + 1)^2*(2*c^2 - 2*c + 1)^3*x*y^2
+ (c - 1)^8*(3*c^2 - 3*c + 1)^2*(2*c^2 - 2*c + 1)^4*(12*c^4 - 42*c^3 + 57*c^2 - 33*c + 7)*y^2
+ 2*(c - 1)^6*(9*c^2 - 10*c + 3)*(3*c^2 - 3*c + 1)^2*(2*c^2 - 2*c + 1)^4*x^2*y
+ (2*c - 1)^2*(c - 1)^8*(3*c^2 - 3*c + 1)^4*(2*c^2 - 2*c + 1)^5*x*y
- (2*c - 1)^2*(c - 1)^13*(3*c^2 - 3*c + 1)^4*(2*c^2 - 2*c + 1)^6*y
+ (2*c - 1)^2*(c - 1)^10*(3*c^2 - 3*c + 1)^4*(2*c^2 - 2*c + 1)^6*x^2);
self._s = ("y^4 \n" +
"+ (6*c^2 - 8*c + 3)*(2*c^2 - 2*c + 1)*x*y^3 \n" +
"+ (c - 1)^4*(2*c^2 - 2*c + 1)^2*(36*c^4 - 90*c^3 + 96*c^2 - 49*c + 10)*y^3 \n" +
"+ 3*(c - 1)^2*(5*c^2 - 6*c + 2)*(2*c^2 - 2*c + 1)^2*x^2*y^2 \n" +
"+ (c - 1)^4*(14*c^2 - 16*c + 5)*(3*c^2 - 3*c + 1)^2*(2*c^2 - 2*c + 1)^3*x*y^2 \n" +
"+ (c - 1)^8*(3*c^2 - 3*c + 1)^2*(2*c^2 - 2*c + 1)^4*(12*c^4 - 42*c^3 + 57*c^2 - 33*c + 7)*y^2 \n" +
"+ 2*(c - 1)^6*(9*c^2 - 10*c + 3)*(3*c^2 - 3*c + 1)^2*(2*c^2 - 2*c + 1)^4*x^2*y \n" +
"+ (2*c - 1)^2*(c - 1)^8*(3*c^2 - 3*c + 1)^4*(2*c^2 - 2*c + 1)^5*x*y \n" +
"- (2*c - 1)^2*(c - 1)^13*(3*c^2 - 3*c + 1)^4*(2*c^2 - 2*c + 1)^6*y \n" +
"+ (2*c - 1)^2*(c - 1)^10*(3*c^2 - 3*c + 1)^4*(2*c^2 - 2*c + 1)^6*x^2");
return;
# Assert im Konstruktor erzeugen, damit kein Objekt für nicht unterstützte n-Werte erzeugt wird
assert false;
# Gibt die String-Repräsentation des Objektes zurück.
def __repr__(self):
str = "Weil function info object of elliptic curve with %i-torsion point (0,0):\n- "%(self._n)
str = str + self._E.__repr__() + "\n"
str = str + "- %i-torsion point P = "%(self._n) + self._P.__repr__() + "\n"
str = str + "- Weil function f(x,y) = " + self._f.__repr__() + "\n"
return str
# Gibt die Weil-Funktion als wohlformatierten String auf der Konsole aus.
def PrintNiceString(self):
print self._s
# Gibt den Wert der Weilfunktion f im Punkt (xval, yval) zurück.
def f(self, xval, yval):
return self._f(x = xval, y = yval)
# Gibt f(-(x,y)) zurück, wobei mit -(x,y) das Inverse bzgl. der Addition auf der elliptischen Kurve gemeint ist.
def GetMinusWeilFunction(self):
a1, a2, a3, a4, a6 = self._E.a_invariants();
return self.f(x, -a1 * x - a3 - y);
# Reduziert den Ausdruck <term> im Quotienten-Ring modulo der Gleichung der elliptischen Kurve.
# @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(self, term):
# Koeffizienten der elliptischen Kurve bestimmen
a1, a2, a3, a4, a6 = self._E.a_invariants();
# Polynomring, Ideal und Quotientenring bestimmen
Q.<x,y,c> = QQ['x','y','c'];
J = Q.ideal(y^2 + a1*x*y + a3*y - x^3 - a2*x^2 - a4*x - a6);
R.<X,Y,C> = QuotientRing(Q, J);
# In <term> 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});
# Statische Methode.
# Ersetzt in <term> alle Vorkommnisse von <search> durch (<repl>).
# Unterstützt werden die Variablen x, y und c.
def ReplaceSubTerm(self, 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})
# Drückt x^m für m >= 3 durch kleinere x-Potenzen aus.
def XPower(self, m):
a1, a2, a3, a4, a6 = self._E.a_invariants();
if m <= 2:
return x^m
if m >= 3:
return expand(x^(m-3) * (y^2 + a1*x*y + a3*y - a2*x^2 - a4*x - a6))
# Ersetzt in <term> alle Vorkommnisse von <x^m> durch Potenzen x^2, x, 1.
# Parameter term: Polynom in x, y und c.
def ReplaceHighXMonoms(self, term):
# Polynomring als "Monom-Wörterbuch"
R.<x,y,c> = sage.rings.polynomial.multi_polynomial_ring.MPolynomialRing_polydict(QQ, 3, order='lex');
# term als Polynom in x, y, c interpreteren
poly = R(term)
# Maximalgrad von x im Polynom <term> bestimmen.
max_x_degree = poly.degree(x)
# Iteration von max_x_degree bis 3
for e in range(max_x_degree, 2, -1):
# Ersetze x^e durch kleinere Potenzen.
term = self.ReplaceSubTerm(term, x^e, self.XPower(e));
# Ausmultiplizieren, damit im nächsten Schritt x^(e-1) ersetzt werden kann.
term = expand(term);
return term;
# Drückt y^m für m >= 2 durch kleinere y-Potenzen aus, FALLS x = 1 gilt.
def YPowerX1(self, m):
a1, a2, a3, a4, a6 = self._E.a_invariants();
if m < 2:
return y^m
if m >= 2:
return expand(y^(m-2) * ((-a1-a3)* y + a2 + a4 + a6 + 1)) # Achtung: x ist 1 !!!
# Ersetzt in <term> alle Vorkommnisse von <y^m> durch Potenzen y^2, y, 1,
# FALLS x = 1 gilt.
# Parameter term: Polynom in y und c.
def ReplaceHighYMonomsX1(self, term):
# Polynomring als "Monom-Wörterbuch"
R.<x,y,c> = sage.rings.polynomial.multi_polynomial_ring.MPolynomialRing_polydict(QQ, 3, order='lex');
# term als Polynom in y, c interpreteren
poly = R(term)
# Maximalgrad von y im Polynom <term> bestimmen.
max_y_degree = poly.degree(y)
# Iteration von max_y_degree bis 2
for e in range(max_y_degree, 1, -1):
# Ersetze y^e durch kleinere Potenzen.
term = self.ReplaceSubTerm(term, y^e, self.YPowerX1(e));
# Ausmultiplizieren, damit im nächsten Schritt y^(e-1) ersetzt werden kann.
term = expand(term);
return term;
# Bestimmt die komplexe Zahl alpha^n, die die folgende Eigenschaft besitzt:
# f((x,y)) * f(-(x,y)) = alpha^n * x^n, wobei mit -(x,y) das Inverse bzgl. der Addition
# auf der elliptischen Kurve gemeint ist.
# -> Wenn alles korrekt ist, muss alpha = (-1)^n gelten.
#
def CalcAlpha(self):
# Minus-Funktion g der Weil-Funktion berechnen, g(Q) = f(-Q) für alle Q in E
g(x, y) = self.GetMinusWeilFunction()
# Produkt aus Weil-Funktion und Minus-Funktion bestimmen
p(x, y) = self.f(x, y) * g(x, y)
# alpha^n berechnen, indem spezieller Wert x = 1 eingesetzt wird
alpha = p(x = 1, y = y)
# ausmultiplizieren
alpha = expand(alpha)
# Potenzen y^m für m >= 2 ersetzen durch Linearkombination aus 1 und y;
# die Variable y muss aus dem Term komplett herausfallen, so dass eine Konstante übrig bleibt
alpha = self.ReplaceHighYMonomsX1(alpha)
# Den Wert (alpha^n) zurückgeben.
return alpha
# Überprüft die Integrität des Objektes:
# - Ist P wirklich ein n-Torsionspunkt von E ?
# - Gilt wirklich alpha = (-1)^n ?
# - Gilt wirklich f((x,y)) * f(-(x,y)) = alpha * x^n ?
def Verify(self):
# 1. Ist P wirklich ein n-Torsionspunkt von E ?
assert (self._n * self._P == self._E(0));
# 2. Gilt wirklich alpha^n = (-1)^n ?
alpha = self.CalcAlpha();
assert (alpha == (-1)^self._n);
# 3. Gilt wirklich f((x,y)) * f(-(x,y)) = alpha^n * x^n ?
# 3a. Minus-Funktion g der Weil-Funktion berechnen, g(Q) = f(-Q) für alle Q in E
g(x, y) = self.GetMinusWeilFunction()
# 3b. Produkt aus Weil-Funktion und Minus-Funktion bestimmen
p(x, y) = self.f(x, y) * g(x, y)
#3c. Differenz p - (alpha^n * x^n) im Ideal reduzieren
diff = p(x,y)-(alpha * x^self._n)
diff = self.ReduceInQuotientRing(diff)
#3d. Differenz muss Null sein
assert (diff == 0);
# 4. Alles OK.
return true;