# Ticket #3674: ell_rational_field.py

File ell_rational_field.py, 160.8 KB (added by , 9 years ago) |
---|

Line | |
---|---|

1 | """ |

2 | Elliptic curves over the rational numbers |

3 | |

4 | AUTHORS: |

5 | -- William Stein (2005): first version |

6 | -- William Stein (2006-02-26): fixed Lseries_extended which didn't work |

7 | because of changes elsewhere in SAGE. |

8 | -- David Harvey (2006-09): Added padic_E2, padic_sigma, padic_height, |

9 | padic_regulator methods. |

10 | -- David Harvey (2007-02): reworked padic-height related code |

11 | -- Christian Wuthrich (2007): added padic sha computation |

12 | -- David Roe (2007-9): moved sha, l-series and p-adic functionality to separate files. |

13 | -- John Cremona (2008-01) |

14 | -- Tobias Nagel & Michael Mardaus (2008-07): added integral_points |

15 | """ |

16 | |

17 | #***************************************************************************** |

18 | # Copyright (C) 2005,2006,2007 William Stein <wstein@gmail.com> |

19 | # |

20 | # Distributed under the terms of the GNU General Public License (GPL) |

21 | # |

22 | # This code is distributed in the hope that it will be useful, |

23 | # but WITHOUT ANY WARRANTY; without even the implied warranty of |

24 | # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |

25 | # General Public License for more details. |

26 | # |

27 | # The full text of the GPL is available at: |

28 | # |

29 | # http://www.gnu.org/licenses/ |

30 | #***************************************************************************** |

31 | |

32 | import ell_point |

33 | import formal_group |

34 | import rational_torsion |

35 | from ell_generic import EllipticCurve_generic |

36 | from ell_number_field import EllipticCurve_number_field |

37 | |

38 | import sage.groups.all |

39 | import sage.rings.arith as arith |

40 | import sage.rings.all as rings |

41 | import sage.rings.number_field.number_field as number_field |

42 | import sage.misc.misc as misc |

43 | from sage.misc.all import verbose |

44 | import sage.functions.constants as constants |

45 | import sage.modular.modform.constructor |

46 | import sage.modular.modform.element |

47 | from sage.misc.functional import log |

48 | from sage.rings.padics.factory import Zp, Qp |

49 | |

50 | # Use some interval arithmetic to guarantee correctness. We assume |

51 | # that alpha is computed to the precision of a float. |

52 | # IR = rings.RIF |

53 | #from sage.rings.interval import IntervalRing; IR = IntervalRing() |

54 | |

55 | import sage.matrix.all as matrix |

56 | import sage.databases.cremona |

57 | from sage.libs.pari.all import pari |

58 | import sage.functions.transcendental as transcendental |

59 | import math |

60 | from sage.calculus.calculus import sqrt, arcsin, floor, ceil |

61 | import sage.libs.mwrank.all as mwrank |

62 | import constructor |

63 | from sage.interfaces.all import gp |

64 | |

65 | import ell_modular_symbols |

66 | import padic_lseries |

67 | import padics |

68 | |

69 | from lseries_ell import Lseries_ell |

70 | |

71 | import mod5family |

72 | |

73 | from sage.rings.all import ( |

74 | PowerSeriesRing, LaurentSeriesRing, O, |

75 | infinity as oo, |

76 | Integer, |

77 | Integers, |

78 | IntegerRing, RealField, |

79 | ComplexField, RationalField) |

80 | |

81 | import gp_cremona |

82 | import sea |

83 | |

84 | from gp_simon import simon_two_descent |

85 | |

86 | import ell_tate_curve |

87 | |

88 | factor = arith.factor |

89 | exp = math.exp |

90 | mul = misc.mul |

91 | next_prime = arith.next_prime |

92 | kronecker_symbol = arith.kronecker_symbol |

93 | |

94 | Q = RationalField() |

95 | C = ComplexField() |

96 | R = RealField() |

97 | Z = IntegerRing() |

98 | IR = rings.RealIntervalField(20) |

99 | |

100 | _MAX_HEIGHT=21 |

101 | |

102 | # CMJ is a dict of pairs (j,D) where j is a rational CM j-invariant |

103 | # and D is the corresponding quadratic discriminant |

104 | |

105 | CMJ={ 0: -3, 54000: -12, -12288000: -27, 1728: -4, 287496: -16, |

106 | -3375: -7, 16581375: -28, 8000: -8, -32768: -11, -884736: -19, |

107 | -884736000: -43, -147197952000: -67, -262537412640768000: -163} |

108 | |

109 | |

110 | |

111 | class EllipticCurve_rational_field(EllipticCurve_number_field): |

112 | """ |

113 | Elliptic curve over the Rational Field. |

114 | """ |

115 | def __init__(self, ainvs, extra=None): |

116 | if extra != None: # possibility of two arguments (the first would be the field) |

117 | ainvs = extra |

118 | if isinstance(ainvs, str): |

119 | label = ainvs |

120 | X = sage.databases.cremona.CremonaDatabase()[label] |

121 | EllipticCurve_number_field.__init__(self, Q, X.a_invariants()) |

122 | for attr in ['rank', 'torsion_order', 'cremona_label', 'conductor', |

123 | 'modular_degree', 'gens', 'regulator']: |

124 | s = "_EllipticCurve_rational_field__"+attr |

125 | if hasattr(X,s): |

126 | setattr(self, s, getattr(X, s)) |

127 | return |

128 | EllipticCurve_number_field.__init__(self, Q, ainvs) |

129 | self.__np = {} |

130 | self.__gens = {} |

131 | self.__rank = {} |

132 | self.__regulator = {} |

133 | if self.base_ring() != Q: |

134 | raise TypeError, "Base field (=%s) must be the Rational Field."%self.base_ring() |

135 | |

136 | def _set_rank(self, r): |

137 | """ |

138 | Internal function to set the cached rank of this elliptic curve to r. |

139 | |

140 | WARNING: No checking is done! Not intended for use by users. |

141 | |

142 | EXAMPLES: |

143 | sage: E=EllipticCurve('37a1') |

144 | sage: E._set_rank(99) # bogus value -- not checked |

145 | sage: E.rank() # returns bogus cached value |

146 | 99 |

147 | sage: E.gens() # causes actual rank to be computed |

148 | [(0 : -1 : 1)] |

149 | sage: E.rank() # the correct rank |

150 | 1 |

151 | """ |

152 | self.__rank = {} |

153 | self.__rank[True] = Integer(r) |

154 | |

155 | def _set_torsion_order(self, t): |

156 | """ |

157 | Internal function to set the cached torsion order of this elliptic curve to t. |

158 | |

159 | WARNING: No checking is done! Not intended for use by users. |

160 | |

161 | EXAMPLES: |

162 | sage: E=EllipticCurve('37a1') |

163 | sage: E._set_torsion_order(99) # bogus value -- not checked |

164 | sage: E.torsion_order() # returns bogus cached value |

165 | 99 |

166 | sage: T = E.torsion_subgroup() # causes actual torsion to be computed |

167 | sage: E.torsion_order() # the correct value |

168 | 1 |

169 | """ |

170 | self.__torsion_order = Integer(t) |

171 | |

172 | def _set_cremona_label(self, L): |

173 | """ |

174 | Internal function to set the cached label of this elliptic curve to L. |

175 | |

176 | WARNING: No checking is done! Not intended for use by users. |

177 | |

178 | EXAMPLES: |

179 | sage: E=EllipticCurve('37a1') |

180 | sage: E._set_cremona_label('bogus') |

181 | sage: E.label() |

182 | 'bogus' |

183 | sage: E.database_curve().label() |

184 | '37a1' |

185 | sage: E.label() # no change |

186 | 'bogus' |

187 | sage: E._set_cremona_label(E.database_curve().label()) |

188 | sage: E.label() # now it is correct |

189 | '37a1' |

190 | """ |

191 | self.__cremona_label = L |

192 | |

193 | def _set_conductor(self, N): |

194 | """ |

195 | Internal function to set the cached conductor of this elliptic curve to N. |

196 | |

197 | WARNING: No checking is done! Not intended for use by users. |

198 | Setting to the wrong value will cause strange problems (see |

199 | examples). |

200 | |

201 | EXAMPLES: |

202 | sage: E=EllipticCurve('37a1') |

203 | sage: E._set_conductor(99) # bogus value -- not checked |

204 | sage: E.conductor() # returns bogus cached value |

205 | 99 |

206 | |

207 | This will not work since the conductor is used when searching |

208 | the database: |

209 | sage: E._set_conductor(E.database_curve().conductor()) |

210 | Traceback (most recent call last): |

211 | ... |

212 | RuntimeError: Elliptic curve ... not in the database. |

213 | sage: E._set_conductor(EllipticCurve(E.a_invariants()).database_curve().conductor()) |

214 | sage: E.conductor() # returns correct value |

215 | 37 |

216 | """ |

217 | self.__conductor_pari = Integer(N) |

218 | |

219 | def _set_modular_degree(self, deg): |

220 | """ |

221 | Internal function to set the cached modular degree of this elliptic curve to deg. |

222 | |

223 | WARNING: No checking is done! |

224 | |

225 | EXAMPLES: |

226 | sage: E=EllipticCurve('5077a1') |

227 | sage: E.modular_degree() |

228 | 1984 |

229 | sage: E._set_modular_degree(123456789) |

230 | sage: E.modular_degree() |

231 | 123456789 |

232 | sage: E._set_modular_degree(E.database_curve().modular_degree()) |

233 | sage: E.modular_degree() |

234 | 1984 |

235 | """ |

236 | self.__modular_degree = Integer(deg) |

237 | |

238 | def _set_gens(self, gens): |

239 | """ |

240 | Internal function to set the cached generators of this elliptic curve to gens. |

241 | |

242 | WARNING: No checking is done! |

243 | |

244 | EXAMPLES: |

245 | sage: E=EllipticCurve('5077a1') |

246 | sage: E.rank() |

247 | 3 |

248 | sage: E.gens() # random |

249 | [(-2 : 3 : 1), (-7/4 : 25/8 : 1), (1 : -1 : 1)] |

250 | sage: E._set_gens([]) # bogus list |

251 | sage: E.rank() # unchanged |

252 | 3 |

253 | sage: E._set_gens(E.database_curve().gens()) |

254 | sage: E.gens() |

255 | [(-2 : 3 : 1), (-7/4 : 25/8 : 1), (1 : -1 : 1)] |

256 | """ |

257 | self.__gens = {} |

258 | self.__gens[True] = [self.point(x, check=True) for x in gens] |

259 | self.__gens[True].sort() |

260 | |

261 | def is_integral(self): |

262 | """ |

263 | Returns True if this elliptic curve has integral coefficients (in Z) |

264 | |

265 | EXAMPLES: |

266 | sage: E=EllipticCurve(QQ,[1,1]); E |

267 | Elliptic Curve defined by y^2 = x^3 + x +1 over Rational Field |

268 | sage: E.is_integral() |

269 | True |

270 | sage: E2=E.change_weierstrass_model(2,0,0,0); E2 |

271 | Elliptic Curve defined by y^2 = x^3 + 1/16*x + 1/64 over Rational Field |

272 | sage: E2.is_integral() |

273 | False |

274 | """ |

275 | try: |

276 | return self.__is_integral |

277 | except AttributeError: |

278 | one = Integer(1) |

279 | self.__is_integral = bool(misc.mul([x.denominator() == 1 for x in self.ainvs()])) |

280 | return self.__is_integral |

281 | |

282 | |

283 | def mwrank(self, options=''): |

284 | """ |

285 | Run Cremona's mwrank program on this elliptic curve and |

286 | return the result as a string. |

287 | |

288 | INPUT: |

289 | options -- string; passed when starting mwrank. The format is |

290 | q p<precision> v<verbosity> b<hlim_q> x<naux> c<hlim_c> l t o s d>] |

291 | |

292 | OUTPUT: |

293 | string -- output of mwrank on this curve |

294 | |

295 | NOTE: The output is a raw string and completely illegible |

296 | using automatic display, so it is recommended to use print |

297 | for legible outout |

298 | |

299 | EXAMPLES: |

300 | sage: E = EllipticCurve('37a1') |

301 | sage: E.mwrank() #random |

302 | ... |

303 | sage: print E.mwrank() |

304 | Curve [0,0,1,-1,0] : Basic pair: I=48, J=-432 |

305 | disc=255744 |

306 | ... |

307 | Generator 1 is [0:-1:1]; height 0.05111... |

308 | <BLANKLINE> |

309 | Regulator = 0.05111... |

310 | <BLANKLINE> |

311 | The rank and full Mordell-Weil basis have been determined unconditionally. |

312 | ... |

313 | |

314 | Options to mwrank can be passed: |

315 | sage: E = EllipticCurve([0,0,0,877,0]) |

316 | |

317 | Run mwrank with 'verbose' flag set to 0 but list generators if found |

318 | sage: print E.mwrank('-v0 -l') |

319 | Curve [0,0,0,877,0] : 0 <= rank <= 1 |

320 | Regulator = 1 |

321 | |

322 | Run mwrank again, this time with a higher bound for point |

323 | searching on homogeneous spaces: |

324 | |

325 | sage: print E.mwrank('-v0 -l -b11') |

326 | Curve [0,0,0,877,0] : Rank = 1 |

327 | Generator 1 is [29604565304828237474403861024284371796799791624792913256602210:-256256267988926809388776834045513089648669153204356603464786949:490078023219787588959802933995928925096061616470779979261000]; height 95.980371987964 |

328 | Regulator = 95.980371987964 |

329 | """ |

330 | if options == "": |

331 | from sage.interfaces.all import mwrank |

332 | else: |

333 | from sage.interfaces.all import Mwrank |

334 | mwrank = Mwrank(options=options) |

335 | return mwrank(self.a_invariants()) |

336 | |

337 | def conductor(self, algorithm="pari"): |

338 | """ |

339 | Returns the conductor of the elliptic curve. |

340 | |

341 | INPUT: |

342 | algorithm -- str, (default: "pari") |

343 | "pari" -- use the PARI C-library ellglobalred |

344 | implementation of Tate's algorithm |

345 | "mwrank" -- use Cremona's mwrank implementation of |

346 | Tate's algorithm; can be faster if the |

347 | curve has integer coefficients (TODO: |

348 | limited to small conductor until mwrank |

349 | gets integer factorization) |

350 | "gp" -- use the GP interpreter. |

351 | "all" -- use both implementations, verify that the |

352 | results are the same (or raise an error), |

353 | and output the common value. |

354 | |

355 | EXAMPLE: |

356 | sage: E = EllipticCurve([1, -1, 1, -29372, -1932937]) |

357 | sage: E.conductor(algorithm="pari") |

358 | 3006 |

359 | sage: E.conductor(algorithm="mwrank") |

360 | 3006 |

361 | sage: E.conductor(algorithm="gp") |

362 | 3006 |

363 | sage: E.conductor(algorithm="all") |

364 | 3006 |

365 | |

366 | NOTE: The conductor computed using each algorithm is cached separately. |

367 | Thus calling E.conductor("pari"), then E.conductor("mwrank") and |

368 | getting the same result checks that both systems compute the same answer. |

369 | """ |

370 | |

371 | if algorithm == "pari": |

372 | try: |

373 | return self.__conductor_pari |

374 | except AttributeError: |

375 | self.__conductor_pari = Integer(self.pari_mincurve().ellglobalred()[0]) |

376 | return self.__conductor_pari |

377 | |

378 | elif algorithm == "gp": |

379 | try: |

380 | return self.__conductor_gp |

381 | except AttributeError: |

382 | self.__conductor_gp = Integer(gp.eval('ellglobalred(ellinit(%s,0))[1]'%self.a_invariants())) |

383 | return self.__conductor_gp |

384 | |

385 | elif algorithm == "mwrank": |

386 | try: |

387 | return self.__conductor_mwrank |

388 | except AttributeError: |

389 | if self.is_integral(): |

390 | self.__conductor_mwrank = Integer(self.mwrank_curve().conductor()) |

391 | else: |

392 | self.__conductor_mwrank = Integer(self.minimal_model().mwrank_curve().conductor()) |

393 | return self.__conductor_mwrank |

394 | |

395 | elif algorithm == "all": |

396 | N1 = self.conductor("pari") |

397 | N2 = self.conductor("mwrank") |

398 | N3 = self.conductor("gp") |

399 | if N1 != N2 or N2 != N3: |

400 | raise ArithmeticError, "Pari, mwrank and gp compute different conductors (%s,%s,%s) for %s"%( |

401 | N1, N2, N3, self) |

402 | return N1 |

403 | else: |

404 | raise RuntimeError, "algorithm '%s' is not known."%algorithm |

405 | |

406 | #################################################################### |

407 | # Access to PARI curves related to this curve. |

408 | #################################################################### |

409 | |

410 | def pari_curve(self, prec = None, factor = 1): |

411 | """ |

412 | Return the PARI curve corresponding to this elliptic curve. |

413 | |

414 | INPUT: |

415 | prec -- The precision of quantities calculated for the |

416 | returned curve (in decimal digits). if None, defaults |

417 | to factor * the precision of the largest cached curve |

418 | (or 10 if none yet computed) |

419 | factor -- the factor to increase the precision over the |

420 | maximum previously computed precision. Only used if |

421 | prec (which gives an explicit precision) is None. |

422 | |

423 | EXAMPLES: |

424 | sage: E = EllipticCurve([0, 0,1,-1,0]) |

425 | sage: e = E.pari_curve() |

426 | sage: type(e) |

427 | <type 'sage.libs.pari.gen.gen'> |

428 | sage: e.type() |

429 | 't_VEC' |

430 | sage: e.ellan(10) |

431 | [1, -2, -3, 2, -2, 6, -1, 0, 6, 4] |

432 | |

433 | sage: E = EllipticCurve(RationalField(), ['1/3', '2/3']) |

434 | sage: e = E.pari_curve() |

435 | sage: e.type() |

436 | 't_VEC' |

437 | sage: e[:5] |

438 | [0, 0, 0, 1/3, 2/3] |

439 | """ |

440 | if prec is None: |

441 | try: |

442 | L = self._pari_curve.keys() |

443 | L.sort() |

444 | if factor == 1: |

445 | return self._pari_curve[L[-1]] |

446 | else: |

447 | prec = int(factor * L[-1]) |

448 | self._pari_curve[prec] = pari(self.a_invariants()).ellinit(precision=prec) |

449 | return self._pari_curve[prec] |

450 | except AttributeError: |

451 | pass |

452 | try: |

453 | return self._pari_curve[prec] |

454 | except AttributeError: |

455 | prec = 10 |

456 | self._pari_curve = {} |

457 | except KeyError: |

458 | pass |

459 | self._pari_curve[prec] = pari(self.a_invariants()).ellinit(precision=prec) |

460 | return self._pari_curve[prec] |

461 | |

462 | def pari_mincurve(self, prec = None): |

463 | """ |

464 | Return the PARI curve corresponding to a minimal model |

465 | for this elliptic curve. |

466 | |

467 | INPUT: |

468 | prec -- The precision of quantities calculated for the |

469 | returned curve (in decimal digits). if None, defaults |

470 | to the precision of the largest cached curve (or 10 if |

471 | none yet computed) |

472 | |

473 | EXAMPLES: |

474 | sage: E = EllipticCurve(RationalField(), ['1/3', '2/3']) |

475 | sage: e = E.pari_mincurve() |

476 | sage: e[:5] |

477 | [0, 0, 0, 27, 486] |

478 | sage: E.conductor() |

479 | 47232 |

480 | sage: e.ellglobalred() |

481 | [47232, [1, 0, 0, 0], 2] |

482 | """ |

483 | if prec is None: |

484 | try: |

485 | L = self._pari_mincurve.keys() |

486 | L.sort() |

487 | return self._pari_mincurve[L[len(L) - 1]] |

488 | except AttributeError: |

489 | pass |

490 | try: |

491 | return self._pari_mincurve[prec] |

492 | except AttributeError: |

493 | prec = 10 |

494 | self._pari_mincurve = {} |

495 | except KeyError: |

496 | pass |

497 | e = self.pari_curve(prec) |

498 | mc, change = e.ellminimalmodel() |

499 | self._pari_mincurve[prec] = mc |

500 | # self.__min_transform = change |

501 | return mc |

502 | |

503 | def database_curve(self): |

504 | """ |

505 | Return the curve in the elliptic curve database isomorphic to |

506 | this curve, if possible. Otherwise raise a RuntimeError |

507 | exception. |

508 | |

509 | EXAMPLES: |

510 | sage: E = EllipticCurve([0,1,2,3,4]) |

511 | sage: E.database_curve() |

512 | Elliptic Curve defined by y^2 = x^3 + x^2 + 3*x + 5 over Rational Field |

513 | |

514 | NOTES: The model of the curve in the database can be different |

515 | than the Weierstrass model for this curve, e.g., |

516 | database models are always minimal. |

517 | """ |

518 | try: |

519 | return self.__database_curve |

520 | except AttributeError: |

521 | misc.verbose("Looking up %s in the database."%self) |

522 | D = sage.databases.cremona.CremonaDatabase() |

523 | ainvs = self.minimal_model().ainvs() |

524 | try: |

525 | self.__database_curve = D.elliptic_curve_from_ainvs(self.conductor(), ainvs) |

526 | except RuntimeError: |

527 | raise RuntimeError, "Elliptic curve %s not in the database."%self |

528 | return self.__database_curve |

529 | |

530 | |

531 | def Np(self, p): |

532 | """ |

533 | The number of points on E modulo p, where p is a prime, not |

534 | necessarily of good reduction. (When p is a bad prime, also |

535 | counts the singular point.) |

536 | |

537 | EXAMPLES: |

538 | sage: E = EllipticCurve([0, -1, 1, -10, -20]) |

539 | sage: E.Np(2) |

540 | 5 |

541 | sage: E.Np(3) |

542 | 5 |

543 | sage: E.conductor() |

544 | 11 |

545 | sage: E.Np(11) |

546 | 11 |

547 | """ |

548 | if self.conductor() % p == 0: |

549 | return p + 1 - self.ap(p) |

550 | #raise ArithmeticError, "p (=%s) must be a prime of good reduction"%p |

551 | if p < 1125899906842624: # TODO: choose more wisely? |

552 | return p+1 - self.ap(p) |

553 | else: |

554 | return self.sea(p) |

555 | |

556 | def sea(self, p, early_abort=False): |

557 | r""" |

558 | Return the number of points on $E$ over $\F_p$ computed using |

559 | the SEA algorithm, as implemented in PARI by Christophe Doche |

560 | and Sylvain Duquesne. |

561 | |

562 | INPUT: |

563 | p -- a prime number |

564 | early_abort -- bool (default: Falst); if True an early abort technique |

565 | is used and the computation is interrupted as soon |

566 | as a small divisor of the order is detected. |

567 | |

568 | \note{As of 2006-02-02 this function does not work on |

569 | Microsoft Windows under Cygwin (though it works under |

570 | vmware of course).} |

571 | |

572 | EXAMPLES: |

573 | sage: E = EllipticCurve('37a') |

574 | sage: E.sea(next_prime(10^30)) |

575 | 1000000000000001426441464441649 |

576 | """ |

577 | p = rings.Integer(p) |

578 | return sea.ellsea(self.minimal_model().a_invariants(), p, early_abort=early_abort) |

579 | |

580 | #def __pari_double_prec(self): |

581 | # EllipticCurve_number_field._EllipticCurve__pari_double_prec(self) |

582 | # try: |

583 | # del self._pari_mincurve |

584 | # except AttributeError: |

585 | # pass |

586 | |

587 | #################################################################### |

588 | # Access to mwrank |

589 | #################################################################### |

590 | def mwrank_curve(self, verbose=False): |

