Attachment 'three_descent.sage'
Download 1 R.<x> = QQ[]
2 K.<w> = NumberField(x^2 - x + 1)
3
4 def norm_in_K(x):
5 if x in ZZ:
6 return x*x
7 else:
8 return x.norm()
9
10 def get_primary_prime(p):
11 if p%3 == 2:
12 fp = p
13 elif p == 3:
14 fp = 1 + w
15 else:
16 R.<x> = GF(p)[]
17 fp = (x^2 - x + 1).roots()[1][0].lift()
18 fp = find_gcd(p, fp-w)
19 fp = K(fp)
20 for u in [1, -1, w, -w, w^2, -w^2]:
21 if (fp*u)[0] % 3 == 2 and (fp*u)[1] % 3 == 0:
22 return fp*u
23 assert False
24
25 def is_in_ring_of_int(a):
26 return a in a.parent().maximal_order()
27
28 def find_gcd(f1, f2):
29 g1 = find_primary(f1)
30 g2 = find_primary(f2)
31 if is_in_ring_of_int(g1/g2): return g2
32 if is_in_ring_of_int(g2/g1): return g1
33 if norm_in_K(g1) >= norm_in_K(g2):
34 return find_gcd(g2, find_primary(g1-g2))
35 else:
36 return find_gcd(g1, find_primary(g2-g1))
37
38 def find_primary(inp):
39 nm = norm_in_K(inp)
40 if nm == 1:
41 return inp
42 while nm/3 in ZZ:
43 inp = inp/(1+w)
44 nm = nm/3
45 if is_in_ring_of_int((inp - 1)/3): return inp
46 if is_in_ring_of_int((-inp - 1)/3): return -inp
47 if is_in_ring_of_int((inp*w - 1)/3): return inp*w
48 if is_in_ring_of_int((-inp*w-1)/3): return -inp*w
49 if is_in_ring_of_int((inp*w^2-1)/3): return inp*w^2
50 if is_in_ring_of_int((-inp*w^2-1)/3): return -inp*w^2
51 return inp
52
53
54 def p_reduction(p, ell, n):
55 # reduction by p, if p|l, p^3|n, then return reduction(l/p,n/p^3)
56 e = min(ell.valuation(p), n.valuation(p)//3)
57 if e == 0:
58 return ell, n
59 else:
60 return ell//p**e, n//p**(3*e)
61
62 def mconcat(v):
63 return "".join(map(str, v))
64
65 def eqn_print(coef, vars, is_first=False):
66 if coef == 0:
67 return ""
68 if coef < 0:
69 prefix = "-"
70 coef = -coef
71 else:
72 prefix = "" if is_first else "+"
73 if vars == "":
74 return prefix + str(coef)
75 elif coef == 1:
76 return prefix + vars
77 else:
78 return prefix + str(coef) + "*" + vars
79
80 def cubic_residue(p, q):
81 """
82 sage: for p in primes(20, 50):
83 ... print p
84 ... for a in range(1, p):
85 ... if cubic_residue(p, a^2) != 2*cubic_residue(p, a):
86 ... print "\t", a, cubic_residue(p, a^2), 2*cubic_residue(p, a)
87 ... assert cubic_residue(p, a^3) == 0
88 23
89 29
90 31
91 37
92 41
93 43
94 47
95
96 sage: for p in primes(10, 50):
97 ... print p, get_primary_prime(p) == gp.GetPrimaryPrime(p)
98 ... for a in range(1,p):
99 ... if cubic_residue(p, a) != gp.CubicResidue(p, a):
100 ... print "\t", a, cubic_residue(p, a), gp.CubicResidue(p, a)
101
102 """
103 if is_prime(p) and p%3 == 2 and q%p != 0:
104 return mod(0, 3)
105 if is_prime(p) and p%3 == 1 and q%p != 0:
106 cr = mod(q, p)^((p-1)/3)
107 if cr == 1:
108 return mod(0, 3)
109 # pp=pa+pb*w^2 ( w^2 = (-1+sqrt(-3))/2 )
110 pp = get_primary_prime(p)
111 pa = pp[0] + pp[1]
112 pb = pp[1]
113 if cr == mod(-pa/pb, p):
114 return mod(1, 3)
115 else:
116 return mod(2, 3)
117 elif p == 9:
118 mq = mod(q, p)
119 if mq in [1,8]:
120 return mod(0, 3)
121 elif mq in [2, 7]:
122 return mod(2, 3)
123 else:
124 return mod(1, 3)
125 else:
126 raise ValueError, "p = %s, q = %s" % (p,q)
127
128 def three_isogeny_descent(A1, A3, upto=10, verbose = -1):
129 #local(da1,da3,cf,fda3,nda3,a1,a3,fac,p,p1,pv,qv,pred,mpow,clist,rlist,m9a1,m9a3,twist,RET,KERNEL,cr,ci,tmp,rcn,rrn,mapto)
130
131 rrn = rcn = 0
132
133 if A3 * (A1^3 - 27*A3) == 0:
134 raise ValueError, "Singular curve in ThreeIsogenyDescent"
135 t = cputime()
136
137 da1 = A1.denominator()
138 da3 = A3.denominator()
139
140 if A1 == 1 or A3==1: # Changed to upper case. ?
141 cf = 1
142 else:
143 fda3 = factor(da3)
144 nda3 = 1
145 for k in range(1, len(fda3)+1):
146 nda3 *= fda3[k][0]**(floor((fda3[k][1] - 1)/3) + 1)
147 cf = lcm(da1,nda3)
148
149 a1 = A1*cf
150 a3 = A3*cf^3
151
152 fac = factor(a3.abs())
153
154 clist=[]
155 rlist=[]
156
157 pred=[a1,a3]
158
159 for i in range(len(fac)):
160 p = fac[i][0]
161 pred = p_reduction(p, pred[0], pred[1])
162
163 mpow = pred[1].valuation(p)
164 if mpow > 0:
165 clist.append((p, mpow))
166
167 if verbose>0:
168 print "clist: ", clist
169 print "reduced [a1,a3]: ", pred
170
171 fac = factor((pred[0]**3-27*pred[1]).abs())
172
173 for i in range(len(fac)):
174 p = fac[i][0]
175 if mod(p, 3) == 1:
176 p1 = get_primary_prime(p)
177 L = p1.galois_conjugates(p1.parent())
178 pconj = L[0] if L[0] != p1 else L[1]
179 rlist.append((p, p1/pconj))
180
181 if verbose > 0:
182 print "rlist: ", rlist
183
184 m9a1 = mod(pred[0], 9)
185 m9a3 = mod(pred[1], 9);
186
187 if (m9a1==0 and m9a3==0) \
188 or (m9a1==0 and m9a3==1) \
189 or (m9a1==0 and m9a3==8) \
190 or (m9a1==3 and m9a3==6) \
191 or (m9a1==6 and m9a3==3) \
192 or (m9a1==3 and m9a3==8) \
193 or (m9a1==6 and m9a3==1) \
194 or (m9a1==3 and m9a3==1 and mod(a1-a3, 27) != 11) \
195 or (m9a1==6 and m9a3==8 and mod(a1-a3, 27) != 16):
196 rlist.append((9,1-w))
197
198 RET = matrix(GF(3), len(rlist), len(clist))
199
200 for i in range(len(rlist)):
201 for j in range(len(clist)):
202 pv = rlist[i]
203 qv = clist[j]
204
205 if mod(pv[0], qv[0]) != 0:
206 RET[i, j] = cubic_residue(pv[0], qv[0])
207 else:
208 RET[i, j] = 2*qv[1]*cubic_residue(pv[0], pred[1]/qv[0]**qv[1])
209 if verbose>1:
210 print "p = %s, q = %s, CubicResidue = %s"%(pv[1], qv[1], RET[i,j])
211
212 if verbose>0:
213 print "Selmer matrix:"
214 print RET
215
216 # when a1^3 - 27*a3 has no prime divisors, we use the 1 by len(clist) zero matrix to get homogeneous spaces
217 if (pred[0]**3 - 27*pred[1]).abs() == 1:
218 RET=matrix(GF(3), 1,len(clist))
219
220 KERNEL = RET.right_kernel().basis_matrix()
221
222 count = 0
223 rcn = KERNEL.nrows()
224
225 for idx in GF(3)**rcn:
226 if count > upto:
227 break
228 count += 1
229 i = 0
230 while i < rcn and idx[i] == 1:
231 i += 1
232
233
234 if i < rcn and idx[i] == 0:
235 kerelmt = (idx - vector([1]*rcn)) * KERNEL
236 twist = 1
237 for j in range(rcn):
238 twist *= clist[j][0]**(kerelmt[j].lift())
239 twist = K(twist)
240 print "=== Curve Over Rational Field Defined By ==="
241 print "-u^3", eqn_print(twist^2, "v^3"), eqn_print(A1*twist, "u*v"), eqn_print(A3*twist, ""), ","
242 print " the Map to the Curve Defined by y^2", eqn_print(A1, "x*y"), eqn_print(A3, "y"), "-x^3"
243 print " (u,v) -> (u*v,", "0" if twist == 0 else eqn_print(twist, "v^3", 1), ")"
244 print ""
245
246 # when a3 has no prime divisors, we use the 1 by len(rlist) zero matrix to get homogeneous spaces
247 if pred[1].abs() == 1:
248 RET = matrix(GF(3), len(rlist), 1)
249 if (pred[0]**3 - 27*pred[1]).abs() != 1:
250 KERNEL=RET.transpose().right_kernel().basis_matrix()
251
252 count=0;
253 rrn=KERNEL.nrows()
254
255 for idx in GF(3)**rrn:
256 if count > upto:
257 break
258 count += 1
259 i = 0
260 while i < rrn and idx[i]==1:
261 i += 1
262
263 if i < rrn and idx[i] == 0:
264 kerelmt = (idx - vector([1]*rrn)) * KERNEL
265 twist = 1
266 for j in range(rrn):
267 twist *= rlist[j][1]**(kerelmt[j].lift())
268 twist = K(twist)
269
270 print("=== Curve Over Rational Field Defined By ");
271
272 # twist = n + mw where w^2 - w + 1 = 0, so let w = (1+\sqrt{-3})/2
273 # then cr = n + m/2, ci = m/2
274
275 cr = twist[0] + twist[1]/2
276 ci = twist[1]/2
277
278 tmp = -(A1**3 - 27*A3)
279 print eqn_print(ci, "(u^3-9*u*v^2)", 1), eqn_print(3*cr, "(u^2*v-v^3)"), eqn_print(3/2*A1, "(u^2+3*v^2)"), eqn_print(3/2*tmp,""), ","
280 print " the Map to the Curve Defined By y^2 - x^3 + 27/4*(", eqn_print(A1, "x", 1), eqn_print(tmp, "", 1) if A1 == 0 else 0, ")^2"
281 mapto = eqn_print(cr, "(u^3-9*u*v^2)", 1) + eqn_print(-9*ci, "(u^2*v-v^3)", 1 if cr==0 else 0)
282 print " (u,v) -> (u^2+3*v^2,", mapto if mapto != "" else "0", ")"
283 print ""
284
285 if verbose > 0:
286 print "Time Elapsed: ", cputime(t)
287 print "Upper Bound of Rank: ", rcn + rrn - 1
288 print "Selmer matrix: ", RET
289 return rcn+rrn-1
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.