| 1 | r""" |
| 2 | Reduction of binary forms and hyperelliptic curve equations |
| 3 | |
| 4 | This file contains functions that reduce the coefficients of a binary form or |
| 5 | hyperelliptic curve by integral transformations, leaving the discriminant |
| 6 | intact. |
| 7 | |
| 8 | AUTHORS: |
| 9 | |
| 10 | - Florian Bouyer (2011 - 09) : covariant and stoll_cremona_reduction |
| 11 | - Marco Streng (2012 - 09) : Hilbert fundamental domains in some cases |
| 12 | |
| 13 | """ |
| 14 | #***************************************************************************** |
| 15 | # Copyright (C) 2011,2012,2013 Florian Bouyer <f.j.s.c.bouyer@gmail.com> |
| 16 | # and Marco Streng <marco.streng@gmail.com> |
| 17 | # |
| 18 | # Distributed under the terms of the GNU General Public License (GPL) |
| 19 | # as published by the Free Software Foundation; either version 2 of |
| 20 | # the License, or (at your option) any later version. |
| 21 | # http://www.gnu.org/licenses/ |
| 22 | #***************************************************************************** |
| 23 | |
| 24 | from sage.structure.sage_object import SageObject |
| 25 | from sage.matrix.all import Matrix |
| 26 | from sage.schemes.plane_conics.constructor import Conic |
| 27 | from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing |
| 28 | from sage.rings.arith import gcd, primes_first_n, factor |
| 29 | from sage.rings.all import (ZZ, QQ, CC, ComplexField, RealField, |
| 30 | is_NumberField, FiniteField) |
| 31 | from sage.schemes.hyperelliptic_curves.constructor import HyperellipticCurve |
| 32 | from sage.misc.misc import prod |
| 33 | from sage.matrix.constructor import identity_matrix |
| 34 | |
| 35 | def stoll_cremona_reduction(f, degree=None, algorithm='default', |
| 36 | precision=53, group='gl', transformation=False): |
| 37 | r""" |
| 38 | Reduces the polynomial f using the Stoll-Cremona method |
| 39 | |
| 40 | INPUT: |
| 41 | |
| 42 | - ``f`` - a univariate polynomial |
| 43 | |
| 44 | - ``degree`` - a positive integer (default: None). Interpret ``f`` as a |
| 45 | homogeneous polynomial of degree ``degree`` in two variables, i.e., |
| 46 | as `F(X,Y) = Y^{degree} f(X/Y)`. If None, use the degree of `f`. |
| 47 | |
| 48 | - ``algorithm`` - ('default' or 'magma'). If set to 'magma', instead of |
| 49 | using :func:`covariant_z0` to find the covariant `z_0`, the algorithm |
| 50 | uses Magma to compute the 'better' covariant `z` (see [SC2002] for more |
| 51 | details) (requires Magma to be installed) |
| 52 | |
| 53 | - ``precision`` - a positive integer (default=53). The precision to use |
| 54 | for the covariants. |
| 55 | |
| 56 | - ``group`` (default='sl') -- a string: 'sl' for SL_2(O_K), 'gl' for |
| 57 | GL_2(O_K), 'plus' for the subgroup GL_2(O_K)^+ of elements of totally |
| 58 | positive determinant in GL_2(O_K). |
| 59 | |
| 60 | - ``transformation`` (default: False). Whether to also output a matrix |
| 61 | (a,b;c,d) such that the output g satisfies |
| 62 | g(z) = (cz+d)^degree * f((az+b)/(cz+d)). |
| 63 | |
| 64 | OUTPUT: |
| 65 | |
| 66 | A polynomial whose associated covariant point in the upper half plane |
| 67 | is in a fundamental domain and that is `SL_2(ZZ)`-equivalent to f |
| 68 | |
| 69 | EXAMPLES: |
| 70 | |
| 71 | We check with Stoll and Cremona's examples:: |
| 72 | |
| 73 | sage: from sage.schemes.hyperelliptic_curves.stoll_cremona import stoll_cremona_reduction |
| 74 | sage: P.<x> = QQ[] |
| 75 | sage: f = 19*x^8 - 262*x^7 + 1507*x^6 - 4784*x^5 + 9202*x^4 - 10962*x^3 + 7844*x^2 - 3040*x + 475 |
| 76 | sage: stoll_cremona_reduction(f) |
| 77 | -x^8 + 6*x^7 - 7*x^6 - 12*x^5 + 27*x^4 - 4*x^3 - 19*x^2 + 10*x - 5 |
| 78 | |
| 79 | sage: Q.<y> = QQ[] |
| 80 | sage: f = y^6 + 30*y^5 + 371*y^4 + 2422*y^3 + 8813*y^2 + 16968*y + 13524 |
| 81 | sage: stoll_cremona_reduction(f) |
| 82 | y^6 - 4*y^4 + 2*y^3 + 8*y^2 - 12*y + 9 |
| 83 | |
| 84 | An example over a real quadratic number field:: |
| 85 | |
| 86 | sage: K.<a> = NumberField(x^2-5) |
| 87 | sage: P.<x> = K[] |
| 88 | sage: f = (x+a)^6+7; f |
| 89 | x^6 + 6*a*x^5 + 75*x^4 + 100*a*x^3 + 375*x^2 + 150*a*x + 132 |
| 90 | sage: stoll_cremona_reduction(f) |
| 91 | x^6 + 7 |
| 92 | |
| 93 | The chosen degree is relevant, in the following example, the Igusa-Clebsch |
| 94 | invariants show that the reduction for degree 5 is not equivalent for |
| 95 | degree 6:: |
| 96 | |
| 97 | sage: Q.<x> = QQ[] |
| 98 | sage: f = 4*x^5 + 26*x^4 + 88*x^3 + 144*x^2 + 112*x + 32 |
| 99 | sage: C = HyperellipticCurve(f) |
| 100 | sage: i = C.igusa_clebsch_invariants() |
| 101 | sage: g = stoll_cremona_reduction(f); g |
| 102 | 2*x^5 + 4*x^4 + 4*x^3 + 24*x^2 - 6*x + 4 |
| 103 | |
| 104 | sage: D = HyperellipticCurve(g) |
| 105 | sage: j = D.igusa_clebsch_invariants() |
| 106 | sage: i[0]^5/i[3] == j[0]^5/j[3] |
| 107 | False |
| 108 | |
| 109 | sage: stoll_cremona_reduction(f, degree=6) |
| 110 | Traceback (most recent call last): |
| 111 | ... |
| 112 | NotImplementedError: Please implement covariant_z0 for binary forms with roots at infinity |
| 113 | |
| 114 | sage: h = stoll_cremona_reduction(f, degree=6, algorithm='magma'); h # optional - magma |
| 115 | 4*x^5 + 6*x^4 + 24*x^3 - 4*x^2 + 4*x - 2 |
| 116 | sage: D = HyperellipticCurve(h) # optional - magma |
| 117 | sage: j = D.igusa_clebsch_invariants() # optional - magma |
| 118 | sage: i[0]^5/i[3] == j[0]^5/j[3] |
| 119 | True |
| 120 | |
| 121 | |
| 122 | With Magma, we can get a better reduction, as the covariant `z` is |
| 123 | implemented there, rather than just the covariant `z0`:: |
| 124 | |
| 125 | sage: P.<x> = QQ[] |
| 126 | sage: f = good example, maybe from [SC2002] |
| 127 | sage: stoll_cremona_reduction(f) |
| 128 | |
| 129 | sage: stoll_cremona_reduction(f, algorithm='magma') # optional - magma |
| 130 | |
| 131 | An example with ``transformation=True``:: |
| 132 | |
| 133 | sage: g = 32*x^6 + 688*x^5 + 6144*x^4 + 29176*x^3 + 77714*x^2 + 110104*x + 64830 |
| 134 | sage: h, T = stoll_cremona_reduction(g, algorithm='magma', transformation=True) |
| 135 | sage: T |
| 136 | [ 3 4] |
| 137 | [-1 -1] |
| 138 | sage: g((3*x+4)/(-x-1))*(-x-1)^6 == h |
| 139 | True |
| 140 | |
| 141 | ALGORITHM: |
| 142 | |
| 143 | This is the algorithm from [SC2002]. |
| 144 | |
| 145 | If the base field is `QQ`, then the algorithm is exactly the same. For |
| 146 | real quadratic fields, we implemented a reduction in the Hilbert upper |
| 147 | half space HH^2. If the real quadratic field has class number one, then |
| 148 | this region is the standard fundamental domain. For more general number |
| 149 | fields, a more general reduction algorithm for the covariant point is |
| 150 | needed. |
| 151 | |
| 152 | REFERENCES |
| 153 | |
| 154 | [SC2002] M. Stoll and J. Cremona "On the reduction theory of binary forms" |
| 155 | J. Reine Angew. Math. 565 (2003) pp. 79-99 |
| 156 | """ |
| 157 | # s = [] |
| 158 | K = f.base_ring() |
| 159 | P = f.parent() |
| 160 | x = P.gen() |
| 161 | h = _HilbertFundDomData(K, embs=precision, group=group) |
| 162 | if degree is None: |
| 163 | degree = f.degree() |
| 164 | R = RealField(precision) |
| 165 | C = ComplexField(precision) |
| 166 | T = Matrix([[1,0],[0,1]]) |
| 167 | previous_height = R(-2) |
| 168 | height = R(-1) |
| 169 | one_minus_epsilon = R(1) - R(2)**(-5) |
| 170 | if algorithm == 'default': |
| 171 | while previous_height < height * one_minus_epsilon: |
| 172 | previous_height = height |
| 173 | z = covariant_z0(f, precision=precision, degree=degree) |
| 174 | M, z = h.step(z) |
| 175 | [[a,b],[c,d]] = M |
| 176 | f = P(f((d*x-b)/(-c*x+a)) * (-c*x+a)**degree) |
| 177 | # s.append((f,z,covariant_z0(f))) |
| 178 | T = T * Matrix([[d, -b], [-c, a]]) |
| 179 | height = prod([t.imag() for t in covariant_z0(f)]).abs() |
| 180 | print f |
| 181 | print height |
| 182 | print previous_height |
| 183 | if height < previous_height * one_minus_epsilon: |
| 184 | raise RuntimeError, "Possibly too low precision in stoll_cremona_reduction" |
| 185 | elif algorithm == 'magma': |
| 186 | from sage.interfaces.all import magma |
| 187 | z = 0; zp = 1 |
| 188 | while z != zp: |
| 189 | Phi = K.embeddings(R) |
| 190 | z = [] |
| 191 | for emb in Phi: |
| 192 | gcoeffs = [emb(f[k]) for k in range(degree+1)] |
| 193 | P = R['x'] |
| 194 | x = P.gen() |
| 195 | g = sum([gcoeffs[indx]*x**indx for indx in range(degree+1)]) |
| 196 | zi = magma(magma.Covariant(degree, magma(g))).sage() |
| 197 | z.append(C(zi)) |
| 198 | M,zp = h.step(z) |
| 199 | [a,b],[c,d] = M |
| 200 | P = K['x'] |
| 201 | x = P.gen() |
| 202 | f = f.parent()(f((d*x-b)/(-c*x+a)) * (-c*x+a)**degree) |
| 203 | T = T * Matrix([[d, -b], [-c, a]]) |
| 204 | else: |
| 205 | raise ValueError("Unkown algorithm '%s'" % algorithm) |
| 206 | if transformation: |
| 207 | return f, T |
| 208 | return f |
| 209 | |
| 210 | |
| 211 | def covariant_z0(f, precision=53, degree=None): |
| 212 | r""" |
| 213 | Returns the covariant `z_0` of f in H (upper half plane) from [SC2002]. |
| 214 | |
| 215 | INPUT: |
| 216 | |
| 217 | - ``f`` - a polynomial |
| 218 | |
| 219 | - ``precision`` - An integer (Default = 53) What precision for real and |
| 220 | complex numbers should the algorithm use |
| 221 | |
| 222 | - ``degree`` - the degree to homogenize with |
| 223 | |
| 224 | OUTPUT: |
| 225 | |
| 226 | A point z in the upper half plane associated to f. If f is defined over |
| 227 | a totally real number field of degree g, then the output is a sequence |
| 228 | of length g. |
| 229 | |
| 230 | EXAMPLES:: |
| 231 | |
| 232 | sage: z0 = sage.schemes.hyperelliptic_curves.reduction.covariant_z0 |
| 233 | sage: f = x^6 + 3*x^3 +4*x^2+ 3 |
| 234 | sage: z0(f) |
| 235 | -0.0478810014556224 + 1.20469810839717*I |
| 236 | sage: -z0(f)^-1 |
| 237 | |
| 238 | sage: z0(f.reverse()) |
| 239 | |
| 240 | sage: z0(f)-1 |
| 241 | |
| 242 | sage: z0(f(x+1)) |
| 243 | |
| 244 | Comparing the different answers using different precision :: |
| 245 | |
| 246 | sage: z0 = sage.schemes.hyperelliptic_curves.reduction.covariant_z0 |
| 247 | sage: P.<x> = QQ[] |
| 248 | sage: f = x^8 + 24*x^7 + 3*x^5 +4*x^2+ 3 |
| 249 | sage: z0(f) |
| 250 | -0.00832825511768667 + 0.835323919995269*I |
| 251 | sage: z0(f,200) |
| 252 | -0.0083282551176867045891624031810609648523736031152551979546157 + 0.83532391999526880130877046174257268976426032236105109987080*I |
| 253 | |
| 254 | If f is not in `\QQ` then ``covariant_z0`` returns a list where every entry |
| 255 | is a z associated to one embeddings of f in `\RR`:: |
| 256 | |
| 257 | sage: z0 = sage.schemes.hyperelliptic_curves.reduction.covariant_z0 |
| 258 | sage: X = var('X') |
| 259 | sage: k = NumberField(X^2-41,'a') |
| 260 | sage: P.<x> = k[] |
| 261 | sage: a = k.an_element() |
| 262 | sage: f = x^8 + 24*x^7 + 3*x^5 +4*x^2+ 3 |
| 263 | sage: z0(f) |
| 264 | [-0.00832825511768667 + 0.835323919995269*I, -0.00832825511768667 + 0.835323919995269*I] |
| 265 | sage: f = x^6 + (a+1)*x^5 + (5*a+2)*x^2 + (a/2+1/2)*x + 5*a |
| 266 | [0.167978514759694 + 1.68636006242584*I, -0.125921323440448 + 1.64323412027265*I] |
| 267 | |
| 268 | Even though Stoll and Cremona allow roots of higher multiplicity, and roots |
| 269 | at infinity (i.e. degree > f.degree()), this is currently not fully |
| 270 | supported:: |
| 271 | |
| 272 | sage: z0 = sage.schemes.hyperelliptic_curves.reduction.covariant_z0 |
| 273 | sage: P.<x> = RR[] |
| 274 | sage: i = (x-3)^2 * (x^2 + 2)^2 |
| 275 | sage: z0(i) |
| 276 | Traceback (most recent call last): |
| 277 | ... |
| 278 | ValueError: Cannot convert NaN or infinity to Pari float |
| 279 | |
| 280 | sage: j = x^8 + 3 |
| 281 | sage: z0(j, degree=9) |
| 282 | Traceback (most recent call last): |
| 283 | ... |
| 284 | NotImplementedError |
| 285 | |
| 286 | ALGORITHM: |
| 287 | |
| 288 | This algorithm is taken straight from [SC2002] page 6 |
| 289 | |
| 290 | Let n be the degree of the polynomial f. For each embeddings of f into |
| 291 | `\RR` let `\alpha_j` be the roots of f. Define a quadratic polynomial |
| 292 | `Q_0` by `Q_0(x) = \sum_{j=1}^{n} |
| 293 | (x-\alpha_j)(x-\bar{\alpha_j})/|f'(\alpha_j)|^{2/(n-2)}` |
| 294 | Finally take `z_0` to be the root with positive imaginary part of `Q_0` |
| 295 | |
| 296 | REFERENCE: |
| 297 | |
| 298 | [SC2002] M. Stoll and J. Cremona "On the reduction theory of binary forms" |
| 299 | J. Reine Angew. Math. 565 (2003) pp. 79-99 |
| 300 | """ |
| 301 | C_out = ComplexField(precision) |
| 302 | |
| 303 | internal_precision = 6*precision # TODO: make an informed choice |
| 304 | C = ComplexField(internal_precision) |
| 305 | R = RealField(internal_precision) |
| 306 | |
| 307 | |
| 308 | if degree is None: |
| 309 | degree = f.degree() |
| 310 | elif degree != f.degree(): |
| 311 | raise NotImplementedError("Please implement covariant_z0 for "\ |
| 312 | "binary forms with roots at infinity") |
| 313 | |
| 314 | k = f.base_ring() |
| 315 | x = R['x'].gen() |
| 316 | if k is QQ or is_NumberField(k): |
| 317 | Phi = k.embeddings(C) |
| 318 | else: |
| 319 | Phi = [C.convert_map_from(k)] |
| 320 | I = C.gen() |
| 321 | |
| 322 | # Actually, the roots should be taken with multiplicities, |
| 323 | # but this does not matter, as the relevant coefficients |
| 324 | # are divided by infinity anyway. |
| 325 | |
| 326 | z_list = [] |
| 327 | for emb in Phi: |
| 328 | gcoeffs = [emb(coeff) for coeff in f.coeffs()] |
| 329 | g = sum([gcoeffs[indx] * x**indx for indx in range(degree+1)]) |
| 330 | roots = g.roots(C, multiplicities=False) |
| 331 | gdiv = g.derivative() |
| 332 | q0 = R(0) |
| 333 | for ai in roots: |
| 334 | q0 += ((x - ai) * (x - ai.conjugate()) / |
| 335 | (abs(gdiv(ai))**(2 / (degree - 2)))) |
| 336 | z1,z2 = q0.roots(C, multiplicities=False) |
| 337 | if z1.imag_part() > 0: |
| 338 | z0 = z1 |
| 339 | else: |
| 340 | z0 = z2 |
| 341 | z_list.append(z0) |
| 342 | |
| 343 | return [C_out(z) for z in z_list] |
| 344 | |
| 345 | |
| 346 | def _HilbertFundDomData(K, embs=None, group='sl'): |
| 347 | """ |
| 348 | Create an object containing the data for a fundamental domain |
| 349 | for the action of SL_2(O_K) on (CC\setminus RR)^g for a totally real |
| 350 | field K of degree g. |
| 351 | |
| 352 | Warning: assumes deg(K) <= 2 and cl(K) = 1. |
| 353 | |
| 354 | Input: |
| 355 | |
| 356 | - ``K`` -- A totally real number field. |
| 357 | |
| 358 | - ``embs`` (default=None) -- None, or the complete set |
| 359 | of embeddings of K into a real field that should be used, |
| 360 | or the precision to use for creating such a set of embeddings. |
| 361 | |
| 362 | - ``group`` (default='sl') -- a string: 'sl' for SL_2(O_K), |
| 363 | 'gl' for GL_2(O_K), 'plus' for |
| 364 | the subgroup GL_2(O_K)^+ of elements of totally positive |
| 365 | determinant in GL_2(O_K). |
| 366 | |
| 367 | EXAMPLES:: |
| 368 | |
| 369 | sage: from sage.schemes.hyperelliptic_curves.stoll_cremona import _HilbertFundDomData |
| 370 | sage: K = QuadraticField(41,'a') |
| 371 | sage: K(K.unit_group().gen(1)).norm() |
| 372 | -1 |
| 373 | sage: z = [0.01*CC.gen()+23.34,0.02+0.03*CC.gen()] |
| 374 | sage: h = _HilbertFundDomData(K, group='sl') |
| 375 | sage: r = h.reduce(z); r |
| 376 | ( |
| 377 | [ 9/2*a - 11/2 -125/2*a + 801/2] |
| 378 | |
| 379 | [ -3/2*a - 17/2 2*a - 13], [1.03241154280086 + |
| 380 | 15.5033759282630*I, -1.33856972357549 + 0.0496716993976175*I] |
| 381 | ) |
| 382 | sage: prod([e.imag() for e in r[1]]), max([log(abs(e.imag())) for e in r[1]]) |
| 383 | (0.770079028756939, 2.74105780203325) |
| 384 | sage: h = _HilbertFundDomData(K, group='gl') |
| 385 | sage: r = h.reduce(z); r |
| 386 | ( |
| 387 | [15/2*a + 99/2 5/2*a - 37/2] |
| 388 | |
| 389 | [-3/2*a - 17/2 2*a - 13], [1.83799715597054 + 0.242181137053011*I, |
| 390 | -2.13249730501399 - 3.17976469236058*I] |
| 391 | ) |
| 392 | sage: prod([e.imag() for e in r[1]]), max([log(abs(e.imag())) for e in r[1]]) |
| 393 | (-0.770079028756901, 1.15680719794121) |
| 394 | """ |
| 395 | if K.degree() == 1: |
| 396 | return _HilbertFundDomData_classical(embs) |
| 397 | if K.degree() == 2: |
| 398 | return _HilbertFundDomData_quadratic(K, embs, group) |
| 399 | return _HilbertFundDomData_generic(K, embs, group) |
| 400 | |
| 401 | |
| 402 | class _HilbertFundDomData_generic(SageObject): |
| 403 | """ |
| 404 | An object containing the data for a fundamental domain |
| 405 | for the action of SL_2(O_K) on HH^g for a totally real |
| 406 | field K of degree g. |
| 407 | """ |
| 408 | |
| 409 | _K = None |
| 410 | _embs = None |
| 411 | _RR = None |
| 412 | |
| 413 | def __init__(self, K, embs, group): |
| 414 | """ |
| 415 | Create the object. |
| 416 | """ |
| 417 | self._K = K |
| 418 | if embs == None: |
| 419 | embs = K.embeddings(RR) |
| 420 | self._RR = RR |
| 421 | elif embs in ZZ: |
| 422 | self._RR = RealField(embs) |
| 423 | embs = K.embeddings(self._RR) |
| 424 | else: |
| 425 | self._RR = embs[0].codomain() |
| 426 | self._embs = embs |
| 427 | |
| 428 | def to_additive(self, x): |
| 429 | """ |
| 430 | Move x in RR^g to an additive fundamental domain. |
| 431 | |
| 432 | Returns an pair (a, x') with a in O_K and x' = x-a in the fundamental domain. |
| 433 | """ |
| 434 | raise NotImplementedError, "Sorry reduction to a Hilbert fundamental domain is not implemented for g>2" |
| 435 | |
| 436 | def step(self, z, bound_on_c=None): |
| 437 | """ |
| 438 | Do one reduction round. |
| 439 | |
| 440 | Returns a pair (M, z') with z' = M(z) nearer the fundamental domain if possible. |
| 441 | """ |
| 442 | M1 = self.up(z, bound_on_c=bound_on_c) |
| 443 | z = self.act(M1, z) |
| 444 | M2 = self.to_multiplicative([e.imag() for e in z])[0] |
| 445 | z = self.act(M2, z) |
| 446 | t = self.to_additive([e.real() for e in z])[0] |
| 447 | M3 = Matrix([[1, -t], [0, 1]]) |
| 448 | z = self.act(M3, z) |
| 449 | return M3*M2*M1, z |
| 450 | |
| 451 | def reduce(self, z, max_rounds = 0, bound_on_c=None): |
| 452 | """ |
| 453 | Do reduction rounds until reduced. |
| 454 | |
| 455 | Returns a pair (M, z') with z' = M(z) in the fundamental domain. |
| 456 | |
| 457 | max_rounds = 0: go on until finishes, risking an infinite loop from bugs or rounding errors |
| 458 | max_rounds > 0: do at most that number of rounds |
| 459 | """ |
| 460 | M = Matrix([[1,0],[0,1]]) |
| 461 | rounds = 0 |
| 462 | while true: |
| 463 | rounds = rounds + 1 |
| 464 | M1, z = self.step(z, bound_on_c=bound_on_c) |
| 465 | M = M1 * M |
| 466 | if M1 == 1 or rounds == max_rounds: |
| 467 | return M, z |
| 468 | |
| 469 | def act(self, M, z): |
| 470 | """ |
| 471 | Return M(z) |
| 472 | """ |
| 473 | a = M[0,0] |
| 474 | b = M[0,1] |
| 475 | c = M[1,0] |
| 476 | d = M[1,1] |
| 477 | embs = self._embs |
| 478 | return [((embs[k](a)*z[k] + embs[k](b))/(embs[k](c)*z[k] + embs[k](d))) for k in range(len(embs))] |
| 479 | |
| 480 | def _up_cd_lll(self, z): |
| 481 | """ |
| 482 | As _up_cd, but instead of making prod(Im(z)) maximal |
| 483 | and log(Im(z)) in a fundamental domain, try to make min(Im(z)) large |
| 484 | using LLL. |
| 485 | """ |
| 486 | # The new imaginary part Im(z) / |cz+d|^2 |
| 487 | # In particular, max(|cz+d|^2/Im(z)) must be as small as possible, |
| 488 | # and less than 1 / mim(Im(z)). |
| 489 | # Let v = ((cz+d)/sqrt(Im(z)))_i. |
| 490 | # Then we want |v|_infty to be as small as possible and less than |
| 491 | # 1 / min(sqrt(Im(z))). |
| 492 | # Note, g*|v|_2 >= g*|v|_infty^2 >= |v|_2, so it is necessary that |
| 493 | # |v|_2 < g/(min(Im(z))) and sufficient that |v|_2 < 1/(min(Im(z))) |
| 494 | g = self._K.degree() |
| 495 | embs = self._embs |
| 496 | m = 2**(z[0].parent().precision()) |
| 497 | O = self._K.maximal_order() |
| 498 | vectors = ( [ [(embs[i](b)) for i in range(g)] + |
| 499 | [0 for i in range(g)] for b in O.basis()] + |
| 500 | [ [(embs[i](b)*z[i].real()) for i in range(g)] + |
| 501 | [ (embs[i](b)*z[i].imag()) for i in range(g)] |
| 502 | for b in O.basis()]) |
| 503 | id = identity_matrix(2*g) |
| 504 | vectors = [id[j].list() + |
| 505 | [(m*vectors[j][i]/(z[i%g].imag().abs().sqrt())).round() |
| 506 | for i in range(len(vectors[j]))] for j in range(2*g)] |
| 507 | # c = [] |
| 508 | # for i in range(len(vectors)): |
| 509 | # for j in range(i, len(vectors)): |
| 510 | # c.append(sum([vectors[i][k]*vectors[j][k] for k in range(2*g)])) |
| 511 | # Q = QuadraticForm(ZZ, 2*g, c) |
| 512 | # l = Q.short_vector_list_up_to_length(g/min([z.imag().abs()])) |
| 513 | # for x in l[1:]: |
| 514 | M = Matrix(vectors).LLL() |
| 515 | d = sum([M[0][i]*O.basis()[i] for i in range(g)]) |
| 516 | c = sum([M[0][g+i]*O.basis()[i] for i in range(g)]) |
| 517 | B = c * O + d * O |
| 518 | if B == 1*O: |
| 519 | return c, d |
| 520 | if B.is_principal(): |
| 521 | g = B.gens_reduced()[0] |
| 522 | return c/g, d/g |
| 523 | # TODO: if B is not principal, find another candidate pair c, d |
| 524 | |
| 525 | def up(self, z, bound_on_c=None): |
| 526 | """ |
| 527 | Move z in HH^g, i.e., increase the norm of its imaginary part, |
| 528 | if possible. |
| 529 | |
| 530 | Returns either None or a matrix M in SL_2(O_K) that such that |
| 531 | Mz has a larger norm of its imaginary part. |
| 532 | |
| 533 | Assumes K has class number one. |
| 534 | |
| 535 | EXAMPLES:: |
| 536 | |
| 537 | sage: from sage.schemes.hyperelliptic_curves.stoll_cremona import _HilbertFundDomData |
| 538 | sage: h = _HilbertFundDomData(QuadraticField(15)) |
| 539 | sage: h.up([123+CC.0*0.05, 456+CC.0*0.0234]) |
| 540 | [ 0 -1] |
| 541 | [ 1 -43*a - 290] |
| 542 | sage: h.act(_, [123+CC.0*0.05, 456+CC.0*0.0234]) |
| 543 | [2.14072851429831 + 0.231823024326897*I, 1.85425166669223 + 0.0806070738044104*I] |
| 544 | h.reduce([123+CC.0*0.05, 456+CC.0*0.0234]) |
| 545 | ( |
| 546 | [ 2*a + 3 -708*a - 2158] |
| 547 | [ -2 86*a + 579], [-0.0404815976687969 + 3.15206362011738*I, -0.618912376653362 + 2.90579848705465*I] |
| 548 | ) |
| 549 | """ |
| 550 | cd = self._up_cd(z, bound_on_c=bound_on_c) |
| 551 | if cd == None: |
| 552 | return Matrix([[1,0],[0,1]]) |
| 553 | c, d = cd |
| 554 | a = d.inverse_mod(c) |
| 555 | b = (a*d-1)/c |
| 556 | return Matrix([[a,b],[c,d]]) |
| 557 | |
| 558 | def _up_cd(self, z): |
| 559 | """ |
| 560 | Move z up in HH^g, i.e., increase the norm of its imaginary part, |
| 561 | if possible. Actually, not implemented in this base class: just calls |
| 562 | _up_cd_lll instead. |
| 563 | """ |
| 564 | return self._up_cd_lll(z) |
| 565 | |
| 566 | |
| 567 | class _HilbertFundDomData_quadratic(_HilbertFundDomData_generic): |
| 568 | """ |
| 569 | As _HilbertFundDomData, but only for quadratic fields. |
| 570 | """ |
| 571 | |
| 572 | _eps = None |
| 573 | _epsmat = None |
| 574 | _epslog = None |
| 575 | _RR = None |
| 576 | |
| 577 | def __init__(self, K, embs, group): |
| 578 | """ |
| 579 | Create the object. |
| 580 | """ |
| 581 | _HilbertFundDomData_generic.__init__(self, K, embs, group) |
| 582 | eps = K(K.unit_group().gens()[1]) |
| 583 | if group == 'sl' or (group == 'plus' and eps.norm() == -1): |
| 584 | self._epsmat = Matrix([[eps, 0], [0, eps**-1]]) |
| 585 | eps = eps**2 |
| 586 | elif group in ['plus', 'gl']: |
| 587 | if self._embs[0](eps) < 0: |
| 588 | eps = -eps |
| 589 | self._epsmat = Matrix([[eps, 0], [0, 1]]) |
| 590 | else: |
| 591 | raise ValueError, "Unknown group in _HilbertFundDomData: %s" %group |
| 592 | self._eps = eps |
| 593 | self._epslog = self._RR(self._embs[1](eps).abs().log()) |
| 594 | |
| 595 | def to_additive(self, x): |
| 596 | """ |
| 597 | Move x in RR^g to an additive fundamental domain. |
| 598 | |
| 599 | Returns an pair (a, x') with a in O_K and x' = x-a in the fundamental domain. |
| 600 | |
| 601 | Examples:: |
| 602 | |
| 603 | sage: from sage.schemes.hyperelliptic_curves.stoll_cremona import _HilbertFundDomData |
| 604 | sage: h = _HilbertFundDomData(QuadraticField(5)) |
| 605 | sage: h.to_additive([3.2,7.9]) |
| 606 | (a + 6, [-0.563932022500210, -0.336067977499789]) |
| 607 | """ |
| 608 | embs = self._embs |
| 609 | K = self._K |
| 610 | om = K.maximal_order().gens()[1] |
| 611 | y = ((x[1]-x[0]) / (embs[1](om) - embs[0](om))).round() |
| 612 | a = y * om |
| 613 | x = [x[0] - y*embs[0](om), x[1] - y*embs[1](om)] |
| 614 | y = ((x[1]+x[0]) / 2).round() |
| 615 | a = a + y |
| 616 | x = [x[0] - y, x[1] - y] |
| 617 | return a, x |
| 618 | |
| 619 | def to_multiplicative(self, y): |
| 620 | """ |
| 621 | Move y in RR_+^g to a multiplicative fundamental domain. |
| 622 | |
| 623 | Returns a pair (M, y') with M a diagonal matrix and M*y = y' and y' in a fundamental domain. |
| 624 | |
| 625 | Examples:: |
| 626 | |
| 627 | sage: from sage.schemes.hyperelliptic_curves.stoll_cremona import _HilbertFundDomData |
| 628 | sage: h = _HilbertFundDomData(QuadraticField(15)) |
| 629 | sage: h.to_multiplicative([1,45678]) |
| 630 | ( |
| 631 | [-a + 4 0] |
| 632 | [ 0 a + 4], [61.9838667696635, 736.933695500892] |
| 633 | ) |
| 634 | sage: h.act(_[0], [CC.0, CC.0*45678]) |
| 635 | [61.9838667696594*I, 736.933695500891*I] |
| 636 | sage: h = _HilbertFundDomData(QuadraticField(15), sl=False) |
| 637 | sage: h.to_multiplicative([1,45678]) |
| 638 | ( |
| 639 | [-63*a + 244 0] |
| 640 | [ 0 1], [487.997950815860, 93.6028520695257] |
| 641 | ) |
| 642 | sage: h.act(_[0], [CC.0, CC.0*45678]) |
| 643 | [487.997950811067*I, 93.6028520686063*I] |
| 644 | """ |
| 645 | embs = self._embs |
| 646 | K = self._K |
| 647 | eps = self._eps |
| 648 | l = self._epslog |
| 649 | e = ((y[1].abs().log() - y[0].abs().log()) / (2*l)).round() |
| 650 | a = eps**e |
| 651 | y = [y[0]/embs[0](a), y[1]/embs[1](a)] |
| 652 | return self._epsmat**-e, y |
| 653 | |
| 654 | |
| 655 | def _up_cd(self, z, bound_on_c=None): |
| 656 | """ |
| 657 | Move z in HH^g, i.e., increase the norm of its imaginary part, |
| 658 | if possible. |
| 659 | |
| 660 | Returns either None or a matrix M in SL_2(O_K) that such that |
| 661 | Mz has a larger norm of its imaginary part. |
| 662 | |
| 663 | Assumes K has class number one. |
| 664 | |
| 665 | bound_on_c None means try a small one first, double each time until a proven one is reached |
| 666 | """ |
| 667 | embs = self._embs |
| 668 | cd = self._up_cd_lll(z) |
| 669 | K = self._K |
| 670 | O = K.maximal_order() |
| 671 | if not cd is None: |
| 672 | (c, d) = cd |
| 673 | cz_plus_d = [embs[i](c)*z[i] + embs[i](d) for i in range(2)] |
| 674 | if abs(prod(cz_plus_d)) < 1: |
| 675 | assert c*O+d*O == 1*O |
| 676 | return c, d |
| 677 | # Assuming class number 1 and the above did not work, we are pretty |
| 678 | # high up in the fundamental domain, so the search ranges below are not |
| 679 | # that large. |
| 680 | # TODO: two different notions of reduced are tested, so something |
| 681 | # may be missed here. In that case, we may be high up after all. |
| 682 | |
| 683 | om = O.gens()[1] |
| 684 | |
| 685 | # M = (a, b; c, d) we need norm(cz+d)<1 |
| 686 | # and we want such c and d only up to units. |
| 687 | # We have 1 > norm(|cz+d|) >= |norm(c)|*norm(im(z)), |
| 688 | # so |norm(c)| < norm(im(z))^-1. |
| 689 | # We can find all such c up to units by enumerating ideals |
| 690 | # and finding generators. |
| 691 | # Given c, we want a coprime d with norm(|cx+d|)^2 minimal. |
| 692 | x = [z[0].real(), z[1].real()] |
| 693 | y = [z[0].imag(), z[1].imag()] |
| 694 | norm_im_z = self._RR(abs(prod(y))) |
| 695 | bound = ZZ((norm_im_z**-1).floor()) |
| 696 | # if bound_on_c is None: |
| 697 | # bound_to_try = 1 |
| 698 | # while bound_to_try < 2*bound: |
| 699 | # m = self._up_cd(z, bound_on_c = bound_to_try) |
| 700 | # if m != None: |
| 701 | # return m |
| 702 | # bound_to_try = bound_to_try * 2 |
| 703 | # return None |
| 704 | if (not bound_on_c is None) and bound > bound_on_c: |
| 705 | bound = bound_on_c |
| 706 | ids = K.ideals_of_bdd_norm(bound) |
| 707 | for n in ids: |
| 708 | for C in ids[n]: |
| 709 | if C.is_principal(): |
| 710 | eps = K(K.unit_group().gens()[1]) |
| 711 | c = C.gens_reduced()[0] |
| 712 | c_times_y = [embs[0](c)*y[0], embs[1](c)*y[1]] |
| 713 | e = ((self._RR(c_times_y[0]/c_times_y[1]).abs().log()) / |
| 714 | (2*(self._RR(embs[0](eps)).abs().log()))) |
| 715 | # print c, e |
| 716 | c = c / eps ** (e.round()) |
| 717 | # print c |
| 718 | |
| 719 | # Next, we try tofind d with |norm(cz+d)|<1. |
| 720 | # We start by moving to an additive fundamental domain, |
| 721 | # and do a and thorough search if that does not work |
| 722 | # (which it usually does not). |
| 723 | c_times_x = [embs[0](c)*x[0], embs[1](c)*x[1]] |
| 724 | d0, xp = self.to_additive(c_times_x) |
| 725 | d0 = -d0 |
| 726 | # now xp = c*x + d_0 is in the fundamental domain |
| 727 | cz_plus_d = [embs[0](c)*z[0] + embs[0](d0), embs[1](c)*z[1] + embs[1](d0)] |
| 728 | # Actually, the following check already uses class number != 1, because |
| 729 | # otherwise the check is supposed to be < Norm(c*O+d0*O) or something. |
| 730 | |
| 731 | if abs(prod(cz_plus_d)) < 1: |
| 732 | if c*O+d0*O == 1*O: |
| 733 | return c, d0 |
| 734 | # In the class number > 1 case, we may want to do something |
| 735 | # if the test c*O+d0*O == 1*O fails. |
| 736 | |
| 737 | # We have 1 > norm(|cz+d|^2) = |
| 738 | |
| 739 | # |norm(c)|^2*norm(im(z))^2 |
| 740 | # + |c_1|^2*im(z_1)^2 *|c_0x_0+d_0|^2 |
| 741 | # + |c_0|^2*im(z_0)^2 *|c_1x_1+d_1|^2 |
| 742 | # + norm(|cx+d|)^2, |
| 743 | # define x_0' = c_0x_0+d_0 |
| 744 | # and x_1' = c_1x_1+d_1, then they have the following |
| 745 | # absolute value bounds: |
| 746 | # |
| 747 | bound_on_x0p = ((1 - n**2 * norm_im_z**2) |
| 748 | / (embs[1](c)**2 * y[1]**2)).sqrt() |
| 749 | bound_on_x1p = ((1 - n**2 * norm_im_z**2) |
| 750 | / (embs[0](c)**2 * y[0]**2)).sqrt() |
| 751 | d0_interval = [-embs[0](c) * x[0] - bound_on_x0p, |
| 752 | -embs[0](c) * x[0] + bound_on_x0p] |
| 753 | d1_interval = [-embs[1](c) * x[1] - bound_on_x1p, |
| 754 | -embs[1](c) * x[1] + bound_on_x1p] |
| 755 | # Now write d = e + f*om = (g + f*delta)/2, |
| 756 | # where g = 2*e + f*trace(om), |
| 757 | # and om = (trace(om) + delta) / 2. |
| 758 | delta = 2*om - om.trace() |
| 759 | assert embs[0](delta) < 0 |
| 760 | assert embs[1](delta) > 0 |
| 761 | delta_num = embs[1](delta) |
| 762 | # Then g = d0 + d1 is in the following interval: |
| 763 | g_interval = [d0_interval[i] + d1_interval[i] |
| 764 | for i in [0,1]] |
| 765 | g_start = sum([-embs[0](c) * x[i] for i in [0,1]]).round() |
| 766 | g_bound = (bound_on_x0p + bound_on_x1p + 3).floor() |
| 767 | t = om.trace() |
| 768 | for i in xrange(g_bound): |
| 769 | # TODO: this range can be way too big, something needs |
| 770 | # to be done about that. |
| 771 | if t == 0 or not ((g_start + i) % 2): |
| 772 | for s in ([-1,1] if i else [1]): |
| 773 | g = g_start + s * i |
| 774 | # now f = (-g + 2*d)/delta, and we know |
| 775 | # bounds on each embedding: |
| 776 | f0_lower = (-g + 2*d0_interval[1]) / -delta_num |
| 777 | f1_lower = (-g + 2*d1_interval[0]) / delta_num |
| 778 | f0_upper = (-g + 2*d0_interval[0]) / -delta_num |
| 779 | f1_upper = (-g + 2*d1_interval[1]) / delta_num |
| 780 | f_lower = max([f0_lower, f1_lower]).ceil() |
| 781 | f_upper = min([f0_upper, f1_upper]).floor() |
| 782 | for f in range(f_lower, f_upper+1): |
| 783 | if not ((g - f*t) % 2): |
| 784 | d = (g+f*delta)/2 |
| 785 | assert d in O |
| 786 | cz_plus_d = [embs[0](c)*z[0] + |
| 787 | embs[0](d), |
| 788 | embs[1](c)*z[1] + |
| 789 | embs[1](d)] |
| 790 | if abs(prod(cz_plus_d)) < 1: |
| 791 | if c*O+d*O == 1*O: |
| 792 | return c, d |
| 793 | # In the class number > 1 case, |
| 794 | # we may want to do something |
| 795 | # if the test c*O+d0*O == 1*O |
| 796 | # fails. For now, the code just |
| 797 | # skips c,d in that case. |
| 798 | # |
| 799 | # # if we write d = e + f*om, then d0+d1 = 2*e + f*tr(om) |
| 800 | # # is in the following interval: |
| 801 | # |
| 802 | # d0-d1 = f*(om_0 - om_1) is in the following interval |
| 803 | # interval = [d0_interval[0]-d1_interval[1], d0_interval[1]-d1_interval[0]] |
| 804 | # f_interval = [interval[i] / (embs[0](om)-embs[1](om)) for i in [0,1]] |
| 805 | # f_interval.sort() |
| 806 | # for f in range(int(self._RR(f_interval[0]).ceil()), int(self._RR(f_interval[1]).floor()+1)): |
| 807 | # e_lower = self._RR(max([d0_interval[0]-f*embs[0](om), d1_interval[0]-f*embs[1](om)])).ceil() |
| 808 | # e_upper = self._RR(min([d0_interval[1]-f*embs[0](om), d1_interval[1]-f*embs[1](om)])).floor() |
| 809 | # for e in range(int(e_lower), int(e_upper+1)): |
| 810 | # d = e + f*om |
| 811 | # cz_plus_d = [embs[0](c)*z[0] + embs[0](d), embs[1](c)*z[1] + embs[1](d)] |
| 812 | # if abs(prod(cz_plus_d)) < 1: |
| 813 | # if c*O+d*O == 1*O: |
| 814 | # return c, d |
| 815 | # # In the class number > 1 case, we may want to do something |
| 816 | # # if the test c*O+d0*O == 1*O fails. |
| 817 | |
| 818 | |
| 819 | class _HilbertFundDomData_classical(_HilbertFundDomData_generic): |
| 820 | """ |
| 821 | As _HilbertFundDomData, but only for the classical upper half |
| 822 | plane and the action of SL_2(ZZ). |
| 823 | """ |
| 824 | |
| 825 | _eps = None |
| 826 | |
| 827 | def __init__(self, embs=None): |
| 828 | """ |
| 829 | Create the object. |
| 830 | """ |
| 831 | self._K = QQ |
| 832 | |
| 833 | def to_additive(self, x): |
| 834 | """ |
| 835 | Move x in RR^1 to an additive fundamental domain. |
| 836 | """ |
| 837 | a = x[0].round() |
| 838 | x = [x[0]-a] |
| 839 | return a, x |
| 840 | |
| 841 | def to_multiplicative(self, y): |
| 842 | """ |
| 843 | Returns (identitymatrix(2), y). |
| 844 | """ |
| 845 | return Matrix([[1,0],[0,1]]), y |
| 846 | |
| 847 | |
| 848 | def up(self, z): |
| 849 | """ |
| 850 | Move z up in HH, i.e., increase its imaginary part, |
| 851 | if possible by the transformation z -> -1/z |
| 852 | |
| 853 | Returns either 1 or a matrix M in SL_2(ZZ) such that |
| 854 | Mz has a larger imaginary part. |
| 855 | |
| 856 | Assumes the real part of x is between -1/2 and 1/2 |
| 857 | """ |
| 858 | if z[0].abs() < 1: |
| 859 | return Matrix(QQ, [[0,-1],[1,0]]) |
| 860 | return Matrix(QQ, [[1,0],[0,1]]) |
| 861 | |
| 862 | def step(self, z): |
| 863 | """ |
| 864 | Do one reduction round. |
| 865 | |
| 866 | Returns a pair (M, z') with z' = M(z) nearer the fundamental domain if possible. |
| 867 | """ |
| 868 | M1 = self.up(z) |
| 869 | z = self.act(M1, z) |
| 870 | t = self.to_additive([e.real() for e in z])[0] |
| 871 | M3 = Matrix([[1, -t], [0, 1]]) |
| 872 | z = self.act(M3, z) |
| 873 | return M3*M1, z |
| 874 | |
| 875 | |
| 876 | def act(self, M, z): |
| 877 | """ |
| 878 | Return M(z) |
| 879 | """ |
| 880 | a = M[0,0] |
| 881 | b = M[0,1] |
| 882 | c = M[1,0] |
| 883 | d = M[1,1] |
| 884 | z = z[0] |
| 885 | return [(a*z+b)/(c*z+d)] |
| 886 | |
| 887 | |
| 888 | |