591 | """ |

592 | Construct an mwrank_EllipticCurve from this elliptic curve |

593 | |

594 | The resulting mwrank_EllipticCurve has available methods from |

595 | John Cremona's eclib library. |

596 | |

597 | EXAMPLES: |

598 | sage: E=EllipticCurve('11a1') |

599 | sage: EE=E.mwrank_curve() |

600 | sage: EE |

601 | y^2+ y = x^3 - x^2 - 10*x - 20 |

602 | sage: type(EE) |

603 | <class 'sage.libs.mwrank.interface.mwrank_EllipticCurve'> |

604 | sage: EE.isogeny_class() |

605 | ([[0, -1, 1, -10, -20], [0, -1, 1, -7820, -263580], [0, -1, 1, 0, 0]], |

606 | [[0, 5, 5], [5, 0, 0], [5, 0, 0]]) |

607 | """ |

608 | try: |

609 | return self.__mwrank_curve |

610 | except AttributeError: |

611 | pass |

612 | self.__mwrank_curve = mwrank.mwrank_EllipticCurve( |

613 | self.ainvs(), verbose=verbose) |

614 | return self.__mwrank_curve |

615 | |

616 | def two_descent(self, verbose=True, |

617 | selmer_only = False, |

618 | first_limit = 20, |

619 | second_limit = 8, |

620 | n_aux = -1, |

621 | second_descent = 1): |

622 | """ |

623 | Compute 2-descent data for this curve. |

624 | |

625 | INPUT: |

626 | verbose -- (default: True) print what mwrank is doing |

627 | If False, *no output* is |

628 | selmer_only -- (default: False) selmer_only switch |

629 | first_limit -- (default: 20) firstlim is bound on |x|+|z| |

630 | second_limit-- (default: 8) secondlim is bound on log max {|x|,|z| }, |

631 | i.e. logarithmic |

632 | n_aux -- (default: -1) n_aux only relevant for general |

633 | 2-descent when 2-torsion trivial; n_aux=-1 causes default |

634 | to be used (depends on method) |

635 | second_descent -- (default: True) second_descent only relevant for |

636 | descent via 2-isogeny |

637 | OUTPUT: |

638 | |

639 | Nothing -- nothing is returned (though much is printed |

640 | unless verbose=False) |

641 | |

642 | EXAMPLES: |

643 | sage: E=EllipticCurve('37a1') |

644 | sage: E.two_descent(verbose=False) # no output |

645 | """ |

646 | self.mwrank_curve().two_descent(verbose, selmer_only, |

647 | first_limit, second_limit, |

648 | n_aux, second_descent) |

649 | |

650 | |

651 | #################################################################### |

652 | # Etc. |

653 | #################################################################### |

654 | |

655 | def aplist(self, n, python_ints=False): |

656 | r""" |

657 | The Fourier coefficients $a_p$ of the modular form attached to |

658 | this elliptic curve, for all primes $p\leq n$. |

659 | |

660 | INPUT: |

661 | n -- integer |

662 | python_ints -- bool (default: False); if True return a list of |

663 | Python ints instead of SAGE integers. |

664 | |

665 | OUTPUT: |

666 | -- list of integers |

667 | |

668 | EXAMPLES: |

669 | sage: e = EllipticCurve('37a') |

670 | sage: e.aplist(1) |

671 | [] |

672 | sage: e.aplist(2) |

673 | [-2] |

674 | sage: e.aplist(10) |

675 | [-2, -3, -2, -1] |

676 | sage: v = e.aplist(13); v |

677 | [-2, -3, -2, -1, -5, -2] |

678 | sage: type(v[0]) |

679 | <type 'sage.rings.integer.Integer'> |

680 | sage: type(e.aplist(13, python_ints=True)[0]) |

681 | <type 'int'> |

682 | """ |

683 | # How is this result dependant on the real precision in pari? At all? |

684 | e = self.pari_mincurve() |

685 | v = e.ellaplist(n, python_ints=True) |

686 | if python_ints: |

687 | return v |

688 | else: |

689 | return [Integer(a) for a in v] |

690 | |

691 | |

692 | |

693 | def anlist(self, n, python_ints=False): |

694 | """ |

695 | The Fourier coefficients up to and including $a_n$ of the |

696 | modular form attached to this elliptic curve. The ith element |

697 | of the return list is a[i]. |

698 | |

699 | INPUT: |

700 | n -- integer |

701 | python_ints -- bool (default: False); if True return a list of |

702 | Python ints instead of SAGE integers. |

703 | |

704 | OUTPUT: |

705 | -- list of integers |

706 | |

707 | EXAMPLES: |

708 | sage: E = EllipticCurve([0, -1, 1, -10, -20]) |

709 | sage: E.anlist(3) |

710 | [0, 1, -2, -1] |

711 | |

712 | sage: E = EllipticCurve([0,1]) |

713 | sage: E.anlist(20) |

714 | [0, 1, 0, 0, 0, 0, 0, -4, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 8, 0] |

715 | """ |

716 | n = int(n) |

717 | # How is this result dependent on the real precision in Pari? At all? |

718 | e = self.pari_mincurve() |

719 | if n >= 2147483648: |

720 | raise RuntimeError, "anlist: n (=%s) must be < 2147483648."%n |

721 | |

722 | v = [0] + e.ellan(n, python_ints=True) |

723 | if not python_ints: |

724 | v = [Integer(x) for x in v] |

725 | return v |

726 | |

727 | |

728 | # There is some overheard associated with coercing the PARI |

729 | # list back to Python, but it's not bad. It's better to do it |

730 | # this way instead of trying to eval the whole list, since the |

731 | # int conversion is done very sensibly. NOTE: This would fail |

732 | # if a_n won't fit in a C int, i.e., is bigger than |

733 | # 2147483648; however, we wouldn't realistically compute |

734 | # anlist for n that large anyways. |

735 | # |

736 | # Some relevant timings: |

737 | # |

738 | # E <--> [0, 1, 1, -2, 0] 389A |

739 | # E = EllipticCurve([0, 1, 1, -2, 0]); // SAGE or MAGMA |

740 | # e = E.pari_mincurve() |

741 | # f = ellinit([0,1,1,-2,0]); |

742 | # |

743 | # Computation Time (1.6Ghz Pentium-4m laptop) |

744 | # time v:=TracesOfFrobenius(E,10000); // MAGMA 0.120 |

745 | # gettime;v=ellan(f,10000);gettime/1000 0.046 |

746 | # time v=e.ellan (10000) 0.04 |

747 | # time v=E.anlist(10000) 0.07 |

748 | |

749 | # time v:=TracesOfFrobenius(E,100000); // MAGMA 1.620 |

750 | # gettime;v=ellan(f,100000);gettime/1000 0.676 |

751 | # time v=e.ellan (100000) 0.7 |

752 | # time v=E.anlist(100000) 0.83 |

753 | |

754 | # time v:=TracesOfFrobenius(E,1000000); // MAGMA 20.850 |

755 | # gettime;v=ellan(f,1000000);gettime/1000 9.238 |

756 | # time v=e.ellan (1000000) 9.61 |

757 | # time v=E.anlist(1000000) 10.95 (13.171 in cygwin vmware) |

758 | |

759 | # time v:=TracesOfFrobenius(E,10000000); //MAGMA 257.850 |

760 | # gettime;v=ellan(f,10000000);gettime/1000 FAILS no matter how many allocatemem()'s!! |

761 | # time v=e.ellan (10000000) 139.37 |

762 | # time v=E.anlist(10000000) 136.32 |

763 | # |

764 | # The last SAGE comp retries with stack size 40MB, |

765 | # 80MB, 160MB, and succeeds last time. It's very interesting that this |

766 | # last computation is *not* possible in GP, but works in py_pari! |

767 | # |

768 | |

769 | def q_expansion(self, prec): |

770 | r""" |

771 | Return the $q$-expansion to precision prec of the newform |

772 | attached to this elliptic curve. |

773 | |

774 | INPUT: |

775 | prec -- an integer |

776 | |

777 | OUTPUT: |

778 | a power series (in th evariable 'q') |

779 | |

780 | NOTE: If you want the output to be a modular form and not just |

781 | a $q$-expansion, use \code{self.modular_form()}. |

782 | |

783 | EXAMPLES: |

784 | sage: E=EllipticCurve('37a1') |

785 | sage: E.q_expansion(20) |

786 | q - 2*q^2 - 3*q^3 + 2*q^4 - 2*q^5 + 6*q^6 - q^7 + 6*q^9 + 4*q^10 - 5*q^11 - 6*q^12 - 2*q^13 + 2*q^14 + 6*q^15 - 4*q^16 - 12*q^18 + O(q^20) |

787 | |

788 | """ |

789 | return PowerSeriesRing(Q, 'q')(self.anlist(prec), prec, check=True) |

790 | |

791 | def modular_form(self): |

792 | r""" |

793 | Return the cuspidal modular form associated to this elliptic curve. |

794 | |

795 | EXAMPLES: |

796 | sage: E = EllipticCurve('37a') |

797 | sage: f = E.modular_form() |

798 | sage: f |

799 | q - 2*q^2 - 3*q^3 + 2*q^4 - 2*q^5 + O(q^6) |

800 | |

801 | If you need to see more terms in the $q$-expansion: |

802 | sage: f.q_expansion(20) |

803 | q - 2*q^2 - 3*q^3 + 2*q^4 - 2*q^5 + 6*q^6 - q^7 + 6*q^9 + 4*q^10 - 5*q^11 - 6*q^12 - 2*q^13 + 2*q^14 + 6*q^15 - 4*q^16 - 12*q^18 + O(q^20) |

804 | |

805 | |

806 | NOTE: If you just want the $q$-expansion, use |

807 | \code{self.q_expansion(prec)}. |

808 | """ |

809 | try: |

810 | return self.__modular_form |

811 | except AttributeError: |

812 | M = sage.modular.modform.constructor.ModularForms(self.conductor(),weight=2) |

813 | f = sage.modular.modform.element.ModularFormElement_elliptic_curve(M, self) |

814 | self.__modular_form = f |

815 | return f |

816 | |

817 | def modular_symbol_space(self, sign=1, base_ring=Q, bound=None): |

818 | r""" |

819 | Return the space of cuspidal modular symbols associated to |

820 | this elliptic curve, with given sign and base ring. |

821 | |

822 | INPUT: |

823 | sign -- 0, -1, or 1 |

824 | base_ring -- a ring |

825 | |

826 | EXAMPLES: |

827 | sage: f = EllipticCurve('37b') |

828 | sage: f.modular_symbol_space() |

829 | Modular Symbols subspace of dimension 1 of Modular Symbols space of dimension 3 for Gamma_0(37) of weight 2 with sign 1 over Rational Field |

830 | sage: f.modular_symbol_space(-1) |

831 | Modular Symbols subspace of dimension 1 of Modular Symbols space of dimension 2 for Gamma_0(37) of weight 2 with sign -1 over Rational Field |

832 | sage: f.modular_symbol_space(0, bound=3) |

833 | Modular Symbols subspace of dimension 2 of Modular Symbols space of dimension 5 for Gamma_0(37) of weight 2 with sign 0 over Rational Field |

834 | |

835 | NOTE: If you just want the $q$-expansion, use |

836 | \code{self.q_expansion(prec)}. |

837 | """ |

838 | typ = (sign, base_ring) |

839 | try: |

840 | return self.__modular_symbol_space[typ] |

841 | except AttributeError: |

842 | self.__modular_symbol_space = {} |

843 | except KeyError: |

844 | pass |

845 | M = ell_modular_symbols.modular_symbol_space(self, sign, base_ring, bound=bound) |

846 | self.__modular_symbol_space[typ] = M |

847 | return M |

848 | |

849 | def modular_symbol(self, sign=1, normalize=True): |

850 | r""" |

851 | Return the modular symbol associated to this elliptic curve, |

852 | with given sign and base ring. This is the map that sends r/s |

853 | to a fixed multiple of 2*pi*I*f(z)dz from oo to r/s, |

854 | normalized so that all values of this map take values in QQ. |

855 | |

856 | If sign=1, the normalization is such that the p-adic |

857 | L-function associated to this modular symbol is correct. |

858 | I.e., the normalization is the same as for the integral period |

859 | mapping divided by 2. |

860 | |

861 | INPUT: |

862 | sign -- -1, or 1 |

863 | base_ring -- a ring |

864 | normalize -- (default: True); if True, the modular symbol |

865 | is correctly normalized (up to possibly a factor of |

866 | -1 or 2). If False, the modular symbol is almost certainly |

867 | not correctly normalized, i.e., all values will be a |

868 | fixed scalar multiple of what they should be. But |

869 | the initial computation of the modular symbol is |

870 | much faster, though evaluation of it after computing |

871 | it won't be any faster. |

872 | |

873 | EXAMPLES: |

874 | sage: E=EllipticCurve('37a1') |

875 | sage: M=E.modular_symbol(); M |

876 | Modular symbol with sign 1 over Rational Field attached to Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field |

877 | sage: M(1/2) |

878 | 0 |

879 | sage: M(1/5) |

880 | 1/2 |

881 | """ |

882 | typ = (sign, normalize) |

883 | try: |

884 | return self.__modular_symbol[typ] |

885 | except AttributeError: |

886 | self.__modular_symbol = {} |

887 | except KeyError: |

888 | pass |

889 | M = ell_modular_symbols.ModularSymbol(self, sign, normalize) |

890 | self.__modular_symbol[typ] = M |

891 | return M |

892 | |

893 | padic_lseries = padics.padic_lseries |

894 | |

895 | def newform(self): |

896 | r""" |

897 | Same as \code{self.modular_form()}. |

898 | |

899 | EXAMPLES: |

900 | sage: E=EllipticCurve('37a1') |

901 | sage: E.newform() |

902 | q - 2*q^2 - 3*q^3 + 2*q^4 - 2*q^5 + O(q^6) |

903 | sage: E.newform() == E.modular_form() |

904 | True |

905 | """ |

906 | return self.modular_form() |

907 | |

908 | def q_eigenform(self, prec): |

909 | r""" |

910 | Synonym for \code{self.q_expansion(prec)}. |

911 | |

912 | EXAMPLES: |

913 | sage: E=EllipticCurve('37a1') |

914 | sage: E.q_eigenform(10) |

915 | q - 2*q^2 - 3*q^3 + 2*q^4 - 2*q^5 + 6*q^6 - q^7 + 6*q^9 + O(q^10) |

916 | sage: E.q_eigenform(10) == E.q_expansion(10) |

917 | True |

918 | """ |

919 | return self.q_expansion(prec) |

920 | |

921 | def analytic_rank(self, algorithm="cremona"): |

922 | r""" |

923 | Return an integer that is \emph{probably} the analytic rank of |

924 | this elliptic curve. |

925 | |

926 | INPUT: |

927 | algorithm: |

928 | -- 'cremona' (default) -- Use the Buhler-Gross algorithm |

929 | as implemented in GP by Tom Womack and John Cremona, |

930 | who note that their implementation is practical for |

931 | any rank and conductor $\leq 10^{10}$ in 10 minutes. |

932 | |

933 | -- 'sympow' --use Watkins's program sympow |

934 | |

935 | -- 'rubinstein' -- use Rubinstein's L-function C++ program lcalc. |

936 | |

937 | -- 'magma' -- use MAGMA |

938 | |

939 | -- 'all' -- compute with all other free algorithms, check that the |

940 | answers agree, and return the common answer. |

941 | |

942 | \note{If the curve is loaded from the large Cremona database, |

943 | then the modular degree is taken from the database.} |

944 | |

945 | Of the three above, probably Rubinstein's is the most |

946 | efficient (in some limited testing I've done). |

947 | |

948 | \note{It is an open problem to \emph{prove} that \emph{any} |

949 | particular elliptic curve has analytic rank $\geq 4$.} |

950 | |

951 | EXAMPLES: |

952 | sage: E = EllipticCurve('389a') |

953 | sage: E.analytic_rank(algorithm='cremona') |

954 | 2 |

955 | sage: E.analytic_rank(algorithm='rubinstein') |

956 | 2 |

957 | sage: E.analytic_rank(algorithm='sympow') |

958 | 2 |

959 | sage: E.analytic_rank(algorithm='magma') # optional |

960 | 2 |

961 | sage: E.analytic_rank(algorithm='all') |

962 | 2 |

963 | """ |

964 | if algorithm == 'cremona': |

965 | return rings.Integer(gp_cremona.ellanalyticrank(self.minimal_model().a_invariants())) |

966 | elif algorithm == 'rubinstein': |

967 | from sage.lfunctions.lcalc import lcalc |

968 | return lcalc.analytic_rank(L=self) |

969 | elif algorithm == 'sympow': |

970 | from sage.lfunctions.sympow import sympow |

971 | return sympow.analytic_rank(self)[0] |

972 | elif algorithm == 'magma': |

973 | return rings.Integer(self._magma_().AnalyticRank()) |

974 | elif algorithm == 'all': |

975 | S = list(set([self.analytic_rank('cremona'), |

976 | self.analytic_rank('rubinstein'), self.analytic_rank('sympow')])) |

977 | if len(S) != 1: |

978 | raise RuntimeError, "Bug in analytic rank; algorithms don't agree! (E=%s)"%E |

979 | return S[0] |

980 | else: |

981 | raise ValueError, "algorithm %s not defined"%algorithm |

982 | |

983 | def p_isogenous_curves(self, p=None): |

984 | r""" |

985 | Return a list of pairs $(p, L)$ where $p$ is a prime and $L$ |

986 | is a list of the elliptic curves over $\Q$ that are |

987 | $p$-isogenous to this elliptic curve. |

988 | |

989 | INPUT: |

990 | p -- prime or None (default: None); if a prime, returns |

991 | a list of the p-isogenous curves. Otherwise, returns |

992 | a list of all prime-degree isogenous curves sorted |

993 | by isogeny degree. |

994 | |

995 | This is implemented using Cremona's GP script \code{allisog.gp}. |

996 | |

997 | EXAMPLES: |

998 | sage: E = EllipticCurve([0,-1,0,-24649,1355209]) |

999 | sage: E.p_isogenous_curves() |

1000 | [(2, [Elliptic Curve defined by y^2 = x^3 - x^2 - 91809*x - 9215775 over Rational Field, Elliptic Curve defined by y^2 = x^3 - x^2 - 383809*x + 91648033 over Rational Field, Elliptic Curve defined by y^2 = x^3 - x^2 + 1996*x + 102894 over Rational Field])] |

1001 | |

1002 | The isogeny class of the curve 11a2 has three curves in it. |

1003 | But \code{p_isogenous_curves} only returns one curves, since |

1004 | there is only one curve $5$-isogenous to 11a2. |

1005 | sage: E = EllipticCurve('11a2') |

1006 | sage: E.p_isogenous_curves() |

1007 | [(5, [Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field])] |

1008 | sage: E.p_isogenous_curves(5) |

1009 | [Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field] |

1010 | sage: E.p_isogenous_curves(3) |

1011 | [] |

1012 | |

1013 | In contrast, the curve 11a1 admits two $5$-isogenies: |

1014 | sage: E = EllipticCurve('11a1') |

1015 | sage: E.p_isogenous_curves(5) |

1016 | [Elliptic Curve defined by y^2 + y = x^3 - x^2 - 7820*x - 263580 over Rational Field, |

1017 | Elliptic Curve defined by y^2 + y = x^3 - x^2 over Rational Field] |

1018 | """ |

1019 | if p is None: |

1020 | X = eval(gp_cremona.allisog(self.minimal_model().a_invariants())) |

1021 | Y = [(p, [constructor.EllipticCurve(ainvs) for ainvs in L]) for p, L in X] |

1022 | Y.sort() |

1023 | return Y |

1024 | else: |

1025 | X = eval(gp_cremona.p_isog(self.minimal_model().a_invariants(), p)) |

1026 | Y = [constructor.EllipticCurve(ainvs) for ainvs in X] |

1027 | Y.sort() |

1028 | return Y |

1029 | |

1030 | def simon_two_descent(self, verbose=0, lim1=5, lim3=50, limtriv=10, maxprob=20, limbigprime=30): |

1031 | r""" |

1032 | Given a curve with no 2-torsion, computes (probably) the rank |

1033 | of the Mordell-Weil group, with certainty the rank of the |

1034 | 2-Selmer group, and a list of independent points on the curve. |

1035 | |

1036 | INPUT: |

1037 | verbose -- integer, 0,1,2,3; (default: 0), the verbosity level |

1038 | lim1 -- (default: 5) limite des points triviaux sur les quartiques |

1039 | lim3 -- (default: 50) limite des points sur les quartiques ELS |

1040 | limtriv -- (default: 10) limite des points triviaux sur la |

1041 | courbe elliptique |

1042 | maxprob -- (default: 20) |

1043 | limbigprime -- (default: 30) to distinguish between small and large prime |

1044 | numbers. Use probabilistic tests for large |

1045 | primes. If 0, don't any probabilistic tests. |

1046 | |

1047 | OUTPUT: |

1048 | integer -- "probably" the rank of self |

1049 | integer -- the 2-rank of the Selmer group |

1050 | list -- list of independent points on the curve. |

1051 | |

1052 | IMPLEMENTATION: Uses {\bf Denis Simon's} GP/PARI scripts from |

1053 | \url{http://www.math.unicaen.fr/~simon/} |

1054 | |

1055 | EXAMPLES: |

1056 | These computations use pseudo-random numbers, so we set the |

1057 | seed for reproducible testing. |

1058 | |

1059 | We compute the ranks of the curves of lowest known conductor up to rank $8$. |

1060 | Amazingly, each of these computations finishes almost instantly! |

1061 | |

1062 | sage: E = EllipticCurve('11a1') |

1063 | sage: set_random_seed(0) |

1064 | sage: E.simon_two_descent() |

1065 | (0, 0, []) |

1066 | sage: E = EllipticCurve('37a1') |

1067 | sage: set_random_seed(0) |

1068 | sage: E.simon_two_descent() |

1069 | (1, 1, [(0 : 0 : 1)]) |

1070 | sage: E = EllipticCurve('389a1') |

1071 | sage: set_random_seed(0) |

1072 | sage: E.simon_two_descent() |

1073 | (2, 2, [(1 : 0 : 1), (-11/9 : -55/27 : 1)]) |

1074 | sage: E = EllipticCurve('5077a1') |

1075 | sage: set_random_seed(0) |

1076 | sage: E.simon_two_descent() |

1077 | (3, 3, [(1 : 0 : 1), (2 : -1 : 1), (0 : 2 : 1)]) |

1078 | |

1079 | |

1080 | In this example Simon's program does not find any points, though |

1081 | it does correctly compute the rank of the 2-Selmer group. |

1082 | sage: E = EllipticCurve([1, -1, 0, -751055859, -7922219731979]) # long (0.6 seconds) |

1083 | sage: set_random_seed(0) |

1084 | sage: E.simon_two_descent () |

1085 | (1, 1, []) |

1086 | |

1087 | The rest of these entries were taken from Tom Womack's page |

1088 | \url{http://tom.womack.net/maths/conductors.htm} |

1089 | |

1090 | sage: E = EllipticCurve([1, -1, 0, -79, 289]) |

1091 | sage: set_random_seed(0) |

1092 | sage: E.simon_two_descent() |

1093 | (4, 4, [(6 : -5 : 1), (4 : 3 : 1), (5 : -3 : 1), (8 : -15 : 1)]) |

1094 | sage: E = EllipticCurve([0, 0, 1, -79, 342]) |

1095 | sage: set_random_seed(0) |

1096 | sage: E.simon_two_descent() |

1097 | (5, 5, [(5 : 8 : 1), (10 : 23 : 1), (3 : 11 : 1), (4 : -10 : 1), (0 : 18 : 1)]) |

1098 | sage: E = EllipticCurve([1, 1, 0, -2582, 48720]) |

1099 | sage: set_random_seed(0) |

1100 | sage: r, s, G = E.simon_two_descent(); r,s |

1101 | (6, 6) |

1102 | sage: E = EllipticCurve([0, 0, 0, -10012, 346900]) |

1103 | sage: set_random_seed(0) |

1104 | sage: r, s, G = E.simon_two_descent(); r,s |

1105 | (7, 7) |

1106 | sage: E = EllipticCurve([0, 0, 1, -23737, 960366]) |

1107 | sage: set_random_seed(0) |

1108 | sage: r, s, G = E.simon_two_descent(); r,s |

1109 | (8, 8) |

1110 | """ |

