test

%hide D2 = load(DATA+'D2.sobj') D5 = load(DATA+'D5.sobj') D7 = load(DATA+'D7.sobj') D13 = load(DATA+'D13.sobj') Ds = [ D2, D5, D7, D13 ] def inject_D(D): from sage.misc.misc import inject_variable for k, v in D.items(): inject_variable(k, v) # assert H.base_ring() is k # assert K.order(basis).is_maximal() if False: # remains to give some ideas of which variables are injected H = D['H'] K = D['K'] Kr = D['Kr'] k = D['k'] k_over_Kr = D['k_over_Kr'] k_over_Kr_into_k = D['k_over_Kr_into_k'] k_into_k_over_Kr = D['k_into_k_over_Kr'] k_into_C = D['k_into_C'] # for basis basis = D['basis'] # for eta Kr_into_Ls = D['Kr_into_Ls'] L = D['L'] K_into_L = D['K_into_L'] def my_red(a, P): r"""Reduce a number field element modulo a degree one prime of special form.""" assert a in P.number_field() p, alpha = P.gens() alpha = alpha.polynomial() d = alpha.degree() if d == 1: Kbar = GF(p) alpha_in_Kbar = Kbar(-alpha(0)) else: Kbar = GF(p^d, 'a', modulus=alpha) alpha_in_Kbar = Kbar.gen() pp = a.polynomial() return pp.numerator()(alpha_in_Kbar)/Kbar(pp.denominator()) def my_mod(H, P): r"""Reduce a hyperelliptic curve over a number field modulo a degree one prime of special form.""" L = [ my_red(a, P) for a in H.hyperelliptic_polynomials()[0] ] return HyperellipticCurve( L[0].parent()['x']( L ) ) def my_below(P): x = ZZ['x'].gen() Kr_gen_in_k = k_over_Kr_into_k( k_over_Kr(Kr.gen()) ) r = my_red( Kr_gen_in_k, P ) if r.parent() == 1: pol = x - r Pbelow = Kr.ideal( [ P.gens()[0], pol(Kr.gen()) ] ) else: # pol = r.polynomial().change_ring(ZZ) Pbelow = k_into_k_over_Kr(P).relative_norm() # XXX improvements? return Pbelow def is_large_prime(p): return p > 5*10^3 def is_large_prime_ideal(P): (p, alpha) = P.gens() return is_large_prime(ZZ(p)) or alpha.polynomial().degree() > 1 chi = lambda tt: QQ(DirichletGroup(4).gen()(tt)) def good_primes(k): Ps_deg1 = k.primes_of_degree_one_list(10)[3:6] Ps_large = [] p = 2^30 while p < 2^30 + 2^12 and len(Ps_large) < 3: p = next_probable_prime(p) rts = k.polynomial().roots(GF(p), multiplicities=False) if len(rts) > 0: P = k.ideal([ p, k.gen() - ZZ(rts[0]) ]) Ps_large.append(P) Ps_deg2 = [] for p in prime_range(150, 180): # a non-degree one prime, in two generator form facs = [ fac for fac, e in GF(p)['t'](k.polynomial()).factor() if fac.degree() > 1 ] if facs: gg = choice(facs) P = k.ideal([ p, gg.change_ring(ZZ) ]) Ps_deg2.append( P ) return (Ps_deg1, Ps_deg2, Ps_large) def my_cm_eta(prime_of_Kr): pp = prime_of_Kr # an ideal of Kr eta_pp_in_L = prod( emb(pp) for emb in Kr_into_Ls ) L_over_K = L.relativize(K_into_L, L.variable_name()) L_over_K_into_L, L_into_L_over_K = L_over_K.structure() J = L_over_K.ideal([ L_into_L_over_K(g) for g in eta_pp_in_L.gens() ]) JJ = J.ideal_below() return JJ # -*- sage -*- def _frobenius_coefficient_bounds(self): # Bound from Weil conjectures, via Kedlaya's paper: # | a_i | <= (2g choose i) * q**(i/2) p = self.base_ring().characteristic() # XXX higher degree sqrtp = RR(p).sqrt() g = self.genus() Bs = [] for i in range(2*g + 1): B = 2 * binomial(2*g, i) * sqrtp**i Bs.append(ZZ(ceil(B))) Bs.reverse() return Bs def _frobenius_coefficient_bound(self): assert self.base_ring().is_finite() p = self.base_ring().characteristic() # XXX higher degree Bs = self._frobenius_coefficient_bounds() M = ZZ(max(Bs)) B = M.exact_log(p) if p**B < M: B += 1 assert p**B >= M # B += 2*self.genus() return B # XXX def frobenius_matrix(self, N=None): from sage.schemes.hyperelliptic_curves.hypellfrob import hypellfrob p = self.base_ring().characteristic() f, h = self.hyperelliptic_polynomials() if h != 0: # need y^2 = f(x) raise NotImplementedError sign = 1 if not f.is_monic(): # at this time we need a monic f c = f.leading_coefficient() f = f / c if c.is_square(): # solutions of $y^2 = c * f(x)$ correspond naturally to solutions of # $(sqrt(c) y)^2 = f(x)$ pass else: # we'll count points on a twist and then correct by untwisting... sign = -1 if not f.is_monic(): raise NotImplementedError if N is None: N = self._frobenius_coefficient_bound() matrix_of_frobenius = hypellfrob(p, N, f) matrix_of_frobenius = sign * matrix_of_frobenius return matrix_of_frobenius def cnt_wo_g(f, kbar): if kbar is None: kbar = Hbar.base_ring() cnt = 1 # pt at \infty if f.degree() % 2 == 0: cnt = 2 for x in kbar: c = f(x) if c == 0: cnt += 1 elif c.is_square(): cnt += 2 return cnt def cnt_w_g(f, g, kbar): if kbar is None: kbar = Hbar.base_ring() cnt = 1 # pt at \infty # if f.degree() % 2 == 0: # cnt = 2 y = kbar['y'].gen() for x in kbar: print x p = y**2 + g(x)*y - f(x) cnt += len( p.roots() ) return cnt def cnt(Hbar, kbar=None): if kbar is None: kbar = Hbar.base_ring() f, g = Hbar.hyperelliptic_polynomials() if g == 0: return cnt_wo_g(f, kbar) else: return cnt_w_g(f, g, kbar) # cnt = 1 # pt at \infty # for x in kbar: # c = f(x) # if c == 0: # cnt += 1 # elif c.is_square(): # cnt += 2 # return cnt def frobenius_polynomial_r(Hbar): g = Hbar.genus() kbar = Hbar.base_ring() assert kbar.is_finite() q = kbar.order() e = SFAElementary(QQ) p = SFAPower(QQ) P = [0] * (g+1) P[1] = q + 1 - cnt(Hbar) for i in range(2, g+1): # N_k = q^k + 1 - \sum_{i=1}^{2 g} \alpha_i^k . kbar2 = GF(ZZ(kbar.order())**i, 'a') P[i] = q**i + 1 - cnt(Hbar, kbar2) A = [0] * (2*g+1) A[0] = 1 for i in range(1, g+1): f = p(e([i])) # express coefficient of frob poly in terms of power basis A[i] = sum( coeff*prod(P[mon] for mon in mons) for mons, coeff in f.monomial_coefficients().items()) for i in range(0, g+1): # a_{2 g - i} = q^{g - i} a_i for i = 0, \ldots, g. A[2*g - i] = q**(g - i) * A[i] return QQ['x'](ZZ['x'](A).reverse()(-ZZ['x'].gen())) def frobenius_polynomial(self, N=None): if N is None: N = self._frobenius_coefficient_bound() try: p = self.frobenius_matrix(N=N).charpoly() except ValueError: return frobenius_polynomial_r(self) try: Bs = self._frobenius_coefficient_bounds() cs = [] for c, B in zip(p.coeffs(), Bs): d = ZZ(c) if d > B: d = -ZZ(-c) assert abs(d) <= B cs.append(d) return ZZ['x'](cs) except: return self.frobenius_polynomial(2*N) # otherwise increase precision def frobenius_order(self): f = self.frobenius_polynomial().factor()[0][0] return ZZ.extension(f, names='phi') def frobenius(self): return self.frobenius_order().fraction_field().gen() def _reduce(self, P): return self.change_ring(self.base_ring().residue_field(P)) sage.schemes.hyperelliptic_curves.hyperelliptic_generic.HyperellipticCurve_generic._frobenius_coefficient_bounds = _frobenius_coefficient_bounds sage.schemes.hyperelliptic_curves.hyperelliptic_generic.HyperellipticCurve_generic._frobenius_coefficient_bound = _frobenius_coefficient_bound sage.schemes.hyperelliptic_curves.hyperelliptic_generic.HyperellipticCurve_generic.frobenius_matrix = frobenius_matrix sage.schemes.hyperelliptic_curves.hyperelliptic_generic.HyperellipticCurve_generic.frobenius_polynomial = frobenius_polynomial sage.schemes.hyperelliptic_curves.hyperelliptic_generic.HyperellipticCurve_generic.frobenius_order = frobenius_order sage.schemes.hyperelliptic_curves.hyperelliptic_generic.HyperellipticCurve_generic.frobenius = frobenius sage.schemes.hyperelliptic_curves.hyperelliptic_generic.HyperellipticCurve_generic.mod = _reduce # sage: x = QQ['x'].0 # sage: f = x^7 + 143*x^6 + 41*x^5 + 58*x^4 + 61*x^3 + 52*x^2 + 164*x + 64 # sage: p = 211 # sage: HyperellipticCurve(GF(211)['x']('x^7 + 143*x^6 + 41*x^5 + 58*x^4 + 61*x^3 + 52*x^2 + 164*x + 64')).frobenius_polynomial() def short_weierstrass_model(self, complete=True): r""" Return a Weierstrass model of the form $y^2 = f(x)$. If \code{complete} is \code{True}, depress the resulting polynomial $f(x)$ so that if $d = \deg(f)$ it has no term $x^{d-1}$. This requires that the characteristic of the ground field does not divide $\deg f$. EXAMPLES: sage: x = polygen(QQ) sage: H = HyperellipticCurve(x^7 + x^6 + 2, x); H Hyperelliptic Curve over Rational Field defined by y^2 + x*y = x^7 + x^6 + 2 sage: H.short_weierstrass_model() Hyperelliptic Curve over Rational Field defined by y^2 = x^7 - 3/7*x^5 + 10/49*x^4 - 15/343*x^3 + 2449/9604*x^2 - 5/16807*x + 1647092/823543 sage: Hbar = H.change_ring(GF(2)) sage: Hbar.short_weierstrass_model() Traceback (most recent call last): ... ValueError: short_weierstrass_model(): no short model for Hyperelliptic Curve over Finite Field of size 7 defined by y^2 + x*y = x^7 + x^6 + 2 (characteristic is 2) sage: Hbar = H.change_ring(GF(7)) sage: Hbar.short_weierstrass_model() Traceback (most recent call last): ... ValueError: short_weierstrass_model(): no short model for Hyperelliptic Curve over Finite Field of size 7 defined by y^2 + x*y = x^7 + x^6 + 2 (characteristic is 7) sage: Hbar.short_weierstrass_model(complete=False) Hyperelliptic Curve over Finite Field of size 7 defined by y^2 = x^7 + x^6 + 2*x^2 + 2 """ f, h = self.hyperelliptic_polynomials() K = self.base_ring() if K.characteristic() == 2: raise ValueError, "short_weierstrass_model(): no short model for %s (characteristic is %s)" % (self, K.characteristic()) new_f = f + h*h/4 # always complete the square if complete: d = ZZ(f.degree()) if ZZ(K.characteristic()).divides(d): raise ValueError, "short_weierstrass_model(): no short model for %s (characteristic is %s)" % (self, K.characteristic()) x = f.parent().gen() new_f = new_f(x - (f[d-1]/d)) # complete the $d - 1$ power return HyperellipticCurve(new_f) def quadratic_twist(self, d=None): r""" XXX characteristics, d == 0 """ if d is None: d = self.base_ring().quadratic_nonresidue() f, _ = self.short_weierstrass_model().hyperelliptic_polynomials() n = f.degree() fp = f.parent()([ f[i]*d**(n-i) for i in range(0, n+1) ]) return HyperellipticCurve(fp, 0) def apply_morphism(self, phi): f, h = self.hyperelliptic_polynomials() P = f.parent().change_ring(phi.codomain()) fp = P([ phi(c) for c in f.coeffs() ]) hp = P([ phi(c) for c in h.coeffs() ]) return HyperellipticCurve(fp, hp) sage.schemes.hyperelliptic_curves.hyperelliptic_generic.HyperellipticCurve_generic.short_weierstrass_model = short_weierstrass_model sage.schemes.hyperelliptic_curves.hyperelliptic_generic.HyperellipticCurve_generic.quadratic_twist = quadratic_twist sage.schemes.hyperelliptic_curves.hyperelliptic_generic.HyperellipticCurve_generic.apply_morphism = apply_morphism def random_points_on_curve(C, numpts=10): pts = [] K = C.base_ring() while len(pts) < numpts: x = K.random_element() try: pt = C.lift_x(x) except ValueError: continue pts.append(pt) return pts def random_points_on_jacobian(J, numpts=10): C = J.curve() pts = random_points_on_curve(C) g = C.genus() divs = [] for ptsg in Subsets(pts, g): div = sum(J(pt) for pt in list(ptsg)) divs.append(div) if len(divs) > numpts: break return divs def random_points_on_jacobian_of_curve(C, numpts=10): J = C.jacobian()(C.base_ring()) return random_points_on_jacobian(J, numpts) def is_multiple_of_order(self, N, numpts=10): r""" Return True if N is a multiple of the order of this curve. """ N = ZZ(N) pts = random_points_on_jacobian(self, numpts) return all((N*pt).is_zero() for pt in pts) sage.schemes.hyperelliptic_curves.jacobian_homset.JacobianHomset_divisor_classes.is_multiple_of_order = is_multiple_of_order def is_multiple_of_order_over(self, R, N, numpts=3): return self.change_ring(R).jacobian().point_set(R).is_multiple_of_order(N, numpts=numpts) def is_multiple_of_order_2(self, N, numpts=3): R = self.base_ring() assert R.is_finite() return self.jacobian().point_set(R).is_multiple_of_order(N, numpts=numpts) sage.schemes.hyperelliptic_curves.hyperelliptic_generic.HyperellipticCurve_generic.is_multiple_of_order_over = is_multiple_of_order_over sage.schemes.hyperelliptic_curves.hyperelliptic_generic.HyperellipticCurve_generic.is_multiple_of_order = is_multiple_of_order_2 def magma_init(self, magma): """ Internal function. Returns a string to initialize this elliptic curve in the Magma subsystem. EXAMPLES: sage: E = EllipticCurve(QQ, [1, 1]); E Elliptic Curve defined by y^2 = x^3 + x +1 over Rational Field sage: E._magma_init_(magma) # optional - magma 'EllipticCurve([0,0,0,1,1])' sage: magma(E) # optional - magma Elliptic Curve defined by y^2 = x^3 + x + 1 over Rational Field sage: magma(E)._sage_() # optional - magma Elliptic Curve defined by y^2 = x^3 + x +1 over Rational Field sage: E == magma(E)._sage_() # optional - magma True sage: H = HyperellipticCurve(QQ['x']('x^5 + 1'), QQ['x']('x')); H Hyperelliptic Curve over Rational Field defined by y^2 = x^5 + 1 sage: magma(H) # optional - magma Hyperelliptic Curve defined by y^2 + x*y = x^5 + 1 over Rational Field sage: magma(H)._sage_() # optional - magma Hyperelliptic Curve over Rational Field defined by y^2 + x*y = x^5 + 1 sage: H == magma(H)._sage_() # optional - magma True """ return 'HyperellipticCurve(%s)' % ', '.join([x._magma_init_(magma) for x in self.hyperelliptic_polynomials()]) sage.schemes.hyperelliptic_curves.hyperelliptic_generic.HyperellipticCurve_generic._magma_init_ = magma_init def items_for(list_of_dicts): L = [] for D in list_of_dicts: K = D['K'] L.append( (D, 'CM by %s (%s, class number %s)' % (K.polynomial(), 'Cyclic Galois' if K.is_galois() else 'Dihedral non-Galois', K.class_number())) ) return L def cap(x, bnd=1e-10): if isinstance(x.parent(), sage.rings.polynomial.multi_polynomial_ring_generic.MPolynomialRing_generic): s = 0 for m in x.monomials(): c = cap(x.monomial_coefficient(m), bnd=bnd) if c != 0: s += m*c return s if isinstance(x.parent(), sage.rings.polynomial.polynomial_ring.PolynomialRing_commutative): return x.parent()( [ cap(c, bnd=bnd) for c in x ] ) if x.parent() is CDF: r, i = x.real(), x.imag() elif hasattr(x, 'apply_map'): return x.apply_map(lambda u: cap(u, bnd=bnd)) else: r, i = list(x) if abs(r) < bnd: r = 0 if abs(i) < bnd: i = 0 if x.parent() is CDF: return CDF(r, i) return x.parent()([r, i]) def test_one(P, Pbelow): if True: Hp = my_mod(H, P) etaPbelow = my_cm_eta(Pbelow) (fr, ) = etaPbelow.gens_reduced() lamb = choice([-1, 1])*fr # little extra randomness if not is_large_prime_ideal(P): f1 = Hp.frobenius_polynomial() else: f1 = lamb.charpoly() if not Hp.is_multiple_of_order( f1(1) ): f1 = f1(-f1.parent().gen()) assert not Hp.is_multiple_of_order( f1(-1) ) assert Hp.is_multiple_of_order( f1(1) ) # print f1.factor() # check that my theorem is working :) v = matrix(map(vector, basis)).solve_left(vector(lamb)) a, b, c, d = v epsilon = chi((c + d) % 4) f2 = (epsilon * lamb).charpoly() good = (f2 == f1) if good: html(r"<table style='background-color:LightGreen'>") else: html(r"<table style='background-color:LightPink'>") if True: html("<tr>") html(r"<th>$\text{degree of }\mathcal{P}$</th>") html(r"<td>$%s$</td>" % P.gens()[1].polynomial().degree()) html("</tr>") html("<tr>") AN = P.absolute_norm() html(r"<th>$\mathrm{N}_{k/\mathbf{Q}}(\mathcal{P})$</th>") html('<td>$%s = %s$</td>' % (latex(AN), latex(AN.factor()))) html("</tr>") html("<tr>") html(r"<th>$\mathrm{N}_{k/K'}(\mathcal{P})$</th>") html('<td>$%s$</td>' % latex(Pbelow)) html("</tr>") html("<tr>") html(r"<th>$\lambda = \eta(\mathrm{N}_{k/K'}(\mathcal{P}))$</th>") html(r"<td>$%s</td>$" % latex(K.ideal(lamb))) html("</tr>") html("<tr>") html(r"<th>$\lambda = a z + b z w' + c w + d$</th>") html(r"<td>$(%s) \cdot z + (%s) \cdot z w' + (%s) \cdot w + (%s)$</td>" % (a, b, c, d)) html("</tr>") html("<tr>") html(r"<th>$\epsilon$</th>") html(r"<td>$%s$</td>" % epsilon) html("</tr>") html("<tr>") html(r"<th>$N = \mathrm{N}_{K/\mathbf{Q}}(1 - \epsilon \cdot \lambda)$</th>") html(r"<td>$%s$</td>" % f2(1)) html("</tr>") html("<tr>") html(r"<th>$H \pmod{\mathcal{P}}$</th>") html(r"<td>$%s$</td>" % latex(Hp.hyperelliptic_polynomials()[0])) html("</tr>") html("<tr>") html(r"<th>$\# \mathrm{Jac}(H)(O_k/\mathcal{P})$</th>") html(r"<td>$%s$</td>" % f1(1)) html("</tr>") html(r"</table>") def show_test_summary(): html("<table>") html("<tr><td>") html("$K = %s$ with $%s \mapsto %s \in \mathbf{C}$" % (latex(K), latex(K.gen()), latex(cap(CDF(K_into_QQbar(K.gen())))))) html("</td></tr>") html("<tr><td>") K0 = K.subfields(2)[0][0] html("$K_0 = %s$ with $D_{K_0} = %s$" % (latex(K0), K0.disc())) html("</td></tr>") html("<tr><td>") html("$O_K = (%s) O_{K_0} + O_{K_0}$" % latex(basis[0])) html("</td></tr>") html("<tr><td>") html("this choice of $z$ leads to $z_0 \sim %s \in \mathcal{H}_2$" % latex(z_0.change_ring(CDF).apply_map(cap))) html("</td></tr>") html("<tr><td>") html("$k = %s$" % latex(k.change_names('b'))) html("</td></tr>") html("<tr><td>") html("with $%s \mapsto %s \in \mathbf{C}$" % (latex(k.change_names('b').gen()), latex(cap(CDF(k_into_QQbar(k.gen())))))) html("</td></tr>") html("<tr><td>") html("$[k : \mathbf{Q}] = %s$ and $[k : K'] = %s$" % (latex(k.degree()), latex(k.degree() // Kr.degree()))) html("</td></tr>") html("<tr><td>") html("$H \sim %s$" % latex(CDF['x'](map(cap, map(CDF, map(k_into_QQbar, H.hyperelliptic_polynomials()[0])))))) html("</td></tr>") @interact def _(example_dict = selector(items_for(Ds), label="Example"), primes = input_box(default = '', label="Rational primes equal or after", type=str), degrees = input_box(default = '1', label="of degrees", type=str), ): inject_D(example_dict) degrees = set(map(abs, vector(ZZ, degrees.split(',')))) Ps1, Ps2, Psl = good_primes(k) html("<p>Some rational primes below small degree one primes are:") html(', '.join([ "%s" % P.gens()[0] for P in Ps1 ]) ) html("</p>") html("<p>Some rational primes below small degree two or larger primes are:") html(', '.join([ "%s" % P.absolute_norm().factor()[0][0] for P in Ps2 ]) ) html("</p>") html("<p>Some rational primes below large degree one primes are:") html(', '.join([ "%s" % P.gens()[0] for P in Psl ]) ) html("</p>") html("<table style='border:1px solid black'>") try: ps = vector(ZZ, primes.split(',')) except: ps = [] for p in ps: if not p: continue if not is_pseudoprime(p): p = next_probable_prime(p) html(r"<tr>") html(r"<td>p = $%s$</td>" % p) # a prime, in two generator form facs = [ fac for fac, e in GF(p)['t'](k.polynomial()).factor() if fac.degree() in degrees ] if not facs: html(r"<td style='background-color:red'>There is no prime above $%s$ with degree in $%s$</td>" % (p, latex(vector(sorted(degrees))).replace("(", r"\{").replace(")", r"\}"))) html("</tr>") # close row continue afac = choice(facs) P = k.ideal([ p, afac.change_ring(ZZ) ]) try: Pbelow = my_below(P) except ZeroDivisionError: html(r"<td style='background-color:red'><td>Please try a prime larger than %s</td>" % p) html("</tr>") # close row continue html("<td>") test_one(P, Pbelow) html("</td>") html("</tr>") # close row html('</table>') show_test_summary()