descent_two_isogeny.pyx

"""
Descent on elliptic curves over QQ with a 2-isogeny.
"""

#*****************************************************************************
#        Copyright (C) 2009 Robert L. Miller <[email protected]>
#
# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
#                         http://www.gnu.org/licenses/
#*****************************************************************************

from sage.rings.all import ZZ
from sage.rings.polynomial.polynomial_ring import polygen
cdef object x_ZZ = polygen(ZZ)
from sage.rings.polynomial.real_roots import real_roots
from sage.rings.arith import prime_divisors
from sage.misc.all import walltime, cputime
from sage.all import ntl

from sage.rings.integer cimport Integer

include "sage/ext/cdefs.pxi"
include "sage/ext/interrupt.pxi"
include "sage/libs/flint/fmpz_poly.pxi"

from sage.libs.flint.nmod_poly cimport *, nmod_poly_t
from sage.libs.flint.ulong_extras cimport *, n_factor_t
from sage.libs.ratpoints cimport ratpoints_mpz_exists_only

cdef int N_RES_CLASSES_BSD = 10

cdef unsigned long ui0 = <unsigned long>0
cdef unsigned long ui1 = <unsigned long>1
cdef unsigned long ui2 = <unsigned long>2
cdef unsigned long ui3 = <unsigned long>3
cdef unsigned long ui4 = <unsigned long>4
cdef unsigned long ui8 = <unsigned long>8

cdef unsigned long valuation(mpz_t a, mpz_t p):
    """
    Return the number of times p divides a.
    """
    cdef mpz_t aa
    cdef unsigned long v
    mpz_init(aa)
    v = mpz_remove(aa,a,p)
    mpz_clear(aa)
    return v

def test_valuation(a, p):
    """
    Doctest function for cdef long valuation(mpz_t, mpz_t).

    EXAMPLE::

        sage: from sage.schemes.elliptic_curves.descent_two_isogeny import test_valuation as tv
        sage: for i in [1..20]:
        ...    print '%10s'%factor(i), tv(i,2), tv(i,3), tv(i,5)
                 1 0 0 0
                 2 1 0 0
                 3 0 1 0
               2^2 2 0 0
                 5 0 0 1
             2 * 3 1 1 0
                 7 0 0 0
               2^3 3 0 0
               3^2 0 2 0
             2 * 5 1 0 1
                11 0 0 0
           2^2 * 3 2 1 0
                13 0 0 0
             2 * 7 1 0 0
             3 * 5 0 1 1
               2^4 4 0 0
                17 0 0 0
           2 * 3^2 1 2 0
                19 0 0 0
           2^2 * 5 2 0 1

    """
    cdef Integer A = Integer(a)
    cdef Integer P = Integer(p)
    return valuation(A.value, P.value)

cdef int padic_square(mpz_t a, mpz_t p):
    """
    Test if a is a p-adic square.
    """
    cdef unsigned long v
    cdef mpz_t aa
    cdef int result

    if mpz_sgn(a) == 0: return 1

    v = valuation(a,p)
    if v&(ui1): return 0

    mpz_init_set(aa,a)
    while v:
        v -= 1
        mpz_divexact(aa, aa, p)
    if mpz_cmp_ui(p, ui2)==0:
        result = ( mpz_fdiv_ui(aa,ui8)==ui1 )
    else:
        result = ( mpz_legendre(aa,p)==1 )
    mpz_clear(aa)
    return result

def test_padic_square(a, p):
    """
    Doctest function for cdef int padic_square(mpz_t, unsigned long).

    EXAMPLE::

        sage: from sage.schemes.elliptic_curves.descent_two_isogeny import test_padic_square as ps
        sage: for i in [1..300]:
        ...    for p in prime_range(100):
        ...        if not Qp(p)(i).is_square()==bool(ps(i,p)):
        ...            print i, p

    """
    cdef Integer A = Integer(a)
    cdef Integer P = Integer(p)
    return padic_square(A.value, P.value)

cdef int lemma6(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e,
                mpz_t x, mpz_t p, unsigned long nu):
    """
    Implements Lemma 6 of BSD's "Notes on elliptic curves, I" for odd p.

    Returns -1 for insoluble, 0 for undecided, +1 for soluble.
    """
    cdef mpz_t g_of_x, g_prime_of_x
    cdef unsigned long lambd, mu
    cdef int result = -1

    mpz_init(g_of_x)
    mpz_mul(g_of_x, a, x)
    mpz_add(g_of_x, g_of_x, b)
    mpz_mul(g_of_x, g_of_x, x)
    mpz_add(g_of_x, g_of_x, c)
    mpz_mul(g_of_x, g_of_x, x)
    mpz_add(g_of_x, g_of_x, d)
    mpz_mul(g_of_x, g_of_x, x)
    mpz_add(g_of_x, g_of_x, e)

    if padic_square(g_of_x, p):
        mpz_clear(g_of_x)
        return +1 # soluble

    mpz_init_set(g_prime_of_x, x)
    mpz_mul(g_prime_of_x, a, x)
    mpz_mul_ui(g_prime_of_x, g_prime_of_x, ui4)
    mpz_addmul_ui(g_prime_of_x, b, ui3)
    mpz_mul(g_prime_of_x, g_prime_of_x, x)
    mpz_addmul_ui(g_prime_of_x, c, ui2)
    mpz_mul(g_prime_of_x, g_prime_of_x, x)
    mpz_add(g_prime_of_x, g_prime_of_x, d)

    lambd = valuation(g_of_x, p)
    if mpz_sgn(g_prime_of_x)==0:
        if lambd >= 2*nu: result = 0 # undecided
    else:
        mu = valuation(g_prime_of_x, p)
        if lambd > 2*mu: result = +1 # soluble
        elif lambd >= 2*nu and mu >= nu: result = 0 # undecided

    mpz_clear(g_prime_of_x)
    mpz_clear(g_of_x)
    return result