1111 | t = simon_two_descent(self, verbose=verbose, lim1=lim1, lim3=lim3, limtriv=limtriv, |

1112 | maxprob=maxprob, limbigprime=limbigprime) |

1113 | prob_rank = rings.Integer(t[0]) |

1114 | two_selmer_rank = rings.Integer(t[1]) |

1115 | prob_gens = [self(P) for P in t[2]] |

1116 | return prob_rank, two_selmer_rank, prob_gens |

1117 | |

1118 | two_descent_simon = simon_two_descent |

1119 | |

1120 | def three_selmer_rank(self, bound=0, method=2): |

1121 | r""" |

1122 | Return the 3-selmer rank of this elliptic curve, computed |

1123 | using Magma. |

1124 | |

1125 | This is not implemented for all curves; a NotImplementedError |

1126 | exception is raised when this function is called on curves |

1127 | for which 3-descent isn't implemented. |

1128 | |

1129 | \note{Use a slightly modified version of Michael Stoll's MAGMA |

1130 | file \code{3descent.m}. You must have Magma to use this |

1131 | function.} |

1132 | |

1133 | EXAMPLES: |

1134 | sage: EllipticCurve('37a').three_selmer_rank() # optional & long -- Magma |

1135 | 1 |

1136 | |

1137 | sage: EllipticCurve('14a1').three_selmer_rank() # optional |

1138 | Traceback (most recent call last): |

1139 | ... |

1140 | NotImplementedError: Currently, only the case with irreducible phi3 is implemented. |

1141 | """ |

1142 | import magma_3descent |

1143 | try: |

1144 | return magma_3descent.three_selmer_rank(self, bound, method) |

1145 | except RuntimeError, msg: |

1146 | msg = str(msg) |

1147 | i = msg.rfind(':') |

1148 | raise NotImplementedError, msg[i+1:] |

1149 | |

1150 | |

1151 | def rank(self, use_database=False, verbose=False, |

1152 | only_use_mwrank=True, |

1153 | algorithm='mwrank_shell', |

1154 | proof=None): |

1155 | """ |

1156 | Return the rank of this elliptic curve, assuming no conjectures. |

1157 | |

1158 | If we fail to provably compute the rank, raises a RuntimeError |

1159 | exception. |

1160 | |

1161 | INPUT: |

1162 | use_database (bool) -- (default: False), if True, try to |

1163 | look up the regulator in the Cremona database. |

1164 | verbose -- (default: None), if specified changes the |

1165 | verbosity of mwrank computations. |

1166 | algorithm -- 'mwrank_shell' -- call mwrank shell command |

1167 | -- 'mwrank_lib' -- call mwrank c library |

1168 | only_use_mwrank -- (default: True) if False try using |

1169 | analytic rank methods first. |

1170 | proof -- bool or None (default: None, see proof.elliptic_curve or |

1171 | sage.structure.proof). Note that results |

1172 | obtained from databases are considered proof = True |

1173 | |

1174 | OUTPUT: |

1175 | rank (int) -- the rank of the elliptic curve. |

1176 | |

1177 | IMPLEMENTATION: Uses L-functions, mwrank, and databases. |

1178 | |

1179 | EXAMPLES: |

1180 | sage: EllipticCurve('11a').rank() |

1181 | 0 |

1182 | sage: EllipticCurve('37a').rank() |

1183 | 1 |

1184 | sage: EllipticCurve('389a').rank() |

1185 | 2 |

1186 | sage: EllipticCurve('5077a').rank() |

1187 | 3 |

1188 | sage: EllipticCurve([1, -1, 0, -79, 289]).rank() # long time. This will use the default proof behavior of True. |

1189 | 4 |

1190 | sage: EllipticCurve([0, 0, 1, -79, 342]).rank(proof=False) # long time -- but under a minute |

1191 | 5 |

1192 | sage: EllipticCurve([0, 0, 1, -79, 342]).simon_two_descent()[0] # much faster -- almost instant. |

1193 | 5 |

1194 | |

1195 | Examples with denominators in defining equations: |

1196 | sage: E = EllipticCurve( [0, 0, 0, 0, -675/4]) |

1197 | sage: E.rank() |

1198 | 0 |

1199 | sage: E = EllipticCurve( [0, 0, 1/2, 0, -1/5]) |

1200 | sage: E.rank() |

1201 | 1 |

1202 | sage: E.minimal_model().rank() |

1203 | 1 |

1204 | """ |

1205 | if proof is None: |

1206 | from sage.structure.proof.proof import get_flag |

1207 | proof = get_flag(proof, "elliptic_curve") |

1208 | else: |

1209 | proof = bool(proof) |

1210 | try: |

1211 | return self.__rank[proof] |

1212 | except KeyError: |

1213 | if proof is False and self.__rank.has_key(True): |

1214 | return self.__rank[True] |

1215 | if use_database: |

1216 | try: |

1217 | self.__rank[True] = self.database_curve().rank() |

1218 | return self.__rank[True] |

1219 | except (AttributeError, RuntimeError): |

1220 | pass |

1221 | if not only_use_mwrank: |

1222 | N = self.conductor() |

1223 | prec = int(4*float(sqrt(N))) + 10 |

1224 | if self.root_number() == 1: |

1225 | L, err = self.lseries().at1(prec) |

1226 | if abs(L) > err + R(0.0001): # definitely doesn't vanish |

1227 | misc.verbose("rank 0 because L(E,1)=%s"%L) |

1228 | self.__rank[proof] = 0 |

1229 | return self.__rank[proof] |

1230 | else: |

1231 | Lprime, err = self.lseries().deriv_at1(prec) |

1232 | if abs(Lprime) > err + R(0.0001): # definitely doesn't vanish |

1233 | misc.verbose("rank 1 because L'(E,1)=%s"%Lprime) |

1234 | self.__rank[proof] = 1 |

1235 | return self.__rank[proof] |

1236 | |

1237 | if algorithm == 'mwrank_lib': |

1238 | misc.verbose("using mwrank lib") |

1239 | C = self.mwrank_curve() |

1240 | C.set_verbose(verbose) |

1241 | r = C.rank() |

1242 | if not C.certain(): |

1243 | del self.__mwrank_curve |

1244 | raise RuntimeError, "Unable to compute the rank with certainty (lower bound=%s). This could be because Sha(E/Q)[2] is nontrivial."%C.rank() + "\nTrying calling something like two_descent(second_limit=13) on the curve then trying this command again. You could also try rank with only_use_mwrank=False." |

1245 | self.__rank[proof] = r |

1246 | elif algorithm == 'mwrank_shell': |

1247 | misc.verbose("using mwrank shell") |

1248 | X = self.mwrank() |

1249 | if not 'The rank and full Mordell-Weil basis have been determined unconditionally' in X: |

1250 | if proof: |

1251 | raise RuntimeError, '%s\nRank not provably correct.'%X |

1252 | else: |

1253 | misc.verbose("Warning -- rank not provably correct", level=1) |

1254 | elif proof is False: |

1255 | proof = True #since we actually provably found the rank |

1256 | i = X.find('Rank = ') |

1257 | assert i != -1 |

1258 | j = i + X[i:].find('\n') |

1259 | self.__rank[proof] = Integer(X[i+7:j]) |

1260 | return self.__rank[proof] |

1261 | |

1262 | def gens(self, verbose=False, rank1_search=10, |

1263 | algorithm='mwrank_shell', |

1264 | only_use_mwrank=True, |

1265 | proof = None): |

1266 | """ |

1267 | Compute and return generators for the Mordell-Weil group E(Q) |

1268 | *modulo* torsion. |

1269 | |

1270 | HINT: If you would like to control the height bounds used in |

1271 | the 2-descent, first call the two_descent function with those |

1272 | height bounds. However that function, while it displays a lot |

1273 | of output, returns no values. |

1274 | |

1275 | TODO: (1) Right now this function assumes that the input curve |

1276 | is in minimal Weierstrass form. This restriction will be |

1277 | removed in the future. This function raises a |

1278 | NotImplementedError if a non-minimal curve is given as input. |

1279 | |

1280 | (2) Allow passing of command-line parameters to mwrank. |

1281 | |

1282 | WARNING: If the program fails to give a provably correct |

1283 | result, it prints a warning message, but does not raise an |

1284 | exception. Use the gens_certain command to find out if |

1285 | this warning message was printed. |

1286 | |

1287 | INPUT: |

1288 | verbose -- (default: None), if specified changes the |

1289 | verbosity of mwrank computations. |

1290 | rank1_search -- (default: 16), if the curve has analytic |

1291 | rank 1, try to find a generator by a direct |

1292 | search up to this logarithmic height. If this |

1293 | fails the usual mwrank procedure is called. |

1294 | algorithm -- 'mwrank_shell' (default) -- call mwrank shell command |

1295 | -- 'mwrank_lib' -- call mwrank c library |

1296 | only_use_mwrank -- bool (default True) if false, attempts to |

1297 | first use more naive, natively implemented methods. |

1298 | proof -- bool or None (default None, see proof.elliptic_curve or |

1299 | sage.structure.proof). |

1300 | OUTPUT: |

1301 | generators -- List of generators for the Mordell-Weil group. |

1302 | |

1303 | IMPLEMENTATION: Uses Cremona's mwrank C library. |

1304 | |

1305 | EXAMPLES: |

1306 | sage: E = EllipticCurve('389a') |

1307 | sage: E.gens() # random output |

1308 | [(-1 : 1 : 1), (0 : 0 : 1)] |

1309 | """ |

1310 | if proof is None: |

1311 | from sage.structure.proof.proof import get_flag |

1312 | proof = get_flag(proof, "elliptic_curve") |

1313 | else: |

1314 | proof = bool(proof) |

1315 | try: |

1316 | return list(self.__gens[proof]) # return copy so not changed |

1317 | except KeyError: |

1318 | if proof is False and self.__gens.has_key(True): |

1319 | return self.__gens[True] |

1320 | if self.conductor() > 10**7: |

1321 | only_use_mwrank = True |

1322 | |

1323 | if not only_use_mwrank: |

1324 | try: |

1325 | misc.verbose("Trying to compute rank.") |

1326 | r = self.rank(only_use_mwrank = False) |

1327 | misc.verbose("Got r = %s."%r) |

1328 | if r == 0: |

1329 | misc.verbose("Rank = 0, so done.") |

1330 | self.__gens[True] = [] |

1331 | self.__regulator[True] = R(1) |

1332 | return self.__gens[True] |

1333 | if r == 1 and rank1_search: |

1334 | misc.verbose("Rank = 1, so using direct search.") |

1335 | h = 6 |

1336 | while h <= rank1_search: |

1337 | misc.verbose("Trying direct search up to height %s"%h) |

1338 | G = self.point_search(h, verbose) |

1339 | G = [P for P in G if P.order() == oo] |

1340 | if len(G) > 0: |

1341 | misc.verbose("Direct search succeeded.") |

1342 | G, _, reg = self.saturation(G, verbose=verbose) |

1343 | misc.verbose("Computed saturation.") |

1344 | self.__gens[True] = G |

1345 | self.__regulator[True] = reg |

1346 | return self.__gens[True] |

1347 | h += 2 |

1348 | misc.verbose("Direct search FAILED.") |

1349 | except RuntimeError: |

1350 | pass |

1351 | # end if (not_use_mwrank) |

1352 | if not self.is_integral(): |

1353 | raise NotImplementedError, "gens via mwrank only implemented for curves with integer coefficients." |

1354 | if algorithm == "mwrank_lib": |

1355 | misc.verbose("Calling mwrank C++ library.") |

1356 | C = self.mwrank_curve(verbose) |

1357 | if not (verbose is None): |

1358 | C.set_verbose(verbose) |

1359 | G = C.gens() |

1360 | if proof is True and C.certain() is False: |

1361 | del self.__mwrank_curve |

1362 | raise RuntimeError, "Unable to compute the rank, hence generators, with certainty (lower bound=%s). This could be because Sha(E/Q)[2] is nontrivial."%C.rank() + \ |

1363 | "\nTrying calling something like two_descent(second_limit=13) on the curve then trying this command again." |

1364 | else: |

1365 | proof = C.certain() |

1366 | else: |

1367 | X = self.mwrank() |

1368 | misc.verbose("Calling mwrank shell.") |

1369 | if not 'The rank and full Mordell-Weil basis have been determined unconditionally' in X: |

1370 | msg = 'Generators not provably computed.' |

1371 | if proof: |

1372 | raise RuntimeError, '%s\n%s'%(X,msg) |

1373 | else: |

1374 | misc.verbose("Warning -- %s"%msg, level=1) |

1375 | elif proof is False: |

1376 | proof = True |

1377 | G = [] |

1378 | i = X.find('Generator ') |

1379 | while i != -1: |

1380 | j = i + X[i:].find(';') |

1381 | k = i + X[i:].find('[') |

1382 | G.append(eval(X[k:j].replace(':',','))) |

1383 | X = X[j:] |

1384 | i = X.find('Generator ') |

1385 | i = X.find('Regulator = ') |

1386 | j = i + X[i:].find('\n') |

1387 | self.__regulator[proof] = R(X[i+len('Regulator = '):j]) |

1388 | #### |

1389 | self.__gens[proof] = [self.point(x, check=True) for x in G] |

1390 | self.__gens[proof].sort() |

1391 | self.__rank[proof] = len(self.__gens[proof]) |

1392 | return self.__gens[proof] |

1393 | |

1394 | def gens_certain(self): |

1395 | """ |

1396 | Return True if the generators have been proven correct. |

1397 | |

1398 | EXAMPLES: |

1399 | sage: E=EllipticCurve('37a1') |

1400 | sage: E.gens() |

1401 | [(0 : -1 : 1)] |

1402 | sage: E.gens_certain() |

1403 | True |

1404 | """ |

1405 | return self.__gens.has_key(True) |

1406 | |

1407 | def ngens(self, proof = None): |

1408 | """ |

1409 | Return the number of generators of this elliptic curve. |

1410 | |

1411 | NOTE: See gens() for further documentation. The function |

1412 | ngens() calls gens() if not already done, but only with |

1413 | default parameters. Better results may be obtained by calling |

1414 | mwrank() with carefully chosen parameters. |

1415 | |

1416 | EXAMPLES: |

1417 | sage: E=EllipticCurve('37a1') |

1418 | sage: E.ngens() |

1419 | 1 |

1420 | |

1421 | TO DO: This example should not cause a run-time error. |

1422 | sage: E=EllipticCurve([0,0,0,877,0]) |

1423 | sage: # E.ngens() ######## causes run-time error |

1424 | |

1425 | sage: print E.mwrank('-v0 -b12 -l') |

1426 | Curve [0,0,0,877,0] : Rank = 1 |

1427 | Generator 1 is [29604565304828237474403861024284371796799791624792913256602210:-256256267988926809388776834045513089648669153204356603464786949:490078023219787588959802933995928925096061616470779979261000]; height 95.980371987964 |

1428 | Regulator = 95.980... |

1429 | |

1430 | |

1431 | """ |

1432 | return len(self.gens(proof = proof)) |

1433 | |

1434 | def regulator(self, use_database=True, verbose=None, proof=None): |

1435 | """ |

1436 | Returns the regulator of this curve, which must be defined |

1437 | over Q. |

1438 | |

1439 | INPUT: |

1440 | use_database -- bool (default: False), if True, try to |

1441 | look up the regulator in the Cremona database. |

1442 | verbose -- (default: None), if specified changes the |

1443 | verbosity of mwrank computations. |

1444 | proof -- bool or None (default: None, see proof.[tab] or |

1445 | sage.structure.proof). Note that results from |

1446 | databases are considered proof = True |

1447 | |

1448 | EXAMPLES: |

1449 | sage: E = EllipticCurve([0, 0, 1, -1, 0]) |

1450 | sage: E.regulator() # long time (1 second) |

1451 | 0.0511114082399688 |

1452 | sage: EllipticCurve('11a').regulator() |

1453 | 1.00000000000000 |

1454 | sage: EllipticCurve('37a').regulator() |

1455 | 0.0511114082399688 |

1456 | sage: EllipticCurve('389a').regulator() |

1457 | 0.152460177943144 |

1458 | sage: EllipticCurve('5077a').regulator() # random low order bit |

1459 | 0.417143558758385 |

1460 | sage: EllipticCurve([1, -1, 0, -79, 289]).regulator() # long time (seconds) |

1461 | 1.50434488827528 |

1462 | sage: EllipticCurve([0, 0, 1, -79, 342]).regulator(proof=False) # long time (seconds) |

1463 | 14.7905275701310 |

1464 | """ |

1465 | if proof is None: |

1466 | from sage.structure.proof.proof import get_flag |

1467 | proof = get_flag(proof, "elliptic_curve") |

1468 | else: |

1469 | proof = bool(proof) |

1470 | try: |

1471 | return self.__regulator[proof] |

1472 | except KeyError: |

1473 | if proof is False and self.__regulator.has_key(True): |

1474 | return self.__regulator[True] |

1475 | if use_database: |

1476 | try: |

1477 | self.__regulator[True] = R(self.database_curve().db_extra[3]) |

1478 | return self.__regulator[True] |

1479 | except (AttributeError, RuntimeError): |

1480 | pass |

1481 | G = self.gens(proof=proof) |

1482 | try: # in some cases self.gens() efficiently computes regulator. |

1483 | return self.__regulator[proof] |

1484 | except KeyError: |

1485 | if proof is False and self.__regulator.has_key(True): |

1486 | return self.__regulator[True] |

1487 | C = self.mwrank_curve() |

1488 | reg = R(C.regulator()) |

1489 | if proof is True and not C.certain(): |

1490 | raise RuntimeError, "Unable to compute the rank, hence regulator, with certainty (lower bound=%s)."%C.rank() |

1491 | proof = C.certain() |

1492 | self.__regulator[proof] = reg |

1493 | return self.__regulator[proof] |

1494 | |

1495 | def saturation(self, points, verbose=False, max_prime=0, odd_primes_only=False): |

1496 | """ |

1497 | Given a list of rational points on E, compute the saturation |

1498 | in E(Q) of the subgroup they generate. |

1499 | |

1500 | INPUT: |

1501 | points (list) -- list of points on E |

1502 | verbose (bool) -- (default: False), if True, give verbose output |

1503 | max_prime (int) -- (default: 0), saturation is performed |

1504 | for all primes up to max_prime. If max_prime==0, |

1505 | perform saturation at *all* primes, i.e., compute |

1506 | the true saturation. |

1507 | odd_primes_only (bool) -- only do saturation at odd primes |

1508 | |

1509 | OUTPUT: |

1510 | saturation (list) -- points that form a basis for the saturation |

1511 | index (int) -- the index of the group generated by points |

1512 | in their saturation |

1513 | regulator (float) -- regulator of saturated points. |

1514 | |

1515 | IMPLEMENTATION: |

1516 | Uses Cremona's mwrank package. With max_prime=0, we call |

1517 | mwrank with successively larger prime bounds until the |

1518 | full saturation is provably found. The results of |

1519 | saturation at the previous primes is stored in each case, |

1520 | so this should be reasonably fast. |

1521 | |

1522 | EXAMPLES: |

1523 | sage: E=EllipticCurve('37a1') |

1524 | sage: P=E.gens()[0] |

1525 | sage: Q=5*P; Q |

1526 | (1/4 : -3/8 : 1) |

1527 | sage: E.saturation([Q]) |

1528 | ([(0 : -1 : 1)], '5', 0.0511114075779915) |

1529 | """ |

1530 | if not isinstance(points, list): |

1531 | raise TypeError, "points (=%s) must be a list."%points |

1532 | |

1533 | v = [] |

1534 | for P in points: |

1535 | if not isinstance(P, ell_point.EllipticCurvePoint_field): |

1536 | P = self(P) |

1537 | elif P.curve() != self: |

1538 | raise ArithmeticError, "point (=%s) must be %s."%(P,self) |

1539 | x, y = P.xy() |

1540 | d = x.denominator().lcm(y.denominator()) |

1541 | v.append((x*d, y*d, d)) |

1542 | c = self.mwrank_curve() |

1543 | mw = mwrank.mwrank_MordellWeil(c, verbose) |

1544 | mw.process(v) |

1545 | if max_prime == 0: |

1546 | repeat_until_saturated = True |

1547 | max_prime = 97 |

1548 | while True: |

1549 | ok, index, unsat = mw.saturate(max_prime=max_prime, odd_primes_only = odd_primes_only) |

1550 | reg = mw.regulator() |

1551 | if not ok and repeat_until_saturated: |

1552 | max_prime = arith.next_prime(max_prime + 100) |

1553 | ok, index, unsat = mw.saturate(max_prime=max_prime, odd_primes_only = odd_primes_only) |

1554 | reg = mw.regulator() |

1555 | else: |

1556 | break |

1557 | sat = mw.points() |

1558 | sat = [self(P) for P in sat] |

1559 | return sat, index, R(reg) |

1560 | |

1561 | def CPS_height_bound(self): |

1562 | """ |

1563 | Return the Cremona-Prickett-Siksek height bound. This is a |

1564 | floating point number B such that if P is a point on the curve, |

1565 | then the naive logarithmetic height of P is off from the |

1566 | canonical height by at most B. |

1567 | |

1568 | EXAMPLES: |

1569 | sage: E = EllipticCurve("11a") |

1570 | sage: E.CPS_height_bound() |

1571 | 2.8774743273580445 |

1572 | sage: E = EllipticCurve("5077a") |

1573 | sage: E.CPS_height_bound() |

1574 | 0.0 |

1575 | sage: E = EllipticCurve([1,2,3,4,1]) |

1576 | sage: E.CPS_height_bound() |

1577 | Traceback (most recent call last): |

1578 | ... |

1579 | RuntimeError: curve must be minimal. |

1580 | sage: F = E.quadratic_twist(-19) |

1581 | sage: F |

1582 | Elliptic Curve defined by y^2 + x*y + y = x^3 - x^2 + 1376*x - 130 over Rational Field |

1583 | sage: F.CPS_height_bound() |

1584 | 0.65551583769728516 |

1585 | |

1586 | IMPLEMENTATION: |

1587 | Call the corresponding mwrank C++ library function. |

1588 | """ |

1589 | if not self.is_minimal(): |

1590 | raise RuntimeError, "curve must be minimal." |

1591 | return self.mwrank_curve().CPS_height_bound() |

1592 | |

1593 | |

1594 | def silverman_height_bound(self): |

1595 | r""" Return the Silverman height bound. This is a positive real |

1596 | (floating point) number B such that for all rational points |

1597 | $P$ on the curve, |

1598 | $$ |

1599 | h(P) \le \hat{h}(P) + B |

1600 | $$ |

1601 | where h(P) is the logarithmic height of $P$ and $\hat{h}(P)$ |

1602 | is the canonical height. |

1603 | |

1604 | Note that the CPS_height_bound is often better (i.e. smaller) |

1605 | than the Silverman bound. |

1606 | |

1607 | EXAMPLES: |

1608 | sage: E=EllipticCurve('37a1') |

1609 | sage: E.silverman_height_bound() |

1610 | 4.8254007581809182 |

1611 | sage: E.CPS_height_bound() |

1612 | 0.16397076103046915 |

1613 | """ |

1614 | return self.mwrank_curve().silverman_bound() |

1615 | |

1616 | |

1617 | def point_search(self, height_limit, verbose=True): |

