Ticket #6583: trac_6583-rebase.patch

File trac_6583-rebase.patch, 48.9 KB (added by was, 10 years ago)

  • module_list.py 

    # HG changeset patch
# User rlm@rlmiller.org
# Date 1260318825 28800
# Node ID 59766f3a2186ad62598fc28e04f1d05c8de431b8
# Parent  18a4ec53e18a7dd51639dc0ea4f3bd9541430e32
trac 6583 -- Two-descent by a two-isogeny

diff --git a/module_list.py b/module_list.py
    a b  
    446446   
    447447    Extension('sage.libs.ratpoints',
    448448              sources = ["sage/libs/ratpoints.pyx"],
    449               #depends = [SAGE_ROOT + 'local/include/ratpoints.h'],
     449              depends = [SAGE_ROOT + '/local/include/ratpoints.h'],
    450450              libraries = ["ratpoints", "gmp"]),
    451451   
    452452    Extension('sage.libs.singular.singular',
    … …  
    12761276    ##
    12771277    ################################
    12781278
     1279    Extension('sage.schemes.elliptic_curves.descent_two_isogeny',
     1280              sources = ['sage/schemes/elliptic_curves/descent_two_isogeny.pyx'],
     1281              extra_compile_args=["-std=c99"],
     1282              depends = [SAGE_ROOT + '/local/include/ratpoints.h',
     1283                         SAGE_ROOT + '/local/include/gmp.h',
     1284                         SAGE_ROOT + '/local/include/FLINT/flint.h'],
     1285              include_dirs = [SAGE_ROOT+'/local/include/FLINT/'],
     1286              libraries = ['flint', 'gmp', 'ratpoints']),
     1287
    12791288    Extension('sage.schemes.hyperelliptic_curves.hypellfrob',
    12801289              sources = ['sage/schemes/hyperelliptic_curves/hypellfrob.pyx',
    12811290                         'sage/schemes/hyperelliptic_curves/hypellfrob/hypellfrob.cpp',

  • sage/libs/flint/fmpz_poly.pxi 

    diff --git a/sage/libs/flint/fmpz_poly.pxi b/sage/libs/flint/fmpz_poly.pxi
    a b  
    2525    void fmpz_poly_truncate(fmpz_poly_t poly, unsigned long length)
    2626
    2727    void fmpz_poly_set(fmpz_poly_t result, fmpz_poly_t poly)
     28    void fmpz_poly_zero(fmpz_poly_t poly)
    2829    void fmpz_poly_set_coeff_si(fmpz_poly_t poly, unsigned long n, long x)
    2930    void fmpz_poly_set_coeff_ui(fmpz_poly_t poly, unsigned long n, \
    3031            unsigned long x)
    … …  
    104105    unsigned long fmpz_poly_resultant_bound(fmpz_poly_t a, fmpz_poly_t b)
    105106    void fmpz_poly_resultant(fmpz_t r, fmpz_poly_t a, fmpz_poly_t b)
    106107
     108    void fmpz_poly_invmod(fmpz_t d, fmpz_poly_t H, fmpz_poly_t poly1, fmpz_poly_t poly2)
     109    void fmpz_poly_derivative(fmpz_poly_t der, fmpz_poly_t poly)
     110    void fmpz_poly_evaluate(fmpz_t output, fmpz_poly_t poly, fmpz_t val)
     111    void fmpz_poly_compose(fmpz_poly_t output, fmpz_poly_t f, fmpz_poly_t g)
     112    void fmpz_poly_scalar_div_ui(fmpz_poly_t output, fmpz_poly_t poly, unsigned long x)
     113
    107114    unsigned long fmpz_poly_max_limbs(fmpz_poly_t poly)

  • sage/libs/flint/long_extras.pxd 

    diff --git a/sage/libs/flint/long_extras.pxd b/sage/libs/flint/long_extras.pxd
    a b  
    9393    cdef unsigned long z_primitive_root(unsigned long p)
    9494
    9595    cdef unsigned long z_primitive_root_precomp(unsigned long p, double p_inv)
     96
     97    ctypedef struct factor_t:
     98        int num
     99        unsigned long p[15]
     100        unsigned long exp[15]
     101
     102    cdef unsigned long z_factor(factor_t *, unsigned long, int)
     103
     104

  • sage/libs/flint/zmod_poly.pxd 

    diff --git a/sage/libs/flint/zmod_poly.pxd b/sage/libs/flint/zmod_poly.pxd
    a b  
    1515
    1616    ctypedef zmod_poly_struct* zmod_poly_t
    1717
     18    ctypedef struct zmod_poly_factor_struct:
     19        zmod_poly_t *factors
     20        unsigned long *exponents
     21        unsigned long alloc
     22        unsigned long num_factors
     23
     24    ctypedef zmod_poly_factor_struct* zmod_poly_factor_t
     25
    1826    cdef void zmod_poly_init(zmod_poly_t poly, unsigned long p)
    1927    cdef void zmod_poly_init_precomp(zmod_poly_t poly, unsigned long p, double p_inv)
    2028    cdef void zmod_poly_init2(zmod_poly_t poly, unsigned long p, unsigned long alloc)
    … …  
    242250    cdef int zmod_poly_gcd_invert(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2)
    243251    cdef void zmod_poly_xgcd(zmod_poly_t res, zmod_poly_t s, zmod_poly_t t, zmod_poly_t poly1, zmod_poly_t poly2)
    244252
     253
     254
    245255    # Composition / evaluation
    246256
    247257    cdef unsigned long zmod_poly_evaluate(zmod_poly_t, unsigned long)
    … …  
    263273    cdef void zmod_poly_factor_square_free(zmod_poly_factor_t, zmod_poly_t)
    264274    cdef void zmod_poly_factor(zmod_poly_factor_t, zmod_poly_t)
    265275
    266     # FLINT 1.1 will have:
    267     # cdef void zmod_poly_powmod(zmod_poly_t res, zmod_poly_t pol, long exp, zmod_poly_t f)
     276    #
     277    # Differentiation
     278    #
     279 
     280    cdef void zmod_poly_derivative(zmod_poly_t res, zmod_poly_t poly)
     281   
     282    #
     283    # Arithmetic modulo a polynomial
     284    #
     285   
     286    cdef void zmod_poly_mulmod(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2, zmod_poly_t f)
     287    cdef void zmod_poly_powmod(zmod_poly_t res,zmod_poly_t pol, long exp, zmod_poly_t f)
     288
     289    #
     290    # Polynomial factorization
     291    #
     292
     293    cdef void zmod_poly_factor_init(zmod_poly_factor_t fac)
     294    cdef void zmod_poly_factor_clear(zmod_poly_factor_t fac)
     295    cdef void zmod_poly_factor_add(zmod_poly_factor_t fac, zmod_poly_t poly)
     296    cdef void zmod_poly_factor_concat(zmod_poly_factor_t res, zmod_poly_factor_t fac)
     297    cdef void zmod_poly_factor_print(zmod_poly_factor_t fac)
     298    cdef void zmod_poly_factor_pow(zmod_poly_factor_t fac, unsigned long exp)
     299    cdef void zmod_poly_factor_square_free(zmod_poly_factor_t res, zmod_poly_t f)
     300    cdef void zmod_poly_factor_berlekamp(zmod_poly_factor_t factors, zmod_poly_t f)
     301    cdef unsigned long zmod_poly_factor(zmod_poly_factor_t result, zmod_poly_t input)
     302    cdef int zmod_poly_isirreducible(zmod_poly_t f)
     303

  • sage/libs/pari/decl.pxi 

    diff --git a/sage/libs/pari/decl.pxi b/sage/libs/pari/decl.pxi
    a b  
    10861086    int     opgs2(int (*f)(GEN, GEN), GEN y, long s)
    10871087
    10881088    long    ZX_valuation(GEN x, GEN *Z)
     1089    GEN     ZX_mul(GEN x, GEN y)
    10891090    GEN     cgetimag()
    10901091    GEN     cgetp(GEN x)
    10911092    GEN     cvtop(GEN x, GEN p, long l)
    … …  
    11801181    GEN     gfrac(GEN x)
    11811182    GEN     gge(GEN x, GEN y)
    11821183    GEN     ggprecision(GEN x)
     1184    GEN     ggrando(GEN x, long n)
    11831185    GEN     ggt(GEN x, GEN y)
    11841186    GEN     gimag(GEN x)
    11851187    GEN     ginv(GEN x)
    … …  
    11871189    GEN     glt(GEN x, GEN y)
    11881190    GEN     gmod(GEN x, GEN y)
    11891191    GEN     gmodulcp(GEN x,GEN y)
     1192    GEN     gmodgs(GEN x, long y)
    11901193    GEN     gmodulo(GEN x,GEN y)
    11911194    GEN     gmodulsg(long x, GEN y)
    11921195    GEN     gmodulss(long x, long y)
    … …  
    12041207    GEN     gshift(GEN x, long n)
    12051208    GEN     gsubst(GEN x, long v, GEN y)
    12061209    GEN     gsubstpol(GEN x, GEN v, GEN y)
     1210    GEN     gsubstvec(GEN x, GEN v, GEN y)
    12071211    GEN     gtocol(GEN x)
    12081212    GEN     gtopoly(GEN x, long v)
    12091213    GEN     gtopolyrev(GEN x, long v)
    … …  
    12211225    int     iscomplex(GEN x)
    12221226    int     isexactzero(GEN g)
    12231227    int     isexactzeroscalar(GEN g)
     1228    int     isinexact(GEN x)
    12241229    int     isinexactreal(GEN x)
    12251230    long    isint(GEN n, long *ptk)
     1231    int     issmall(GEN n, long *ptk)
    12261232    int     ismonome(GEN x)
    12271233    GEN     lift(GEN x)
    12281234    GEN     lift0(GEN x,long v)
    12291235    GEN     lift_intern0(GEN x,long v)
    12301236    GEN     truncr(GEN x)
     1237    GEN     mkcoln(long n, ...)
     1238    GEN     mkintn(long n, ...)
     1239    GEN     mkpoln(long n, ...)
     1240    GEN     mkvecn(long n, ...)
    12311241    GEN     mulmat_real(GEN x, GEN y)
    12321242    GEN     numer(GEN x)
    12331243    GEN     padicappr(GEN f, GEN a)
    … …  
    12501260    GEN     scalarser(GEN x, long v, long prec)
    12511261    GEN     simplify(GEN x)
    12521262    GEN     simplify_i(GEN x)
     1263    GEN     tayl(GEN x, long v, long precdl)
     1264    GEN     toser_i(GEN x)
    12531265    GEN     truecoeff(GEN x, long n)
    12541266    GEN     trunc0(GEN x, GEN *pte)
    12551267    GEN     u2toi(ulong a, ulong b)
    … …  
    19631975    long *   int_precW(long * xp)
    19641976    long *   int_nextW(long * xp)
    19651977
     1978   
     1979

  • new file sage/libs/ratpoints.pxd 

    diff --git a/sage/libs/ratpoints.pxd b/sage/libs/ratpoints.pxd
new file mode 100644
    - +  
     1
     2include '../ext/cdefs.pxi'
     3include '../ext/stdsage.pxi'
     4from sage.rings.integer cimport Integer
     5
     6cdef extern from "ratpoints.h":
     7    long RATPOINTS_MAX_DEGREE
     8    long RATPOINTS_ARRAY_SIZE
     9    long RATPOINTS_DEFAULT_SP1
     10    long RATPOINTS_DEFAULT_SP2
     11    long RATPOINTS_DEFAULT_NUM_PRIMES
     12    long RATPOINTS_DEFAULT_MAX_FORBIDDEN
     13    long RATPOINTS_DEFAULT_STURM
     14    long RATPOINTS_NON_SQUAREFREE
     15    long RATPOINTS_BAD_ARGS
     16
     17    # for args flags:
     18    long RATPOINTS_NO_CHECK # when set, do not check whether the surviving
     19                            # x-coordinates give rise to rational points
     20    long RATPOINTS_NO_Y # when set, only list x coordinates instead of actual points
     21    long RATPOINTS_NO_REVERSE # when set, do not modify the mpz_t array
     22    long RATPOINTS_NO_JACOBI # when set, prevent use of Jacobi symbol test
     23    long RATPOINTS_VERBOSE # when set, print some output on what ratpoints is doing
     24    # define RATPOINTS_FLAGS_INPUT_MASK \
     25    # (RATPOINTS_NO_CHECK | RATPOINTS_NO_Y | RATPOINTS_NO_REVERSE | \
     26    #  RATPOINTS_NO_JACOBI | RATPOINTS_VERBOSE)
     27
     28
     29    ctypedef struct ratpoints_interval:
     30        double low
     31        double up
     32    ctypedef struct ratpoints_args:
     33        mpz_t *cof
     34        long degree
     35        long height
     36        ratpoints_interval *domain
     37        long num_inter
     38        long b_low
     39        long b_high
     40        long sp1
     41        long sp2
     42        long array_size
     43        long sturm
     44        long num_primes
     45        long max_forbidden
     46        unsigned int flags
     47        # from here: private data
     48        # mpz_t *work
     49        # void *se_buffer
     50        # void *se_next
     51        # void *ba_buffer
     52        # void *ba_next
     53        # int *int_buffer
     54        # int *int_next
     55        # void *sieve_list
     56    long find_points(ratpoints_args*, int proc(long, long, mpz_t, void*, int*), void*)
     57    void find_points_init(ratpoints_args*)
     58    long find_points_work(ratpoints_args*, int proc(long, long, mpz_t, void*, int*), void*)
     59    void find_points_clear(ratpoints_args*)
     60
     61ctypedef struct point_list:
     62    long *xes
     63    mpz_t *ys
     64    long *zs
     65    long array_size
     66    long num_points
     67    long max_num_points
     68
     69ctypedef struct info_struct_exists_only:
     70    int verbose
     71
     72cdef bint ratpoints_mpz_exists_only(mpz_t *, long, int, bint)
     73
     74
     75
     76

  • sage/libs/ratpoints.pyx 

    diff --git a/sage/libs/ratpoints.pyx b/sage/libs/ratpoints.pyx
    a b  
    33
    44"""
    55
    6 include '../ext/cdefs.pxi'
    7 include '../ext/stdsage.pxi'
    8 from sage.rings.integer cimport Integer
    9 cdef extern from "ratpoints.h":
    10     long RATPOINTS_MAX_DEGREE
    11     long RATPOINTS_ARRAY_SIZE
    12     long RATPOINTS_DEFAULT_SP1
    13     long RATPOINTS_DEFAULT_SP2
    14     long RATPOINTS_DEFAULT_NUM_PRIMES
    15     long RATPOINTS_DEFAULT_MAX_FORBIDDEN
    16     long RATPOINTS_DEFAULT_STURM
    17     long RATPOINTS_NON_SQUAREFREE
    18     long RATPOINTS_BAD_ARGS
    19 
    20     # for args flags:
    21     long RATPOINTS_NO_CHECK # when set, do not check whether the surviving
    22                             # x-coordinates give rise to rational points
    23     long RATPOINTS_NO_Y # when set, only list x coordinates instead of actual points
    24     long RATPOINTS_NO_REVERSE # when set, do not modify the mpz_t array
    25     long RATPOINTS_NO_JACOBI # when set, prevent use of Jacobi symbol test
    26     long RATPOINTS_VERBOSE # when set, print some output on what ratpoints is doing
    27     # define RATPOINTS_FLAGS_INPUT_MASK \
    28     # (RATPOINTS_NO_CHECK | RATPOINTS_NO_Y | RATPOINTS_NO_REVERSE | \
    29     #  RATPOINTS_NO_JACOBI | RATPOINTS_VERBOSE)
    30 
    31 
    32     ctypedef struct ratpoints_interval:
    33         double low
    34         double up
    35     ctypedef struct ratpoints_args:
    36         mpz_t *cof
    37         long degree
    38         long height
    39         ratpoints_interval *domain
    40         long num_inter
    41         long b_low
    42         long b_high
    43         long sp1
    44         long sp2
    45         long array_size
    46         long sturm
    47         long num_primes
    48         long max_forbidden
    49         unsigned int flags
    50         # from here: private data
    51         # mpz_t *work
    52         # void *se_buffer
    53         # void *se_next
    54         # void *ba_buffer
    55         # void *ba_next
    56         # int *int_buffer
    57         # int *int_next
    58         # void *sieve_list
    59     long find_points(ratpoints_args*, int proc(long, long, mpz_t, void*, int*), void*)
    60     void find_points_init(ratpoints_args*)
    61     long find_points_work(ratpoints_args*, int proc(long, long, mpz_t, void*, int*), void*)
    62     void find_points_clear(ratpoints_args*)
    63 
    64 ctypedef struct point_list:
    65     long *xes
    66     mpz_t *ys
    67     long *zs
    68     long array_size
    69     long num_points
    70     long max_num_points
    71 
    726cdef int process(long x, long z, mpz_t y, void *info0, int *quit):
    737    # ratpoints calls this function when it finds a point [x : y : z]
    748    # info0 is the pointer passed to ratpoints originally
    … …  
    217151   
    218152    return L
    219153
     154cdef int process_exists_only(long x, long z, mpz_t y, void *info0, int *quit):
     155    cdef info_struct_exists_only *info_s = <info_struct_exists_only *>info0
     156    cdef Integer YY
     157    if info_s.verbose:
     158        YY = Integer(0); mpz_set(YY.value, y)
     159        print 'Found point [ %d : %d : %d ], quitting'%(x,YY,z)
     160    quit[0] = -1
     161    return 1
    220162
     163cdef bint ratpoints_mpz_exists_only(mpz_t *coeffs, long H, int degree, bint verbose):
     164    """
     165    Direct call to ratpoints to search for existence only.
     166   
     167    WARNING - The coefficient array will be modified by ratpoints.
     168    """
     169    cdef ratpoints_args args
     170    cdef info_struct_exists_only info_s
     171    cdef long total, verby = ~0 if verbose else 0
     172    info_s.verbose = verbose
     173    assert degree <= RATPOINTS_MAX_DEGREE
     174    args.degree = degree
     175    args.cof = coeffs
     176    args.domain = <ratpoints_interval *> sage_malloc(2*args.degree * sizeof(ratpoints_interval))
     177    args.height = H
     178    args.num_inter = 0
     179    args.b_low = 1
     180    args.b_high = H
     181    args.sp1 = RATPOINTS_DEFAULT_SP1
     182    args.sp2 = RATPOINTS_DEFAULT_SP2
     183    args.array_size = RATPOINTS_ARRAY_SIZE
     184    args.sturm = RATPOINTS_DEFAULT_STURM
     185    args.num_primes = RATPOINTS_DEFAULT_NUM_PRIMES
     186    args.max_forbidden = RATPOINTS_DEFAULT_MAX_FORBIDDEN
     187    args.flags = (RATPOINTS_VERBOSE & verby)
     188    total = find_points(&args, process_exists_only, <void *>(&info_s))
     189    sage_free(args.domain)
     190    if total == RATPOINTS_NON_SQUAREFREE:
     191        raise RuntimeError('Polynomial must be square-free')
     192    if total == RATPOINTS_BAD_ARGS:
     193        raise RuntimeError('Bad arguments to ratpoints')
     194    return (total > 0)
     195
     196
     197
     198

  • new file sage/schemes/elliptic_curves/descent_two_isogeny.pyx 

    diff --git a/sage/schemes/elliptic_curves/descent_two_isogeny.pyx b/sage/schemes/elliptic_curves/descent_two_isogeny.pyx
new file mode 100644
    - +  
     1"""
     2Descent 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
     12from sage.rings.all import ZZ
     13from sage.rings.polynomial.polynomial_ring import polygen
     14cdef object x_ZZ = polygen(ZZ)
     15from sage.rings.polynomial.real_roots import real_roots
     16from sage.rings.arith import prime_divisors
     17from sage.misc.all import walltime, cputime
     18from sage.libs.pari.gen import pari
     19
     20from sage.rings.integer cimport Integer
     21
     22include "../../ext/cdefs.pxi"
     23include "../../libs/flint/fmpz_poly.pxi"
     24
     25from sage.libs.flint.zmod_poly cimport *, zmod_poly_t
     26from sage.libs.flint.long_extras cimport *, factor_t
     27from sage.libs.ratpoints cimport ratpoints_mpz_exists_only
     28
     29cdef int N_RES_CLASSES_BSD = 10
     30
     31cdef unsigned long ui0 = <unsigned long>0
     32cdef unsigned long ui1 = <unsigned long>1
     33cdef unsigned long ui2 = <unsigned long>2
     34cdef unsigned long ui3 = <unsigned long>3
     35cdef unsigned long ui4 = <unsigned long>4
     36cdef unsigned long ui8 = <unsigned long>8
     37
     38cdef 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
     49def 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
     84cdef 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
     108def 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
     125cdef 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
     171cdef 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
     232cdef 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
     266cdef 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
     503cdef 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
     542cdef 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
     557cdef 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
     580def 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
     609cdef 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
     641def 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
     663cdef 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
     792def 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
     865def 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