cdef int lemma7(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e,
                mpz_t x, mpz_t p, unsigned long nu):
    """
    Implements Lemma 7 of BSD's "Notes on elliptic curves, I" for p=2.

    Returns -1 for insoluble, 0 for undecided, +1 for soluble.
    """
    cdef mpz_t g_of_x, g_prime_of_x, g_of_x_odd_part
    cdef unsigned long lambd, mu, g_of_x_odd_part_mod_4
    cdef int result = -1

    mpz_init(g_of_x)
    mpz_mul(g_of_x, a, x)
    mpz_add(g_of_x, g_of_x, b)
    mpz_mul(g_of_x, g_of_x, x)
    mpz_add(g_of_x, g_of_x, c)
    mpz_mul(g_of_x, g_of_x, x)
    mpz_add(g_of_x, g_of_x, d)
    mpz_mul(g_of_x, g_of_x, x)
    mpz_add(g_of_x, g_of_x, e)

    if padic_square(g_of_x, p):
        mpz_clear(g_of_x)
        return +1 # soluble

    mpz_init_set(g_prime_of_x, x)
    mpz_mul(g_prime_of_x, a, x)
    mpz_mul_ui(g_prime_of_x, g_prime_of_x, ui4)
    mpz_addmul_ui(g_prime_of_x, b, ui3)
    mpz_mul(g_prime_of_x, g_prime_of_x, x)
    mpz_addmul_ui(g_prime_of_x, c, ui2)
    mpz_mul(g_prime_of_x, g_prime_of_x, x)
    mpz_add(g_prime_of_x, g_prime_of_x, d)

    lambd = valuation(g_of_x, p)
    mpz_init_set(g_of_x_odd_part, g_of_x)
    while mpz_even_p(g_of_x_odd_part):
        mpz_divexact_ui(g_of_x_odd_part, g_of_x_odd_part, ui2)
    g_of_x_odd_part_mod_4 = mpz_fdiv_ui(g_of_x_odd_part, ui4)
    if mpz_sgn(g_prime_of_x)==0:
        if lambd >= 2*nu: result = 0 # undecided
        elif lambd == 2*nu-2 and g_of_x_odd_part_mod_4==1:
            result = 0 # undecided
    else:
        mu = valuation(g_prime_of_x, p)
        if lambd > 2*mu: result = +1 # soluble
        elif nu > mu:
            if lambd >= mu+nu: result = +1 # soluble
            elif lambd+1 == mu+nu and lambd&ui1==0:
                result = +1 # soluble
            elif lambd+2 == mu+nu and lambd&ui1==0 and g_of_x_odd_part_mod_4==1:
                result = +1 # soluble
        else: # nu <= mu
            if lambd >= 2*nu: result = 0 # undecided
            elif lambd+2 == 2*nu and g_of_x_odd_part_mod_4==1:
                result = 0 # undecided

    mpz_clear(g_prime_of_x)
    mpz_clear(g_of_x)
    return result

cdef int Zp_soluble_BSD(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e,
                        mpz_t x_k, mpz_t p, unsigned long k):
    """
    Uses the approach of BSD's "Notes on elliptic curves, I" to test for
    solubility of y^2 == ax^4 + bx^3 + cx^2 + dx + e over Zp, with
    x=x_k (mod p^k).
    """
    # returns solubility of y^2 = ax^4 + bx^3 + cx^2 + dx + e
    # in Zp with x=x_k (mod p^k)
    cdef int code
    cdef unsigned long t
    cdef mpz_t s

    if mpz_cmp_ui(p, ui2)==0:
        code = lemma7(a,b,c,d,e,x_k,p,k)
    else:
        code = lemma6(a,b,c,d,e,x_k,p,k)
    if code == 1:
        return 1
    if code == -1:
        return 0

    # now code == 0
    t = 0
    mpz_init(s)
    while code == 0 and mpz_cmp_ui(p, t) > 0 and t < N_RES_CLASSES_BSD:
        mpz_pow_ui(s, p, k)
        mpz_mul_ui(s, s, t)
        mpz_add(s, s, x_k)
        code = Zp_soluble_BSD(a,b,c,d,e,s,p,k+1)
        t += 1
    mpz_clear(s)
    return code