1618 | """ |

1619 | Search for points on a curve up to an input bound on the naive |

1620 | logarithmic height. |

1621 | |

1622 | INPUT: |

1623 | |

1624 | height_limit (float) -- bound on naive height (at most 21, |

1625 | or mwrank overflows -- see below) |

1626 | |

1627 | verbose (bool) -- (default: True) |

1628 | |

1629 | If True, report on each point as found |

1630 | together with linear relations between |

1631 | the points found and the saturation |

1632 | process. |

1633 | |

1634 | If False, just return the result. |

1635 | |

1636 | OUTPUT: points (list) -- list of independent points which |

1637 | generate the subgroup of the Mordell-Weil group generated |

1638 | by the points found and then p-saturated for p<20. |

1639 | |

1640 | WARNING: |

1641 | height_limit is logarithmic, so increasing by 1 will cause |

1642 | the running time to increase by a factor of approximately |

1643 | 4.5 (=exp(1.5)). The limit of 21 is to prevent overflow, |

1644 | but in any case using height_limit=20 takes rather a long |

1645 | time! |

1646 | |

1647 | IMPLEMENTATION: Uses Cremona's mwrank package. At the heart of |

1648 | this function is Cremona's port of Stoll's ratpoints program |

1649 | (version 1.4). |

1650 | |

1651 | EXAMPLES: |

1652 | sage: E=EllipticCurve('389a1') |

1653 | sage: E.point_search(5, verbose=False) |

1654 | [(0 : -1 : 1), (-1 : 1 : 1)] |

1655 | |

1656 | Increasing the height_limit takes longer, but finds no more points: |

1657 | sage: E.point_search(10, verbose=False) |

1658 | [(0 : -1 : 1), (-1 : 1 : 1)] |

1659 | |

1660 | In fact this curve has rank 2 so no more than 2 points will |

1661 | ever be output, but we are not using this fact. |

1662 | |

1663 | sage: E.saturation(_) |

1664 | ([(0 : -1 : 1), (-1 : 1 : 1)], '1', 0.152460172772408) |

1665 | |

1666 | What this shows is that if the rank is 2 then the points |

1667 | listed do generate the Mordell-Weil group (mod torsion). |

1668 | Finally, |

1669 | |

1670 | sage: E.rank() |

1671 | 2 |

1672 | |

1673 | """ |

1674 | height_limit = float(height_limit) |

1675 | if height_limit > _MAX_HEIGHT: |

1676 | raise OverflowError, "height_limit (=%s) must be at most %s."%(height_limit,_MAX_HEIGHT) |

1677 | c = self.mwrank_curve() |

1678 | mw = mwrank.mwrank_MordellWeil(c, verbose) |

1679 | mw.search(height_limit, verbose=verbose) |

1680 | v = mw.points() |

1681 | return [self(P) for P in v] |

1682 | |

1683 | def two_torsion_rank(self): |

1684 | r""" |

1685 | Return the dimension of the 2-torsion subgroup of $E(\Q)$. |

1686 | |

1687 | This will be 0, 1 or 2. |

1688 | |

1689 | NOTE: As s side-effect of calling this function, the full |

1690 | torsion subgroup of the curve is computed (if not already |

1691 | cached). A simpler implementation of this function would be |

1692 | possible (by counting the roots of the 2-division polynomial), |

1693 | but the full torsion subgroup computation is not expensive. |

1694 | |

1695 | EXAMPLES: |

1696 | sage: EllipticCurve('11a1').two_torsion_rank() |

1697 | 0 |

1698 | sage: EllipticCurve('14a1').two_torsion_rank() |

1699 | 1 |

1700 | sage: EllipticCurve('15a1').two_torsion_rank() |

1701 | 2 |

1702 | """ |

1703 | A = self.torsion_subgroup().invariants() |

1704 | if len(A) == 2: |

1705 | return rings.Integer(2) |

1706 | elif len(A) == 1 and A[0] % 2 == 0: |

1707 | return rings.Integer(1) |

1708 | else: |

1709 | return rings.Integer(0) |

1710 | |

1711 | def selmer_rank_bound(self): |

1712 | """ |

1713 | Bound on the rank of the curve, computed using the |

1714 | 2-selmer group. This is the rank of the curve |

1715 | minus the rank of the 2-torsion, minus a number |

1716 | determined by whatever mwrank was able to determine |

1717 | related to Sha[2]. Thus in many cases, this is |

1718 | the actual rank of the curve. |

1719 | |

1720 | EXAMPLE: |

1721 | The following is the curve 960D1, which has rank 0, |

1722 | but Sha of order 4. |

1723 | sage: E = EllipticCurve([0, -1, 0, -900, -10098]) |

1724 | sage: E.selmer_rank_bound() |

1725 | 0 |

1726 | |

1727 | It gives 0 instead of 2, because it knows Sha is nontrivial. |

1728 | In contrast, for the curve 571A, also with rank 0 and Sha |

1729 | of order 4, we get a worse bound: |

1730 | sage: E = EllipticCurve([0, -1, 1, -929, -10595]) |

1731 | sage: E.selmer_rank_bound() |

1732 | 2 |

1733 | sage: E.rank(only_use_mwrank=False) # uses L-function |

1734 | 0 |

1735 | """ |

1736 | try: |

1737 | return self.__selmer_rank_bound |

1738 | except AttributeError: |

1739 | C = self.mwrank_curve() |

1740 | self.__selmer_rank_bound = C.selmer_rank_bound() |

1741 | return self.__selmer_rank_bound |

1742 | |

1743 | |

1744 | def an(self, n): |

1745 | """ |

1746 | The n-th Fourier coefficient of the modular form corresponding |

1747 | to this elliptic curve, where n is a positive integer. |

1748 | |

1749 | EXAMPLES: |

1750 | sage: E=EllipticCurve('37a1') |

1751 | sage: [E.an(n) for n in range(20) if n>0] |

1752 | [1, -2, -3, 2, -2, 6, -1, 0, 6, 4, -5, -6, -2, 2, 6, -4, 0, -12, 0] |

1753 | """ |

1754 | return Integer(self.pari_mincurve().ellak(n)) |

1755 | |

1756 | def ap(self, p): |

1757 | """ |

1758 | The p-th Fourier coefficient of the modular form corresponding |

1759 | to this elliptic curve, where p is prime. |

1760 | |

1761 | EXAMPLES: |

1762 | sage: E=EllipticCurve('37a1') |

1763 | sage: [E.ap(p) for p in prime_range(50)] |

1764 | [-2, -3, -2, -1, -5, -2, 0, 0, 2, 6, -4, -1, -9, 2, -9] |

1765 | |

1766 | """ |

1767 | if not arith.is_prime(p): |

1768 | raise ArithmeticError, "p must be prime" |

1769 | return Integer(self.pari_mincurve().ellap(p)) |

1770 | |

1771 | def quadratic_twist(self, D): |

1772 | """ |

1773 | Return the quadratic twist of this elliptic curve by D. |

1774 | |

1775 | D must be a nonzero rational number. |

1776 | |

1777 | NOTE: This function overrides the generic quadratic_twist() |

1778 | function for elliptic curves, returning a minimal model |

1779 | |

1780 | EXAMPLES: |

1781 | sage: E=EllipticCurve('37a1') |

1782 | sage: E2=E.quadratic_twist(2); E2 |

1783 | Elliptic Curve defined by y^2 = x^3 - 4*x + 2 over Rational Field |

1784 | sage: E2.conductor() |

1785 | 2368 |

1786 | sage: E2.quadratic_twist(2) == E |

1787 | True |

1788 | """ |

1789 | return EllipticCurve_number_field.quadratic_twist(self, D).minimal_model() |

1790 | |

1791 | def minimal_model(self): |

1792 | r""" |

1793 | Return the unique minimal Weierstrass equation for this |

1794 | elliptic curve. This is the model with minimal discriminant |

1795 | and $a_1,a_2,a_3 \in \{0,\pm 1\}$. |

1796 | |

1797 | EXAMPLES: |

1798 | sage: E=EllipticCurve([10,100,1000,10000,1000000]) |

1799 | sage: E.minimal_model() |

1800 | Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 + x +1 over Rational Field |

1801 | """ |

1802 | try: |

1803 | return self.__minimal_model |

1804 | except AttributeError: |

1805 | F = self.pari_mincurve() |

1806 | self.__minimal_model = EllipticCurve_rational_field([Q(F[i]) for i in range(5)]) |

1807 | return self.__minimal_model |

1808 | |

1809 | def is_minimal(self): |

1810 | r""" |

1811 | Return True iff this elliptic curve is a reduced minimal model. |

1812 | |

1813 | the unique minimal Weierstrass equation for this |

1814 | elliptic curve. This is the model with minimal discriminant |

1815 | and $a_1,a_2,a_3 \in \{0,\pm 1\}$. |

1816 | |

1817 | TO DO: This is not very efficient since it just computes the |

1818 | minimal model and compares. A better implementation using the |

1819 | Kraus conditions would be preferable. |

1820 | |

1821 | EXAMPLES: |

1822 | sage: E=EllipticCurve([10,100,1000,10000,1000000]) |

1823 | sage: E.is_minimal() |

1824 | False |

1825 | sage: E=E.minimal_model() |

1826 | sage: E.is_minimal() |

1827 | True |

1828 | """ |

1829 | return self.ainvs() == self.minimal_model().ainvs() |

1830 | |

1831 | # Duplicate! |

1832 | # |

1833 | # def is_integral(self): |

1834 | # for n in self.ainvs(): |

1835 | # if n.denominator() != 1: |

1836 | # return False |

1837 | # return True |

1838 | |

1839 | def kodaira_type(self, p): |

1840 | """ |

1841 | Local Kodaira type of the elliptic curve at $p$. |

1842 | |

1843 | INPUT: |

1844 | -- p, an integral prime |

1845 | OUTPUT: |

1846 | -- the kodaira type of this elliptic curve at p, as a KodairaSymbol. |

1847 | |

1848 | EXAMPLES: |

1849 | sage: E = EllipticCurve('124a') |

1850 | sage: E.kodaira_type(2) |

1851 | IV |

1852 | """ |

1853 | if not arith.is_prime(p): |

1854 | raise ArithmeticError, "p must be prime" |

1855 | try: |

1856 | self.__kodaira_type |

1857 | except AttributeError: |

1858 | self.__kodaira_type = {} |

1859 | self.__tamagawa_number = {} |

1860 | if not self.__kodaira_type.has_key(p): |

1861 | v = self.pari_mincurve().elllocalred(p) |

1862 | from kodaira_symbol import KodairaSymbol |

1863 | self.__kodaira_type[p] = KodairaSymbol(v[1]) |

1864 | self.__tamagawa_number[p] = Integer(v[3]) |

1865 | return self.__kodaira_type[p] |

1866 | |

1867 | def tamagawa_number(self, p): |

1868 | """ |

1869 | The Tamagawa number of the elliptic curve at $p$. |

1870 | |

1871 | EXAMPLES: |

1872 | sage: E = EllipticCurve('11a') |

1873 | sage: E.tamagawa_number(11) |

1874 | 5 |

1875 | sage: E = EllipticCurve('37b') |

1876 | sage: E.tamagawa_number(37) |

1877 | 3 |

1878 | """ |

1879 | if not arith.is_prime(p): |

1880 | raise ArithmeticError, "p must be prime" |

1881 | try: |

1882 | return self.__tamagawa_number[p] |

1883 | except (AttributeError, KeyError): |

1884 | self.kodaira_type(p) |

1885 | return self.__tamagawa_number[p] |

1886 | |

1887 | def tamagawa_numbers(self): |

1888 | """ |

1889 | Return a list of all Tamagawa numbers for all prime divisors of |

1890 | the conductor (in order). |

1891 | |

1892 | EXAMPLES: |

1893 | sage: e = EllipticCurve('30a1') |

1894 | sage: e.tamagawa_numbers() |

1895 | [2, 3, 1] |

1896 | sage: vector(e.tamagawa_numbers()) |

1897 | (2, 3, 1) |

1898 | """ |

1899 | return [self.tamagawa_number(p) for p in arith.prime_divisors(self.conductor())] |

1900 | |

1901 | def tamagawa_product(self): |

1902 | """ |

1903 | Returns the product of the Tamagawa numbers. |

1904 | |

1905 | EXAMPLES: |

1906 | sage: E = EllipticCurve('54a') |

1907 | sage: E.tamagawa_product () |

1908 | 3 |

1909 | """ |

1910 | try: |

1911 | return self.__tamagawa_product |

1912 | except AttributeError: |

1913 | self.__tamagawa_product = self.pari_mincurve().ellglobalred()[2].python() |

1914 | return self.__tamagawa_product |

1915 | |

1916 | def real_components(self): |

1917 | """ |

1918 | Returns 1 if there is 1 real component and 2 if there are 2. |

1919 | |

1920 | EXAMPLES: |

1921 | sage: E = EllipticCurve('37a') |

1922 | sage: E.real_components () |

1923 | 2 |

1924 | sage: E = EllipticCurve('37b') |

1925 | sage: E.real_components () |

1926 | 2 |

1927 | sage: E = EllipticCurve('11a') |

1928 | sage: E.real_components () |

1929 | 1 |

1930 | """ |

1931 | invs = self.short_weierstrass_model().ainvs() |

1932 | x = rings.polygen(self.base_ring()) |

1933 | f = x**3 + invs[3]*x + invs[4] |

1934 | if f.discriminant() > 0: |

1935 | return 2 |

1936 | else: |

1937 | return 1 |

1938 | |

1939 | def period_lattice(self): |

1940 | r""" |

1941 | Returns the period lattice of the elliptic curve. |

1942 | |

1943 | EXAMPLES: |

1944 | sage: E = EllipticCurve('37a') |

1945 | sage: E.period_lattice() |

1946 | Period lattice associated to Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field |

1947 | """ |

1948 | from sage.schemes.elliptic_curves.period_lattice import PeriodLattice_ell |

1949 | return PeriodLattice_ell(self) |

1950 | |

1951 | def lseries(self): |

1952 | """ |

1953 | Returns the L-series of this elliptic curve. |

1954 | |

1955 | Further documentation is available for the functions which |

1956 | apply to the L-series. |

1957 | |

1958 | EXAMPLES: |

1959 | sage: E=EllipticCurve('37a1') |

1960 | sage: E.lseries() |

1961 | Complex L-series of the Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field |

1962 | """ |

1963 | try: |

1964 | return self.__lseries |

1965 | except AttributeError: |

1966 | from lseries_ell import Lseries_ell |

1967 | self.__lseries = Lseries_ell(self) |

1968 | return self.__lseries |

1969 | |

1970 | def Lambda(self, s, prec): |

1971 | r""" |

1972 | Returns the value of the Lambda-series of the elliptic curve E |

1973 | at s, where s can be any complex number. |

1974 | |

1975 | IMPLEMENTATION: Fairly *slow* computation using the definitions |

1976 | and implemented in Python. |

1977 | |

1978 | Uses prec terms of the power series. |

1979 | |

1980 | EXAMPLES: |

1981 | sage: E = EllipticCurve('389a') |

1982 | sage: E.Lambda(1.4+0.5*I, 50) |

1983 | -0.354172680517... + 0.874518681720...*I |

1984 | """ |

1985 | s = C(s) |

1986 | N = self.conductor() |

1987 | pi = R(constants.pi) |

1988 | Gamma = transcendental.gamma |

1989 | Gamma_inc = transcendental.gamma_inc |

1990 | a = self.anlist(prec) |

1991 | eps = self.root_number() |

1992 | sqrtN = float(N.sqrt()) |

1993 | def _F(n, t): |

1994 | return Gamma_inc(t+1, 2*pi*n/sqrtN) * C(sqrtN/(2*pi*n))**(t+1) |

1995 | return sum([a[n]*(_F(n,s-1) + eps*_F(n,1-s)) for n in xrange(1,prec+1)]) |

1996 | |

1997 | def is_local_integral_model(self,*p): |

1998 | r""" |

1999 | Tests if self is integral at the prime $p$, or at all the |

2000 | primes if $p$ is a list or tuple of primes |

2001 | |

2002 | EXAMPLES: |

2003 | sage: E=EllipticCurve([1/2,1/5,1/5,1/5,1/5]) |

2004 | sage: [E.is_local_integral_model(p) for p in (2,3,5)] |

2005 | [False, True, False] |

2006 | sage: E.is_local_integral_model(2,3,5) |

2007 | False |

2008 | sage: Eint2=E.local_integral_model(2) |

2009 | sage: Eint2.is_local_integral_model(2) |

2010 | True |

2011 | """ |

2012 | if len(p)==1: p=p[0] |

2013 | if isinstance(p,(tuple,list)): |

2014 | return misc.forall(p, lambda x : self.is_local_integral_model(x))[0] |

2015 | assert p.is_prime(), "p must be prime in is_local_integral_model()" |

2016 | return misc.forall(self.ainvs(), lambda x : x.valuation(p) >= 0)[0] |

2017 | |

2018 | def local_integral_model(self,p): |

2019 | r""" |

2020 | Return a model of self which is integral at the prime $p$ |

2021 | |

2022 | EXAMPLES: |

2023 | sage: E=EllipticCurve([0, 0, 1/216, -7/1296, 1/7776]) |

2024 | sage: E.local_integral_model(2) |

2025 | Elliptic Curve defined by y^2 + 1/27*y = x^3 - 7/81*x + 2/243 over Rational Field |

2026 | sage: E.local_integral_model(3) |

2027 | Elliptic Curve defined by y^2 + 1/8*y = x^3 - 7/16*x + 3/32 over Rational Field |

2028 | sage: E.local_integral_model(2).local_integral_model(3) == EllipticCurve('5077a1') |

2029 | True |

2030 | """ |

2031 | assert p.is_prime(), "p must be prime in local_integral_model()" |

2032 | ai = self.a_invariants() |

2033 | e = min([(ai[i].valuation(p)/[1,2,3,4,6][i]) for i in range(5)]).floor() |

2034 | return constructor.EllipticCurve([ai[i]/p**(e*[1,2,3,4,6][i]) for i in range(5)]) |

2035 | |

2036 | def is_global_integral_model(self): |

2037 | r""" |

2038 | Return true iff self is integral at all primes |

2039 | |

2040 | EXAMPLES: |

2041 | sage: E=EllipticCurve([1/2,1/5,1/5,1/5,1/5]) |

2042 | sage: E.is_global_integral_model() |

2043 | False |

2044 | sage: Emin=E.global_integral_model() |

2045 | sage: Emin.is_global_integral_model() |

2046 | True |

2047 | """ |

2048 | return self.is_integral() |

2049 | |

2050 | def global_integral_model(self): |

2051 | r""" |

2052 | Return a model of self which is integral at all primes |

2053 | |

2054 | EXAMPLES: |

2055 | sage: E = EllipticCurve([0, 0, 1/216, -7/1296, 1/7776]) |

2056 | sage: F = E.global_integral_model(); F |

2057 | Elliptic Curve defined by y^2 + y = x^3 - 7*x + 6 over Rational Field |

2058 | sage: F == EllipticCurve('5077a1') |

2059 | True |

2060 | """ |

2061 | ai = self.a_invariants() |

2062 | for a in ai: |

2063 | if not a.is_integral(): |

2064 | for p, _ in a.denom().factor(): |

2065 | e = min([(ai[i].valuation(p)/[1,2,3,4,6][i]) for i in range(5)]).floor() |

2066 | ai = [ai[i]/p**(e*[1,2,3,4,6][i]) for i in range(5)] |

2067 | for z in ai: |

2068 | assert z.denominator() == 1, "bug in global_integral_model: %s" % ai |

2069 | return constructor.EllipticCurve(ai) |

2070 | |

2071 | def integral_model(self): |

2072 | r""" |

2073 | Return a weierstrass model, $F$, of self with integral coefficients, |

2074 | along with a morphism $\phi$ of points on self to points on $F$. |

2075 | |

2076 | EXAMPLES: |

2077 | sage: E = EllipticCurve([1/2,0,0,5,1/3]) |

2078 | sage: F, phi = E.integral_model() |

2079 | sage: F |

2080 | Elliptic Curve defined by y^2 + 3*x*y = x^3 + 6480*x + 15552 over Rational Field |

2081 | sage: phi |

2082 | Generic morphism: |

2083 | From: Abelian group of points on Elliptic Curve defined by y^2 + 1/2*x*y = x^3 + 5*x + 1/3 over Rational Field |

2084 | To: Abelian group of points on Elliptic Curve defined by y^2 + 3*x*y = x^3 + 6480*x + 15552 over Rational Field |

2085 | Via: (u,r,s,t) = (1/6, 0, 0, 0) |

2086 | sage: P = E([4/9,41/27]) |

2087 | sage: phi(P) |

2088 | (16 : 328 : 1) |

2089 | sage: phi(P) in F |

2090 | True |

2091 | """ |

2092 | F = self.global_integral_model() |

2093 | return F, self.isomorphism_to(F) |

2094 | |

2095 | def integral_weierstrass_model(self): |

2096 | r""" |

2097 | Return a model of the form $y^2 = x^3 + a*x + b$ for this curve with $a,b\in\Z$. |

2098 | |

2099 | EXAMPLES: |

2100 | sage: E = EllipticCurve('17a1') |

2101 | sage: E.integral_weierstrass_model() |

2102 | Elliptic Curve defined by y^2 = x^3 - 11*x - 890 over Rational Field |

2103 | """ |

2104 | F = self.minimal_model().short_weierstrass_model() |

2105 | _,_,_,A,B = F.ainvs() |

2106 | for p in [2,3]: |

2107 | e=min(A.valuation(p)/4,B.valuation(p)/6).floor() |

2108 | A /= Integer(p**(4*e)) |

2109 | B /= Integer(p**(6*e)) |

2110 | return constructor.EllipticCurve([A,B]) |

2111 | |

2112 | def modular_degree(self, algorithm='sympow'): |

2113 | r""" |

2114 | Return the modular degree of this elliptic curve. |

2115 | |

2116 | The result is cached. Subsequence calls, even with a |

2117 | different algorithm, just returned the cached result. |

2118 | |

2119 | INPUT: |

2120 | algorithm -- string: |

2121 | 'sympow' -- (default) use Mark Watkin's (newer) C program sympow |

2122 | 'magma' -- requires that MAGMA be installed (also implemented |

2123 | by Mark Watkins) |

2124 | |

2125 | \note{On 64-bit computers ec does not work, so \sage uses |

2126 | sympow even if ec is selected on a 64-bit computer.} |

2127 | |

2128 | The correctness of this function when called with algorithm "ec" |

2129 | is subject to the following three hypothesis: |

2130 | |

2131 | \begin{itemize} |

2132 | |

2133 | \item Manin's conjecture: the Manin constant is 1 |

2134 | |

2135 | \item Steven's conjecture: the $X_1(N)$-optimal quotient |

2136 | is the curve with minimal Faltings height. (This is proved |

2137 | in most cases.) |

2138 | |

2139 | \item The modular degree fits in a machine double, so it |

2140 | better be less than about 50-some bits. (If you use sympow |

2141 | this contraint does not apply.) |

2142 | |

2143 | \end{itemize} |

2144 | |

2145 | Moreover for all algorithms, computing a certain value of an |

2146 | $L$-function ``uses a heuristic method that discerns when the |

2147 | real-number approximation to the modular degree is within epsilon |

2148 | [=0.01 for algorithm="sympow"] of the same integer for 3 |

2149 | consecutive trials (which occur maybe every 25000 coefficients |

2150 | or so). Probably it could just round at some point. For |

2151 | rigour, you would need to bound the tail by assuming |

2152 | (essentially) that all the $a_n$ are as large as possible, but |

2153 | in practise they exhibit significant (square root) |

2154 | cancellation. One difficulty is that it doesn't do the sum in |

2155 | 1-2-3-4 order; it uses 1-2-4-8---3-6-12-24--9-18-- (Euler |

2156 | product style) instead, and so you have to guess ahead of time |

2157 | at what point to curtail this expansion.'' (Quote from an |

2158 | email of Mark Watkins.) |

2159 | |

2160 | \note{If the curve is loaded from the large Cremona database, |

2161 | then the modular degree is taken from the database.} |

2162 | |

2163 | EXAMPLES: |

2164 | sage: E = EllipticCurve([0, -1, 1, -10, -20]) |

2165 | sage: E |

2166 | Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field |

2167 | sage: E.modular_degree() |

2168 | 1 |

2169 | sage: E = EllipticCurve('5077a') |

2170 | sage: E.modular_degree() |

2171 | 1984 |

2172 | sage: factor(1984) |

2173 | 2^6 * 31 |

2174 | |

2175 | sage: EllipticCurve([0, 0, 1, -7, 6]).modular_degree() |

2176 | 1984 |

2177 | sage: EllipticCurve([0, 0, 1, -7, 6]).modular_degree(algorithm='sympow') |

2178 | 1984 |

2179 | sage: EllipticCurve([0, 0, 1, -7, 6]).modular_degree(algorithm='magma') # optional |

2180 | 1984 |

2181 | |

2182 | We compute the modular degree of the curve with rank 4 having smallest |

2183 | (known) conductor: |

2184 | |

2185 | sage: E = EllipticCurve([1, -1, 0, -79, 289]) |

2186 | sage: factor(E.conductor()) # conductor is 234446 |

2187 | 2 * 117223 |

2188 | sage: factor(E.modular_degree()) |

2189 | 2^7 * 2617 |

2190 | """ |

