Ticket #6583: trac_6583_large_primes.patch

File trac_6583_large_primes.patch, 12.7 KB (added by Robert Miller, 13 years ago)

  • sage/schemes/elliptic_curves/descent_two_isogeny.pyx 

    # HG changeset patch
# User Robert L. Miller <rlm@rlmiller.org>
# Date 1253307292 25200
# Node ID 32deafdfe447b6382167b38fe671672aafa59e54
# Parent  9fb66d817bf332aaf44297050575a08442d82509
Use the Python-based NTL interface to handle factoring modulo large primes.

diff -r 9fb66d817bf3 -r 32deafdfe447 sage/schemes/elliptic_curves/descent_two_isogeny.pyx
    a b  
    1616from sage.rings.arith import prime_divisors
    1717from sage.misc.all import walltime, cputime
    1818from sage.libs.pari.gen import pari
     19from sage.all import ntl
    1920
    2021from sage.rings.integer cimport Integer
    2122
    … …  
    500501            mpz_clear(ee)
    501502        return result
    502503
     504cdef bint Zp_soluble_siksek_large_p(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e, mpz_t pp,
     505                                    fmpz_poly_t f1, fmpz_poly_t linear):
     506    """
     507    Uses the approach of Algorithm 5.3.1 of Siksek's thesis to test for
     508    solubility of y^2 == ax^4 + bx^3 + cx^2 + dx + e over Zp.
     509    """   
     510    cdef unsigned long v_min, v
     511    cdef mpz_t roots[4]
     512    cdef int i, j, has_roots, has_single_roots
     513    cdef bint result
     514
     515    cdef mpz_t aa, bb, cc, dd, ee
     516    cdef mpz_t aaa, bbb, ccc, ddd, eee
     517    cdef mpz_t qq, rr, ss, tt
     518    cdef Integer A,B,C,D,E,P
     519
     520    # Step 0: divide out all common p from the quartic
     521    v_min = valuation(a, pp)
     522    if mpz_cmp_ui(b, ui0) != 0:
     523        v = valuation(b, pp)
     524        if v < v_min: v_min = v
     525    if mpz_cmp_ui(c, ui0) != 0:
     526        v = valuation(c, pp)
     527        if v < v_min: v_min = v
     528    if mpz_cmp_ui(d, ui0) != 0:
     529        v = valuation(d, pp)
     530        if v < v_min: v_min = v
     531    if mpz_cmp_ui(e, ui0) != 0:
     532        v = valuation(e, pp)
     533        if v < v_min: v_min = v
     534    for 0 <= v < v_min:
     535        mpz_divexact(a, a, pp)
     536        mpz_divexact(b, b, pp)
     537        mpz_divexact(c, c, pp)
     538        mpz_divexact(d, d, pp)
     539        mpz_divexact(e, e, pp)
     540
     541    if not v_min%2:
     542        # Step I in Alg. 5.3.1 of Siksek's thesis
     543        A = Integer(0); B = Integer(0); C = Integer(0); D = Integer(0); E = Integer(0); P = Integer(0)
     544        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)
     545        f = ntl.ZZ_pX([E,D,C,B,A], P)
     546        f /= ntl.ZZ_pX([A], P) # now f is monic, and we are done with A,B,C,D,E
     547        mpz_set(qq, A.value) # qq is the leading coefficient of the polynomial
     548        f_factzn = f.factor()
     549        result = 0
     550        for factor, exponent in f_factzn:
     551            if exponent&1:
     552                result = 1
     553                break
     554        if result == 0 and mpz_legendre(qq, pp) == 1:
     555            result = 1
     556        if result:
     557            return 1
     558
     559        f = ntl.ZZ_pX([1], P)
     560        for factor, exponent in f_factzn:
     561            for j from 0 <= j < (exponent/2):
     562                f *= factor
     563
     564        f /= f.leading_coefficient()
     565        f_factzn = f.factor()
     566
     567        has_roots = 0
     568        j = 0
     569        for factor, exponent in f_factzn:
     570            if factor.degree() == 1:
     571                has_roots = 1
     572                A = P - Integer(factor[0])
     573                mpz_set(roots[j], A.value)
     574                j += 1
     575        if not has_roots:
     576            return 0
     577       
     578        i = f.degree()
     579        mpz_init(aaa)
     580        mpz_init(bbb)
     581        mpz_init(ccc)
     582        mpz_init(ddd)
     583        mpz_init(eee)
     584
     585        if i == 0: # g == 1
     586            mpz_set(aaa, a)
     587            mpz_set(bbb, b)
     588            mpz_set(ccc, c)
     589            mpz_set(ddd, d)
     590            mpz_sub(eee, e, qq)
     591        elif i == 1: # g == x + rr
     592            mpz_set(aaa, a)
     593            mpz_set(bbb, b)
     594            mpz_sub(ccc, c, qq)
     595            A = Integer(f[0])
     596            mpz_set(rr, A.value)
     597            mpz_mul(ss, rr, qq)
     598            mpz_set(ddd,d)
     599            mpz_sub(ddd, ddd, ss)
     600            mpz_sub(ddd, ddd, ss)
     601            mpz_set(eee,e)
     602            mpz_mul(ss, ss, rr)
     603            mpz_sub(eee, eee, ss)
     604            mpz_divexact(ss, ss, rr)
     605        elif i == 2: # g == x^2 + rr*x + ss
     606            mpz_sub(aaa, a, qq)
     607            A = Integer(f[1])
     608            mpz_set(rr, A.value)
     609            mpz_init(tt)
     610            mpz_mul(tt, rr, qq)
     611            mpz_set(bbb,b)
     612            mpz_submul_ui(bbb, tt, ui2)
     613            mpz_set(ccc,c)
     614            mpz_submul(ccc, tt, rr)
     615            A = Integer(f[0])
     616            mpz_set(ss, A.value)           
     617            mpz_mul(tt, ss, qq)
     618            mpz_set(eee,e)
     619            mpz_submul(eee, tt, ss)
     620            mpz_mul_ui(tt, tt, ui2)
     621            mpz_sub(ccc, ccc, tt)
     622            mpz_set(ddd,d)
     623            mpz_submul(ddd, tt, rr)
     624            mpz_clear(tt)
     625        mpz_divexact(aaa, aaa, pp)
     626        mpz_divexact(bbb, bbb, pp)
     627        mpz_divexact(ccc, ccc, pp)
     628        mpz_divexact(ddd, ddd, pp)
     629        mpz_divexact(eee, eee, pp)
     630        # now aaa,bbb,ccc,ddd,eee represents h(x)
     631
     632        result = 0
     633        mpz_init(tt)
     634        for i from 0 <= i < j:
     635            mpz_mul(tt, aaa, roots[i])
     636            mpz_add(tt, tt, bbb)
     637            mpz_mul(tt, tt, roots[i])
     638            mpz_add(tt, tt, ccc)
     639            mpz_mul(tt, tt, roots[i])
     640            mpz_add(tt, tt, ddd)
     641            mpz_mul(tt, tt, roots[i])
     642            mpz_add(tt, tt, eee)
     643            # tt == h(r) mod p
     644            mpz_mod(tt, tt, pp)
     645            if mpz_sgn(tt) == 0:
     646                fmpz_poly_zero(f1)
     647                fmpz_poly_zero(linear)
     648                fmpz_poly_set_coeff_mpz(f1, 0, e)
     649                fmpz_poly_set_coeff_mpz(f1, 1, d)
     650                fmpz_poly_set_coeff_mpz(f1, 2, c)
     651                fmpz_poly_set_coeff_mpz(f1, 3, b)
     652                fmpz_poly_set_coeff_mpz(f1, 4, a)
     653                fmpz_poly_set_coeff_mpz(linear, 0, roots[i])
     654                fmpz_poly_set_coeff_mpz(linear, 1, pp)
     655                fmpz_poly_compose(f1, f1, linear)
     656                fmpz_poly_scalar_div_mpz(f1, f1, pp)
     657                fmpz_poly_scalar_div_mpz(f1, f1, pp)
     658
     659                mpz_init(aa)
     660                mpz_init(bb)
     661                mpz_init(cc)
     662                mpz_init(dd)
     663                mpz_init(ee)
     664                fmpz_poly_get_coeff_mpz(aa, f1, 4)
     665                fmpz_poly_get_coeff_mpz(bb, f1, 3)
     666                fmpz_poly_get_coeff_mpz(cc, f1, 2)
     667                fmpz_poly_get_coeff_mpz(dd, f1, 1)
     668                fmpz_poly_get_coeff_mpz(ee, f1, 0)
     669                result = Zp_soluble_siksek_large_p(aa, bb, cc, dd, ee, pp, f1, linear)
     670                mpz_clear(aa)
     671                mpz_clear(bb)
     672                mpz_clear(cc)
     673                mpz_clear(dd)
     674                mpz_clear(ee)
     675                if result == 1:
     676                    break
     677        mpz_clear(aaa)
     678        mpz_clear(bbb)
     679        mpz_clear(ccc)
     680        mpz_clear(ddd)
     681        mpz_clear(eee)
     682        mpz_clear(tt)
     683        return result
     684    else:
     685        # Step II in Alg. 5.3.1 of Siksek's thesis
     686        A = Integer(0); B = Integer(0); C = Integer(0); D = Integer(0); E = Integer(0); P = Integer(0)
     687        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)
     688        f = ntl.ZZ_pX([E,D,C,B,A], P)
     689        f /= ntl.ZZ_pX([A], P) # now f is monic
     690        f_factzn = f.factor()
     691
     692        has_roots = 0
     693        has_single_roots = 0
     694        j = 0
     695        for factor, exponent in f_factzn:
     696            if factor.degree() == 1:
     697                has_roots = 1
     698                if exponent == 1:
     699                    has_single_roots = 1
     700                    break
     701                A = P - Integer(factor[0])
     702                mpz_set(roots[j], A.value)
     703                j += 1
     704
     705        if not has_roots: return 0
     706        if has_single_roots: return 1
     707
     708        result = 0
     709        if j > 0:
     710            mpz_init(aa)
     711            mpz_init(bb)
     712            mpz_init(cc)
     713            mpz_init(dd)
     714            mpz_init(ee)       
     715        for i from 0 <= i < j:
     716            fmpz_poly_zero(f1)
     717            fmpz_poly_zero(linear)
     718            fmpz_poly_set_coeff_mpz(f1, 0, e)
     719            fmpz_poly_set_coeff_mpz(f1, 1, d)
     720            fmpz_poly_set_coeff_mpz(f1, 2, c)
     721            fmpz_poly_set_coeff_mpz(f1, 3, b)
     722            fmpz_poly_set_coeff_mpz(f1, 4, a)
     723            fmpz_poly_set_coeff_mpz(linear, 0, roots[i])
     724            fmpz_poly_set_coeff_mpz(linear, 1, pp)
     725            fmpz_poly_compose(f1, f1, linear)
     726            fmpz_poly_scalar_div_mpz(f1, f1, pp)
     727            fmpz_poly_get_coeff_mpz(aa, f1, 4)
     728            fmpz_poly_get_coeff_mpz(bb, f1, 3)
     729            fmpz_poly_get_coeff_mpz(cc, f1, 2)
     730            fmpz_poly_get_coeff_mpz(dd, f1, 1)
     731            fmpz_poly_get_coeff_mpz(ee, f1, 0)
     732            result = Zp_soluble_siksek_large_p(aa, bb, cc, dd, ee, pp, f1, linear)
     733            if result == 1:
     734                break
     735        if j > 0:
     736            mpz_clear(aa)
     737            mpz_clear(bb)
     738            mpz_clear(cc)
     739            mpz_clear(dd)
     740            mpz_clear(ee)
     741        return result
     742
    503743cdef bint Qp_soluble_siksek(mpz_t A, mpz_t B, mpz_t C, mpz_t D, mpz_t E,
    504744                            mpz_t p, unsigned long P,
    505745                            zmod_poly_factor_t f_factzn, fmpz_poly_t f1,
    … …  
    539779    zmod_poly_clear(f)
    540780    return result
    541781
     782cdef bint Qp_soluble_siksek_large_p(mpz_t A, mpz_t B, mpz_t C, mpz_t D, mpz_t E,
     783                                    mpz_t p, fmpz_poly_t f1, fmpz_poly_t linear):
     784    """
     785    Uses Samir Siksek's thesis results to determine whether the quartic is
     786    locally soluble at p, when p is bigger than wordsize, and we can't use
     787    FLINT.
     788    """
     789    cdef int result = 0
     790    cdef mpz_t a,b,c,d,e
     791
     792    mpz_init_set(a,A)
     793    mpz_init_set(b,B)
     794    mpz_init_set(c,C)
     795    mpz_init_set(d,D)
     796    mpz_init_set(e,E)
     797
     798    if Zp_soluble_siksek_large_p(a,b,c,d,e,p,f1,linear):
     799        result = 1
     800    else:
     801        mpz_set(a,A)
     802        mpz_set(b,B)
     803        mpz_set(c,C)
     804        mpz_set(d,D)
     805        mpz_set(e,E)
     806        if Zp_soluble_siksek_large_p(e,d,c,b,a,p,f1,linear):
     807            result = 1
     808
     809    mpz_clear(a)
     810    mpz_clear(b)
     811    mpz_clear(c)
     812    mpz_clear(d)
     813    mpz_clear(e)
     814    return result
     815
    542816cdef bint Qp_soluble_BSD(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e, mpz_t p):
    543817    """
    544818    Uses the original test of Birch and Swinnerton-Dyer to test for local
    … …  
    563837    cdef unsigned long pp
    564838    cdef fmpz_poly_t f1, linear
    565839    cdef zmod_poly_factor_t f_factzn
    566     fmpz_poly_init(f1)
    567     fmpz_poly_init(linear)
    568     zmod_poly_factor_init(f_factzn)
    569840    bsd_sol = Qp_soluble_BSD(a, b, c, d, e, p)
    570841    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)
     842        fmpz_poly_init(f1)
     843        fmpz_poly_init(linear)
     844        if mpz_fits_ulong_p(p):
     845            zmod_poly_factor_init(f_factzn)
     846            pp = mpz_get_ui(p)
     847            sik_sol = Qp_soluble_siksek(a,b,c,d,e,p,pp,f_factzn,f1,linear)
     848            zmod_poly_factor_clear(f_factzn)
     849        else:
     850            sik_sol = Qp_soluble_siksek_large_p(a,b,c,d,e,p,f1,linear)
     851        fmpz_poly_clear(f1)
     852        fmpz_poly_clear(linear)
    573853    else:
    574854        sik_sol = bsd_sol
    575     fmpz_poly_clear(f1)
    576     fmpz_poly_clear(linear)
    577     zmod_poly_factor_clear(f_factzn)
    578855    return sik_sol
    579856
    580857def test_qpls(a,b,c,d,e,p):
    … …  
    632909    Delta = f.discriminant()
    633910    for p in prime_divisors(Delta):
    634911        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.')
    637912        if not Qp_soluble(a,b,c,d,e,p.value): return 0
    638913
    639914    return 1
    … …  
    8201095        sage: log(n1,2) + log(n1_prime,2) - 2 # the rank
    8211096        3
    8221097
    823     ::
     1098    Using the verbosity option::
    8241099
    8251100        sage: E = EllipticCurve('14a')
    8261101        sage: two_descent_by_two_isogeny(E, verbosity=1)
    … …  
    8361111        0 <= rank of E(Q) = rank of E'(Q) <= 0
    8371112        (2, 2, 2, 2)
    8381113
     1114    Handling curves whose discriminants involve larger than wordsize primes::
     1115   
     1116        sage: E = EllipticCurve('14a')
     1117        sage: E = E.quadratic_twist(next_prime(10^20))
     1118        sage: E
     1119        Elliptic Curve defined by y^2 = x^3 + x^2 + 716666666666666667225666666666666666775672*x - 391925925925925926384240370370370370549019837037037037060249356 over Rational Field
     1120        sage: E.discriminant().factor()
     1121        -1 * 2^18 * 7^3 * 100000000000000000039^6
     1122        sage: log(100000000000000000039.0, 2.0)
     1123        66.438...
     1124        sage: n1, n2, n1_prime, n2_prime = two_descent_by_two_isogeny(E)
     1125        sage: log(n1,2) + log(n1_prime,2) - 2 # the rank
     1126        0
    8391127
    8401128    """
    8411129    cdef Integer a1, a2, a3, a4, a6, s2, s4, s6