cdef bint Zp_soluble_siksek(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e,
                            mpz_t pp, unsigned long pp_ui,
                            nmod_poly_factor_t f_factzn, nmod_poly_t f,
                            fmpz_poly_t f1, fmpz_poly_t linear):
    """
    Uses the approach of Algorithm 5.3.1 of Siksek's thesis to test for
    solubility of y^2 == ax^4 + bx^3 + cx^2 + dx + e over Zp.
    """
    cdef unsigned long v_min, v
    cdef unsigned long roots[4]
    cdef int i, j, has_roots, has_single_roots
    cdef bint result

    cdef mpz_t aa, bb, cc, dd, ee
    cdef mpz_t aaa, bbb, ccc, ddd, eee
    cdef unsigned long qq
    cdef unsigned long rr, ss
    cdef mpz_t tt

    # Step 0: divide out all common p from the quartic
    v_min = valuation(a, pp)
    if mpz_cmp_ui(b, ui0) != 0:
        v = valuation(b, pp)
        if v < v_min: v_min = v
    if mpz_cmp_ui(c, ui0) != 0:
        v = valuation(c, pp)
        if v < v_min: v_min = v
    if mpz_cmp_ui(d, ui0) != 0:
        v = valuation(d, pp)
        if v < v_min: v_min = v
    if mpz_cmp_ui(e, ui0) != 0:
        v = valuation(e, pp)
        if v < v_min: v_min = v
    for 0 <= v < v_min:
        mpz_divexact(a, a, pp)
        mpz_divexact(b, b, pp)
        mpz_divexact(c, c, pp)
        mpz_divexact(d, d, pp)
        mpz_divexact(e, e, pp)

    if not v_min%2:
        # Step I in Alg. 5.3.1 of Siksek's thesis
        nmod_poly_set_coeff_ui(f, 0, mpz_fdiv_ui(e, pp_ui))
        nmod_poly_set_coeff_ui(f, 1, mpz_fdiv_ui(d, pp_ui))
        nmod_poly_set_coeff_ui(f, 2, mpz_fdiv_ui(c, pp_ui))
        nmod_poly_set_coeff_ui(f, 3, mpz_fdiv_ui(b, pp_ui))
        nmod_poly_set_coeff_ui(f, 4, mpz_fdiv_ui(a, pp_ui))

        result = 0
        (<nmod_poly_factor_struct *>f_factzn)[0].num = 0 # reset data struct
        qq = nmod_poly_factor(f_factzn, f)
        for i from 0 <= i < f_factzn.num:
            if f_factzn.exp[i]&1:
                result = 1
                break
        if result == 0 and n_jacobi(qq, pp_ui) == 1:
            result = 1
        if result:
            return 1

        nmod_poly_zero(f)
        nmod_poly_set_coeff_ui(f, 0, ui1)
        for i from 0 <= i < f_factzn.num:
            for j from 0 <= j < (f_factzn.exp[i]>>1):
                nmod_poly_mul(f, f, &f_factzn.p[i])

        (<nmod_poly_factor_struct *>f_factzn)[0].num = 0 # reset data struct
        nmod_poly_factor(f_factzn, f)
        has_roots = 0
        j = 0
        for i from 0 <= i < f_factzn.num:
            if nmod_poly_degree(&f_factzn.p[i]) == 1 and 0 != nmod_poly_get_coeff_ui(&f_factzn.p[i], 1):
                has_roots = 1
                roots[j] = pp_ui - nmod_poly_get_coeff_ui(&f_factzn.p[i], 0)
                j += 1
        if not has_roots:
            return 0

        i = nmod_poly_degree(f)
        mpz_init(aaa)
        mpz_init(bbb)
        mpz_init(ccc)
        mpz_init(ddd)
        mpz_init(eee)

        if i == 0: # g == 1
            mpz_set(aaa, a)
            mpz_set(bbb, b)
            mpz_set(ccc, c)
            mpz_set(ddd, d)
            mpz_sub_ui(eee, e, qq)
        elif i == 1: # g == x + rr
            mpz_set(aaa, a)
            mpz_set(bbb, b)
            mpz_sub_ui(ccc, c, qq)
            rr = nmod_poly_get_coeff_ui(f, 0)
            ss = rr*qq
            mpz_set(ddd,d)
            mpz_sub_ui(ddd, ddd, ss*ui2)
            mpz_set(eee,e)
            mpz_sub_ui(eee, eee, ss*rr)
        elif i == 2: # g == x^2 + rr*x + ss
            mpz_sub_ui(aaa, a, qq)
            rr = nmod_poly_get_coeff_ui(f, 1)
            mpz_init(tt)
            mpz_set_ui(tt, rr*qq)
            mpz_set(bbb,b)
            mpz_submul_ui(bbb, tt, ui2)
            mpz_set(ccc,c)
            mpz_submul_ui(ccc, tt, rr)
            ss = nmod_poly_get_coeff_ui(f, 0)
            mpz_set_ui(tt, ss*qq)
            mpz_set(eee,e)
            mpz_submul_ui(eee, tt, ss)
            mpz_mul_ui(tt, tt, ui2)
            mpz_sub(ccc, ccc, tt)
            mpz_set(ddd,d)
            mpz_submul_ui(ddd, tt, rr)
            mpz_clear(tt)
        mpz_divexact(aaa, aaa, pp)
        mpz_divexact(bbb, bbb, pp)
        mpz_divexact(ccc, ccc, pp)
        mpz_divexact(ddd, ddd, pp)
        mpz_divexact(eee, eee, pp)
        # now aaa,bbb,ccc,ddd,eee represents h(x)

        result = 0
        mpz_init(tt)
        for i from 0 <= i < j:
            mpz_mul_ui(tt, aaa, roots[i])
            mpz_add(tt, tt, bbb)
            mpz_mul_ui(tt, tt, roots[i])
            mpz_add(tt, tt, ccc)
            mpz_mul_ui(tt, tt, roots[i])
            mpz_add(tt, tt, ddd)
            mpz_mul_ui(tt, tt, roots[i])
            mpz_add(tt, tt, eee)
            # tt == h(r) mod p
            mpz_mod(tt, tt, pp)
            if mpz_sgn(tt) == 0:
                fmpz_poly_zero(f1)
                fmpz_poly_zero(linear)
                fmpz_poly_set_coeff_mpz(f1, 0, e)
                fmpz_poly_set_coeff_mpz(f1, 1, d)
                fmpz_poly_set_coeff_mpz(f1, 2, c)
                fmpz_poly_set_coeff_mpz(f1, 3, b)
                fmpz_poly_set_coeff_mpz(f1, 4, a)
                fmpz_poly_set_coeff_ui(linear, 0, roots[i])
                fmpz_poly_set_coeff_mpz(linear, 1, pp)
                fmpz_poly_compose(f1, f1, linear)
                fmpz_poly_scalar_fdiv_ui(f1, f1, pp_ui)
                fmpz_poly_scalar_fdiv_ui(f1, f1, pp_ui)
                mpz_init(aa)
                mpz_init(bb)
                mpz_init(cc)
                mpz_init(dd)
                mpz_init(ee)
                fmpz_poly_get_coeff_mpz(aa, f1, 4)
                fmpz_poly_get_coeff_mpz(bb, f1, 3)
                fmpz_poly_get_coeff_mpz(cc, f1, 2)
                fmpz_poly_get_coeff_mpz(dd, f1, 1)
                fmpz_poly_get_coeff_mpz(ee, f1, 0)
                result = Zp_soluble_siksek(aa, bb, cc, dd, ee, pp, pp_ui, f_factzn, f, f1, linear)
                mpz_clear(aa)
                mpz_clear(bb)
                mpz_clear(cc)
                mpz_clear(dd)
                mpz_clear(ee)
                if result == 1:
                    break
        mpz_clear(aaa)
        mpz_clear(bbb)
        mpz_clear(ccc)
        mpz_clear(ddd)
        mpz_clear(eee)
        mpz_clear(tt)
        return result
    else:
        # Step II in Alg. 5.3.1 of Siksek's thesis
        nmod_poly_set_coeff_ui(f, 0, mpz_fdiv_ui(e, pp_ui))
        nmod_poly_set_coeff_ui(f, 1, mpz_fdiv_ui(d, pp_ui))
        nmod_poly_set_coeff_ui(f, 2, mpz_fdiv_ui(c, pp_ui))
        nmod_poly_set_coeff_ui(f, 3, mpz_fdiv_ui(b, pp_ui))
        nmod_poly_set_coeff_ui(f, 4, mpz_fdiv_ui(a, pp_ui))
        (<nmod_poly_factor_struct *>f_factzn)[0].num = 0 # reset data struct
        nmod_poly_factor(f_factzn, f)
        has_roots = 0
        has_single_roots = 0
        j = 0
        for i from 0 <= i < f_factzn.num:
            if nmod_poly_degree(&f_factzn.p[i]) == 1 and 0 != nmod_poly_get_coeff_ui(&f_factzn.p[i], 1):
                has_roots = 1
                if f_factzn.exp[i] == 1:
                    has_single_roots = 1
                    break
                roots[j] = pp_ui - nmod_poly_get_coeff_ui(&f_factzn.p[i], 0)
                j += 1

        if not has_roots: return 0
        if has_single_roots: return 1

        result = 0
        if j > 0:
            mpz_init(aa)
            mpz_init(bb)
            mpz_init(cc)
            mpz_init(dd)
            mpz_init(ee)
        for i from 0 <= i < j:
            fmpz_poly_zero(f1)
            fmpz_poly_zero(linear)
            fmpz_poly_set_coeff_mpz(f1, 0, e)
            fmpz_poly_set_coeff_mpz(f1, 1, d)
            fmpz_poly_set_coeff_mpz(f1, 2, c)
            fmpz_poly_set_coeff_mpz(f1, 3, b)
            fmpz_poly_set_coeff_mpz(f1, 4, a)
            fmpz_poly_set_coeff_ui(linear, 0, roots[i])
            fmpz_poly_set_coeff_mpz(linear, 1, pp)
            fmpz_poly_compose(f1, f1, linear)
            fmpz_poly_scalar_fdiv_ui(f1, f1, pp_ui)
            fmpz_poly_get_coeff_mpz(aa, f1, 4)
            fmpz_poly_get_coeff_mpz(bb, f1, 3)
            fmpz_poly_get_coeff_mpz(cc, f1, 2)
            fmpz_poly_get_coeff_mpz(dd, f1, 1)
            fmpz_poly_get_coeff_mpz(ee, f1, 0)
            result = Zp_soluble_siksek(aa, bb, cc, dd, ee, pp, pp_ui, f_factzn, f, f1, linear)
            if result == 1:
                break
        if j > 0:
            mpz_clear(aa)
            mpz_clear(bb)
            mpz_clear(cc)
            mpz_clear(dd)
            mpz_clear(ee)
        return result