2191 | try: |

2192 | return self.__modular_degree |

2193 | |

2194 | except AttributeError: |

2195 | if algorithm == 'sympow': |

2196 | from sage.lfunctions.all import sympow |

2197 | m = sympow.modular_degree(self) |

2198 | elif algorithm == 'magma': |

2199 | m = rings.Integer(self._magma_().ModularDegree()) |

2200 | else: |

2201 | raise ValueError, "unknown algorithm %s"%algorithm |

2202 | self.__modular_degree = m |

2203 | return m |

2204 | |

2205 | def modular_parametrization(self): |

2206 | r""" |

2207 | Computes and returns the modular parametrization of this elliptic curve. |

2208 | |

2209 | The curve is converted to a minimal model. |

2210 | |

2211 | OUTPUT: A list of two Larent series [X(x),Y(x)] of degrees -2, |

2212 | -3 respectively, which satisfy the equation of the (minimal |

2213 | model of the) elliptic curve. The are modular functions on |

2214 | $\Gamma_0(N)$ where $N$ is the conductor. |

2215 | |

2216 | X.deriv()/(2*Y+a1*X+a3) should equal f(q)dq/q where f is |

2217 | self.q_expansion(). |

2218 | |

2219 | EXAMPLES: |

2220 | sage: E=EllipticCurve('389a1') |

2221 | sage: X,Y=E.modular_parametrization() |

2222 | sage: X |

2223 | q^-2 + 2*q^-1 + 4 + 7*q + 13*q^2 + 18*q^3 + 31*q^4 + 49*q^5 + 74*q^6 + 111*q^7 + 173*q^8 + 251*q^9 + 379*q^10 + 560*q^11 + 824*q^12 + 1199*q^13 + 1773*q^14 + 2365*q^15 + 3463*q^16 + 4508*q^17 + O(q^18) |

2224 | sage: Y |

2225 | -q^-3 - 3*q^-2 - 8*q^-1 - 17 - 33*q - 61*q^2 - 110*q^3 - 186*q^4 - 320*q^5 - 528*q^6 - 861*q^7 - 1383*q^8 - 2218*q^9 - 3472*q^10 - 5451*q^11 - 8447*q^12 - 13020*q^13 - 20083*q^14 - 29512*q^15 - 39682*q^16 + O(q^17) |

2226 | |

2227 | The following should give 0, but only approximately: |

2228 | sage: q = X.parent().gen() |

2229 | sage: E.defining_polynomial()(X,Y,1) + O(q^11) == 0 |

2230 | True |

2231 | |

2232 | Note that below we have to change variable from x to q |

2233 | sage: a1,_,a3,_,_=E.a_invariants() |

2234 | sage: f=E.q_expansion(17) |

2235 | sage: q=f.parent().gen() |

2236 | sage: f/q == (X.derivative()/(2*Y+a1*X+a3)) |

2237 | True |

2238 | |

2239 | """ |

2240 | R = LaurentSeriesRing(RationalField(),'q') |

2241 | XY = self.pari_mincurve().elltaniyama() |

2242 | return [1/R(1/XY[0]),1/R(1/XY[1])] |

2243 | |

2244 | def cremona_label(self, space=False): |

2245 | """ |

2246 | Return the Cremona label associated to (the minimal model) of |

2247 | this curve, if it is known. If not, raise a RuntimeError |

2248 | exception. |

2249 | |

2250 | EXAMPLES: |

2251 | sage: E=EllipticCurve('389a1') |

2252 | sage: E.cremona_label() |

2253 | '389a1' |

2254 | |

2255 | The default database only contains conductors up to 10000, so |

2256 | any curve with conductor greater than that will cause an error |

2257 | to be raised: |

2258 | |

2259 | sage: E=EllipticCurve([1,2,3,4,5]); |

2260 | sage: E.conductor() |

2261 | 10351 |

2262 | sage: E.cremona_label() |

2263 | Traceback (most recent call last): |

2264 | ... |

2265 | RuntimeError: Cremona label not known for Elliptic Curve defined by y^2 + x*y + 3*y = x^3 + 2*x^2 + 4*x + 5 over Rational Field. |

2266 | """ |

2267 | try: |

2268 | if not space: |

2269 | return self.__cremona_label.replace(' ','') |

2270 | return self.__cremona_label |

2271 | except AttributeError: |

2272 | try: |

2273 | X = self.database_curve() |

2274 | except RuntimeError: |

2275 | raise RuntimeError, "Cremona label not known for %s."%self |

2276 | self.__cremona_label = X.__cremona_label |

2277 | return self.cremona_label(space) |

2278 | |

2279 | label = cremona_label |

2280 | |

2281 | def torsion_order(self): |

2282 | """ |

2283 | Return the order of the torsion subgroup. |

2284 | |

2285 | EXAMPLES: |

2286 | sage: e = EllipticCurve('11a') |

2287 | sage: e.torsion_order() |

2288 | 5 |

2289 | sage: type(e.torsion_order()) |

2290 | <type 'sage.rings.integer.Integer'> |

2291 | sage: e = EllipticCurve([1,2,3,4,5]) |

2292 | sage: e.torsion_order() |

2293 | 1 |

2294 | sage: type(e.torsion_order()) |

2295 | <type 'sage.rings.integer.Integer'> |

2296 | """ |

2297 | try: |

2298 | return self.__torsion_order |

2299 | except AttributeError: |

2300 | self.__torsion_order = self.torsion_subgroup().order() |

2301 | return self.__torsion_order |

2302 | |

2303 | def torsion_subgroup(self, flag=0): |

2304 | """ |

2305 | Returns the torsion subgroup of this elliptic curve. |

2306 | |

2307 | INPUT: |

2308 | flag -- (default: 0) chooses PARI algorithm: |

2309 | flag = 0: uses Doud algorithm |

2310 | flag = 1: uses Lutz-Nagell algorithm |

2311 | |

2312 | OUTPUT: |

2313 | The EllipticCurveTorsionSubgroup instance associated to this elliptic curve. |

2314 | |

2315 | EXAMPLES: |

2316 | sage: EllipticCurve('11a').torsion_subgroup() |

2317 | Torsion Subgroup isomorphic to Multiplicative Abelian Group isomorphic to C5 associated to the Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field |

2318 | sage: EllipticCurve('37b').torsion_subgroup() |

2319 | Torsion Subgroup isomorphic to Multiplicative Abelian Group isomorphic to C3 associated to the Elliptic Curve defined by y^2 + y = x^3 + x^2 - 23*x - 50 over Rational Field |

2320 | |

2321 | sage: e = EllipticCurve([-1386747,368636886]);e |

2322 | Elliptic Curve defined by y^2 = x^3 - 1386747*x + 368636886 over Rational Field |

2323 | sage: G = e.torsion_subgroup(); G |

2324 | Torsion Subgroup isomorphic to Multiplicative Abelian |

2325 | Group isomorphic to C8 x C2 associated to the Elliptic |

2326 | Curve defined by y^2 = x^3 - 1386747*x + 368636886 over |

2327 | Rational Field |

2328 | sage: G.0 |

2329 | (1227 : 22680 : 1) |

2330 | sage: G.1 |

2331 | (282 : 0 : 1) |

2332 | sage: list(G) |

2333 | [1, P1, P0, P0*P1, P0^2, P0^2*P1, P0^3, P0^3*P1, P0^4, P0^4*P1, P0^5, P0^5*P1, P0^6, P0^6*P1, P0^7, P0^7*P1] |

2334 | """ |

2335 | try: |

2336 | return self.__torsion_subgroup |

2337 | except AttributeError: |

2338 | self.__torsion_subgroup = rational_torsion.EllipticCurveTorsionSubgroup(self, flag) |

2339 | self.__torsion_order = self.__torsion_subgroup.order() |

2340 | return self.__torsion_subgroup |

2341 | |

2342 | ## def newform_eval(self, z, prec): |

2343 | ## """ |

2344 | ## The value of the newform attached to this elliptic curve at |

2345 | ## the point z in the complex upper half plane, computed using |

2346 | ## prec terms of the power series expansion. Note that the power |

2347 | ## series need not converge well near the real axis. |

2348 | ## """ |

2349 | ## raise NotImplementedError |

2350 | |

2351 | def root_number(self): |

2352 | """ |

2353 | Returns the root number of this elliptic curve. |

2354 | |

2355 | This is 1 if the order of vanishing of the L-function L(E,s) |

2356 | at 1 is even, and -1 if it is odd. |

2357 | |

2358 | EXAMPLES: |

2359 | sage: EllipticCurve('11a1').root_number() |

2360 | 1 |

2361 | sage: EllipticCurve('37a1').root_number() |

2362 | -1 |

2363 | sage: EllipticCurve('389a1').root_number() |

2364 | 1 |

2365 | """ |

2366 | try: |

2367 | return self.__root_number |

2368 | except AttributeError: |

2369 | self.__root_number = int(self.pari_mincurve().ellrootno()) |

2370 | return self.__root_number |

2371 | |

2372 | |

2373 | def has_cm(self): |

2374 | """ |

2375 | Returns True iff this elliptic curve has Complex Multiplication |

2376 | |

2377 | EXAMPLES: |

2378 | sage: E=EllipticCurve('37a1') |

2379 | sage: E.has_cm() |

2380 | False |

2381 | sage: E=EllipticCurve('32a1') |

2382 | sage: E.has_cm() |

2383 | True |

2384 | sage: E.j_invariant() |

2385 | 1728 |

2386 | """ |

2387 | |

2388 | return CMJ.has_key(self.j_invariant()) |

2389 | |

2390 | def cm_discriminant(self): |

2391 | """ |

2392 | Returns the associated quadratic discriminant if this elliptic |

2393 | curve has Complex Multiplication. |

2394 | |

2395 | A ValueError is raised if the curve does not have CM (see the |

2396 | function has_cm()) |

2397 | |

2398 | EXAMPLES: |

2399 | sage: E=EllipticCurve('32a1') |

2400 | sage: E.cm_discriminant() |

2401 | -4 |

2402 | sage: E=EllipticCurve('121b1') |

2403 | sage: E.cm_discriminant() |

2404 | -11 |

2405 | sage: E=EllipticCurve('37a1') |

2406 | sage: E.cm_discriminant() |

2407 | Traceback (most recent call last): |

2408 | ... |

2409 | ValueError: Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field does not have CM |

2410 | """ |

2411 | |

2412 | try: |

2413 | return CMJ[self.j_invariant()] |

2414 | except KeyError: |

2415 | raise ValueError, "%s does not have CM"%self |

2416 | |

2417 | |

2418 | def quadratic_twist(self, D): |

2419 | """ |

2420 | Return the global minimal model of the quadratic twist of this |

2421 | curve by D. |

2422 | |

2423 | EXAMPLES: |

2424 | sage: E=EllipticCurve('37a1') |

2425 | sage: E7=E.quadratic_twist(7); E7 |

2426 | Elliptic Curve defined by y^2 = x^3 - 784*x + 5488 over Rational Field |

2427 | sage: E7.conductor() |

2428 | 29008 |

2429 | sage: E7.quadratic_twist(7) == E |

2430 | True |

2431 | """ |

2432 | return EllipticCurve_number_field.quadratic_twist(self, D).minimal_model() |

2433 | |

2434 | |

2435 | ########################################################## |

2436 | # Isogeny class (currently just uses Cremona database.) |

2437 | ########################################################## |

2438 | def isogeny_class(self, algorithm="mwrank", verbose=False): |

2439 | r""" |

2440 | Return all curves over $\Q$ in the isogeny class of this |

2441 | elliptic curve. |

2442 | |

2443 | INPUT: |

2444 | algorithm -- string: |

2445 | "mwrank" -- (default) use the mwrank C++ library |

2446 | "database" -- use the Cremona database (only works if |

2447 | curve is isomorphic to a curve in the database) |

2448 | |

2449 | OUTPUT: |

2450 | Returns the sorted list of the curves isogenous to self. |

2451 | If algorithm is "mwrank", also returns the isogeny matrix (otherwise |

2452 | returns None as second return value). |

2453 | |

2454 | \note{The result is \emph{not} provably correct, in the sense |

2455 | that when the numbers are huge isogenies could be missed |

2456 | because of precision issues.} |

2457 | |

2458 | \note{The ordering depends on which algorithm is used.} |

2459 | |

2460 | EXAMPLES: |

2461 | sage: I, A = EllipticCurve('37b').isogeny_class('mwrank') |

2462 | sage: I # randomly ordered |

2463 | [Elliptic Curve defined by y^2 + y = x^3 + x^2 - 23*x - 50 over Rational Field, |

2464 | Elliptic Curve defined by y^2 + y = x^3 + x^2 - 1873*x - 31833 over Rational Field, |

2465 | Elliptic Curve defined by y^2 + y = x^3 + x^2 - 3*x +1 over Rational Field] |

2466 | sage: A |

2467 | [0 3 3] |

2468 | [3 0 0] |

2469 | [3 0 0] |

2470 | |

2471 | sage: I, _ = EllipticCurve('37b').isogeny_class('database'); I |

2472 | [Elliptic Curve defined by y^2 + y = x^3 + x^2 - 1873*x - 31833 over Rational Field, |

2473 | Elliptic Curve defined by y^2 + y = x^3 + x^2 - 23*x - 50 over Rational Field, |

2474 | Elliptic Curve defined by y^2 + y = x^3 + x^2 - 3*x +1 over Rational Field] |

2475 | |

2476 | This is an example of a curve with a $37$-isogeny: |

2477 | sage: E = EllipticCurve([1,1,1,-8,6]) |

2478 | sage: E.isogeny_class () |

2479 | ([Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 - 8*x + 6 over Rational Field, |

2480 | Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 - 208083*x - 36621194 over Rational Field], |

2481 | [ 0 37] |

2482 | [37 0]) |

2483 | |

2484 | This curve had numerous $2$-isogenies: |

2485 | sage: e=EllipticCurve([1,0,0,-39,90]) |

2486 | sage: e.isogeny_class () |

2487 | ([Elliptic Curve defined by y^2 + x*y = x^3 - 39*x + 90 over Rational Field, |

2488 | Elliptic Curve defined by y^2 + x*y = x^3 - 4*x -1 over Rational Field, |

2489 | Elliptic Curve defined by y^2 + x*y = x^3 + x over Rational Field, |

2490 | Elliptic Curve defined by y^2 + x*y = x^3 - 49*x - 136 over Rational Field, |

2491 | Elliptic Curve defined by y^2 + x*y = x^3 - 34*x - 217 over Rational Field, |

2492 | Elliptic Curve defined by y^2 + x*y = x^3 - 784*x - 8515 over Rational Field], |

2493 | [0 2 0 0 0 0] |

2494 | [2 0 2 2 0 0] |

2495 | [0 2 0 0 0 0] |

2496 | [0 2 0 0 2 2] |

2497 | [0 0 0 2 0 0] |

2498 | [0 0 0 2 0 0]) |

2499 | |

2500 | See \url{http://modular.ucsd.edu/Tables/nature/} for more interesting |

2501 | examples of isogeny structures. |

2502 | """ |

2503 | #if algorithm == "gp": |

2504 | |

2505 | # return sum([L for _, L in self.isogenous_curves(algorithm="gp")], [self]) |

2506 | |

2507 | if algorithm == "mwrank": |

2508 | try: |

2509 | E = self.mwrank_curve() |

2510 | except ValueError: |

2511 | E = self.minimal_model().mwrank_curve() |

2512 | I, A = E.isogeny_class(verbose=verbose) |

2513 | mat = matrix.MatrixSpace(rings.IntegerRing(), len(A))(A) |

2514 | I = [constructor.EllipticCurve(ainvs) for ainvs in I] |

2515 | return I, mat |

2516 | |

2517 | elif algorithm == "database": |

2518 | |

2519 | try: |

2520 | label = self.cremona_label(space=False) |

2521 | except RuntimeError: |

2522 | raise RuntimeError, "unable to to find %s in the database"%self |

2523 | db = sage.databases.cremona.CremonaDatabase() |

2524 | I = db.isogeny_class(label) |

2525 | I.sort() |

2526 | return I, None |

2527 | |

2528 | else: |

2529 | |

2530 | raise ValueError, "unknown algorithm '%s'%"%algorithm |

2531 | |

2532 | def isogeny_graph(self): |

2533 | r""" |

2534 | Returns a graph representing the isogeny class of this elliptic curve, |

2535 | where the vertices are isogenous curves over $\Q$ and the edges are |

2536 | prime degree isogenies labeled by their degree. |

2537 | |

2538 | EXAMPLES: |

2539 | sage: LL = [] |

2540 | sage: for e in cremona_optimal_curves(range(1, 38)): |

2541 | ... G = e.isogeny_graph() |

2542 | ... already = False |

2543 | ... for H in LL: |

2544 | ... if G.is_isomorphic(H): |

2545 | ... already = True |

2546 | ... break |

2547 | ... if not already: |

2548 | ... LL.append(G) |

2549 | ... |

2550 | sage: graphs_list.show_graphs(LL) |

2551 | |

2552 | sage: E = EllipticCurve('195a') |

2553 | sage: G = E.isogeny_graph() |

2554 | sage: for v in G: print v, G.get_vertex(v) |

2555 | ... |

2556 | 0 Elliptic Curve defined by y^2 + x*y = x^3 - 110*x + 435 over Rational Field |

2557 | 1 Elliptic Curve defined by y^2 + x*y = x^3 - 115*x + 392 over Rational Field |

2558 | 2 Elliptic Curve defined by y^2 + x*y = x^3 + 210*x + 2277 over Rational Field |

2559 | 3 Elliptic Curve defined by y^2 + x*y = x^3 - 520*x - 4225 over Rational Field |

2560 | 4 Elliptic Curve defined by y^2 + x*y = x^3 + 605*x - 19750 over Rational Field |

2561 | 5 Elliptic Curve defined by y^2 + x*y = x^3 - 8125*x - 282568 over Rational Field |

2562 | 6 Elliptic Curve defined by y^2 + x*y = x^3 - 7930*x - 296725 over Rational Field |

2563 | 7 Elliptic Curve defined by y^2 + x*y = x^3 - 130000*x - 18051943 over Rational Field |

2564 | sage: G.plot(edge_labels=True).save('isogeny_graph.png') |

2565 | """ |

2566 | from sage.graphs.graph import Graph |

2567 | L, M = self.isogeny_class() |

2568 | G = Graph(M, format='weighted_adjacency_matrix') |

2569 | d = {} |

2570 | for v in G.vertices(): |

2571 | d[v] = L[v] |

2572 | G.set_vertices(d) |

2573 | return G |

2574 | |

2575 | ########################################################## |

2576 | # Galois Representations |

2577 | ########################################################## |

2578 | |

2579 | def is_reducible(self, p): |

2580 | """ |

2581 | Return True if the mod-p representation attached |

2582 | to E is reducible. |

2583 | |

2584 | INPUT: |

2585 | p -- a prime number |

2586 | |

2587 | NOTE: The answer is cached. |

2588 | |

2589 | EXAMPLES: |

2590 | sage: E = EllipticCurve('121a'); E |

2591 | Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 - 30*x - 76 over Rational Field |

2592 | sage: E.is_reducible(7) |

2593 | False |

2594 | sage: E.is_reducible(11) |

2595 | True |

2596 | sage: EllipticCurve('11a').is_reducible(5) |

2597 | True |

2598 | sage: e = EllipticCurve('11a2') |

2599 | sage: e.is_reducible(5) |

2600 | True |

2601 | sage: e.torsion_order() |

2602 | 1 |

2603 | """ |

2604 | try: |

2605 | return self.__is_reducible[p] |

2606 | except AttributeError: |

2607 | self.__is_reducible = {} |

2608 | except KeyError: |

2609 | pass |

2610 | |

2611 | if not arith.is_prime(p): |

2612 | raise ValueError, 'p (=%s) must be prime'%p |

2613 | # we do is_surjective first, since this is |

2614 | # much easier than computing isogeny_class |

2615 | t, why = self.is_surjective(p) |

2616 | if t == True: |

2617 | self.__is_reducible[p] = False |

2618 | return False # definitely not reducible |

2619 | isogeny_matrix = self.isogeny_class()[ 1 ] |

2620 | v = isogeny_matrix.row(0) # first row |

2621 | for a in v: |

2622 | if a != 0 and a % p == 0: |

2623 | self.__is_reducible[p] = True |

2624 | return True |

2625 | self.__is_reducible[p] = False |

2626 | return False |

2627 | |

2628 | def is_irreducible(self, p): |

2629 | """ |

2630 | Return True if the mod p represenation is irreducible. |

2631 | |

2632 | EXAMPLES: |

2633 | sage: e = EllipticCurve('37b') |

2634 | sage: e.is_irreducible(2) |

2635 | True |

2636 | sage: e.is_irreducible(3) |

2637 | False |

2638 | sage: e.is_reducible(2) |

2639 | False |

2640 | sage: e.is_reducible(3) |

2641 | True |

2642 | """ |

2643 | return not self.is_reducible(p) |

2644 | |

2645 | def is_surjective(self, p, A=1000): |

2646 | """ |

2647 | Return True if the mod-p representation attached to E |

2648 | is surjective, False if it is not, or None if we were |

2649 | unable to determine whether it is or not. |

2650 | |

2651 | NOTE: The answer is cached. |

2652 | |

2653 | INPUT: |

2654 | p -- int (a prime number) |

2655 | A -- int (a bound on the number of a_p to use) |

2656 | |

2657 | OUTPUT: |

2658 | a 2-tuple: |

2659 | -- surjective or (probably) not |

2660 | -- information about what it is if not surjective |

2661 | |

2662 | EXAMPLES: |

2663 | sage: e = EllipticCurve('37b') |

2664 | sage: e.is_surjective(2) |

2665 | (True, None) |

2666 | sage: e.is_surjective(3) |

2667 | (False, '3-torsion') |

2668 | |

2669 | |

2670 | REMARKS: |

2671 | |

2672 | 1. If p >= 5 then the mod-p representation is surjective |

2673 | if and only if the p-adic representation is |

2674 | surjective. When p = 2, 3 there are counterexamples. |

2675 | See a very recent paper of Elkies for more details |

2676 | when p=3. |

2677 | |

2678 | 2. When p <= 3 this function always gives the correct |

2679 | result irregardless of A, since it explicitly |

2680 | determines the p-division polynomial. |

2681 | |

2682 | """ |

2683 | if not arith.is_prime(p): |

2684 | raise TypeError, "p (=%s) must be prime."%p |

