| 1 | """ |
| 2 | Descent on elliptic curves over QQ with a 2-isogeny. |
| 3 | """ |
| 4 | |
| 5 | #***************************************************************************** |
| 6 | # Copyright (C) 2009 Robert L. Miller <rlmillster@gmail.com> |
| 7 | # |
| 8 | # Distributed under the terms of the GNU General Public License (GPL) |
| 9 | # http://www.gnu.org/licenses/ |
| 10 | #***************************************************************************** |
| 11 | |
| 12 | from sage.rings.all import ZZ |
| 13 | from sage.rings.polynomial.polynomial_ring import polygen |
| 14 | cdef object x_ZZ = polygen(ZZ) |
| 15 | from sage.rings.polynomial.real_roots import real_roots |
| 16 | from sage.rings.arith import prime_divisors |
| 17 | from sage.misc.all import walltime, cputime |
| 18 | from sage.libs.pari.gen import pari |
| 19 | |
| 20 | from sage.rings.integer cimport Integer |
| 21 | |
| 22 | include "../../ext/cdefs.pxi" |
| 23 | include "../../libs/flint/fmpz_poly.pxi" |
| 24 | |
| 25 | from sage.libs.flint.zmod_poly cimport *, zmod_poly_t |
| 26 | from sage.libs.flint.long_extras cimport *, factor_t |
| 27 | from sage.libs.ratpoints cimport ratpoints_mpz_exists_only |
| 28 | |
| 29 | cdef int N_RES_CLASSES_BSD = 10 |
| 30 | |
| 31 | cdef unsigned long ui0 = <unsigned long>0 |
| 32 | cdef unsigned long ui1 = <unsigned long>1 |
| 33 | cdef unsigned long ui2 = <unsigned long>2 |
| 34 | cdef unsigned long ui3 = <unsigned long>3 |
| 35 | cdef unsigned long ui4 = <unsigned long>4 |
| 36 | cdef unsigned long ui8 = <unsigned long>8 |
| 37 | |
| 38 | cdef unsigned long valuation(mpz_t a, mpz_t p): |
| 39 | """ |
| 40 | Return the number of times p divides a. |
| 41 | """ |
| 42 | cdef mpz_t aa |
| 43 | cdef unsigned long v |
| 44 | mpz_init(aa) |
| 45 | v = mpz_remove(aa,a,p) |
| 46 | mpz_clear(aa) |
| 47 | return v |
| 48 | |
| 49 | def test_valuation(a, p): |
| 50 | """ |
| 51 | Doctest function for cdef long valuation(mpz_t, mpz_t). |
| 52 | |
| 53 | EXAMPLE:: |
| 54 | |
| 55 | sage: from sage.schemes.elliptic_curves.descent import test_valuation as tv |
| 56 | sage: for i in [1..20]: |
| 57 | ... print '%10s'%factor(i), tv(i,2), tv(i,3), tv(i,5) |
| 58 | 1 0 0 0 |
| 59 | 2 1 0 0 |
| 60 | 3 0 1 0 |
| 61 | 2^2 2 0 0 |
| 62 | 5 0 0 1 |
| 63 | 2 * 3 1 1 0 |
| 64 | 7 0 0 0 |
| 65 | 2^3 3 0 0 |
| 66 | 3^2 0 2 0 |
| 67 | 2 * 5 1 0 1 |
| 68 | 11 0 0 0 |
| 69 | 2^2 * 3 2 1 0 |
| 70 | 13 0 0 0 |
| 71 | 2 * 7 1 0 0 |
| 72 | 3 * 5 0 1 1 |
| 73 | 2^4 4 0 0 |
| 74 | 17 0 0 0 |
| 75 | 2 * 3^2 1 2 0 |
| 76 | 19 0 0 0 |
| 77 | 2^2 * 5 2 0 1 |
| 78 | |
| 79 | """ |
| 80 | cdef Integer A = Integer(a) |
| 81 | cdef Integer P = Integer(p) |
| 82 | return valuation(A.value, P.value) |
| 83 | |
| 84 | cdef int padic_square(mpz_t a, mpz_t p): |
| 85 | """ |
| 86 | Test if a is a p-adic square. |
| 87 | """ |
| 88 | cdef unsigned long v |
| 89 | cdef mpz_t aa |
| 90 | cdef int result |
| 91 | |
| 92 | if mpz_sgn(a) == 0: return 1 |
| 93 | |
| 94 | v = valuation(a,p) |
| 95 | if v&(ui1): return 0 |
| 96 | |
| 97 | mpz_init_set(aa,a) |
| 98 | while v: |
| 99 | v -= 1 |
| 100 | mpz_divexact(aa, aa, p) |
| 101 | if mpz_cmp_ui(p, ui2)==0: |
| 102 | result = ( mpz_fdiv_ui(aa,ui8)==ui1 ) |
| 103 | else: |
| 104 | result = ( mpz_legendre(aa,p)==1 ) |
| 105 | mpz_clear(aa) |
| 106 | return result |
| 107 | |
| 108 | def test_padic_square(a, p): |
| 109 | """ |
| 110 | Doctest function for cdef int padic_square(mpz_t, unsigned long). |
| 111 | |
| 112 | EXAMPLE:: |
| 113 | |
| 114 | sage: from sage.schemes.elliptic_curves.descent_two_isogeny import test_padic_square as ps |
| 115 | sage: for i in [1..300]: |
| 116 | ... for p in prime_range(100): |
| 117 | ... if not Qp(p)(i).is_square()==bool(ps(i,p)): |
| 118 | ... print i, p |
| 119 | |
| 120 | """ |
| 121 | cdef Integer A = Integer(a) |
| 122 | cdef Integer P = Integer(p) |
| 123 | return padic_square(A.value, P.value) |
| 124 | |
| 125 | cdef int lemma6(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e, |
| 126 | mpz_t x, mpz_t p, unsigned long nu): |
| 127 | """ |
| 128 | Implements Lemma 6 of BSD's "Notes on elliptic curves, I" for odd p. |
| 129 | |
| 130 | Returns -1 for insoluble, 0 for undecided, +1 for soluble. |
| 131 | """ |
| 132 | cdef mpz_t g_of_x, g_prime_of_x |
| 133 | cdef unsigned long lambd, mu |
| 134 | cdef int result = -1 |
| 135 | |
| 136 | mpz_init(g_of_x) |
| 137 | mpz_mul(g_of_x, a, x) |
| 138 | mpz_add(g_of_x, g_of_x, b) |
| 139 | mpz_mul(g_of_x, g_of_x, x) |
| 140 | mpz_add(g_of_x, g_of_x, c) |
| 141 | mpz_mul(g_of_x, g_of_x, x) |
| 142 | mpz_add(g_of_x, g_of_x, d) |
| 143 | mpz_mul(g_of_x, g_of_x, x) |
| 144 | mpz_add(g_of_x, g_of_x, e) |
| 145 | |
| 146 | if padic_square(g_of_x, p): |
| 147 | mpz_clear(g_of_x) |
| 148 | return +1 # soluble |
| 149 | |
| 150 | mpz_init_set(g_prime_of_x, x) |
| 151 | mpz_mul(g_prime_of_x, a, x) |
| 152 | mpz_mul_ui(g_prime_of_x, g_prime_of_x, ui4) |
| 153 | mpz_addmul_ui(g_prime_of_x, b, ui3) |
| 154 | mpz_mul(g_prime_of_x, g_prime_of_x, x) |
| 155 | mpz_addmul_ui(g_prime_of_x, c, ui2) |
| 156 | mpz_mul(g_prime_of_x, g_prime_of_x, x) |
| 157 | mpz_add(g_prime_of_x, g_prime_of_x, d) |
| 158 | |
| 159 | lambd = valuation(g_of_x, p) |
| 160 | if mpz_sgn(g_prime_of_x)==0: |
| 161 | if lambd >= 2*nu: result = 0 # undecided |
| 162 | else: |
| 163 | mu = valuation(g_prime_of_x, p) |
| 164 | if lambd > 2*mu: result = +1 # soluble |
| 165 | elif lambd >= 2*nu and mu >= nu: result = 0 # undecided |
| 166 | |
| 167 | mpz_clear(g_prime_of_x) |
| 168 | mpz_clear(g_of_x) |
| 169 | return result |
| 170 | |
| 171 | cdef int lemma7(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e, |
| 172 | mpz_t x, mpz_t p, unsigned long nu): |
| 173 | """ |
| 174 | Implements Lemma 7 of BSD's "Notes on elliptic curves, I" for p=2. |
| 175 | |
| 176 | Returns -1 for insoluble, 0 for undecided, +1 for soluble. |
| 177 | """ |
| 178 | cdef mpz_t g_of_x, g_prime_of_x, g_of_x_odd_part |
| 179 | cdef unsigned long lambd, mu, g_of_x_odd_part_mod_4 |
| 180 | cdef int result = -1 |
| 181 | |
| 182 | mpz_init(g_of_x) |
| 183 | mpz_mul(g_of_x, a, x) |
| 184 | mpz_add(g_of_x, g_of_x, b) |
| 185 | mpz_mul(g_of_x, g_of_x, x) |
| 186 | mpz_add(g_of_x, g_of_x, c) |
| 187 | mpz_mul(g_of_x, g_of_x, x) |
| 188 | mpz_add(g_of_x, g_of_x, d) |
| 189 | mpz_mul(g_of_x, g_of_x, x) |
| 190 | mpz_add(g_of_x, g_of_x, e) |
| 191 | |
| 192 | if padic_square(g_of_x, p): |
| 193 | mpz_clear(g_of_x) |
| 194 | return +1 # soluble |
| 195 | |
| 196 | mpz_init_set(g_prime_of_x, x) |
| 197 | mpz_mul(g_prime_of_x, a, x) |
| 198 | mpz_mul_ui(g_prime_of_x, g_prime_of_x, ui4) |
| 199 | mpz_addmul_ui(g_prime_of_x, b, ui3) |
| 200 | mpz_mul(g_prime_of_x, g_prime_of_x, x) |
| 201 | mpz_addmul_ui(g_prime_of_x, c, ui2) |
| 202 | mpz_mul(g_prime_of_x, g_prime_of_x, x) |
| 203 | mpz_add(g_prime_of_x, g_prime_of_x, d) |
| 204 | |
| 205 | lambd = valuation(g_of_x, p) |
| 206 | mpz_init_set(g_of_x_odd_part, g_of_x) |
| 207 | while mpz_even_p(g_of_x_odd_part): |
| 208 | mpz_divexact_ui(g_of_x_odd_part, g_of_x_odd_part, ui2) |
| 209 | g_of_x_odd_part_mod_4 = mpz_fdiv_ui(g_of_x_odd_part, ui4) |
| 210 | if mpz_sgn(g_prime_of_x)==0: |
| 211 | if lambd >= 2*nu: result = 0 # undecided |
| 212 | elif lambd == 2*nu-2 and g_of_x_odd_part_mod_4==1: |
| 213 | result = 0 # undecided |
| 214 | else: |
| 215 | mu = valuation(g_prime_of_x, p) |
| 216 | if lambd > 2*mu: result = +1 # soluble |
| 217 | elif nu > mu: |
| 218 | if lambd >= mu+nu: result = +1 # soluble |
| 219 | elif lambd+1 == mu+nu and lambd&ui1==0: |
| 220 | result = +1 # soluble |
| 221 | elif lambd+2 == mu+nu and lambd&ui1==0 and g_of_x_odd_part_mod_4==1: |
| 222 | result = +1 # soluble |
| 223 | else: # nu <= mu |
| 224 | if lambd >= 2*nu: result = 0 # undecided |
| 225 | elif lambd+2 == 2*nu and g_of_x_odd_part_mod_4==1: |
| 226 | result = 0 # undecided |
| 227 | |
| 228 | mpz_clear(g_prime_of_x) |
| 229 | mpz_clear(g_of_x) |
| 230 | return result |
| 231 | |
| 232 | cdef int Zp_soluble_BSD(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e, |
| 233 | mpz_t x_k, mpz_t p, unsigned long k): |
| 234 | """ |
| 235 | Uses the approach of BSD's "Notes on elliptic curves, I" to test for |
| 236 | solubility of y^2 == ax^4 + bx^3 + cx^2 + dx + e over Zp, with |
| 237 | x=x_k (mod p^k). |
| 238 | """ |
| 239 | # returns solubility of y^2 = ax^4 + bx^3 + cx^2 + dx + e |
| 240 | # in Zp with x=x_k (mod p^k) |
| 241 | cdef int code |
| 242 | cdef unsigned long t |
| 243 | cdef mpz_t s |
| 244 | |
| 245 | if mpz_cmp_ui(p, ui2)==0: |
| 246 | code = lemma7(a,b,c,d,e,x_k,p,k) |
| 247 | else: |
| 248 | code = lemma6(a,b,c,d,e,x_k,p,k) |
| 249 | if code == 1: |
| 250 | return 1 |
| 251 | if code == -1: |
| 252 | return 0 |
| 253 | |
| 254 | # now code == 0 |
| 255 | t = 0 |
| 256 | mpz_init(s) |
| 257 | while code == 0 and mpz_cmp_ui(p, t) > 0 and t < N_RES_CLASSES_BSD: |
| 258 | mpz_pow_ui(s, p, k) |
| 259 | mpz_mul_ui(s, s, t) |
| 260 | mpz_add(s, s, x_k) |
| 261 | code = Zp_soluble_BSD(a,b,c,d,e,s,p,k+1) |
| 262 | t += 1 |
| 263 | mpz_clear(s) |
| 264 | return code |
| 265 | |
| 266 | cdef bint Zp_soluble_siksek(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e, |
| 267 | mpz_t pp, unsigned long pp_ui, |
| 268 | zmod_poly_factor_t f_factzn, zmod_poly_t f, |
| 269 | fmpz_poly_t f1, fmpz_poly_t linear, |
| 270 | double pp_ui_inv): |
| 271 | """ |
| 272 | Uses the approach of Algorithm 5.3.1 of Siksek's thesis to test for |
| 273 | solubility of y^2 == ax^4 + bx^3 + cx^2 + dx + e over Zp. |
| 274 | """ |
| 275 | cdef unsigned long v_min, v |
| 276 | cdef unsigned long roots[4] |
| 277 | cdef int i, j, has_roots, has_single_roots |
| 278 | cdef bint result |
| 279 | |
| 280 | cdef mpz_t aa, bb, cc, dd, ee |
| 281 | cdef mpz_t aaa, bbb, ccc, ddd, eee |
| 282 | cdef unsigned long qq |
| 283 | cdef unsigned long rr, ss |
| 284 | cdef mpz_t tt |
| 285 | |
| 286 | # Step 0: divide out all common p from the quartic |
| 287 | v_min = valuation(a, pp) |
| 288 | if mpz_cmp_ui(b, ui0) != 0: |
| 289 | v = valuation(b, pp) |
| 290 | if v < v_min: v_min = v |
| 291 | if mpz_cmp_ui(c, ui0) != 0: |
| 292 | v = valuation(c, pp) |
| 293 | if v < v_min: v_min = v |
| 294 | if mpz_cmp_ui(d, ui0) != 0: |
| 295 | v = valuation(d, pp) |
| 296 | if v < v_min: v_min = v |
| 297 | if mpz_cmp_ui(e, ui0) != 0: |
| 298 | v = valuation(e, pp) |
| 299 | if v < v_min: v_min = v |
| 300 | for 0 <= v < v_min: |
| 301 | mpz_divexact(a, a, pp) |
| 302 | mpz_divexact(b, b, pp) |
| 303 | mpz_divexact(c, c, pp) |
| 304 | mpz_divexact(d, d, pp) |
| 305 | mpz_divexact(e, e, pp) |
| 306 | |
| 307 | if not v_min%2: |
| 308 | # Step I in Alg. 5.3.1 of Siksek's thesis |
| 309 | zmod_poly_set_coeff_ui(f, 0, mpz_fdiv_ui(e, pp_ui)) |
| 310 | zmod_poly_set_coeff_ui(f, 1, mpz_fdiv_ui(d, pp_ui)) |
| 311 | zmod_poly_set_coeff_ui(f, 2, mpz_fdiv_ui(c, pp_ui)) |
| 312 | zmod_poly_set_coeff_ui(f, 3, mpz_fdiv_ui(b, pp_ui)) |
| 313 | zmod_poly_set_coeff_ui(f, 4, mpz_fdiv_ui(a, pp_ui)) |
| 314 | |
| 315 | result = 0 |
| 316 | (<zmod_poly_factor_struct *>f_factzn)[0].num_factors = 0 # reset data struct |
| 317 | qq = zmod_poly_factor(f_factzn, f) |
| 318 | for i from 0 <= i < f_factzn.num_factors: |
| 319 | if f_factzn.exponents[i]&1: |
| 320 | result = 1 |
| 321 | break |
| 322 | if result == 0 and z_legendre_precomp(qq, pp_ui, pp_ui_inv) == 1: |
| 323 | result = 1 |
| 324 | if result: |
| 325 | return 1 |
| 326 | |
| 327 | zmod_poly_zero(f) |
| 328 | zmod_poly_set_coeff_ui(f, 0, ui1) |
| 329 | for i from 0 <= i < f_factzn.num_factors: |
| 330 | for j from 0 <= j < (f_factzn.exponents[i]>>1): |
| 331 | zmod_poly_mul(f, f, f_factzn.factors[i]) |
| 332 | |
| 333 | (<zmod_poly_factor_struct *>f_factzn)[0].num_factors = 0 # reset data struct |
| 334 | zmod_poly_factor(f_factzn, f) |
| 335 | has_roots = 0 |
| 336 | j = 0 |
| 337 | for i from 0 <= i < f_factzn.num_factors: |
| 338 | if zmod_poly_degree(f_factzn.factors[i]) == 1 and 0 != zmod_poly_get_coeff_ui(f_factzn.factors[i], 1): |
| 339 | has_roots = 1 |
| 340 | roots[j] = pp_ui - zmod_poly_get_coeff_ui(f_factzn.factors[i], 0) |
| 341 | j += 1 |
| 342 | if not has_roots: |
| 343 | return 0 |
| 344 | |
| 345 | i = zmod_poly_degree(f) |
| 346 | mpz_init(aaa) |
| 347 | mpz_init(bbb) |
| 348 | mpz_init(ccc) |
| 349 | mpz_init(ddd) |
| 350 | mpz_init(eee) |
| 351 | |
| 352 | if i == 0: # g == 1 |
| 353 | mpz_set(aaa, a) |
| 354 | mpz_set(bbb, b) |
| 355 | mpz_set(ccc, c) |
| 356 | mpz_set(ddd, d) |
| 357 | mpz_sub_ui(eee, e, qq) |
| 358 | elif i == 1: # g == x + rr |
| 359 | mpz_set(aaa, a) |
| 360 | mpz_set(bbb, b) |
| 361 | mpz_sub_ui(ccc, c, qq) |
| 362 | rr = zmod_poly_get_coeff_ui(f, 0) |
| 363 | ss = rr*qq |
| 364 | mpz_set(ddd,d) |
| 365 | mpz_sub_ui(ddd, ddd, ss*ui2) |
| 366 | mpz_set(eee,e) |
| 367 | mpz_sub_ui(eee, eee, ss*rr) |
| 368 | elif i == 2: # g == x^2 + rr*x + ss |
| 369 | mpz_sub_ui(aaa, a, qq) |
| 370 | rr = zmod_poly_get_coeff_ui(f, 1) |
| 371 | mpz_init(tt) |
| 372 | mpz_set_ui(tt, rr*qq) |
| 373 | mpz_set(bbb,b) |
| 374 | mpz_submul_ui(bbb, tt, ui2) |
| 375 | mpz_set(ccc,c) |
| 376 | mpz_submul_ui(ccc, tt, rr) |
| 377 | ss = zmod_poly_get_coeff_ui(f, 0) |
| 378 | mpz_set_ui(tt, ss*qq) |
| 379 | mpz_set(eee,e) |
| 380 | mpz_submul_ui(eee, tt, ss) |
| 381 | mpz_mul_ui(tt, tt, ui2) |
| 382 | mpz_sub(ccc, ccc, tt) |
| 383 | mpz_set(ddd,d) |
| 384 | mpz_submul_ui(ddd, tt, rr) |
| 385 | mpz_clear(tt) |
| 386 | mpz_divexact(aaa, aaa, pp) |
| 387 | mpz_divexact(bbb, bbb, pp) |
| 388 | mpz_divexact(ccc, ccc, pp) |
| 389 | mpz_divexact(ddd, ddd, pp) |
| 390 | mpz_divexact(eee, eee, pp) |
| 391 | # now aaa,bbb,ccc,ddd,eee represents h(x) |
| 392 | |
| 393 | result = 0 |
| 394 | mpz_init(tt) |
| 395 | for i from 0 <= i < j: |
| 396 | mpz_mul_ui(tt, aaa, roots[i]) |
| 397 | mpz_add(tt, tt, bbb) |
| 398 | mpz_mul_ui(tt, tt, roots[i]) |
| 399 | mpz_add(tt, tt, ccc) |
| 400 | mpz_mul_ui(tt, tt, roots[i]) |
| 401 | mpz_add(tt, tt, ddd) |
| 402 | mpz_mul_ui(tt, tt, roots[i]) |
| 403 | mpz_add(tt, tt, eee) |
| 404 | # tt == h(r) mod p |
| 405 | mpz_mod(tt, tt, pp) |
| 406 | if mpz_sgn(tt) == 0: |
| 407 | fmpz_poly_zero(f1) |
| 408 | fmpz_poly_zero(linear) |
| 409 | fmpz_poly_set_coeff_mpz(f1, 0, e) |
| 410 | fmpz_poly_set_coeff_mpz(f1, 1, d) |
| 411 | fmpz_poly_set_coeff_mpz(f1, 2, c) |
| 412 | fmpz_poly_set_coeff_mpz(f1, 3, b) |
| 413 | fmpz_poly_set_coeff_mpz(f1, 4, a) |
| 414 | fmpz_poly_set_coeff_ui(linear, 0, roots[i]) |
| 415 | fmpz_poly_set_coeff_mpz(linear, 1, pp) |
| 416 | fmpz_poly_compose(f1, f1, linear) |
| 417 | fmpz_poly_scalar_div_ui(f1, f1, pp_ui) |
| 418 | fmpz_poly_scalar_div_ui(f1, f1, pp_ui) |
| 419 | mpz_init(aa) |
| 420 | mpz_init(bb) |
| 421 | mpz_init(cc) |
| 422 | mpz_init(dd) |
| 423 | mpz_init(ee) |
| 424 | fmpz_poly_get_coeff_mpz(aa, f1, 4) |
| 425 | fmpz_poly_get_coeff_mpz(bb, f1, 3) |
| 426 | fmpz_poly_get_coeff_mpz(cc, f1, 2) |
| 427 | fmpz_poly_get_coeff_mpz(dd, f1, 1) |
| 428 | fmpz_poly_get_coeff_mpz(ee, f1, 0) |
| 429 | result = Zp_soluble_siksek(aa, bb, cc, dd, ee, pp, pp_ui, f_factzn, f, f1, linear, pp_ui_inv) |
| 430 | mpz_clear(aa) |
| 431 | mpz_clear(bb) |
| 432 | mpz_clear(cc) |
| 433 | mpz_clear(dd) |
| 434 | mpz_clear(ee) |
| 435 | if result == 1: |
| 436 | break |
| 437 | mpz_clear(aaa) |
| 438 | mpz_clear(bbb) |
| 439 | mpz_clear(ccc) |
| 440 | mpz_clear(ddd) |
| 441 | mpz_clear(eee) |
| 442 | mpz_clear(tt) |
| 443 | return result |
| 444 | else: |
| 445 | # Step II in Alg. 5.3.1 of Siksek's thesis |
| 446 | zmod_poly_set_coeff_ui(f, 0, mpz_fdiv_ui(e, pp_ui)) |
| 447 | zmod_poly_set_coeff_ui(f, 1, mpz_fdiv_ui(d, pp_ui)) |
| 448 | zmod_poly_set_coeff_ui(f, 2, mpz_fdiv_ui(c, pp_ui)) |
| 449 | zmod_poly_set_coeff_ui(f, 3, mpz_fdiv_ui(b, pp_ui)) |
| 450 | zmod_poly_set_coeff_ui(f, 4, mpz_fdiv_ui(a, pp_ui)) |
| 451 | (<zmod_poly_factor_struct *>f_factzn)[0].num_factors = 0 # reset data struct |
| 452 | zmod_poly_factor(f_factzn, f) |
| 453 | has_roots = 0 |
| 454 | has_single_roots = 0 |
| 455 | j = 0 |
| 456 | for i from 0 <= i < f_factzn.num_factors: |
| 457 | if zmod_poly_degree(f_factzn.factors[i]) == 1 and 0 != zmod_poly_get_coeff_ui(f_factzn.factors[i], 1): |
| 458 | has_roots = 1 |
| 459 | if f_factzn.exponents[i] == 1: |
| 460 | has_single_roots = 1 |
| 461 | break |
| 462 | roots[j] = pp_ui - zmod_poly_get_coeff_ui(f_factzn.factors[i], 0) |
| 463 | j += 1 |
| 464 | |
| 465 | if not has_roots: return 0 |
| 466 | if has_single_roots: return 1 |
| 467 | |
| 468 | result = 0 |
| 469 | if j > 0: |
| 470 | mpz_init(aa) |
| 471 | mpz_init(bb) |
| 472 | mpz_init(cc) |
| 473 | mpz_init(dd) |
| 474 | mpz_init(ee) |
| 475 | for i from 0 <= i < j: |
| 476 | fmpz_poly_zero(f1) |
| 477 | fmpz_poly_zero(linear) |
| 478 | fmpz_poly_set_coeff_mpz(f1, 0, e) |
| 479 | fmpz_poly_set_coeff_mpz(f1, 1, d) |
| 480 | fmpz_poly_set_coeff_mpz(f1, 2, c) |
| 481 | fmpz_poly_set_coeff_mpz(f1, 3, b) |
| 482 | fmpz_poly_set_coeff_mpz(f1, 4, a) |
| 483 | fmpz_poly_set_coeff_ui(linear, 0, roots[i]) |
| 484 | fmpz_poly_set_coeff_mpz(linear, 1, pp) |
| 485 | fmpz_poly_compose(f1, f1, linear) |
| 486 | fmpz_poly_scalar_div_ui(f1, f1, pp_ui) |
| 487 | fmpz_poly_get_coeff_mpz(aa, f1, 4) |
| 488 | fmpz_poly_get_coeff_mpz(bb, f1, 3) |
| 489 | fmpz_poly_get_coeff_mpz(cc, f1, 2) |
| 490 | fmpz_poly_get_coeff_mpz(dd, f1, 1) |
| 491 | fmpz_poly_get_coeff_mpz(ee, f1, 0) |
| 492 | result = Zp_soluble_siksek(aa, bb, cc, dd, ee, pp, pp_ui, f_factzn, f, f1, linear, pp_ui_inv) |
| 493 | if result == 1: |
| 494 | break |
| 495 | if j > 0: |
| 496 | mpz_clear(aa) |
| 497 | mpz_clear(bb) |
| 498 | mpz_clear(cc) |
| 499 | mpz_clear(dd) |
| 500 | mpz_clear(ee) |
| 501 | return result |
| 502 | |
| 503 | cdef bint Qp_soluble_siksek(mpz_t A, mpz_t B, mpz_t C, mpz_t D, mpz_t E, |
| 504 | mpz_t p, unsigned long P, |
| 505 | zmod_poly_factor_t f_factzn, fmpz_poly_t f1, |
| 506 | fmpz_poly_t linear): |
| 507 | """ |
| 508 | Uses Samir Siksek's thesis results to determine whether the quartic is |
| 509 | locally soluble at p. |
| 510 | """ |
| 511 | cdef int result = 0 |
| 512 | cdef mpz_t a,b,c,d,e |
| 513 | cdef zmod_poly_t f |
| 514 | cdef double pp_ui_inv = z_precompute_inverse(P) |
| 515 | zmod_poly_init(f, P) |
| 516 | |
| 517 | mpz_init_set(a,A) |
| 518 | mpz_init_set(b,B) |
| 519 | mpz_init_set(c,C) |
| 520 | mpz_init_set(d,D) |
| 521 | mpz_init_set(e,E) |
| 522 | |
| 523 | if Zp_soluble_siksek(a,b,c,d,e,p,P,f_factzn, f, f1, linear, pp_ui_inv): |
| 524 | result = 1 |
| 525 | else: |
| 526 | mpz_set(a,A) |
| 527 | mpz_set(b,B) |
| 528 | mpz_set(c,C) |
| 529 | mpz_set(d,D) |
| 530 | mpz_set(e,E) |
| 531 | if Zp_soluble_siksek(e,d,c,b,a,p,P,f_factzn, f, f1, linear, pp_ui_inv): |
| 532 | result = 1 |
| 533 | |
| 534 | mpz_clear(a) |
| 535 | mpz_clear(b) |
| 536 | mpz_clear(c) |
| 537 | mpz_clear(d) |
| 538 | mpz_clear(e) |
| 539 | zmod_poly_clear(f) |
| 540 | return result |
| 541 | |
| 542 | cdef bint Qp_soluble_BSD(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e, mpz_t p): |
| 543 | """ |
| 544 | Uses the original test of Birch and Swinnerton-Dyer to test for local |
| 545 | solubility of the quartic at p. |
| 546 | """ |
| 547 | cdef mpz_t zero |
| 548 | cdef int result = 0 |
| 549 | mpz_init_set_ui(zero, ui0) |
| 550 | if Zp_soluble_BSD(a,b,c,d,e,zero,p,0): |
| 551 | result = 1 |
| 552 | elif Zp_soluble_BSD(e,d,c,b,a,zero,p,1): |
| 553 | result = 1 |
| 554 | mpz_clear(zero) |
| 555 | return result |
| 556 | |
| 557 | cdef bint Qp_soluble(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e, mpz_t p): |
| 558 | """ |
| 559 | Try the BSD approach for a few residue classes and if no solution is found, |
| 560 | switch to Siksek to try to prove insolubility. |
| 561 | """ |
| 562 | cdef int bsd_sol, sik_sol |
| 563 | cdef unsigned long pp |
| 564 | cdef fmpz_poly_t f1, linear |
| 565 | cdef zmod_poly_factor_t f_factzn |
| 566 | fmpz_poly_init(f1) |
| 567 | fmpz_poly_init(linear) |
| 568 | zmod_poly_factor_init(f_factzn) |
| 569 | bsd_sol = Qp_soluble_BSD(a, b, c, d, e, p) |
| 570 | if mpz_cmp_ui(p,N_RES_CLASSES_BSD)>0 and not bsd_sol: |
| 571 | pp = mpz_get_ui(p) |
| 572 | sik_sol = Qp_soluble_siksek(a,b,c,d,e,p,pp,f_factzn,f1,linear) |
| 573 | else: |
| 574 | sik_sol = bsd_sol |
| 575 | fmpz_poly_clear(f1) |
| 576 | fmpz_poly_clear(linear) |
| 577 | zmod_poly_factor_clear(f_factzn) |
| 578 | return sik_sol |
| 579 | |
| 580 | def test_qpls(a,b,c,d,e,p): |
| 581 | """ |
| 582 | Testing function for Qp_soluble. |
| 583 | |
| 584 | EXAMPLE: |
| 585 | sage: from sage.schemes.elliptic_curves.descent_two_isogeny import test_qpls as tq |
| 586 | sage: tq(1,2,3,4,5,7) |
| 587 | 1 |
| 588 | |
| 589 | """ |
| 590 | cdef Integer A,B,C,D,E,P |
| 591 | cdef int i, result |
| 592 | cdef mpz_t aa,bb,cc,dd,ee,pp |
| 593 | A=Integer(a); B=Integer(b); C=Integer(c); D=Integer(d); E=Integer(e); P=Integer(p) |
| 594 | mpz_init_set(aa, A.value) |
| 595 | mpz_init_set(bb, B.value) |
| 596 | mpz_init_set(cc, C.value) |
| 597 | mpz_init_set(dd, D.value) |
| 598 | mpz_init_set(ee, E.value) |
| 599 | mpz_init_set(pp, P.value) |
| 600 | result = Qp_soluble(aa, bb, cc, dd, ee, pp) |
| 601 | mpz_clear(aa) |
| 602 | mpz_clear(bb) |
| 603 | mpz_clear(cc) |
| 604 | mpz_clear(dd) |
| 605 | mpz_clear(ee) |
| 606 | mpz_clear(pp) |
| 607 | return result |
| 608 | |
| 609 | cdef int everywhere_locally_soluble(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e) except -1: |
| 610 | """ |
| 611 | Returns whether the quartic has local solutions at all primes p. |
| 612 | """ |
| 613 | cdef Integer A,B,C,D,E,Delta,p |
| 614 | cdef mpz_t mpz_2 |
| 615 | A=Integer(0); B=Integer(0); C=Integer(0); D=Integer(0); E=Integer(0) |
| 616 | mpz_set(A.value, a); mpz_set(B.value, b); mpz_set(C.value, c); mpz_set(D.value, d); mpz_set(E.value, e); |
| 617 | f = (((A*x_ZZ + B)*x_ZZ + C)*x_ZZ + D)*x_ZZ + E |
| 618 | |
| 619 | # RR soluble: |
| 620 | if mpz_sgn(a)!=1: |
| 621 | if len(real_roots(f)) == 0: |
| 622 | return 0 |
| 623 | |
| 624 | # Q2 soluble: |
| 625 | mpz_init_set_ui(mpz_2, ui2) |
| 626 | if not Qp_soluble(a,b,c,d,e,mpz_2): |
| 627 | mpz_clear(mpz_2) |
| 628 | return 0 |
| 629 | mpz_clear(mpz_2) |
| 630 | |
| 631 | # Odd finite primes |
| 632 | Delta = f.discriminant() |
| 633 | for p in prime_divisors(Delta): |
| 634 | if p == 2: continue |
| 635 | if not mpz_fits_ulong_p(p.value): |
| 636 | raise ValueError('Modulus must be word-sized to use FLINT for factoring.') |
| 637 | if not Qp_soluble(a,b,c,d,e,p.value): return 0 |
| 638 | |
| 639 | return 1 |
| 640 | |
| 641 | def test_els(a,b,c,d,e): |
| 642 | """ |
| 643 | Doctest function for cdef int everywhere_locally_soluble(mpz_t, mpz_t, mpz_t, mpz_t, mpz_t). |
| 644 | |
| 645 | EXAMPLE:: |
| 646 | |
| 647 | sage: from sage.schemes.elliptic_curves.descent_two_isogeny import test_els |
| 648 | sage: from sage.libs.ratpoints import ratpoints |
| 649 | sage: for _ in range(1000): |
| 650 | ... a,b,c,d,e = randint(1,1000), randint(1,1000), randint(1,1000), randint(1,1000), randint(1,1000) |
| 651 | ... if len(ratpoints([e,d,c,b,a], 1000)) > 0: |
| 652 | ... try: |
| 653 | ... if not test_els(a,b,c,d,e): |
| 654 | ... print "This never happened", a,b,c,d,e |
| 655 | ... except ValueError: |
| 656 | ... continue |
| 657 | |
| 658 | """ |
| 659 | cdef Integer A,B,C,D,E,Delta |
| 660 | A=Integer(a); B=Integer(b); C=Integer(c); D=Integer(d); E=Integer(e) |
| 661 | return everywhere_locally_soluble(A.value, B.value, C.value, D.value, E.value) |
| 662 | |
| 663 | cdef int count(mpz_t c_mpz, mpz_t d_mpz, mpz_t *p_list, unsigned long p_list_len, |
| 664 | int global_limit_small, int global_limit_large, |
| 665 | int verbosity, bint selmer_only, mpz_t n1, mpz_t n2): |
| 666 | """ |
| 667 | Count the number of els/gls quartic 2-covers of E. |
| 668 | """ |
| 669 | cdef unsigned long n_primes, i |
| 670 | cdef bint found_global_points, els, check_negs, verbose = (verbosity > 4) |
| 671 | cdef Integer a_Int, c_Int, e_Int |
| 672 | cdef mpz_t c_sq_mpz, d_prime_mpz |
| 673 | cdef mpz_t n_divisors, j |
| 674 | cdef mpz_t *coeffs_ratp |
| 675 | |
| 676 | |
| 677 | mpz_init(c_sq_mpz) |
| 678 | mpz_mul(c_sq_mpz, c_mpz, c_mpz) |
| 679 | mpz_init_set(d_prime_mpz, c_sq_mpz) |
| 680 | mpz_submul_ui(d_prime_mpz, d_mpz, ui4) |
| 681 | check_negs = 0 |
| 682 | if mpz_sgn(d_prime_mpz) > 0: |
| 683 | if mpz_sgn(c_mpz) >= 0 or mpz_cmp(c_sq_mpz, d_prime_mpz) <= 0: |
| 684 | check_negs = 1 |
| 685 | mpz_clear(c_sq_mpz) |
| 686 | mpz_clear(d_prime_mpz) |
| 687 | |
| 688 | |
| 689 | # Set up coefficient array, and static variables |
| 690 | cdef mpz_t *coeffs = <mpz_t *> sage_malloc(5 * sizeof(mpz_t)) |
| 691 | for i from 0 <= i <= 4: |
| 692 | mpz_init(coeffs[i]) |
| 693 | mpz_set_ui(coeffs[1], ui0) # |
| 694 | mpz_set(coeffs[2], c_mpz) # These never change |
| 695 | mpz_set_ui(coeffs[3], ui0) # |
| 696 | |
| 697 | if not selmer_only: |
| 698 | # allocate space for ratpoints |
| 699 | coeffs_ratp = <mpz_t *> sage_malloc(5 * sizeof(mpz_t)) |
| 700 | for i from 0 <= i <= 4: |
| 701 | mpz_init(coeffs_ratp[i]) |
| 702 | |
| 703 | # Get prime divisors, and put them in an mpz_t array |
| 704 | # (this block, by setting check_negs, takes care of |
| 705 | # local solubility over RR) |
| 706 | cdef mpz_t *p_div_d_mpz = <mpz_t *> sage_malloc((p_list_len+1) * sizeof(mpz_t)) |
| 707 | n_primes = 0 |
| 708 | for i from 0 <= i < p_list_len: |
| 709 | if mpz_divisible_p(d_mpz, p_list[i]): |
| 710 | mpz_init(p_div_d_mpz[n_primes]) |
| 711 | mpz_set(p_div_d_mpz[n_primes], p_list[i]) |
| 712 | n_primes += 1 |
| 713 | if check_negs: |
| 714 | mpz_init(p_div_d_mpz[n_primes]) |
| 715 | mpz_set_si(p_div_d_mpz[n_primes], -1) |
| 716 | n_primes += 1 |
| 717 | mpz_init_set_ui(n_divisors, ui1) |
| 718 | mpz_mul_2exp(n_divisors, n_divisors, n_primes) |
| 719 | # if verbosity > 3: |
| 720 | # print '\nDivisors of d which may lead to RR-soluble quartics:', p_div_d |
| 721 | |
| 722 | mpz_init_set_ui(j, ui0) |
| 723 | if not selmer_only: |
| 724 | mpz_set_ui(n1, ui0) |
| 725 | mpz_set_ui(n2, ui0) |
| 726 | while mpz_cmp(j, n_divisors) < 0: |
| 727 | mpz_set_ui(coeffs[4], ui1) |
| 728 | for i from 0 <= i < n_primes: |
| 729 | if mpz_tstbit(j, i): |
| 730 | mpz_mul(coeffs[4], coeffs[4], p_div_d_mpz[i]) |
| 731 | if verbosity > 3: |
| 732 | a_Int = Integer(0); mpz_set(a_Int.value, coeffs[4]) |
| 733 | print '\nSquarefree divisor:', a_Int |
| 734 | mpz_divexact(coeffs[0], d_mpz, coeffs[4]) |
| 735 | found_global_points = 0 |
| 736 | if not selmer_only: |
| 737 | if verbose: |
| 738 | print "\nCalling ratpoints for small point search" |
| 739 | for i from 0 <= i <= 4: |
| 740 | mpz_set(coeffs_ratp[i], coeffs[i]) |
| 741 | if ratpoints_mpz_exists_only(coeffs_ratp, global_limit_small, 4, verbose): |
| 742 | if verbosity > 2: |
| 743 | a_Int = Integer(0); mpz_set(a_Int.value, coeffs[4]) |
| 744 | c_Int = Integer(0); mpz_set(c_Int.value, coeffs[2]) |
| 745 | e_Int = Integer(0); mpz_set(e_Int.value, coeffs[0]) |
| 746 | print 'Found small global point, quartic (%d,%d,%d,%d,%d)'%(a_Int,0,c_Int,0,e_Int) |
| 747 | found_global_points = 1 |
| 748 | mpz_add_ui(n1, n1, ui1) |
| 749 | mpz_add_ui(n2, n2, ui1) |
| 750 | if verbose: |
| 751 | print "\nDone calling ratpoints for small point search" |
| 752 | if not found_global_points: |
| 753 | # Test whether the quartic is everywhere locally soluble: |
| 754 | els = 1 |
| 755 | for i from 0 <= i < p_list_len: |
| 756 | if not Qp_soluble(coeffs[4], coeffs[3], coeffs[2], coeffs[1], coeffs[0], p_list[i]): |
| 757 | els = 0 |
| 758 | break |
| 759 | if els: |
| 760 | if verbosity > 2: |
| 761 | a_Int = Integer(0); mpz_set(a_Int.value, coeffs[4]) |
| 762 | c_Int = Integer(0); mpz_set(c_Int.value, coeffs[2]) |
| 763 | e_Int = Integer(0); mpz_set(e_Int.value, coeffs[0]) |
| 764 | print 'ELS without small global points, quartic (%d,%d,%d,%d,%d)'%(a_Int,0,c_Int,0,e_Int) |
| 765 | mpz_add_ui(n2, n2, ui1) |
| 766 | if not selmer_only: |
| 767 | if verbose: |
| 768 | print "\nCalling ratpoints for large point search" |
| 769 | for i from 0 <= i <= 4: |
| 770 | mpz_set(coeffs_ratp[i], coeffs[i]) |
| 771 | if ratpoints_mpz_exists_only(coeffs_ratp, global_limit_large, 4, verbose): |
| 772 | if verbosity > 2: |
| 773 | print ' -- Found large global point.' |
| 774 | mpz_add_ui(n1, n1, ui1) |
| 775 | if verbose: |
| 776 | print "\nDone calling ratpoints for large point search" |
| 777 | mpz_add_ui(j, j, ui1) |
| 778 | if not selmer_only: |
| 779 | for i from 0 <= i <= 4: |
| 780 | mpz_clear(coeffs_ratp[i]) |
| 781 | sage_free(coeffs_ratp) |
| 782 | mpz_clear(j) |
| 783 | for i from 0 <= i < n_primes: |
| 784 | mpz_clear(p_div_d_mpz[i]) |
| 785 | sage_free(p_div_d_mpz) |
| 786 | mpz_clear(n_divisors) |
| 787 | for i from 0 <= i <= 4: |
| 788 | mpz_clear(coeffs[i]) |
| 789 | sage_free(coeffs) |
| 790 | return 0 |
| 791 | |
| 792 | def two_descent_by_two_isogeny(E, |
| 793 | int global_limit_small = 10, |
| 794 | int global_limit_large = 10000, |
| 795 | int verbosity = 0, |
| 796 | bint selmer_only = 0, bint proof = 1): |
| 797 | """ |
| 798 | Given an elliptic curve E with a two-isogeny phi : E --> E' and dual isogeny |
| 799 | phi', runs a two-isogeny descent on E, returning n1, n2, n1' and n2'. Here |
| 800 | n1 is the number of quartic covers found with a rational point, and n2 is |
| 801 | the number which are ELS. |
| 802 | |
| 803 | EXAMPLES:: |
| 804 | |
| 805 | sage: from sage.schemes.elliptic_curves.descent_two_isogeny import two_descent_by_two_isogeny |
| 806 | sage: E = EllipticCurve('14a') |
| 807 | sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny(E) |
| 808 | sage: log(n1,2) + log(n1_prime,2) - 2 # the rank |
| 809 | 0 |
| 810 | sage: E = EllipticCurve('65a') |
| 811 | sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny(E) |
| 812 | sage: log(n1,2) + log(n1_prime,2) - 2 # the rank |
| 813 | 1 |
| 814 | sage: E = EllipticCurve('1088j1') |
| 815 | sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny(E) |
| 816 | sage: log(n1,2) + log(n1_prime,2) - 2 # the rank |
| 817 | 2 |
| 818 | sage: E = EllipticCurve('59450i') |
| 819 | sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny(E) |
| 820 | sage: log(n1,2) + log(n1_prime,2) - 2 # the rank |
| 821 | 3 |
| 822 | |
| 823 | :: |
| 824 | |
| 825 | sage: E = EllipticCurve('14a') |
| 826 | sage: two_descent_by_two_isogeny(E, verbosity=1) |
| 827 | 2-isogeny |
| 828 | Results: |
| 829 | 2 <= #E(Q)/phi'(E'(Q)) <= 2 |
| 830 | 2 <= #E'(Q)/phi(E(Q)) <= 2 |
| 831 | #Sel^(phi')(E'/Q) = 2 |
| 832 | #Sel^(phi)(E/Q) = 2 |
| 833 | 1 <= #Sha(E'/Q)[phi'] <= 1 |
| 834 | 1 <= #Sha(E/Q)[phi] <= 1 |
| 835 | 1 <= #Sha(E/Q)[2], #Sha(E'/Q)[2] <= 1 |
| 836 | 0 <= rank of E(Q) = rank of E'(Q) <= 0 |
| 837 | (2, 2, 2, 2) |
| 838 | |
| 839 | |
| 840 | """ |
| 841 | cdef Integer a1, a2, a3, a4, a6, s2, s4, s6 |
| 842 | cdef Integer c, d, x0 |
| 843 | cdef list x_list |
| 844 | assert E.torsion_order()%2==0, 'Need rational two-torsion for isogeny descent.' |
| 845 | if verbosity > 0: |
| 846 | print '\n2-isogeny' |
| 847 | if verbosity > 1: |
| 848 | print '\nchanging coordinates' |
| 849 | a1 = Integer(E.a1()) |
| 850 | a2 = Integer(E.a2()) |
| 851 | a3 = Integer(E.a3()) |
| 852 | a4 = Integer(E.a4()) |
| 853 | a6 = Integer(E.a6()) |
| 854 | if a1==0 and a3==0: |
| 855 | s2=a2; s4=a4; s6=a6 |
| 856 | else: |
| 857 | s2=a1*a1+4*a2; s4=8*(a1*a3+2*a4); s6=16*(a3*a3+4*a6) |
| 858 | f = ((x_ZZ + s2)*x_ZZ + s4)*x_ZZ + s6 |
| 859 | x_list = f.roots() # over ZZ -- use FLINT directly? |
| 860 | x0 = x_list[0][0] |
| 861 | c = 3*x0+s2; d = (c+s2)*x0+s4 |
| 862 | return two_descent_by_two_isogeny_work(c, d, |
| 863 | global_limit_small, global_limit_large, verbosity, selmer_only, proof) |
| 864 | |
| 865 | def two_descent_by_two_isogeny_work(Integer c, Integer d, |
| 866 | int global_limit_small = 10, int global_limit_large = 10000, |
| 867 | int verbosity = 0, bint selmer_only = 0, bint proof = 1): |
| 868 | """ |
| 869 | Do all the work in doing a two-isogeny descent. |
| 870 | |
| 871 | EXAMPLES:: |
| 872 | |
| 873 | sage: from sage.schemes.elliptic_curves.descent_two_isogeny import two_descent_by_two_isogeny_work |
| 874 | sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny_work(13,128) |
| 875 | sage: log(n1,2) + log(n1_prime,2) - 2 # the rank |
| 876 | 0 |
| 877 | sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny_work(1,-16) |
| 878 | sage: log(n1,2) + log(n1_prime,2) - 2 # the rank |
| 879 | 1 |
| 880 | sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny_work(10,8) |
| 881 | sage: log(n1,2) + log(n1_prime,2) - 2 # the rank |
| 882 | 2 |
| 883 | sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny_work(85,320) |
| 884 | sage: log(n1,2) + log(n1_prime,2) - 2 # the rank |
| 885 | 3 |
| 886 | |
| 887 | """ |
| 888 | cdef mpz_t c_mpz, d_mpz, c_prime_mpz, d_prime_mpz |
| 889 | cdef mpz_t *p_list_mpz |
| 890 | cdef unsigned long i, j, p, p_list_len |
| 891 | cdef Integer P, n1, n2, n1_prime, n2_prime, c_prime, d_prime |
| 892 | cdef object PO |
| 893 | cdef bint found, too_big, d_neg, d_prime_neg |
| 894 | cdef factor_t fact |
| 895 | cdef list primes |
| 896 | mpz_init_set(c_mpz, c.value) # |
| 897 | mpz_init_set(d_mpz, d.value) # |
| 898 | mpz_init(c_prime_mpz) # |
| 899 | mpz_init(d_prime_mpz) # |
| 900 | mpz_mul_si(c_prime_mpz, c_mpz, -2) |
| 901 | mpz_mul(d_prime_mpz, c_mpz, c_mpz) |
| 902 | mpz_submul_ui(d_prime_mpz, d_mpz, ui4) |
| 903 | |
| 904 | d_neg = 0 |
| 905 | d_prime_neg = 0 |
| 906 | if mpz_sgn(d_mpz) < 0: |
| 907 | d_neg = 1 |
| 908 | mpz_neg(d_mpz, d_mpz) |
| 909 | if mpz_sgn(d_prime_mpz) < 0: |
| 910 | d_prime_neg = 1 |
| 911 | mpz_neg(d_prime_mpz, d_prime_mpz) |
| 912 | if mpz_fits_ulong_p(d_mpz) and mpz_fits_ulong_p(d_prime_mpz): |
| 913 | # Factor very quickly using FLINT. |
| 914 | p_list_mpz = <mpz_t *> sage_malloc(20 * sizeof(mpz_t)) |
| 915 | mpz_init_set_ui(p_list_mpz[0], ui2) |
| 916 | p_list_len = 1 |
| 917 | z_factor(&fact, mpz_get_ui(d_mpz), proof) |
| 918 | for i from 0 <= i < fact.num: |
| 919 | p = fact.p[i] |
| 920 | if p != ui2: |
| 921 | mpz_init_set_ui(p_list_mpz[p_list_len], p) |
| 922 | p_list_len += 1 |
| 923 | z_factor(&fact, mpz_get_ui(d_prime_mpz), proof) |
| 924 | for i from 0 <= i < fact.num: |
| 925 | p = fact.p[i] |
| 926 | found = 0 |
| 927 | for j from 0 <= j < p_list_len: |
| 928 | if mpz_cmp_ui(p_list_mpz[j], p)==0: |
| 929 | found = 1 |
| 930 | break |
| 931 | if not found: |
| 932 | mpz_init_set_ui(p_list_mpz[p_list_len], p) |
| 933 | p_list_len += 1 |
| 934 | else: |
| 935 | # Factor more slowly using Pari via Python. |
| 936 | d = Integer(0) |
| 937 | mpz_set(d.value, d_mpz) |
| 938 | primes = list(pari(d).factor()[0]) |
| 939 | d_prime = Integer(0) |
| 940 | mpz_set(d_prime.value, d_prime_mpz) |
| 941 | for PO in pari(d_prime).factor()[0]: |
| 942 | if PO not in primes: |
| 943 | primes.append(PO) |
| 944 | P = Integer(2) |
| 945 | if P not in primes: primes.append(P) |
| 946 | p_list_len = len(primes) |
| 947 | p_list_mpz = <mpz_t *> sage_malloc(p_list_len * sizeof(mpz_t)) |
| 948 | for i from 0 <= i < p_list_len: |
| 949 | P = Integer(primes[i]) |
| 950 | mpz_init_set(p_list_mpz[i], P.value) |
| 951 | if d_neg: |
| 952 | mpz_neg(d_mpz, d_mpz) |
| 953 | if d_prime_neg: |
| 954 | mpz_neg(d_prime_mpz, d_prime_mpz) |
| 955 | |
| 956 | if verbosity > 1: |
| 957 | c_prime = -2*c |
| 958 | d_prime = c*c-4*d |
| 959 | print '\nnew curve is y^2 == x( x^2 + (%d)x + (%d) )'%(int(c),int(d)) |
| 960 | print 'new isogenous curve is' + \ |
| 961 | ' y^2 == x( x^2 + (%d)x + (%d) )'%(int(c_prime),int(d_prime)) |
| 962 | |
| 963 | n1 = Integer(0); n2 = Integer(0) |
| 964 | n1_prime = Integer(0); n2_prime = Integer(0) |
| 965 | count(c.value, d.value, p_list_mpz, p_list_len, |
| 966 | global_limit_small, global_limit_large, verbosity, selmer_only, |
| 967 | n1.value, n2.value) |
| 968 | count(c_prime_mpz, d_prime_mpz, p_list_mpz, p_list_len, |
| 969 | global_limit_small, global_limit_large, verbosity, selmer_only, |
| 970 | n1_prime.value, n2_prime.value) |
| 971 | |
| 972 | for i from 0 <= i < p_list_len: |
| 973 | mpz_clear(p_list_mpz[i]) |
| 974 | sage_free(p_list_mpz) |
| 975 | |
| 976 | if verbosity > 0: |
| 977 | print "\nResults:" |
| 978 | print n1, "<= #E(Q)/phi'(E'(Q)) <=", n2 |
| 979 | print n1_prime, "<= #E'(Q)/phi(E(Q)) <=", n2_prime |
| 980 | print "#Sel^(phi')(E'/Q) =", n2 |
| 981 | print "#Sel^(phi)(E/Q) =", n2_prime |
| 982 | print "1 <= #Sha(E'/Q)[phi'] <=", n2/n1 |
| 983 | print "1 <= #Sha(E/Q)[phi] <=", n2_prime/n1_prime |
| 984 | print "1 <= #Sha(E/Q)[2], #Sha(E'/Q)[2] <=", (n2_prime/n1_prime)*(n2/n1) |
| 985 | a = Integer(n1*n1_prime).log(Integer(2)) |
| 986 | e = Integer(n2*n2_prime).log(Integer(2)) |
| 987 | print a - 2, "<= rank of E(Q) = rank of E'(Q) <=", e - 2 |
| 988 | |
| 989 | return n1, n2, n1_prime, n2_prime |
| 990 | |
| 991 | |