cdef bint Zp_soluble_siksek_large_p(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e, mpz_t pp,
                                    fmpz_poly_t f1, fmpz_poly_t linear):
    """
    Uses the approach of Algorithm 5.3.1 of Siksek's thesis to test for
    solubility of y^2 == ax^4 + bx^3 + cx^2 + dx + e over Zp.
    """
    cdef unsigned long v_min, v
    cdef mpz_t roots[4]
    cdef int i, j, has_roots, has_single_roots
    cdef bint result

    cdef mpz_t aa, bb, cc, dd, ee
    cdef mpz_t aaa, bbb, ccc, ddd, eee
    cdef mpz_t qq, rr, ss, tt
    cdef Integer A,B,C,D,E,P

    # Step 0: divide out all common p from the quartic
    v_min = valuation(a, pp)
    if mpz_cmp_ui(b, ui0) != 0:
        v = valuation(b, pp)
        if v < v_min: v_min = v
    if mpz_cmp_ui(c, ui0) != 0:
        v = valuation(c, pp)
        if v < v_min: v_min = v
    if mpz_cmp_ui(d, ui0) != 0:
        v = valuation(d, pp)
        if v < v_min: v_min = v
    if mpz_cmp_ui(e, ui0) != 0:
        v = valuation(e, pp)
        if v < v_min: v_min = v
    for 0 <= v < v_min:
        mpz_divexact(a, a, pp)
        mpz_divexact(b, b, pp)
        mpz_divexact(c, c, pp)
        mpz_divexact(d, d, pp)
        mpz_divexact(e, e, pp)

    if not v_min%2:
        # Step I in Alg. 5.3.1 of Siksek's thesis
        A = Integer(0); B = Integer(0); C = Integer(0); D = Integer(0); E = Integer(0); P = Integer(0)
        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); mpz_set(P.value, pp)
        f = ntl.ZZ_pX([E,D,C,B,A], P)
        f /= ntl.ZZ_pX([A], P) # now f is monic, and we are done with A,B,C,D,E
        mpz_set(qq, A.value) # qq is the leading coefficient of the polynomial
        f_factzn = f.factor()
        result = 0
        for factor, exponent in f_factzn:
            if exponent&1:
                result = 1
                break
        if result == 0 and mpz_legendre(qq, pp) == 1:
            result = 1
        if result:
            return 1

        f = ntl.ZZ_pX([1], P)
        for factor, exponent in f_factzn:
            for j from 0 <= j < (exponent/2):
                f *= factor

        f /= f.leading_coefficient()
        f_factzn = f.factor()

        has_roots = 0
        j = 0
        for factor, exponent in f_factzn:
            if factor.degree() == 1:
                has_roots = 1
                A = P - Integer(factor[0])
                mpz_set(roots[j], A.value)
                j += 1
        if not has_roots:
            return 0

        i = f.degree()
        mpz_init(aaa)
        mpz_init(bbb)
        mpz_init(ccc)
        mpz_init(ddd)
        mpz_init(eee)

        if i == 0: # g == 1
            mpz_set(aaa, a)
            mpz_set(bbb, b)
            mpz_set(ccc, c)
            mpz_set(ddd, d)
            mpz_sub(eee, e, qq)
        elif i == 1: # g == x + rr
            mpz_set(aaa, a)
            mpz_set(bbb, b)
            mpz_sub(ccc, c, qq)
            A = Integer(f[0])
            mpz_set(rr, A.value)
            mpz_mul(ss, rr, qq)
            mpz_set(ddd,d)
            mpz_sub(ddd, ddd, ss)
            mpz_sub(ddd, ddd, ss)
            mpz_set(eee,e)
            mpz_mul(ss, ss, rr)
            mpz_sub(eee, eee, ss)
            mpz_divexact(ss, ss, rr)
        elif i == 2: # g == x^2 + rr*x + ss
            mpz_sub(aaa, a, qq)
            A = Integer(f[1])
            mpz_set(rr, A.value)
            mpz_init(tt)
            mpz_mul(tt, rr, qq)
            mpz_set(bbb,b)
            mpz_submul_ui(bbb, tt, ui2)
            mpz_set(ccc,c)
            mpz_submul(ccc, tt, rr)
            A = Integer(f[0])
            mpz_set(ss, A.value)
            mpz_mul(tt, ss, qq)
            mpz_set(eee,e)
            mpz_submul(eee, tt, ss)
            mpz_mul_ui(tt, tt, ui2)
            mpz_sub(ccc, ccc, tt)
            mpz_set(ddd,d)
            mpz_submul(ddd, tt, rr)
            mpz_clear(tt)
        mpz_divexact(aaa, aaa, pp)
        mpz_divexact(bbb, bbb, pp)
        mpz_divexact(ccc, ccc, pp)
        mpz_divexact(ddd, ddd, pp)
        mpz_divexact(eee, eee, pp)
        # now aaa,bbb,ccc,ddd,eee represents h(x)

        result = 0
        mpz_init(tt)
        for i from 0 <= i < j:
            mpz_mul(tt, aaa, roots[i])
            mpz_add(tt, tt, bbb)
            mpz_mul(tt, tt, roots[i])
            mpz_add(tt, tt, ccc)
            mpz_mul(tt, tt, roots[i])
            mpz_add(tt, tt, ddd)
            mpz_mul(tt, tt, roots[i])
            mpz_add(tt, tt, eee)
            # tt == h(r) mod p
            mpz_mod(tt, tt, pp)
            if mpz_sgn(tt) == 0:
                fmpz_poly_zero(f1)
                fmpz_poly_zero(linear)
                fmpz_poly_set_coeff_mpz(f1, 0, e)
                fmpz_poly_set_coeff_mpz(f1, 1, d)
                fmpz_poly_set_coeff_mpz(f1, 2, c)
                fmpz_poly_set_coeff_mpz(f1, 3, b)
                fmpz_poly_set_coeff_mpz(f1, 4, a)
                fmpz_poly_set_coeff_mpz(linear, 0, roots[i])
                fmpz_poly_set_coeff_mpz(linear, 1, pp)
                fmpz_poly_compose(f1, f1, linear)
                fmpz_poly_scalar_fdiv_mpz(f1, f1, pp)
                fmpz_poly_scalar_fdiv_mpz(f1, f1, pp)

                mpz_init(aa)
                mpz_init(bb)
                mpz_init(cc)
                mpz_init(dd)
                mpz_init(ee)
                fmpz_poly_get_coeff_mpz(aa, f1, 4)
                fmpz_poly_get_coeff_mpz(bb, f1, 3)
                fmpz_poly_get_coeff_mpz(cc, f1, 2)
                fmpz_poly_get_coeff_mpz(dd, f1, 1)
                fmpz_poly_get_coeff_mpz(ee, f1, 0)
                result = Zp_soluble_siksek_large_p(aa, bb, cc, dd, ee, pp, f1, linear)
                mpz_clear(aa)
                mpz_clear(bb)
                mpz_clear(cc)
                mpz_clear(dd)
                mpz_clear(ee)
                if result == 1:
                    break
        mpz_clear(aaa)
        mpz_clear(bbb)
        mpz_clear(ccc)
        mpz_clear(ddd)
        mpz_clear(eee)
        mpz_clear(tt)
        return result
    else:
        # Step II in Alg. 5.3.1 of Siksek's thesis
        A = Integer(0); B = Integer(0); C = Integer(0); D = Integer(0); E = Integer(0); P = Integer(0)
        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); mpz_set(P.value, pp)
        f = ntl.ZZ_pX([E,D,C,B,A], P)
        f /= ntl.ZZ_pX([A], P) # now f is monic
        f_factzn = f.factor()

        has_roots = 0
        has_single_roots = 0
        j = 0
        for factor, exponent in f_factzn:
            if factor.degree() == 1:
                has_roots = 1
                if exponent == 1:
                    has_single_roots = 1
                    break
                A = P - Integer(factor[0])
                mpz_set(roots[j], A.value)
                j += 1

        if not has_roots: return 0
        if has_single_roots: return 1

        result = 0
        if j > 0:
            mpz_init(aa)
            mpz_init(bb)
            mpz_init(cc)
            mpz_init(dd)
            mpz_init(ee)
        for i from 0 <= i < j:
            fmpz_poly_zero(f1)
            fmpz_poly_zero(linear)
            fmpz_poly_set_coeff_mpz(f1, 0, e)
            fmpz_poly_set_coeff_mpz(f1, 1, d)
            fmpz_poly_set_coeff_mpz(f1, 2, c)
            fmpz_poly_set_coeff_mpz(f1, 3, b)
            fmpz_poly_set_coeff_mpz(f1, 4, a)
            fmpz_poly_set_coeff_mpz(linear, 0, roots[i])
            fmpz_poly_set_coeff_mpz(linear, 1, pp)
            fmpz_poly_compose(f1, f1, linear)
            fmpz_poly_scalar_fdiv_mpz(f1, f1, pp)
            fmpz_poly_get_coeff_mpz(aa, f1, 4)
            fmpz_poly_get_coeff_mpz(bb, f1, 3)
            fmpz_poly_get_coeff_mpz(cc, f1, 2)
            fmpz_poly_get_coeff_mpz(dd, f1, 1)
            fmpz_poly_get_coeff_mpz(ee, f1, 0)
            result = Zp_soluble_siksek_large_p(aa, bb, cc, dd, ee, pp, f1, linear)
            if result == 1:
                break
        if j > 0:
            mpz_clear(aa)
            mpz_clear(bb)
            mpz_clear(cc)
            mpz_clear(dd)
            mpz_clear(ee)
        return result