2685 | A = int(A) |

2686 | key = (p, A) |

2687 | try: |

2688 | return self.__is_surjective[key] |

2689 | except KeyError: |

2690 | pass |

2691 | except AttributeError: |

2692 | self.__is_surjective = {} |

2693 | |

2694 | ans = self._is_surjective(p, A) |

2695 | self.__is_surjective[key] = ans |

2696 | return ans |

2697 | |

2698 | def _is_surjective(self, p, A): |

2699 | T = self.torsion_subgroup().order() |

2700 | if T % p == 0: |

2701 | return False, "%s-torsion"%p |

2702 | |

2703 | if p == 2: |

2704 | # E is isomorphic to [0,b2,0,8*b4,16*b6] |

2705 | b2,b4,b6,b8=self.b_invariants() |

2706 | R = rings.PolynomialRing(self.base_ring(), 'x') |

2707 | x = R.gen() |

2708 | f = x**3 + b2*x**2 + 8*b4*x + 16*b6 |

2709 | if not f.is_irreducible(): |

2710 | return False, '2-torsion' |

2711 | if arith.is_square(f.discriminant()): |

2712 | return False, "A3" |

2713 | return True, None |

2714 | |

2715 | if p == 3: |

2716 | # Algorithm: Let f be the 3-division polynomial, which is |

2717 | # a polynomial of degree 4. Then I claim that this |

2718 | # polynomial has Galois group S_4 if and only if the |

2719 | # representation rhobar_{E,3} is surjective. If the group |

2720 | # is S_4, then S_4 is a quotient of the image of |

2721 | # rhobar_{E,3}. Since S_4 has order 24 and GL_2(F_3) |

2722 | # has order 48, the only possibility we have to consider |

2723 | # is that the image of rhobar is isomorphic to S_4. |

2724 | # But this is not the case because S_4 is not a subgroup |

2725 | # of GL_2(F_3). If it were, it would be normal, since |

2726 | # it would have index 2. But there is a *unique* normal |

2727 | # subgroup of GL_2(F_3) of index 2, namely SL_2(F_3), |

2728 | # and SL_2(F_3) is not isomorphic to S_4 (S_4 has a normal |

2729 | # subgroup of index 2 and SL_2(F_3) does not.) |

2730 | # (What's a simple way to see that SL_2(F_3) is the |

2731 | # unique index-2 normal subgroup? I didn't see an obvious |

2732 | # reason, so just used the NormalSubgroups command in MAGMA |

2733 | # and it output exactly one of index 2.) |

2734 | |

2735 | # Here's Noam Elkies proof for the other direction: |

2736 | |

2737 | #> Let E be an elliptic curve over Q. Is the mod-3 |

2738 | #> representation E[3] surjective if and only if the |

2739 | #> (degree 4) division polynomial has Galois group S_4? I |

2740 | #> can see why the group being S_4 implies the |

2741 | #> representation is surjective, but the converse is not |

2742 | #> clear to me. |

2743 | # I would have thought that this is the easier part: to |

2744 | # say that E[3] is surjective is to say the 3-torsion |

2745 | # field Q(E[3]) has Galois group GL_2(Z/3) over Q. Let |

2746 | # E[3]+ be the subfield fixed by the element -1 of |

2747 | # GL_2(Z/3). Then E[3] has Galois group PGL_2(Z/3), which |

2748 | # is identified with S_4 by its action on the four |

2749 | # 3-element subgroups of E[3]. Each such subgroup is in |

2750 | # turn determined by the x-coordinate shared by its two |

2751 | # nonzero points. So, if E[3] is surjective then any |

2752 | # permutation of those x-coordinates is realized by some |

2753 | # element of Gal(E[3]+/Q). Thus the Galois group of the |

2754 | # division polynomial (whose roots are those |

2755 | # x-coordinates) maps surjectively to S_4, which means it |

2756 | # equals S_4. |

2757 | |

2758 | |

2759 | f = self.division_polynomial(3) |

2760 | if not f.is_irreducible(): |

2761 | return False, "reducible_3-divpoly" |

2762 | n = pari(f).polgalois()[0] |

2763 | if n == 24: |

2764 | return True, None |

2765 | else: |

2766 | return False, "3-divpoly_galgroup_order_%s"%n |

2767 | |

2768 | if self.has_cm(): |

2769 | return False, "CM" |

2770 | an = self.anlist(A) |

2771 | ell = 0 |

2772 | Np = self.conductor() * p |

2773 | signs = [] |

2774 | while True: |

2775 | ell = arith.next_prime(ell) |

2776 | if ell >= A: break |

2777 | if Np % ell != 0: |

2778 | a_ell = an[int(ell)] |

2779 | if a_ell % p != 0: |

2780 | s = arith.kronecker(a_ell**2 - 4*ell, p) |

2781 | #print ell, s |

2782 | if s == 0: continue |

2783 | if not (s in signs): |

2784 | signs.append(s) |

2785 | if len(signs) == 2: |

2786 | return True, None |

2787 | |

2788 | # could do something further here... |

2789 | return False, signs |

2790 | |

2791 | def is_semistable(self): |

2792 | """ |

2793 | Return True iff this elliptic curve is semi-stable at all primes. |

2794 | |

2795 | EXAMPLES: |

2796 | sage: E=EllipticCurve('37a1') |

2797 | sage: E.is_semistable() |

2798 | True |

2799 | sage: E=EllipticCurve('90a1') |

2800 | sage: E.is_semistable() |

2801 | False |

2802 | """ |

2803 | if self.base_ring() != Q: |

2804 | raise NotImplementedError, "is_semistable only implemented for curves over the rational numbers." |

2805 | return self.conductor().is_squarefree() |

2806 | |

2807 | def reducible_primes(self): |

2808 | r""" |

2809 | Returns a list of the primes $p$ such that the mod $p$ |

2810 | representation $\rho_{E,p}$ is reducible. For all other |

2811 | primes the representation is irreducible. |

2812 | |

2813 | NOTE -- this is \emph{not} provably correct in general. |

2814 | See the documentation for \code{self.isogeny_class}. |

2815 | |

2816 | EXAMPLES: |

2817 | sage: E = EllipticCurve('225a') |

2818 | sage: E.reducible_primes() |

2819 | [3] |

2820 | """ |

2821 | try: |

2822 | return self.__reducible_primes |

2823 | except AttributeError: |

2824 | pass |

2825 | C, I = self.isogeny_class(algorithm='mwrank') |

2826 | X = set(I.list()) |

2827 | R = [p for p in X if arith.is_prime(p)] |

2828 | self.__reducible_primes = R |

2829 | return R |

2830 | |

2831 | def non_surjective(self, A=1000): |

2832 | r""" |

2833 | Returns a list of primes p such that the mod-p representation |

2834 | $\rho_{E,p}$ *might* not be surjective (this list usually |

2835 | contains 2, because of shortcomings of the algorithm). If p |

2836 | is not in the returned list, then rho_{E,p} is provably |

2837 | surjective (see A. Cojocaru's paper). If the curve has CM |

2838 | then infinitely many representations are not surjective, so we |

2839 | simply return the sequence [(0,"cm")] and do no further computation. |

2840 | |

2841 | INPUT: |

2842 | A -- an integer |

2843 | OUTPUT: |

2844 | list -- if curve has CM, returns [(0,"cm")]. Otherwise, returns a |

2845 | list of primes where mod-p representation very likely |

2846 | not surjective. At any prime not in this list, |

2847 | the representation is definitely surjective. |

2848 | EXAMPLES: |

2849 | sage: E = EllipticCurve([0, 0, 1, -38, 90]) # 361A |

2850 | sage: E.non_surjective() # CM curve |

2851 | [(0, 'cm')] |

2852 | |

2853 | sage: E = EllipticCurve([0, -1, 1, 0, 0]) # X_1(11) |

2854 | sage: E.non_surjective() |

2855 | [(5, '5-torsion')] |

2856 | |

2857 | sage: E = EllipticCurve([0, 0, 1, -1, 0]) # 37A |

2858 | sage: E.non_surjective() |

2859 | [] |

2860 | |

2861 | sage: E = EllipticCurve([0,-1,1,-2,-1]) # 141C |

2862 | sage: E.non_surjective() |

2863 | [(13, [1])] |

2864 | |

2865 | ALGORITHM: |

2866 | When p<=3 use division polynomials. For 5 <= p <= B, |

2867 | where B is Cojocaru's bound, use the results in Section 2 |

2868 | of Serre's inventiones paper"Sur Les Representations Modulaires Deg |

2869 | Degre 2 de Galqbar Over Q." |

2870 | """ |

2871 | if self.has_cm(): |

2872 | misc.verbose("cm curve") |

2873 | return [(0,"cm")] |

2874 | N = self.conductor() |

2875 | if self.is_semistable(): |

2876 | C = 11 |

2877 | misc.verbose("semistable -- so bound is 11") |

2878 | else: |

2879 | C = 1 + 4*sqrt(6)*int(N)/3 * sqrt(mul([1+1.0/int(p) for p,_ in factor(N)])) |

2880 | misc.verbose("conductor = %s, and bound is %s"%(N,C)) |

2881 | C = 1 + 4*sqrt(6)*int(N)/3 * sqrt(mul([1+1.0/int(p) for p,_ in factor(N)])) |

2882 | misc.verbose("conductor = %s, and bound is %s"%(N,C)) |

2883 | B = [] |

2884 | p = 2 |

2885 | while p <= C: |

2886 | t, v = self.is_surjective(p, A=A) |

2887 | misc.verbose("(%s,%s,%s)"%(p,t,v)) |

2888 | if not t: |

2889 | B.append((p,v)) |

2890 | p = next_prime(p) |

2891 | return B |

2892 | |

2893 | def is_ordinary(self, p, ell=None): |

2894 | """ |

2895 | Return True precisely when the mod-p representation attached |

2896 | to this elliptic curve is ordinary at ell. |

2897 | |

2898 | INPUT: |

2899 | p -- a prime |

2900 | ell - a prime (default: p) |

2901 | |

2902 | OUTPUT: |

2903 | bool |

2904 | |

2905 | EXAMPLES: |

2906 | sage: E=EllipticCurve('37a1') |

2907 | sage: E.is_ordinary(37) |

2908 | True |

2909 | sage: E=EllipticCurve('32a1') |

2910 | sage: E.is_ordinary(2) |

2911 | False |

2912 | sage: [p for p in prime_range(50) if E.is_ordinary(p)] |

2913 | [5, 13, 17, 29, 37, 41] |

2914 | """ |

2915 | if ell is None: |

2916 | ell = p |

2917 | return self.ap(ell) % p != 0 |

2918 | |

2919 | def is_good(self, p, check=True): |

2920 | """ |

2921 | Return True if $p$ is a prime of good reduction for $E$. |

2922 | |

2923 | INPUT: |

2924 | p -- a prime |

2925 | |

2926 | OUTPUT: |

2927 | bool |

2928 | |

2929 | EXAMPLES: |

2930 | sage: e = EllipticCurve('11a') |

2931 | sage: e.is_good(-8) |

2932 | Traceback (most recent call last): |

2933 | ... |

2934 | ValueError: p must be prime |

2935 | sage: e.is_good(-8, check=False) |

2936 | True |

2937 | """ |

2938 | if check: |

2939 | if not arith.is_prime(p): |

2940 | raise ValueError, "p must be prime" |

2941 | return self.conductor() % p != 0 |

2942 | |

2943 | |

2944 | def is_supersingular(self, p, ell=None): |

2945 | """ |

2946 | Return True precisely when p is a prime of good reduction and |

2947 | the mod-p representation attached to this elliptic curve is |

2948 | supersingular at ell. |

2949 | |

2950 | INPUT: |

2951 | p -- a prime |

2952 | ell - a prime (default: p) |

2953 | |

2954 | OUTPUT: |

2955 | bool |

2956 | |

2957 | EXAMPLES: |

2958 | sage: E=EllipticCurve('37a1') |

2959 | sage: E.is_supersingular(37) |

2960 | False |

2961 | sage: E=EllipticCurve('32a1') |

2962 | sage: E.is_supersingular(2) |

2963 | False |

2964 | sage: E.is_supersingular(7) |

2965 | True |

2966 | sage: [p for p in prime_range(50) if E.is_supersingular(p)] |

2967 | [3, 7, 11, 19, 23, 31, 43, 47] |

2968 | """ |

2969 | if ell is None: |

2970 | ell = p |

2971 | return self.is_good(p) and not self.is_ordinary(p, ell) |

2972 | |

2973 | def supersingular_primes(self, B): |

2974 | """ |

2975 | Return a list of all supersingular primes for this elliptic curve |

2976 | up to and possibly including B. |

2977 | |

2978 | EXAMPLES: |

2979 | sage: e = EllipticCurve('11a') |

2980 | sage: e.aplist(20) |

2981 | [-2, -1, 1, -2, 1, 4, -2, 0] |

2982 | sage: e.supersingular_primes(1000) |

2983 | [2, 19, 29, 199, 569, 809] |

2984 | |

2985 | sage: e = EllipticCurve('27a') |

2986 | sage: e.aplist(20) |

2987 | [0, 0, 0, -1, 0, 5, 0, -7] |

2988 | sage: e.supersingular_primes(97) |

2989 | [2, 5, 11, 17, 23, 29, 41, 47, 53, 59, 71, 83, 89] |

2990 | sage: e.ordinary_primes(97) |

2991 | [7, 13, 19, 31, 37, 43, 61, 67, 73, 79, 97] |

2992 | sage: e.supersingular_primes(3) |

2993 | [2] |

2994 | sage: e.supersingular_primes(2) |

2995 | [2] |

2996 | sage: e.supersingular_primes(1) |

2997 | [] |

2998 | """ |

2999 | v = self.aplist(max(B, 3)) |

3000 | P = arith.prime_range(max(B,3)+1) |

3001 | N = self.conductor() |

3002 | return [P[i] for i in [0,1] if P[i] <= B and v[i]%P[i]==0 and N%P[i] != 0] + \ |

3003 | [P[i] for i in range(2,len(v)) if v[i] == 0 and N%P[i] != 0] |

3004 | |

3005 | def ordinary_primes(self, B): |

3006 | """ |

3007 | Return a list of all ordinary primes for this elliptic curve |

3008 | up to and possibly including B. |

3009 | |

3010 | EXAMPLES: |

3011 | sage: e = EllipticCurve('11a') |

3012 | sage: e.aplist(20) |

3013 | [-2, -1, 1, -2, 1, 4, -2, 0] |

3014 | sage: e.ordinary_primes(97) |

3015 | [3, 5, 7, 11, 13, 17, 23, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97] |

3016 | sage: e = EllipticCurve('49a') |

3017 | sage: e.aplist(20) |

3018 | [1, 0, 0, 0, 4, 0, 0, 0] |

3019 | sage: e.supersingular_primes(97) |

3020 | [3, 5, 13, 17, 19, 31, 41, 47, 59, 61, 73, 83, 89, 97] |

3021 | sage: e.ordinary_primes(97) |

3022 | [2, 11, 23, 29, 37, 43, 53, 67, 71, 79] |

3023 | sage: e.ordinary_primes(3) |

3024 | [2] |

3025 | sage: e.ordinary_primes(2) |

3026 | [2] |

3027 | sage: e.ordinary_primes(1) |

3028 | [] |

3029 | """ |

3030 | v = self.aplist(max(B, 3) ) |

3031 | P = arith.prime_range(max(B,3) +1) |

3032 | return [P[i] for i in [0,1] if P[i] <= B and v[i]%P[i]!=0] +\ |

3033 | [P[i] for i in range(2,len(v)) if v[i] != 0] |

3034 | |

3035 | def eval_modular_form(self, points, prec): |

3036 | """ |

3037 | Evaluate the L-series of this elliptic curve at points in CC |

3038 | |

3039 | INPUT: |

3040 | points-- a list of points in the half-plane of convergence |

3041 | prec-- precision |

3042 | |

3043 | OUTPUT: |

3044 | A list of values L(E,s) for s in points |

3045 | |

3046 | NOTE: Better examples are welcome. |

3047 | |

3048 | EXAMPLES: |

3049 | sage: E=EllipticCurve('37a1') |

3050 | sage: E.eval_modular_form([1.5+I,2.0+I,2.5+I],0.000001) |

3051 | [0, 0, 0] |

3052 | """ |

3053 | if not isinstance(points, (list,xrange)): |

3054 | try: |

3055 | points = list(points) |

3056 | except TypeError: |

3057 | return self.eval_modular_form([points],prec) |

3058 | an = self.pari_mincurve().ellan(prec) |

3059 | s = 0 |

3060 | I = pari("I") |

3061 | pi = pari("Pi") |

3062 | c = pari(2)*pi*I |

3063 | ans = [] |

3064 | for z in points: |

3065 | s = pari(0) |

3066 | r0 = (c*z).exp() |

3067 | r = r0 |

3068 | for n in xrange(1,prec): |

3069 | s += an[n-1]*r |

3070 | r *= r0 |

3071 | ans.append(s.python()) |

3072 | return ans |

3073 | |

3074 | def _multiple_of_degree_of_isogeny_to_optimal_curve(self): |

3075 | """ |

3076 | Internal function returning an integer m such that the degree of |

3077 | the isogeny between this curve and the optimal curve in its |

3078 | isogeny class is a divisor of m. |

3079 | |

3080 | EXAMPLES: |

3081 | sage: E=EllipticCurve('11a1') |

3082 | sage: E._multiple_of_degree_of_isogeny_to_optimal_curve() |

3083 | 25 |

3084 | sage: E=EllipticCurve('11a2') |

3085 | sage: E._multiple_of_degree_of_isogeny_to_optimal_curve() |

3086 | 5 |

3087 | sage: E=EllipticCurve('11a3') |

3088 | sage: E._multiple_of_degree_of_isogeny_to_optimal_curve() |

3089 | 5 |

3090 | """ |

3091 | M = self.isogeny_class()[1] |

3092 | return Integer(misc.prod([x for x in M.row(0) if x], 1)) |

3093 | |

3094 | ######################################################################## |

3095 | # Functions related to bounding the order of Sha (provably correctly!) |

3096 | # Heegner points and Kolyvagin's theorem |

3097 | ######################################################################## |

3098 | |

3099 | def sha(self): |

3100 | """ |

3101 | Return an object of class |

3102 | 'sage.schemes.elliptic_curves.sha.Sha' attached to this |

3103 | elliptic curve. |

3104 | |

3105 | This can be used in functions related to bounding the order of |

3106 | Sha (The Tate-Shafarevich group of the curve). |

3107 | |

3108 | EXAMPLES: |

3109 | sage: E=EllipticCurve('37a1') |

3110 | sage: S=E.sha() |

3111 | sage: S |

3112 | <class 'sage.schemes.elliptic_curves.sha.Sha'> |

3113 | sage: S.bound_kolyvagin() |

3114 | ([2], 1) |

3115 | """ |

3116 | try: |

3117 | return self.__sha |

3118 | except AttributeError: |

3119 | from sha import Sha |

3120 | self.__sha = Sha(self) |

3121 | return self.__sha |

3122 | |

3123 | def satisfies_heegner_hypothesis(self, D): |

3124 | """ |

3125 | Returns True precisely when D is a fundamental discriminant |

3126 | that satisfies the Heegner hypothesis for this elliptic curve. |

3127 | |

3128 | EXAMPLES: |

3129 | sage: E = EllipticCurve('11a1') |

3130 | sage: E.satisfies_heegner_hypothesis(-7) |

3131 | True |

3132 | sage: E.satisfies_heegner_hypothesis(-11) |

3133 | False |

3134 | """ |

3135 | if not number_field.is_fundamental_discriminant(D): |

3136 | return False |

3137 | if arith.GCD(D, self.conductor()) != 1: |

3138 | return False |

3139 | for p, _ in factor(self.conductor()): |

3140 | if kronecker_symbol(D,p) != 1: |

3141 | return False |

3142 | return True |

3143 | |

3144 | def heegner_discriminants(self, bound): |

3145 | """ |

3146 | Return the list of self's Heegner discriminants between -1 and -bound. |

3147 | |

3148 | INPUT: |

3149 | bound (int) -- upper bound for -discriminant |

3150 | |

3151 | OUTPUT: |

3152 | The list of Heegner discriminants between -1 and -bound for the given elliptic curve. |

3153 | |

3154 | EXAMPLE: |

3155 | sage: E=EllipticCurve('11a') |

3156 | sage: E.heegner_discriminants(30) |

3157 | [-7, -8, -19, -24] |

3158 | """ |

3159 | return [-D for D in xrange(1,bound) if self.satisfies_heegner_hypothesis(-D)] |

3160 | |

3161 | def heegner_discriminants_list(self, n): |

3162 | """ |

3163 | Return the list of self's first n Heegner discriminants smaller than -5. |

3164 | |

3165 | INPUT: |

3166 | n (int) -- the number of discriminants to compute |

3167 | |

3168 | OUTPUT: |

3169 | The list of the first n Heegner discriminants smaller than -5 for the given elliptic curve. |

3170 | |

3171 | EXAMPLE: |

3172 | sage: E=EllipticCurve('11a') |

3173 | sage: E.heegner_discriminants_list(4) |

3174 | [-7, -8, -19, -24] |

3175 | """ |

3176 | v = [] |

3177 | D = -5 |

3178 | while len(v) < n: |

3179 | while not self.satisfies_heegner_hypothesis(D): |

3180 | D -= 1 |

3181 | v.append(D) |

3182 | D -= 1 |

3183 | return v |

3184 | |

3185 | def heegner_point_height(self, D, prec=2): |

3186 | """ |

3187 | Use the Gross-Zagier formula to compute the Neron-Tate |

3188 | canonical height over K of the Heegner point corresponding to |

3189 | D, as an Interval (since it's computed to some precision using |

3190 | L-functions). |

3191 | |

3192 | INPUT: |

3193 | D (int) -- fundamental discriminant (=/= -3, -4) |

3194 | prec (int) -- (default: 2), use prec*sqrt(N) + 20 terms |

3195 | of L-series in computations, where N is the |

3196 | conductor. |

3197 | |

3198 | OUTPUT: |

3199 | Interval that contains the height of the Heegner point. |

3200 | |

3201 | EXAMPLE: |

3202 | sage: E = EllipticCurve('11a') |

3203 | sage: E.heegner_point_height(-7) |

3204 | [0.22226977 ... 0.22227479] |

3205 | """ |

3206 | |

3207 | if not self.satisfies_heegner_hypothesis(D): |

3208 | raise ArithmeticError, "Discriminant (=%s) must be a fundamental discriminant that satisfies the Heegner hypothesis."%D |

3209 | if D == -3 or D == -4: |

3210 | raise ArithmeticError, "Discriminant (=%s) must not be -3 or -4."%D |

3211 | eps = self.root_number() |

3212 | L1_vanishes = self.lseries().L1_vanishes() |

3213 | if eps == 1 and L1_vanishes: |

3214 | return IR(0) # rank even hence >= 2, so Heegner point is torsion. |

3215 | alpha = R(sqrt(abs(D)))/(2*self.period_lattice().complex_area()) |

3216 | F = self.quadratic_twist(D) |

3217 | E = self |

3218 | k_E = prec*sqrt(E.conductor()) + 20 |

3219 | k_F = prec*sqrt(F.conductor()) + 20 |

3220 | |

3221 | MIN_ERR = R('1e-6') # we assume that regulator and |

3222 | # discriminant, etc., computed to this accuracy (which is easily the case). |

3223 | # this should be made more intelligent / rigorous relative |

3224 | # to the rest of the system. |

3225 | if eps == 1: # E has even rank |

3226 | LF1, err_F = F.lseries().deriv_at1(k_F) |

3227 | LE1, err_E = E.lseries().at1(k_E) |

3228 | err_F = max(err_F, MIN_ERR) |

3229 | err_E = max(err_E, MIN_ERR) |

