divisibility_reductions/jumps.sage


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145

from sage.schemes.elliptic_curves.ell_generic import is_EllipticCurve           
                                                                                
# Pre-tests on E and P                                                          
def suitable( E, P, ell ):                                                      
    if not is_EllipticCurve(E):                                                 
        print "E is not an elliptic curve"                                      
        return False                                                            
    if E.base_field() != QQ:                                                    
        print "E is not defined over Q"                                         
        return False                                                            
    if not P in E:                                                              
        print "P is not in E"                                                   
        return False                                                            
    if P.has_finite_order():                                                    
        print "P has finite order"                                              
        return False                                                            
    if not is_prime(ell):                                                       
        print "ell is not prime"                                                
        return False                                                            
    return True                                                                 
                                                                                
                                                                                
# Stupid auxiliary function. Returns higest power of n dividing m.              
def val( n, m ):                                                                
    ret = 0                                                                     
    while m % (n^(ret+1)) == 0:                                                 
        ret += 1                                                                
    return ret

def test_jump_kl2_label( label, ell ):
    E = EllipticCurve( label )
    if E.rank() < 1:
        print "E has rank 0"
        return
    P = E.gens()[0]
    test_jump_kl2( E, P, ell )

def test_jump_kl2( E, P, ell ):
    if not suitable( E, P, ell ):
        return

    flag = True
    for p in Primes():
        if p > 10^3:
            break
        if p == ell:                                                            
            # print "Skipping", p, "because = ell"                                
            continue                                                            
        if E.discriminant() % p == 0:                                           
            # print "Skipping", p, "because E has bad reduction"                  
            continue                                                            
        if P.reduction(p).order() % ell != 0:                                   
            # print "Skipping", p, "because red of P is infinitely ell-divisible" 
            continue

        # print "Working with prime p =", p

        F = FiniteField(p)

        for i in range(3):   
            # print "Torsion level", i, "..."                                     
                                                                                 
            # We can build the division fields for increasing powers of l       
            # incrementally. To get a division field, we first compute the      
            # splitting field of the division polynomial. We may be off by      
            # a degree 2 extension. If so, by finite field magic we know        
            # exactly which degree 2 extension we need: the unique one!         
            R.<x> = PolynomialRing(F)                                           
            F.<a> = E.reduction(p).division_polynomial(ell^i).splitting_field() 
            E_red = E.reduction(p).base_extend(F)                               
            # extend if necessary                                               
            k = val(ell,E_red.gens()[0].order())                                
            h = val(ell,E_red.order()) - k                                      
            if k < i or h < i:                                                  
                F.<a> = F.extension(2)                                          
            E_red = E_red.base_extend(F)                                        
                                                                                
            P_red = E_red.point(P.reduction(p))

            flag = False
            if P_red.is_divisible_by(ell):
                # print "P ell-divisible in this torsion level"
                flag = True
                break

        if not flag:
            print "Point not divisible mod", p
            print "Stopping here"
            break
    if flag:
        print "*********************************"
        print "*** Candidate counterexample! ***"
        print "*********************************"

# Wrapper                                                
def test_jump_den_label( label, ell ):                                              
    E = EllipticCurve(label)                                                    
    if E.rank() < 1:                                                            
        print "E has rank 0"                                                    
        return                                                                  
    P = E.gens()[0]                                                             
    test_jump_den( E, P, ell ) 

def test_jump_den( E, P, ell ):
    if not suitable( E, P, ell ):
        return

    n_primes = 0
    inf_divisible = 0
    tot_l_part = 0
    tot_non_l_part = 0

    for p in Primes():
        if p > 5*10^3:
            break
        if p == ell or E.discriminant() % p == 0:
            continue

        n_primes += 1

        E_red = E.reduction(p)
        N = E_red.order()                                                       
        k = val(ell,E_red.gens()[0].order())                                    
        h = val(ell,N) - k                                                      
        temp_non_l_part = (N / (ell^(h+k)))-1                                   
        tot_non_l_part += temp_non_l_part                                       
        temp_l_part = N - temp_non_l_part
        tot_l_part += temp_l_part

        if P.reduction(p).order() % ell != 0:
            inf_divisible += 1

    found = RDF( inf_divisible / n_primes )
    expected = RDF( tot_non_l_part / (tot_l_part+tot_non_l_part) )
    print "Found: %0.3f, expected: %0.3f"%(found, expected)
    if abs(found-expected) > 0.05:
        print "Unexpected density! Checking divisibility in reductions..."
        test_jump_kl2( E, P, ell )

def test_from_file( filename ):
    attach(filename)
    for coord in data:
        label = EllipticCurve(coord).label()
        print "Trying curve", label, "with ell =", 3
        test_jump_den_label(label,3)