cdef bint Qp_soluble_siksek(mpz_t A, mpz_t B, mpz_t C, mpz_t D, mpz_t E,
                            mpz_t p, unsigned long P,
                            nmod_poly_factor_t f_factzn, fmpz_poly_t f1,
                            fmpz_poly_t linear):
    """
    Uses Samir Siksek's thesis results to determine whether the quartic is
    locally soluble at p.
    """
    cdef int result = 0
    cdef mpz_t a,b,c,d,e
    cdef nmod_poly_t f
    nmod_poly_init(f, P)

    mpz_init_set(a,A)
    mpz_init_set(b,B)
    mpz_init_set(c,C)
    mpz_init_set(d,D)
    mpz_init_set(e,E)

    if Zp_soluble_siksek(a,b,c,d,e,p,P,f_factzn, f, f1, linear):
        result = 1
    else:
        mpz_set(a,A)
        mpz_set(b,B)
        mpz_set(c,C)
        mpz_set(d,D)
        mpz_set(e,E)
        if Zp_soluble_siksek(e,d,c,b,a,p,P,f_factzn, f, f1, linear):
            result = 1

    mpz_clear(a)
    mpz_clear(b)
    mpz_clear(c)
    mpz_clear(d)
    mpz_clear(e)
    nmod_poly_clear(f)
    return result

cdef bint Qp_soluble_siksek_large_p(mpz_t A, mpz_t B, mpz_t C, mpz_t D, mpz_t E,
                                    mpz_t p, fmpz_poly_t f1, fmpz_poly_t linear):
    """
    Uses Samir Siksek's thesis results to determine whether the quartic is
    locally soluble at p, when p is bigger than wordsize, and we can't use
    FLINT.
    """
    cdef int result = 0
    cdef mpz_t a,b,c,d,e

    mpz_init_set(a,A)
    mpz_init_set(b,B)
    mpz_init_set(c,C)
    mpz_init_set(d,D)
    mpz_init_set(e,E)

    if Zp_soluble_siksek_large_p(a,b,c,d,e,p,f1,linear):
        result = 1
    else:
        mpz_set(a,A)
        mpz_set(b,B)
        mpz_set(c,C)
        mpz_set(d,D)
        mpz_set(e,E)
        if Zp_soluble_siksek_large_p(e,d,c,b,a,p,f1,linear):
            result = 1

    mpz_clear(a)
    mpz_clear(b)
    mpz_clear(c)
    mpz_clear(d)
    mpz_clear(e)
    return result

cdef bint Qp_soluble_BSD(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e, mpz_t p):
    """
    Uses the original test of Birch and Swinnerton-Dyer to test for local
    solubility of the quartic at p.
    """
    cdef mpz_t zero
    cdef int result = 0
    mpz_init_set_ui(zero, ui0)
    if Zp_soluble_BSD(a,b,c,d,e,zero,p,0):
        result = 1
    elif Zp_soluble_BSD(e,d,c,b,a,zero,p,1):
        result = 1
    mpz_clear(zero)
    return result

cdef bint Qp_soluble(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e, mpz_t p):
    """
    Try the BSD approach for a few residue classes and if no solution is found,
    switch to Siksek to try to prove insolubility.
    """
    cdef int bsd_sol, sik_sol
    cdef unsigned long pp
    cdef fmpz_poly_t f1, linear
    cdef nmod_poly_factor_t f_factzn
    bsd_sol = Qp_soluble_BSD(a, b, c, d, e, p)
    if mpz_cmp_ui(p,N_RES_CLASSES_BSD)>0 and not bsd_sol:
        fmpz_poly_init(f1)
        fmpz_poly_init(linear)
        if mpz_fits_ulong_p(p):
            nmod_poly_factor_init(f_factzn)
            pp = mpz_get_ui(p)
            sik_sol = Qp_soluble_siksek(a,b,c,d,e,p,pp,f_factzn,f1,linear)
            nmod_poly_factor_clear(f_factzn)
        else:
            sik_sol = Qp_soluble_siksek_large_p(a,b,c,d,e,p,f1,linear)
        fmpz_poly_clear(f1)
        fmpz_poly_clear(linear)
    else:
        sik_sol = bsd_sol
    return sik_sol