3230 | return IR(alpha-MIN_ERR,alpha+MIN_ERR) * IR(LE1-err_E,LE1+err_E) * IR(LF1-err_F,LF1+err_F) |

3231 | |

3232 | else: # E has odd rank |

3233 | LE1, err_E = E.lseries().deriv_at1(k_E) |

3234 | LF1, err_F = F.lseries().at1(k_F) |

3235 | err_F = max(err_F, MIN_ERR) |

3236 | err_E = max(err_E, MIN_ERR) |

3237 | return IR(alpha-MIN_ERR,alpha+MIN_ERR) * IR(LE1-err_E,LE1+err_E) * IR(LF1-err_F,LF1+err_F) |

3238 | |

3239 | |

3240 | def heegner_index(self, D, min_p=2, prec=5, verbose=False): |

3241 | r""" |

3242 | Return an interval that contains the index of the Heegner |

3243 | point $y_K$ in the group of K-rational points modulo torsion |

3244 | on this elliptic curve, computed using the Gross-Zagier |

3245 | formula and/or a point search, or the index divided by $2$. |

3246 | |

3247 | NOTES: If \code{min_p} is bigger than 2 then the index can be |

3248 | off by any prime less than \code{min_p}. This function |

3249 | returns the index divided by $2$ exactly when $E(\Q)_{/tor}$ |

3250 | has index $2$ in $E(K)_{/tor}$. |

3251 | |

3252 | INPUT: |

3253 | D (int) -- Heegner discriminant |

3254 | min_p (int) -- (default: 2) only rule out primes >= min_p |

3255 | dividing the index. |

3256 | verbose (bool) -- (default: False); print lots of mwrank |

3257 | search status information when computing |

3258 | regulator |

3259 | prec (int) -- (default: 5), use prec*sqrt(N) + 20 terms |

3260 | of L-series in computations, where N is the |

3261 | conductor. |

3262 | |

3263 | OUTPUT: |

3264 | an interval that contains the index |

3265 | |

3266 | EXAMPLES: |

3267 | sage: E = EllipticCurve('11a') |

3268 | sage: E.heegner_discriminants(50) |

3269 | [-7, -8, -19, -24, -35, -39, -40, -43] |

3270 | sage: E.heegner_index(-7) |

3271 | [0.99999332 .. 1.0000077] |

3272 | |

3273 | sage: E = EllipticCurve('37b') |

3274 | sage: E.heegner_discriminants(100) |

3275 | [-3, -4, -7, -11, -40, -47, -67, -71, -83, -84, -95] |

3276 | sage: E.heegner_index(-95) # long time (1 second) |

3277 | [1.9999923 .. 2.0000077] |

3278 | |

3279 | Current discriminants -3 and -4 are not supported: |

3280 | sage: E.heegner_index(-3) |

3281 | Traceback (most recent call last): |

3282 | ... |

3283 | ArithmeticError: Discriminant (=-3) must not be -3 or -4. |

3284 | |

3285 | The curve 681b returns an interval that contains $3/2$. |

3286 | This is because $E(\Q)$ is not saturated in $E(K)$. The |

3287 | true index is $3$: |

3288 | sage: E = EllipticCurve('681b') |

3289 | sage: I = E.heegner_index(-8); I |

3290 | [1.4999942 .. 1.5000058] |

3291 | sage: 2*I |

3292 | [2.9999885 .. 3.0000115] |

3293 | |

3294 | In fact, whenever the returned index has a denominator |

3295 | of $2$, the true index is got by multiplying the returned |

3296 | index by $2$. Unfortunately, this is not an if and only |

3297 | if condition, i.e., sometimes the index must be multiplied |

3298 | by $2$ even though the denominator is not $2$. |

3299 | """ |

3300 | # First compute upper bound on height of Heegner point. |

3301 | tm = misc.verbose("computing heegner point height...") |

3302 | h0 = self.heegner_point_height(D, prec=prec) |

3303 | |

3304 | # We divide by 2 to get the height **over Q** of the |

3305 | # Heegner point on the twist. |

3306 | |

3307 | ht = h0/2 |

3308 | misc.verbose('Height of heegner point = %s'%ht, tm) |

3309 | |

3310 | if self.root_number() == 1: |

3311 | F = self.quadratic_twist(D) |

3312 | else: |

3313 | F = self |

3314 | h = ht.upper() |

3315 | misc.verbose("Heegner height bound = %s"%h) |

3316 | B = F.CPS_height_bound() |

3317 | misc.verbose("CPS bound = %s"%B) |

3318 | c = h/(min_p**2) + B |

3319 | misc.verbose("Search would have to be up to height = %s"%c) |

3320 | |

3321 | if c > _MAX_HEIGHT or F is self: |

3322 | misc.verbose("Doing direct computation of MW group.") |

3323 | reg = F.regulator(verbose=verbose) |

3324 | return self.__adjust_heegner_index(ht/IR(reg)) |

3325 | |

3326 | # Do naive search to eliminate possibility that Heegner point |

3327 | # is divisible by p<min_p, without finding Heegner point. |

3328 | misc.verbose("doing point search") |

3329 | P = F.point_search(c, verbose=verbose) |

3330 | misc.verbose("done with point search") |

3331 | P = [x for x in P if x.order() == oo] |

3332 | if len(P) == 0: |

3333 | return IR(1) |

3334 | |

3335 | misc.verbose("saturating") |

3336 | S, I, reg = F.saturation(P, verbose=verbose) |

3337 | misc.verbose("done saturating") |

3338 | return self.__adjust_heegner_index(ht/IR(reg)) |

3339 | |

3340 | def __adjust_heegner_index(self, a): |

3341 | r""" |

3342 | Take the square root of the interval that contains the Heegner |

3343 | index. |

3344 | |

3345 | EXAMPLES: |

3346 | sage: E = EllipticCurve('11a1') |

3347 | sage: a = RIF(sqrt(2))-1.4142135623730951 |

3348 | sage: E._EllipticCurve_rational_field__adjust_heegner_index(a) |

3349 | [0.0000000... .. 1.490116...e-8] |

3350 | """ |

3351 | if a.lower() < 0: |

3352 | a = IR((0, a.upper())) |

3353 | return a.sqrt() |

3354 | |

3355 | |

3356 | def heegner_index_bound(self, D=0, prec=5, verbose=True, max_height=_MAX_HEIGHT): |

3357 | """ |

3358 | Assume self has rank 0. |

3359 | |

3360 | Return a list v of primes such that if an odd prime p divides |

3361 | the index of the Heegner point in the group of rational |

3362 | points *modulo torsion*, then p is in v. |

3363 | |

3364 | If 0 is in the interval of the height of the Heegner point |

3365 | computed to the given prec, then this function returns v = 0. |

3366 | This does not mean that the Heegner point is torsion, just |

3367 | that it is very likely torsion. |

3368 | |

3369 | If we obtain no information from a search up to max_height, e.g., |

3370 | if the Siksek et al. bound is bigger than max_height, then |

3371 | we return v = -1. |

3372 | |

3373 | INPUT: |

3374 | D (int) -- (deault: 0) Heegner discriminant; if 0, use the |

3375 | first discriminant < -4 that satisfies the Heegner |

3376 | hypothesis |

3377 | verbose (bool) -- (default: True) |

3378 | prec (int) -- (default: 5), use prec*sqrt(N) + 20 terms |

3379 | of L-series in computations, where N is the |

3380 | conductor. |

3381 | max_height (float) -- should be <= 21; bound on logarithmic |

3382 | naive height used in point searches. |

3383 | Make smaller to make this function |

3384 | faster, at the expense of possibly |

3385 | obtaining a worse answer. A good |

3386 | range is between 13 and 21. |

3387 | |

3388 | OUTPUT: |

3389 | v -- list or int (bad primes or 0 or -1) |

3390 | D -- the discriminant that was used (this is useful if D was |

3391 | automatically selected). |

3392 | |

3393 | EXAMPLES: |

3394 | sage: E = EllipticCurve('11a1') |

3395 | sage: E.heegner_index_bound(verbose=False) |

3396 | ([2], -7) |

3397 | """ |

3398 | max_height = min(float(max_height), _MAX_HEIGHT) |

3399 | if self.root_number() != 1: |

3400 | raise RuntimeError, "The rank must be 0." |

3401 | |

3402 | if D == 0: |

3403 | D = -5 |

3404 | while not self.satisfies_heegner_hypothesis(D): |

3405 | D -= 1 |

3406 | |

3407 | # First compute upper bound on Height of Heegner point. |

3408 | ht = self.heegner_point_height(D, prec=prec) |

3409 | if 0 in ht: |

3410 | return 0, D |

3411 | F = self.quadratic_twist(D) |

3412 | h = ht.upper() |

3413 | misc.verbose("Heegner height bound = %s"%h) |

3414 | B = F.CPS_height_bound() |

3415 | misc.verbose("CPS bound = %s"%B) |

3416 | H = h |

3417 | p = 3 |

3418 | while True: |

3419 | c = h/(2*p**2) + B |

3420 | if c < max_height: |

3421 | break |

3422 | if p > 100: |

3423 | break |

3424 | p = next_prime(p) |

3425 | misc.verbose("Using p = %s"%p) |

3426 | |

3427 | if c > max_height: |

3428 | misc.verbose("No information by searching only up to max_height (=%s)."%c) |

3429 | return -1, D |

3430 | |

3431 | misc.verbose("Searching up to height = %s"%c) |

3432 | eps = 10e-5 |

3433 | |

3434 | def _bound(P): |

3435 | """ |

3436 | We will use this function below in two places. It |

3437 | bounds the index using a nontrivial point. |

3438 | """ |

3439 | assert len(P) == 1 |

3440 | |

3441 | S, I, reg = F.saturation(P, verbose=verbose) |

3442 | h = IR(reg-eps,reg+eps) |

3443 | ind2 = ht/(h/2) |

3444 | misc.verbose("index squared = %s"%ind2) |

3445 | ind = ind2.sqrt() |

3446 | misc.verbose("index = %s"%ind) |

3447 | # Compute upper bound on square root of index. |

3448 | if ind.absolute_diameter() < 1: |

3449 | t, i = ind.is_int() |

3450 | if t: # unique integer in interval, so we've found exact index squared. |

3451 | return arith.prime_divisors(i), D |

3452 | raise RuntimeError, "Unable to compute bound for e=%s, D=%s (try increasing precision)"%(self,D) |

3453 | |

3454 | # First try a quick search, in case we get lucky and find |

3455 | # a generator. |

3456 | P = F.point_search(13, verbose=verbose) |

3457 | P = [x for x in P if x.order() == oo] |

3458 | if len(P) > 0: |

3459 | return _bound(P) |

3460 | |

3461 | # Do search to eliminate possibility that Heegner point is |

3462 | # divisible by primes up to p, without finding Heegner point. |

3463 | P = F.point_search(c, verbose=verbose) |

3464 | P = [x for x in P if x.order() == oo] |

3465 | if len(P) == 0: |

3466 | # We've eliminated the possibility of a divisor up to p. |

3467 | return arith.prime_range(3,p), D |

3468 | else: |

3469 | return _bound(P) |

3470 | |

3471 | padic_regulator = padics.padic_regulator |

3472 | |

3473 | padic_height_pairing_matrix = padics.padic_height_pairing_matrix |

3474 | |

3475 | padic_height = padics.padic_height |

3476 | padic_height_via_multiply = padics.padic_height_via_multiply |

3477 | |

3478 | padic_sigma = padics.padic_sigma |

3479 | padic_sigma_truncated = padics.padic_sigma_truncated |

3480 | |

3481 | padic_E2 = padics.padic_E2 |

3482 | |

3483 | matrix_of_frobenius = padics.matrix_of_frobenius |

3484 | |

3485 | # def weierstrass_p(self): |

3486 | # # TODO: add allowing negative valuations for power series |

3487 | # return 1/t**2 + a1/t + rings.frac(1,12)*(a1-8*a2) -a3*t \ |

3488 | # - (a4+a1*a3)*t**2 + O(t**3) |

3489 | |

3490 | |

3491 | def mod5family(self): |

3492 | """ |

3493 | Return the family of all elliptic curves with the same mod-5 |

3494 | representation as self. |

3495 | |

3496 | EXAMPLES: |

3497 | sage: E=EllipticCurve('32a1') |

3498 | sage: E.mod5family() |

3499 | Elliptic Curve defined by y^2 = x^3 + 4*x over Fraction Field of Univariate Polynomial Ring in t over Rational Field |

3500 | """ |

3501 | E = self.short_weierstrass_model() |

3502 | a = E.a4() |

3503 | b = E.a6() |

3504 | return mod5family.mod5family(a,b) |

3505 | |

3506 | def tate_curve(self, p): |

3507 | r""" |

3508 | Creates the Tate Curve over the $p$-adics associated to this elliptic curves. |

3509 | |

3510 | This Tate curve a $p$-adic curve with split multiplicative |

3511 | reduction of the form $y^2+xy=x^3+s_4 x+s_6$ which is |

3512 | isomorphic to the given curve over the algebraic closure of |

3513 | $\QQ_p$. Its points over $\QQ_p$ are isomorphic to |

3514 | $\QQ_p^{\times}/q^{\Z}$ for a certain parameter $q\in\Z_p$. |

3515 | |

3516 | INPUT: |

3517 | |

3518 | p -- a prime where the curve has multiplicative reduction. |

3519 | |

3520 | |

3521 | EXAMPLES: |

3522 | sage: e = EllipticCurve('130a1') |

3523 | sage: e.tate_curve(2) |

3524 | 2-adic Tate curve associated to the Elliptic Curve defined by y^2 + x*y + y = x^3 - 33*x + 68 over Rational Field |

3525 | |

3526 | The input curve must have multiplicative reduction at the prime. |

3527 | sage: e.tate_curve(3) |

3528 | Traceback (most recent call last): |

3529 | ... |

3530 | ValueError: The elliptic curve must have multiplicative reduction at 3 |

3531 | |

3532 | We compute with $p=5$: |

3533 | sage: T = e.tate_curve(5); T |

3534 | 5-adic Tate curve associated to the Elliptic Curve defined by y^2 + x*y + y = x^3 - 33*x + 68 over Rational Field |

3535 | |

3536 | We find the Tate parameter $q$: |

3537 | sage: T.parameter(prec=5) |

3538 | 3*5^3 + 3*5^4 + 2*5^5 + 2*5^6 + 3*5^7 + O(5^8) |

3539 | |

3540 | We compute the $L$-invariant of the curve: |

3541 | sage: T.L_invariant(prec=10) |

3542 | 5^3 + 4*5^4 + 2*5^5 + 2*5^6 + 2*5^7 + 3*5^8 + 5^9 + O(5^10) |

3543 | """ |

3544 | try: |

3545 | return self._tate_curve[p] |

3546 | except AttributeError: |

3547 | self._tate_curve = {} |

3548 | except KeyError: |

3549 | pass |

3550 | |

3551 | Eq = ell_tate_curve.TateCurve(self,p) |

3552 | self._tate_curve[p] = Eq |

3553 | return Eq |

3554 | |

3555 | |

3556 | def integral_points(self, mw_base='auto',tors_points='auto'): |

3557 | """ |

3558 | Computes all integral points on the elliptic curve E which has Mordell-Weil |

3559 | base mw_base and torsions points specified by tors_points. |

3560 | Default 'auto' relies on self.gens() and self.torsion_subgroup().gens() |

3561 | |

3562 | INPUT: |

3563 | self -- EllipticCurve_Rational_Field |

3564 | mw_base -- list of EllipticCurvePoint representing the Mordell-Weil base |

3565 | (default: 'auto' - calls self.gens()) |

3566 | tors_points -- list of EllipticCurvePoint giving all torsion points of E |

3567 | (default: 'auto' - calls self.torsion_subgroup().gens()) |

3568 | |

3569 | OUTPUT: |

3570 | set of all integral points on E |

3571 | |

3572 | HINTS: |

3573 | - The complexity increases exponentially in the rank of curve E. |

3574 | - It can help if you try another Mordell-Weil base, because the |

3575 | computation time depends on this, too. |

3576 | |

3577 | EXAMPLES: |

3578 | A curve of rank 3 with no torsion points |

3579 | |

3580 | sage: E=EllipticCurve([0,0,1,-7,6]) |

3581 | sage: P1=E.point((2,0)); P2=E.point((-1,3)); P3=E.point((4,6)) |

3582 | sage: a=E.integral_points([P1,P2,P3], []); a |

3583 | set([(3 : -4 : 1), (21 : 95 : 1), (-3 : 0 : 1), (21 : -96 : 1), (93 : -897 : 1), (3 : 3 : 1), (-1 : 3 : 1), (8 : 21 : 1), (-3 : -1 : 1), (342 : 6324 : 1), (37 : -225 : 1), (37 : 224 : 1), (11 : -36 : 1), (406 : -8181 : 1), (4 : -7 : 1), (2 : 0 : 1), (52 : 374 : 1), (93 : 896 : 1), (14 : 51 : 1), (14 : -52 : 1), (8 : -22 : 1), (4 : 6 : 1), (-2 : -4 : 1), (816 : -23310 : 1), (11 : 35 : 1), (2 : -1 : 1), (0 : 2 : 1), (1 : -1 : 1), (-2 : 3 : 1), (0 : -3 : 1), (406 : 8180 : 1), (52 : -375 : 1), (816 : 23309 : 1), (-1 : -4 : 1), (1 : 0 : 1), (342 : -6325 : 1)]) |

3584 | |

3585 | It is also possible not to specify mw_base and tors_points but because the |

3586 | Mordell-Weil base found by built-in functions differs this takes longer |

3587 | |

3588 | sage: E=EllipticCurve([0,0,1,-7,6]) |

3589 | sage: a=E.integral_points(); a |

3590 | set([(3 : -4 : 1), (21 : 95 : 1), (-3 : 0 : 1), (21 : -96 : 1), (93 : -897 : 1), (3 : 3 : 1), (-1 : 3 : 1), (4 : 6 : 1), (-3 : -1 : 1), (342 : 6324 : 1), (37 : -225 : 1), (37 : 224 : 1), (11 : -36 : 1), (406 : -8181 : 1), (4 : -7 : 1), (2 : 0 : 1), (52 : 374 : 1), (93 : 896 : 1), (14 : 51 : 1), (14 : -52 : 1), (8 : -22 : 1), (8 : 21 : 1), (-2 : -4 : 1), (816 : -23310 : 1), (11 : 35 : 1), (2 : -1 : 1), (0 : 2 : 1), (1 : -1 : 1), (-2 : 3 : 1), (0 : -3 : 1), (406 : 8180 : 1), (52 : -375 : 1), (816 : 23309 : 1), (-1 : -4 : 1), (1 : 0 : 1), (342 : -6325 : 1)]) |

3591 | |

3592 | |

3593 | Another example with rank 5 and no torsion points |

3594 | |

3595 | sage: E=EllipticCurve([-879984,319138704]) |

3596 | sage: P1=E.point((540,1188)); P2=E.point((576,1836)) |

3597 | sage: P3=E.point((468,3132)); P4=E.point((612,3132)) |

3598 | sage: P5=E.point((432,4428)) |

3599 | sage: a=E.integral_points([P1,P2,P3,P4,P5],[]); a # long time (!) |

3600 | set([(155356 : -61232804 : 1), (2717460 : -4479656508 : 1), (468 : -3132 : 1), (513 : -1647 : 1), (-764 : 23356 : 1), (155356 : 61232804 : 1), (65052 : -16589988 : 1), (792 : 10908 : 1), (2204663452 : -103517423460188 : 1), (36828 : -7065252 : 1), (1732 : 63172 : 1), (-792 : -22788 : 1), (673 : -5633 : 1), (-792 : 22788 : 1), (2592 : 124308 : 1), (-620 : 25028 : 1), (720 : -7668 : 1), (1368 : 40932 : 1), (-620 : -25028 : 1), (540 : -1188 : 1), (1072 : 24652 : 1), (108 : 15012 : 1), (1260648 : -1415437308 : 1), (612 : -3132 : 1), (576 : -1836 : 1), (1072 : -24652 : 1), (1260648 : 1415437308 : 1), (-72 : -19548 : 1), (432 : 4428 : 1), (396 : -5724 : 1), (520 : -1468 : 1), (193320 : 84998268 : 1), (-684 : 24516 : 1), (972 : 19548 : 1), (-900 : 19548 : 1), (-423 : -24813 : 1), (43200 : 8976852 : 1), (9972 : -991548 : 1), (14265 : 1700163 : 1), (78696 : 22074876 : 1), (376 : 6436 : 1), (196992 : 87431508 : 1), (-684 : -24516 : 1), (585 : 2133 : 1), (792 : -10908 : 1), (4680 : 314172 : 1), (-1080 : 3132 : 1), (2448 : -113292 : 1), (1732 : -63172 : 1), (-279 : 23301 : 1), (193320 : -84998268 : 1), (3060 : -162108 : 1), (-72 : 19548 : 1), (172 : -13148 : 1), (376 : -6436 : 1), (2385 : 108567 : 1), (36828 : 7065252 : 1), (114280 : -38631428 : 1), (612 : 3132 : 1), (540 : 1188 : 1), (114280 : 38631428 : 1), (4680 : -314172 : 1), (5940 : -452412 : 1), (2448 : 113292 : 1), (-423 : 24813 : 1), (-864 : 20844 : 1), (-279 : -23301 : 1), (8668 : 802468 : 1), (720 : 7668 : 1), (17856 : 2382804 : 1), (2385 : -108567 : 1), (-1080 : -3132 : 1), (432 : -4428 : 1), (3060 : 162108 : 1), (513 : 1647 : 1), (972 : -19548 : 1), (3537 : 203607 : 1), (368892 : 224051292 : 1), (468 : 3132 : 1), (396 : 5724 : 1), (17856 : -2382804 : 1), (2113 : -88847 : 1), (9972 : 991548 : 1), (14265 : -1700163 : 1), (43200 : -8976852 : 1), (36 : 16956 : 1), (108 : -15012 : 1), (-864 : -20844 : 1), (-900 : -19548 : 1), (2717460 : 4479656508 : 1), (585 : -2133 : 1), (8668 : -802468 : 1), (2204663452 : 103517423460188 : 1), (196992 : -87431508 : 1), (1368 : -40932 : 1), (78696 : -22074876 : 1), (172 : 13148 : 1), (576 : 1836 : 1), (3537 : -203607 : 1), (2592 : -124308 : 1), (520 : 1468 : 1), (2113 : 88847 : 1), (368892 : -224051292 : 1), (65052 : 16589988 : 1), (36 : -16956 : 1), (673 : 5633 : 1), (-764 : -23356 : 1), (5940 : 452412 : 1)]) |

3601 | |

3602 | |

3603 | NOTES: |

3604 | - This function uses the algorithm given in [Co1] and the computation of |

3605 | the complex elliptic logarithm as presented in [Co2]. |

3606 | REFERENCES: |

3607 | -[Co1] Cohen H., Number Theory Vol I: Tools and Diophantine Equations |

3608 | GTM 239, Springer 2007 |

3609 | -[Co2] Cohen H., A Course in Computational Algebraic Number Theory |

3610 | GTM 138, Springer 1996 |

3611 | AUTHORS: |

3612 | - Michael Mardaus (2008-07) |

3613 | - Tobias Nagel (2008-07) |

3614 | """ |

3615 | ############################################################################### |

3616 | # INPUT CHECK ################################################################# |

3617 | if mw_base=='auto': |

3618 | mw_base = self.gens() |

3619 | if tors_points=='auto': |

3620 | tors_points = self.torsion_subgroup().gens() |

3621 | try: |

3622 | len_tors = len(tors_points) |

3623 | except TypeError: |

3624 | raise TypeError, "tors_points must be a list" |

3625 | try: |

3626 | r = len(mw_base) |

3627 | except TypeError: |

3628 | raise TypeError, 'mw_base must be a list' |

3629 | if (r == 0) and (len_tors == 0): |

3630 | raise RuntimeError, 'Both base points and torsions points are not specified' |

3631 | try: # test if all given points are on curve (try not necessary because of no new error) |

3632 | if r < len_tors: |

3633 | for i in range(r): |

3634 | trash = self.point(mw_base[i]) |

3635 | trash = self.point(tors_points[i]) |

3636 | for i in range(r,len_tors): |

3637 | trash = self.point(tors_points[i]) |

3638 | else: |

3639 | for i in range(len_tors): |

3640 | trash = self.point(mw_base[i]) |

3641 | trash = self.point(tors_points[i]) |

3642 | for i in range(len_tors,r): |

3643 | trash = self.point(mw_base[i]) |

3644 | except TypeError: |

3645 | raise |

3646 | #END INPUT-CHECK############################################################### |

3647 | ############################################################################### |

3648 | |

3649 | ############################################################################### |

3650 | # INTERNAL FUNCTIONS ########################################################## |

3651 | |

3652 | ############################## begin ################################ |

3653 | def is_int(a): |

3654 | if floor(a)==ceil(a): |

3655 | return 1 |

3656 | else: |

3657 | return 0 |

3658 | ############################## end ################################ |

3659 | ############################## begin ################################ |

3660 | def height_of_curve(): |

3661 | #Computes height of curve 'self' according to Co1 |

3662 | if j == 0: |

3663 | h_j = 1 |

3664 | else: |

3665 | h_j = max(log(abs(j.numerator()), e), log(abs(j.denominator()), e)) |

3666 | if (g2 != 0) and (g3 != 0): |

3667 | h_gs = max(1, log(abs(g2), e), log(abs(g3), e)) |

3668 | elif g2 == 0: |

3669 | if g3 == 0: |

3670 | return max(1,h_j) |

3671 | h_gs = max(1, log(abs(g3), e)) |

3672 | else: |

3673 | h_gs = max(1, log(abs(g2), e)) |

3674 | return max(1,h_j, h_gs) |

3675 | ############################## end ################################ |

3676 | ############################## begin ################################ |

3677 | def log_plus(x): |

3678 | if x==0: |

3679 | return 1 |

3680 | else: |

3681 | return max(1, log(abs(x), e)) |

3682 | ############################## end ################################ |

3683 | ############################## begin ################################ |

3684 | def height_paaring_matrix(list): |

3685 | #computes the height paaring matrix associated with the bilinear form of the |

3686 | #canonical (Nero-Tate) height function |

3687 | point_height = [] |

3688 | dim=len(list) |

3689 | for i in range(dim): |

3690 | point_height.append(list[i].height()) |

3691 | M = matrix.MatrixSpace(R, dim) |

3692 | mat = M() |

3693 | for j in range(dim): |

3694 | for k in range(j,dim): |

3695 | mat[j,k]=((list[j]+list[k]).height() - point_height[j] - point_height[k])/2 |

3696 | if j!=k: |

3697 | mat[k,j]=mat[j,k] |

3698 | return mat |

3699 | ############################## end ################################ |

3700 | ############################## begin ################################ |

3701 | def min_eigenvalue(M): |

3702 | #computes the minimal eigenvalue of M |

3703 | eigenvalues = [] |

3704 | eigenvalues =M.charpoly().real_roots(); |

3705 | return min(eigenvalues) |

3706 | ############################## end ################################ |

3707 | ############################## begin ################################ |

3708 | def prod(list): |

3709 | #multiplies all elements in list |

3710 | tmp = 1 |

3711 | for i in range(len(list)): |

3712 | tmp = tmp * list[i] |

3713 | return tmp |

3714 | ############################## end ################################ |

3715 | ############################## begin ################################ |

3716 | def extract_realroots(list, eps = 1e-8): |

3717 | #function is used to extract real roots from result of |

3718 | #function roots() listed in list |

3719 | #at the same time cutting off impreciseness of roots() |

3720 | #RETURN: sorted list of roots in RR |

3721 | tmp = [] |

3722 | for i in range(0, len(list)): |

3723 | im = list[i][0].imag() |

3724 | re = list[i][0].real() |

3725 | if ((im < eps) and (im > -eps)): #imaginary part in eps |

3726 | if ((re < eps) and (re > -eps)): #both parts are imprecise (also 0 included) |

3727 | tmp.append(0) |

3728 | else: |

3729 | tmp.append(re) |

3730 | tmp.sort() |

3731 | return tmp |

3732 | ############################## end ################################ |

3733 | ############################## begin ################################ |

3734 | def extract_roots(list, eps = 1e-8): |

3735 | #function is used to extract roots from result of |

3736 | #function roots() listed in list |

3737 | #cutting off impreciseness of roots() |

3738 | #RETURN: sorted list of roots in C |

3739 | tmp = [] |

3740 | for i in range(0, len(list)): |

3741 | im = list[i][0].imag() |

3742 | re = list[i][0].real() |

3743 | if ((im < eps) and (im > -eps)): |

3744 | if ((re < eps) and (re > -eps)): #both parts are imprecise (also 0 included) |

3745 | tmp.append(0 + 0*I) |

3746 | else: |

3747 | tmp.append(re + 0*I) |

3748 | elif ((re < eps) and (re > -eps)): |

3749 | tmp.append(im*I) |

3750 | else: #no impreciseness |

3751 | tmp.append(list[i][0]) |

3752 | return tmp |

3753 | ############################## end ################################ |

3754 | ############################## begin ################################ |

3755 | def in_egg(P,e1,e2): |

3756 | #test if P is on compact component of E |

3757 | #only possible for disc(E)>0, otherwise there is no compact component |

3758 | if (P.xy()[0] >= e1) and (P.xy()[0] <= e2): |

3759 | return True |

3760 | else: |

3761 | return False |

3762 | ############################## end ################################ |

3763 | ############################## begin ################################ |

3764 | def point_preprocessing(list,e1,e2): |

3765 | #Transforms mw_base so that at most one point is on the |

3766 | #compact component of the curve |

3767 | Q = [] |

3768 | mem = -1 |

3769 | for i in range(0,len(list)): |

3770 | if in_egg(list[i],e1,e2): |

3771 | if mem == -1: |

3772 | mem = i |

3773 | else: |

3774 | Q.append(list[i]+list[mem]) |

3775 | mem = i |

3776 | else: |

3777 | Q.append(list[i]) |

3778 | if mem != -1: #add last point, which is not in egg, to Q |

3779 | Q.append(2*list[mem]) |

3780 | return Q |

3781 | ############################## end ################################ |

3782 | ############################## begin ################################ |

3783 | def modified_height(i):#computes modified height if base point i |

3784 | return max(mw_base[i].height(),h_E,c7*mw_base_log[i]**2) |

3785 | ############################## end ################################ |

3786 | ############################## begin ################################ |

3787 | def complex_elliptic_log(E,P): |

3788 | #Algorithm according to Co2 |

3789 | x = P.xy()[0] |

3790 | y = P.xy()[1] |

3791 | #Initialize |

3792 | a1 = E.a1() |

3793 | a3 = E.a3() |

3794 | b2 = E.b2() |

3795 | b4 = E.b4() |

3796 | b6 = E.b6() |

3797 | if disc < 0: |

3798 | #Connected Case |

3799 | PZ = rings.PolynomialRing(C, 'var') |

3800 | var = PZ.gen() |

3801 | tmp = (4*var**3+b2*var**2+2*b4*var+b6).roots() |

3802 | real_roots = extract_realroots(tmp) |

3803 | if len(real_roots) != 1: |

3804 | raise ValueError, ' no or more than one real root despite disc < 0' |

3805 | else: |

3806 | e1 = real_roots[0] |

3807 | beta = sqrt(3*e1**2 + (b2/2)*e1 + b4/2) |

3808 | alpha = 3*e1 + b2/4 |

3809 | a = R(2*sqrt(beta)) |

3810 | b = R(sqrt(alpha + 2*beta)) |

3811 | c = R((x - e1 + beta)/sqrt(x - e1)) |

3812 | while (a - b) > 1e-12: |

3813 | a = R((a + b)/2) |

3814 | b = R(sqrt(a*b)) |

3815 | c = R((c + sqrt(c**2 + b**2 - a**2))/2) |

3816 | if (2*y + a1*x + a3)*((x - e1)**2 - beta**2) < 0: |

3817 | z = R(arcsin(a/c)/a) |

3818 | else: |

3819 | z = R((pi - arcsin(a/c))/a) |

3820 | if (2*y + a1*x + a3) > 0: |

3821 | z = z + pi/a |

3822 | return z |

3823 | else: |

3824 | a2 = E.a2() |

3825 | #Disconnected Case |

3826 | PZ = rings.PolynomialRing(C, 'var') |

3827 | var = PZ.gen() |

3828 | tmp = (4*var**3+b2*var**2+2*b4*var+b6).roots() |

3829 | real_roots = extract_realroots(tmp) |

3830 | if len(real_roots)!=3: |

3831 | raise ValueError, 'Not excatly 3 real roots found despite disc > 0' |

3832 | else: |

3833 | e1 = real_roots[2] |

3834 | e2 = real_roots[1] |

3835 | e3 = real_roots[0] |

3836 | a = R(sqrt(e1 - e3)) |

3837 | b = R(sqrt(e1 - e2)) |

3838 | if x < e1: |

3839 | f = 1 |

3840 | lam = y/(x - e3) |

3841 | x = lam**2 + a1*lam - a2 - x - e3 |

3842 | else: |

3843 | f = 0 |

3844 | c = sqrt(x - e3) |

3845 | while (a-b) > 1e-12: |

3846 | a = R((a + b)/2) |

3847 | b = R(sqrt(a*b)) |

3848 | c = R((c + sqrt(c**2 + b**2 - a**2))/2) |

3849 | #Connected component |

3850 | if (((f == 0) and (2*y + a1*x + a3) < 0) or ((f == 1) and (2*y + a1*x + a3) >= 0)): |

3851 | z = R(arcsin(a/c)/a) |

3852 | else: |

3853 | z = R((pi - arcsin(a/c))/a) |

3854 | if f == 0: |

3855 | return z |

3856 | else: |

3857 | #other component |

3858 | w2=E.period_lattice().basis()[1] |

3859 | return (z + w2/2) |

3860 | ############################## end ################################ |

3861 | ############################## begin ################################ |

3862 | def search_points(lowerb, upperb): |

3863 | # Search of integral points P on E with x(P) between lowerbd |

3864 | # and upperbd, using the explicit defintion of E |

3865 | a = self.a_invariants() |

3866 | for x in range(lowerb,upperb): |

3867 | rs = x**3 + a[1]*x**2 + a[3]*x + a[4] # right hand side of E |

3868 | y_coef = a[0]*x + a[2] #coefficient at y |

3869 | disc = y_coef**2 + 4*rs |

3870 | if disc > 0: |

3871 | sol1 = R((-y_coef + sqrt(disc))/2) |

3872 | sol2 = R((-y_coef - sqrt(disc))/2) |

3873 | if is_int(sol1): #tests if sol1 is integral. |

3874 | int_points.add(self.point((x,sol1))) |

3875 | if is_int(sol2): |

3876 | int_points.add(self.point((x,sol2))) |

3877 | elif disc == 0: |

3878 | sol1 = R((-y_coef + sqrt(disc))/2) |

3879 | if is_int(sol1): #tests if sol1 is integral. |

3880 | int_points.add(self.point((x,sol1))) |

3881 | ############################## end ################################ |

3882 | ############################## begin ################################ |

3883 | def search_remaining_points(): |

3884 | #search integral points on curve E written as linear combination |

3885 | #of n times the mordell-weil base points and torsion points |

3886 | #(n is bounded by H_q, which will be computed at runtime) |

3887 | ###################### begin ######################## |

3888 | #Computes b-adic representation of num |

3889 | #OUTPUT: list with b-adic representation list[0]=base^0 ... list[i]=base^i |

3890 | def badic_repres(num,base): |

3891 | new_num=[] |

3892 | for i in range(r): |

3893 | new_num.append(num%base) |

3894 | num=num//base |

3895 | return new_num |

3896 | ###################### end ######################## |

3897 | OP=self.point((0,1,0)) |

3898 | #qw as tmp memory, for coefficients in linear combination |

3899 | #qw = [] |

3900 | #for go in range(r): |

3901 | # qw.append('init') |

3902 | #end |

3903 | for i in range(1,((2*H_q+1)**r)/2): |

3904 | koeffs = badic_repres(i,2*H_q+1) |

3905 | P = self.point((0,1,0)) |

3906 | for index in range(r): |

3907 | P = P + (koeffs[index]-H_q)*mw_base[index] |

3908 | #qw[index]=koeffs[index]-H_q #s.a. |

3909 | if P!=OP:#P itself could be an integral point (torsion point would be (0:1:0) ) |

3910 | if is_int(P.xy()[0]) and is_int(P.xy()[1]): |

3911 | int_points.add(P) |

3912 | int_points.add(-P) |

3913 | #print P, qw |

3914 | for j in range(len_tors): # P + torsion points |

3915 | P = P + tors_points[j] |

3916 | if P!=OP: |

3917 | if is_int(P.xy()[0]) and is_int(P.xy()[1]): |

3918 | int_points.add(P) |

3919 | int_points.add(-P) |

3920 | ############################## end ################################ |

3921 | |

3922 | # END Internal functions ###################################################### |

3923 | ############################################################################### |

3924 | |

3925 | e = exp(1) |

3926 | I = constants.I |

3927 | pi = R(constants.pi) |

3928 | |

3929 | g2 = R((self.c4()/12)) |

3930 | g3 = R((self.c6()/216)) |

3931 | disc = self.discriminant() |

3932 | j = self.j_invariant() |

3933 | b2 = self.b2() |

3934 | |

3935 | PR = rings.PolynomialRing(C, 'var') |

3936 | var = PR.gen() |

3937 | tmp = (4*var**3-g2*var-g3).roots() |

3938 | if disc > 0: #compact and noncompact component -> 3 roots in RR |

3939 | roots = extract_realroots(tmp) |

3940 | e1 = roots[0] |

3941 | e2 = roots[1] |

3942 | e3 = roots[2] |

3943 | mw_base = point_preprocessing(mw_base,e1,e2) #at most one point in E^{egg} |

3944 | |

3945 | elif disc < 0: #only noncompact component => 1 root in RR (=: e3), 2 roots in C |

3946 | roots = extract_roots(tmp) |

3947 | if roots[0].imag() == 0: |

3948 | e3 = roots[0].real() |

3949 | e2 = roots[1] |

3950 | e1 = roots[2] |

3951 | elif roots[1].imag() == 0: |

3952 | e3 = roots[1].real() |

3953 | e2 = roots[0] |

3954 | e1 = roots[2] |

3955 | else: |

3956 | e3 = roots[2].real() |

3957 | e2 = roots[1] |

3958 | e1 = roots[0] |

3959 | |

3960 | #curve of rank 0 doesn't need the whole algorithm |

3961 | if r == 0: |

3962 | int_points = set([]) |

3963 | if disc > 0: |

3964 | #Points in egg have e2>=x>=e1 |

3965 | search_points(ceil(e1), floor(e2)+1) |

3966 | #Points in noncompact component with x(P)< 2*max(|e1|,|e2|,|e3|) , espec. x(P)>=e3 |

3967 | search_points(ceil(e3), ceil(2*max(abs(e1),abs(e2),abs(e3)))) |

3968 | OP=self.point((0,1,0)) |

3969 | for i in range(len_tors): |

3970 | P_tmp=tors_points[i] |

3971 | while P_tmp != OP: |

3972 | if is_int(P_tmp.xy()[0]) and is_int(P_tmp.xy()[1]): |

3973 | int_points.add(P_tmp) |

3974 | P_tmp = P_tmp + tors_points[i] |

3975 | return list(int_points) |

3976 | # Algorithm presented in [Co1] |

3977 | h_E = height_of_curve() |

3978 | w1 = self.period_lattice().basis()[0] |

3979 | w2 = self.period_lattice().basis()[1] |

3980 | mu = (log((abs(disc)),e) + log_plus(j))/6 + log_plus(b2/12) |

3981 | if b2!=0: |

3982 | mu=mu+log(2,e) |

3983 | |

3984 | c1 = R(exp(mu + 2.14)) |

3985 | c2 = min_eigenvalue(height_paaring_matrix(mw_base)) |

3986 | c3 = (w1**2)*abs(b2)/48 + 8 |

3987 | c5 = R(sqrt(c1*c3)) |

3988 | c7 = abs((3*pi)/((w1**2) * (w1/w2).imag())) |

3989 | |

3990 | mw_base_log = [] #contains \Phi(Q_i) |

3991 | mod_h_list = [] #contains h_m(Q_i) |

3992 | c9_help_list = [] |

3993 | for i in range(0,r): |

3994 | mw_base_log.append(abs(complex_elliptic_log(self,mw_base[i]))) |

3995 | mod_h_list.append(modified_height(i)) |

3996 | c9_help_list.append(sqrt(mod_h_list[i])/mw_base_log[i]) |

3997 | |

3998 | c8 = max(e*h_E,max(mod_h_list)) |

3999 | c9 = e/sqrt(c7) * min(c9_help_list) |

4000 | n=r+1 |

4001 | c10 = R(2 * 10**(8+7*n) * R((2/e)**(2 * n**2)) * (n+1)**(4 * n**2 + 10 * n) * log(c9,e)**(-2*n - 1) * prod(mod_h_list)) |

4002 | |

4003 | top = 128 #arbitrary first upper bound |

4004 | bottom = 0 |

4005 | log_c9=log(c9,e); log_c5=log(c5,e) |

4006 | log_r_top = log(r*(1e1)**top, e) |

4007 | |

4008 | while R(c10*(log_r_top+log_c9)*(log(log_r_top,e)+h_E+log_c9)**(n+1)) > R(c2/2 * ((1e1)**top)**2 - log_c5): |

4009 | #initial bound 'top' too small, upshift of search interval |

4010 | bottom = top |

4011 | top = 2*top |

4012 | while top >= bottom: #binary-search like search for fitting exponent (bound) |

4013 | bound = floor(bottom + (top - bottom)/2) |

4014 | log_r_bound = log(r*(1e1)**bound, e) |

4015 | if R(c10*(log_r_bound+log_c9)*(log(log_r_bound, e)+h_E+log_c9)**(n+1)) > R(c2/2 * ((1e1)**bound)**2 - log_c5): |

4016 | bottom = bound + 1 |

4017 | else: |

4018 | top = bound - 1 |

4019 | |

4020 | H_q = R((1e1)**bound) |

4021 | break_cond = 0 #at least one reduction step |

4022 | #reduction via LLL |

4023 | while break_cond < 0.9: #as long as the improvement of the new bound in comparison to the old is greater than 10% |

4024 | c = R((H_q**n)*10) #c has to be greater than H_q^n |

4025 | M = matrix.MatrixSpace(Z,n) |

4026 | m = M.identity_matrix() |

4027 | for i in range(r): |

4028 | m[i, r] = R(c*mw_base_log[i]).round() |

4029 | m[r,r] = R(c*w1).round() |

4030 | |

4031 | #LLL - implemented in sage - operates on rows not on columns |

4032 | m_LLL = m.LLL() |

4033 | m_gram = m_LLL.gram_schmidt()[0] |

4034 | b1_norm = R(m_LLL.row(0).norm()) |

4035 | |

4036 | #compute constant c1 ~ c1_LLL of Corollary 2.3.17 and hence d(L,0)^2 ~ d_L_0 |

4037 | c1_LLL = -1 |

4038 | for i in range(n): |

4039 | tmp = R(b1_norm/(m_gram.row(i).norm())) |

4040 | if tmp > c1_LLL: |

4041 | c1_LLL = tmp |

4042 | |

4043 | if c1_LLL < 0: |

4044 | raise RuntimeError, 'Unexpected intermediate result. Please try another Mordell-Weil base' |

4045 | |

4046 | d_L_0 = R(b1_norm**2 / c1_LLL) |

4047 | |

4048 | #Reducing of upper bound |

4049 | Q = r * H_q**2 |

4050 | T = (1 + (3/2*r*H_q))/2 |

4051 | if d_L_0 < R(T**2+Q): |

4052 | d_L_0 = 10*(T**2*Q) |

4053 | low_bound = R((sqrt(d_L_0 - Q) - T)/c) |

4054 | |

4055 | #new bound according to low_bound and upper bound [c_5 exp((-c_2*H_q^2)/2)] provided by Corollary 8.7.3 |

4056 | H_q_new = R(sqrt(log(low_bound/c5, e)/(-c2/2))) |

4057 | H_q_new = ceil(H_q_new) |

4058 | break_cond = R(H_q_new/H_q) |

4059 | H_q = H_q_new |

4060 | #END LLL-Reduction loop |

4061 | |

4062 | |

4063 | int_points = set([]) |

4064 | ##for a more detailed view how the function works uncomment |

4065 | ##all print - statements below |

4066 | if disc > 0: |

4067 | ##Points in egg have e2>=x>=e1 |

4068 | search_points(ceil(e1), floor(e2)+1) |

4069 | #print 'Points on compact component \n',list(int_points) |

4070 | |

4071 | ##Points in noncompact component with x(P)< 2*max(|e1|,|e2|,|e3|) , espec. x(P)>=e3 |

4072 | search_points(ceil(e3), ceil(2*max(abs(e1),abs(e2),abs(e3)))) |

4073 | #print 'Points on compact component and non-compact component part 1 \n', list(int_points) |

4074 | |

4075 | #print 'starting search of remainig points using bound ',H_q |

4076 | search_remaining_points() |

4077 | #print 'Integral points:\n',list(int_points) |

4078 | #print 'Total amount of points:',len(int_points) |

4079 | |

4080 | return int_points |

4081 | |

4082 | |

4083 | def cremona_curves(conductors): |

4084 | """ |

4085 | Return iterator over all known curves (in database) with conductor |

4086 | in the list of conductors. |

4087 | |

4088 | EXAMPLES: |

4089 | sage: [(E.label(), E.rank()) for E in cremona_curves(srange(35,40))] |

4090 | [('35a1', 0), |

4091 | ('35a2', 0), |

4092 | ('35a3', 0), |

4093 | ('36a1', 0), |

4094 | ('36a2', 0), |

4095 | ('36a3', 0), |

4096 | ('36a4', 0), |

4097 | ('37a1', 1), |

4098 | ('37b1', 0), |

4099 | ('37b2', 0), |

4100 | ('37b3', 0), |

4101 | ('38a1', 0), |

4102 | ('38a2', 0), |

4103 | ('38a3', 0), |

4104 | ('38b1', 0), |

4105 | ('38b2', 0), |

4106 | ('39a1', 0), |

4107 | ('39a2', 0), |

4108 | ('39a3', 0), |

4109 | ('39a4', 0)] |

4110 | """ |

4111 | if isinstance(conductors, (int,long, rings.RingElement)): |

4112 | conductors = [conductors] |

4113 | return sage.databases.cremona.CremonaDatabase().iter(conductors) |

4114 | |

4115 | def cremona_optimal_curves(conductors): |

4116 | """ |

4117 | Return iterator over all known optimal curves (in database) with |

4118 | conductor in the list of conductors. |

4119 | |

4120 | EXAMPLES: |

4121 | sage: [(E.label(), E.rank()) for E in cremona_optimal_curves(srange(35,40))] |

4122 | [('35a1', 0), |

4123 | ('36a1', 0), |

4124 | ('37a1', 1), |

4125 | ('37b1', 0), |

4126 | ('38a1', 0), |

4127 | ('38b1', 0), |

4128 | ('39a1', 0)] |

4129 | """ |

4130 | if isinstance(conductors, (int,long,rings.RingElement)): |

4131 | conductors = [conductors] |

4132 | return sage.databases.cremona.CremonaDatabase().iter_optimal(conductors) |

4133 |