def test_qpls(a,b,c,d,e,p):
    """
    Testing function for Qp_soluble.

    EXAMPLE:
        sage: from sage.schemes.elliptic_curves.descent_two_isogeny import test_qpls as tq
        sage: tq(1,2,3,4,5,7)
        1

    """
    cdef Integer A,B,C,D,E,P
    cdef int i, result
    cdef mpz_t aa,bb,cc,dd,ee,pp
    A=Integer(a); B=Integer(b); C=Integer(c); D=Integer(d); E=Integer(e); P=Integer(p)
    mpz_init_set(aa, A.value)
    mpz_init_set(bb, B.value)
    mpz_init_set(cc, C.value)
    mpz_init_set(dd, D.value)
    mpz_init_set(ee, E.value)
    mpz_init_set(pp, P.value)
    result = Qp_soluble(aa, bb, cc, dd, ee, pp)
    mpz_clear(aa)
    mpz_clear(bb)
    mpz_clear(cc)
    mpz_clear(dd)
    mpz_clear(ee)
    mpz_clear(pp)
    return result

cdef int everywhere_locally_soluble(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e) except -1:
    """
    Returns whether the quartic has local solutions at all primes p.
    """
    cdef Integer A,B,C,D,E,Delta,p
    cdef mpz_t mpz_2
    A=Integer(0); B=Integer(0); C=Integer(0); D=Integer(0); E=Integer(0)
    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);
    f = (((A*x_ZZ + B)*x_ZZ + C)*x_ZZ + D)*x_ZZ + E

    # RR soluble:
    if mpz_sgn(a)!=1:
        if len(real_roots(f)) == 0:
            return 0

    # Q2 soluble:
    mpz_init_set_ui(mpz_2, ui2)
    if not Qp_soluble(a,b,c,d,e,mpz_2):
        mpz_clear(mpz_2)
        return 0
    mpz_clear(mpz_2)

    # Odd finite primes
    Delta = f.discriminant()
    for p in prime_divisors(Delta):
        if p == 2: continue
        if not Qp_soluble(a,b,c,d,e,p.value): return 0

    return 1

def test_els(a,b,c,d,e):
    """
    Doctest function for cdef int everywhere_locally_soluble(mpz_t, mpz_t, mpz_t, mpz_t, mpz_t).

    EXAMPLE::

        sage: from sage.schemes.elliptic_curves.descent_two_isogeny import test_els
        sage: from sage.libs.ratpoints import ratpoints
        sage: for _ in range(1000):
        ...    a,b,c,d,e = randint(1,1000), randint(1,1000), randint(1,1000), randint(1,1000), randint(1,1000)
        ...    if len(ratpoints([e,d,c,b,a], 1000)) > 0:
        ...        try:
        ...            if not test_els(a,b,c,d,e):
        ...                print "This never happened", a,b,c,d,e
        ...        except ValueError:
        ...            continue

    """
    cdef Integer A,B,C,D,E,Delta
    A=Integer(a); B=Integer(b); C=Integer(c); D=Integer(d); E=Integer(e)
    return everywhere_locally_soluble(A.value, B.value, C.value, D.value, E.value)

cdef int count(mpz_t c_mpz, mpz_t d_mpz, mpz_t *p_list, unsigned long p_list_len,
               int global_limit_small, int global_limit_large,
               int verbosity, bint selmer_only, mpz_t n1, mpz_t n2) except -1:
    """
    Count the number of els/gls quartic 2-covers of E.
    """
    cdef unsigned long n_primes, i
    cdef bint found_global_points, els, check_negs, verbose = (verbosity > 4)
    cdef Integer a_Int, c_Int, e_Int
    cdef mpz_t c_sq_mpz, d_prime_mpz
    cdef mpz_t n_divisors, j
    cdef mpz_t *coeffs_ratp


    mpz_init(c_sq_mpz)
    mpz_mul(c_sq_mpz, c_mpz, c_mpz)
    mpz_init_set(d_prime_mpz, c_sq_mpz)
    mpz_submul_ui(d_prime_mpz, d_mpz, ui4)
    check_negs = 0
    if mpz_sgn(d_prime_mpz) > 0:
        if mpz_sgn(c_mpz) >= 0 or mpz_cmp(c_sq_mpz, d_prime_mpz) <= 0:
            check_negs = 1
    mpz_clear(c_sq_mpz)
    mpz_clear(d_prime_mpz)


    # Set up coefficient array, and static variables
    cdef mpz_t *coeffs = <mpz_t *> sage_malloc(5 * sizeof(mpz_t))
    for i from 0 <= i <= 4:
        mpz_init(coeffs[i])
    mpz_set_ui(coeffs[1], ui0)     #
    mpz_set(coeffs[2], c_mpz)      # These never change
    mpz_set_ui(coeffs[3], ui0)     #

    if not selmer_only:
        # allocate space for ratpoints
        coeffs_ratp = <mpz_t *> sage_malloc(5 * sizeof(mpz_t))
        for i from 0 <= i <= 4:
            mpz_init(coeffs_ratp[i])

    # Get prime divisors, and put them in an mpz_t array
    # (this block, by setting check_negs, takes care of
    # local solubility over RR)
    cdef mpz_t *p_div_d_mpz = <mpz_t *> sage_malloc((p_list_len+1) * sizeof(mpz_t))
    n_primes = 0
    for i from 0 <= i < p_list_len:
        if mpz_divisible_p(d_mpz, p_list[i]):
            mpz_init(p_div_d_mpz[n_primes])
            mpz_set(p_div_d_mpz[n_primes], p_list[i])
            n_primes += 1
    if check_negs:
        mpz_init(p_div_d_mpz[n_primes])
        mpz_set_si(p_div_d_mpz[n_primes], -1)
        n_primes += 1
    mpz_init_set_ui(n_divisors, ui1)
    mpz_mul_2exp(n_divisors, n_divisors, n_primes)
#    if verbosity > 3:
#        print '\nDivisors of d which may lead to RR-soluble quartics:', p_div_d

    mpz_init_set_ui(j, ui0)
    if not selmer_only:
        mpz_set_ui(n1, ui0)
    mpz_set_ui(n2, ui0)
    while mpz_cmp(j, n_divisors) < 0:
        mpz_set_ui(coeffs[4], ui1)
        for i from 0 <= i < n_primes:
            if mpz_tstbit(j, i):
                mpz_mul(coeffs[4], coeffs[4], p_div_d_mpz[i])
        if verbosity > 3:
            a_Int = Integer(0); mpz_set(a_Int.value, coeffs[4])
            print '\nSquarefree divisor:', a_Int
        mpz_divexact(coeffs[0], d_mpz, coeffs[4])
        found_global_points = 0
        if not selmer_only:
            if verbose:
                print "\nCalling ratpoints for small point search"
            for i from 0 <= i <= 4:
                mpz_set(coeffs_ratp[i], coeffs[i])
            sig_on()
            found_global_points = ratpoints_mpz_exists_only(coeffs_ratp, global_limit_small, 4, verbose)
            sig_off()
            if found_global_points:
                if verbosity > 2:
                    a_Int = Integer(0); mpz_set(a_Int.value, coeffs[4])
                    c_Int = Integer(0); mpz_set(c_Int.value, coeffs[2])
                    e_Int = Integer(0); mpz_set(e_Int.value, coeffs[0])
                    print 'Found small global point, quartic (%d,%d,%d,%d,%d)'%(a_Int,0,c_Int,0,e_Int)
                mpz_add_ui(n1, n1, ui1)
                mpz_add_ui(n2, n2, ui1)
            if verbose:
                print "\nDone calling ratpoints for small point search"
        if not found_global_points:
            # Test whether the quartic is everywhere locally soluble:
            els = 1
            for i from 0 <= i < p_list_len:
                if not Qp_soluble(coeffs[4], coeffs[3], coeffs[2], coeffs[1], coeffs[0], p_list[i]):
                    els = 0
                    break
            if els:
                if verbosity > 2:
                    a_Int = Integer(0); mpz_set(a_Int.value, coeffs[4])
                    c_Int = Integer(0); mpz_set(c_Int.value, coeffs[2])
                    e_Int = Integer(0); mpz_set(e_Int.value, coeffs[0])
                    print 'ELS without small global points, quartic (%d,%d,%d,%d,%d)'%(a_Int,0,c_Int,0,e_Int)
                mpz_add_ui(n2, n2, ui1)
                if not selmer_only:
                    if verbose:
                        print "\nCalling ratpoints for large point search"
                    for i from 0 <= i <= 4:
                        mpz_set(coeffs_ratp[i], coeffs[i])
                    sig_on()
                    found_global_points = ratpoints_mpz_exists_only(coeffs_ratp, global_limit_large, 4, verbose)
                    sig_off()
                    if found_global_points:
                        if verbosity > 2:
                            print '  -- Found large global point.'
                        mpz_add_ui(n1, n1, ui1)
                    if verbose:
                        print "\nDone calling ratpoints for large point search"
        mpz_add_ui(j, j, ui1)
    if not selmer_only:
        for i from 0 <= i <= 4:
            mpz_clear(coeffs_ratp[i])
        sage_free(coeffs_ratp)
    mpz_clear(j)
    for i from 0 <= i < n_primes:
        mpz_clear(p_div_d_mpz[i])
    sage_free(p_div_d_mpz)
    mpz_clear(n_divisors)
    for i from 0 <= i <= 4:
        mpz_clear(coeffs[i])
    sage_free(coeffs)
    return 0

def two_descent_by_two_isogeny(E,
                int global_limit_small = 10,
                int global_limit_large = 10000,
                int verbosity = 0,
                bint selmer_only = 0, bint proof = 1):
    """
    Given an elliptic curve E with a two-isogeny phi : E --> E' and dual isogeny
    phi', runs a two-isogeny descent on E, returning n1, n2, n1' and n2'. Here
    n1 is the number of quartic covers found with a rational point, and n2 is
    the number which are ELS.

    EXAMPLES::

        sage: from sage.schemes.elliptic_curves.descent_two_isogeny import two_descent_by_two_isogeny
        sage: E = EllipticCurve('14a')
        sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny(E)
        sage: log(n1,2) + log(n1_prime,2) - 2 # the rank
        0
        sage: E = EllipticCurve('65a')
        sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny(E)
        sage: log(n1,2) + log(n1_prime,2) - 2 # the rank
        1
        sage: x,y = var('x,y')
        sage: E = EllipticCurve(y^2 == x^3 + x^2 - 25*x + 39)
        sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny(E)
        sage: log(n1,2) + log(n1_prime,2) - 2 # the rank
        2
        sage: E = EllipticCurve(y^2 + x*y + y == x^3 - 131*x + 558)
        sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny(E)
        sage: log(n1,2) + log(n1_prime,2) - 2 # the rank
        3

    Using the verbosity option::

        sage: E = EllipticCurve('14a')
        sage: two_descent_by_two_isogeny(E, verbosity=1)
        2-isogeny
        Results:
        2 <= #E(Q)/phi'(E'(Q)) <= 2
        2 <= #E'(Q)/phi(E(Q)) <= 2
        #Sel^(phi')(E'/Q) = 2
        #Sel^(phi)(E/Q) = 2
        1 <= #Sha(E'/Q)[phi'] <= 1
        1 <= #Sha(E/Q)[phi] <= 1
        1 <= #Sha(E/Q)[2], #Sha(E'/Q)[2] <= 1
        0 <= rank of E(Q) = rank of E'(Q) <= 0
        (2, 2, 2, 2)

    Handling curves whose discriminants involve larger than wordsize primes::

        sage: E = EllipticCurve('14a')
        sage: E = E.quadratic_twist(next_prime(10^20))
        sage: E
        Elliptic Curve defined by y^2 = x^3 + x^2 + 716666666666666667225666666666666666775672*x - 391925925925925926384240370370370370549019837037037037060249356 over Rational Field
        sage: E.discriminant().factor()
        -1 * 2^18 * 7^3 * 100000000000000000039^6
        sage: log(100000000000000000039.0, 2.0)
        66.438...
        sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny(E)
        sage: log(n1,2) + log(n1_prime,2) - 2 # the rank
        0

    TESTS:

    Here we contrive an example to demonstrate that a keyboard interrupt
    is caught. Here we let `E` be the smallest optimal curve with two-torsion
    and nontrivial `Sha[2]`. This ensures that the two-descent will be looking
    for rational points which do not exist, and by setting global_limit_large
    to a very high bound, it will still be working when we simulate a ``CTRL-C``::

        sage: from sage.schemes.elliptic_curves.descent_two_isogeny import two_descent_by_two_isogeny
        sage: import sage.tests.interrupt
        sage: E = EllipticCurve('960d'); E
        Elliptic Curve defined by y^2 = x^3 - x^2 - 900*x - 10098 over Rational Field
        sage: E.sha().an()
        4
        sage: try:
        ...     sage.tests.interrupt.interrupt_after_delay(1000)
        ...     two_descent_by_two_isogeny(E, global_limit_large=10^8)
        ... except KeyboardInterrupt:
        ...     print "Caught!"
        Caught!
    """
    cdef Integer a1, a2, a3, a4, a6, s2, s4, s6
    cdef Integer c, d, x0
    cdef list x_list
    assert E.torsion_order()%2==0, 'Need rational two-torsion for isogeny descent.'
    if verbosity > 0:
        print '\n2-isogeny'
        if verbosity > 1:
            print '\nchanging coordinates'
    a1 = Integer(E.a1())
    a2 = Integer(E.a2())
    a3 = Integer(E.a3())
    a4 = Integer(E.a4())
    a6 = Integer(E.a6())
    if a1==0 and a3==0:
        s2=a2; s4=a4; s6=a6
    else:
        s2=a1*a1+4*a2; s4=8*(a1*a3+2*a4); s6=16*(a3*a3+4*a6)
    f = ((x_ZZ + s2)*x_ZZ + s4)*x_ZZ + s6
    x_list = f.roots() # over ZZ -- use FLINT directly?
    x0 = x_list[0][0]
    c = 3*x0+s2;  d = (c+s2)*x0+s4
    return two_descent_by_two_isogeny_work(c, d,
        global_limit_small, global_limit_large, verbosity, selmer_only, proof)

def two_descent_by_two_isogeny_work(Integer c, Integer d,
                int global_limit_small = 10, int global_limit_large = 10000,
                int verbosity = 0, bint selmer_only = 0, bint proof = 1):
    """
    Do all the work in doing a two-isogeny descent.

    EXAMPLES::

        sage: from sage.schemes.elliptic_curves.descent_two_isogeny import two_descent_by_two_isogeny_work
        sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny_work(13,128)
        sage: log(n1,2) + log(n1_prime,2) - 2 # the rank
        0
        sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny_work(1,-16)
        sage: log(n1,2) + log(n1_prime,2) - 2 # the rank
        1
        sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny_work(10,8)
        sage: log(n1,2) + log(n1_prime,2) - 2 # the rank
        2
        sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny_work(85,320)
        sage: log(n1,2) + log(n1_prime,2) - 2 # the rank
        3

    """
    cdef mpz_t c_mpz, d_mpz, c_prime_mpz, d_prime_mpz
    cdef mpz_t *p_list_mpz
    cdef unsigned long i, j, p, p_list_len
    cdef Integer P, n1, n2, n1_prime, n2_prime, c_prime, d_prime
    cdef object PO
    cdef bint found, too_big, d_neg, d_prime_neg
    cdef n_factor_t fact
    cdef list primes
    mpz_init_set(c_mpz, c.value)      #
    mpz_init_set(d_mpz, d.value)      #
    mpz_init(c_prime_mpz)             #
    mpz_init(d_prime_mpz)             #
    mpz_mul_si(c_prime_mpz, c_mpz, -2)
    mpz_mul(d_prime_mpz, c_mpz, c_mpz)
    mpz_submul_ui(d_prime_mpz, d_mpz, ui4)

    d_neg = 0
    d_prime_neg = 0
    if mpz_sgn(d_mpz) < 0:
        d_neg = 1
        mpz_neg(d_mpz, d_mpz)
    if mpz_sgn(d_prime_mpz) < 0:
        d_prime_neg = 1
        mpz_neg(d_prime_mpz, d_prime_mpz)
    if mpz_fits_ulong_p(d_mpz) and mpz_fits_ulong_p(d_prime_mpz):
        # Factor very quickly using FLINT.
        p_list_mpz = <mpz_t *> sage_malloc(20 * sizeof(mpz_t))
        mpz_init_set_ui(p_list_mpz[0], ui2)
        p_list_len = 1
        n_factor_init(&fact)
        n_factor(&fact, mpz_get_ui(d_mpz), proof)
        for i from 0 <= i < fact.num:
            p = fact.p[i]
            if p != ui2:
                mpz_init_set_ui(p_list_mpz[p_list_len], p)
                p_list_len += 1
        n_factor(&fact, mpz_get_ui(d_prime_mpz), proof)
        for i from 0 <= i < fact.num:
            p = fact.p[i]
            found = 0
            for j from 0 <= j < p_list_len:
                if mpz_cmp_ui(p_list_mpz[j], p)==0:
                    found = 1
                    break
            if not found:
                mpz_init_set_ui(p_list_mpz[p_list_len], p)
                p_list_len += 1
    else:
        # Factor more slowly using Pari via Python.
        from sage.libs.pari.all import pari
        d = Integer(0)
        mpz_set(d.value, d_mpz)
        primes = list(pari(d).factor()[0])
        d_prime = Integer(0)
        mpz_set(d_prime.value, d_prime_mpz)
        for PO in pari(d_prime).factor()[0]:
            if PO not in primes:
                primes.append(PO)
        P = Integer(2)
        if P not in primes: primes.append(P)
        p_list_len = len(primes)
        p_list_mpz = <mpz_t *> sage_malloc(p_list_len * sizeof(mpz_t))
        for i from 0 <= i < p_list_len:
            P = Integer(primes[i])
            mpz_init_set(p_list_mpz[i], P.value)
    if d_neg:
        mpz_neg(d_mpz, d_mpz)
    if d_prime_neg:
        mpz_neg(d_prime_mpz, d_prime_mpz)

    if verbosity > 1:
        c_prime = -2*c
        d_prime = c*c-4*d
        print '\nnew curve is y^2 == x( x^2 + (%d)x + (%d) )'%(int(c),int(d))
        print 'new isogenous curve is' + \
               ' y^2 == x( x^2 + (%d)x + (%d) )'%(int(c_prime),int(d_prime))

    n1 = Integer(0); n2 = Integer(0)
    n1_prime = Integer(0); n2_prime = Integer(0)
    count(c.value, d.value, p_list_mpz, p_list_len,
          global_limit_small, global_limit_large, verbosity, selmer_only,
          n1.value, n2.value)
    count(c_prime_mpz, d_prime_mpz, p_list_mpz, p_list_len,
          global_limit_small, global_limit_large, verbosity, selmer_only,
          n1_prime.value, n2_prime.value)

    for i from 0 <= i < p_list_len:
        mpz_clear(p_list_mpz[i])
    sage_free(p_list_mpz)

    if verbosity > 0:
        print "\nResults:"
        print n1, "<= #E(Q)/phi'(E'(Q)) <=", n2
        print n1_prime, "<= #E'(Q)/phi(E(Q)) <=", n2_prime
        print "#Sel^(phi')(E'/Q) =", n2
        print "#Sel^(phi)(E/Q) =", n2_prime
        print "1 <= #Sha(E'/Q)[phi'] <=", n2/n1
        print "1 <= #Sha(E/Q)[phi] <=", n2_prime/n1_prime
        print "1 <= #Sha(E/Q)[2], #Sha(E'/Q)[2] <=", (n2_prime/n1_prime)*(n2/n1)
        a = Integer(n1*n1_prime).log(Integer(2))
        e = Integer(n2*n2_prime).log(Integer(2))
        print a - 2, "<= rank of E(Q) = rank of E'(Q) <=", e - 2

    return n1, n2, n1_prime, n2_prime