Curve_Finding_Methods

%cython ################################################ # This is our code to run the modjon method, # which is a naive search aided by congruence # conditions. We define a class (SievedSearch) # whose __call__ method executes the method. ################################################ from libc.stdint cimport uint8_t, int32_t, uint32_t, int64_t, uint64_t from psage.libs.smalljac.wrapper1 import elliptic_curve_ap from sage.rings.integer cimport Integer from sage.schemes.elliptic_curves.constructor import EllipticCurve from psage.modform.hilbert.sqrt5.sqrt5 import F as K cdef extern from "pari/pari.h": uint32_t Fl_sqrt(uint32_t, uint32_t) ctypedef uint32_t pr_cur[2] ctypedef uint32_t uint_cur[4] ctypedef int32_t int_cur[10] def gal_con(alpha): return alpha[0]+alpha[1]*(1-K.gen(0)) cdef class SievedSearch: cdef uint32_t a11,a12,a31,a32,A1,A2,B1,B2,p,numPrimes,L cdef uint32_t primes[5],lens[5],liftCoef[2] cdef public uint32_t r1,r2 cdef int64_t a21,a22,a41,a42,a61,a62 cdef int_cur var_lifts[2304] cdef int_cur *cur_lift cdef uint32_t q,v,numFirst,numSecond cdef mpz_t s,t,u cdef pr_cur *firstp, *secondp def __init__(self, *args): for i,p in enumerate(args): self.primes[i] = p self.numPrimes = len(args) self.p = 1 for i in range(self.numPrimes): self.p *= self.primes[i] self.q = self.p+1u>>1u self.v = (self.p+1u)/3u if self.p%3u==2u else ((self.p<<1u)+1u)/3u cpdef get_primes(self,uint32_t p): cdef uint32_t g = Fl_sqrt(5u, p)+1u cdef uint32_t w = (<uint64_t>(p+1u)>>1u)*(<uint64_t>g)%p g = p+1u-w if w > g: self.r1=w self.r2=g else: self.r1=g self.r2=w cdef uint_cur *split_prime_combine(self, uint32_t i): cdef uint32_t p = self.primes[i] cdef uint32_t r1,r2,diff_inv self.get_primes(p) r1 = self.r1 r2 = self.r2 diff_inv = self.inv_mod(r1-r2,p) cdef uint_cur *ret ret = <uint_cur *>sage_malloc(sizeof(uint_cur)*self.lens[i]) cdef uint32_t *E, *F, *lift cdef pr_cur *Eptr = self.firstp cdef pr_cur *Fptr cdef uint_cur *liftprt liftptr = ret cdef uint32_t j,k for j in range(self.numFirst): E = Eptr[0] Fptr = self.secondp for k in range(self.numSecond): F = Fptr[0] lift = liftptr[0] self.custom_crt(E[0],F[0],r1,diff_inv,p) lift[0] = self.liftCoef[0] lift[1] = self.liftCoef[1] self.custom_crt(E[1],F[1],r1,diff_inv,p) lift[2] = self.liftCoef[0] lift[3] = self.liftCoef[1] liftptr += 1 Fptr += 1 Eptr += 1 return ret cdef uint_cur *integral_combine(self, uint_cur *first_list, uint_cur *second_list, uint32_t i): cdef uint32_t oldL = self.L self.L *= self.lens[i] cdef uint_cur *ret = <uint_cur *>sage_malloc(sizeof(uint_cur)*self.L) cdef uint32_t j,k,oldpMod,newp if not i: for j in range(self.L): ret[j][0] = second_list[j][0] ret[j][1] = second_list[j][1] ret[j][2] = second_list[j][2] ret[j][3] = second_list[j][3] return ret cdef uint32_t oldp = 1u cdef uint_cur *Eptr,*Fptr,*liftptr cdef uint32_t *E,*F,*lift for j in range(i): oldp *= self.primes[j] newp = oldp*self.primes[i] oldpMod = self.inv_mod(oldp,self.primes[i]) oldpMod *= oldp Eptr = first_list liftptr = ret for j in range(oldL): E = Eptr[0] Fptr = second_list for k in range(self.lens[i]): F = Fptr[0] lift = liftptr[0] lift[0] = self.integral_crt(E[0],F[0],oldpMod,newp) lift[1] = self.integral_crt(E[1],F[1],oldpMod,newp) lift[2] = self.integral_crt(E[2],F[2],oldpMod,newp) lift[3] = self.integral_crt(E[3],F[3],oldpMod,newp) liftptr += 1 Fptr += 1 Eptr += 1 return ret cdef uint32_t integral_crt(self, uint32_t s, uint32_t t, uint32_t pMod, uint32_t newp): return ((newp+t-s)*pMod+s)%newp cdef custom_crt(self, uint32_t s, uint32_t t, uint32_t r1, uint32_t diff_inv, uint32_t p): self.liftCoef[1] = (p+s-t)*diff_inv%p self.liftCoef[0] = (p+s-self.liftCoef[1]*r1%p)%p cdef uint32_t inv_mod(self, uint32_t x, uint32_t p): cdef uint32_t t,u1,u2,v1,v2 u1=x v1=1 v2=(p-p/x)%p u2=v2*x%p while u2>1: t=u1/u2 v1,v2=v2,(p+v1-t*v2%p)%p u1,u2=u2,(p+u1-t*u2%p)%p return v2 cdef pr_cur *curves_in_GF(self, uint32_t which, int32_t ap, uint32_t k): cdef uint32_t p = self.primes[k] cdef pr_cur *curves = <pr_cur *>sage_malloc(sizeof(pr_cur)<<4) cdef uint32_t cur_len = 16u cdef uint32_t ret = 0u cdef uint32_t i,j for i in range(p): for j in range(p): if 64u*i*i*i + 432u*j*j: if elliptic_curve_ap(i, j, p) == ap: if ret == cur_len: cur_len <<= 1u curves = <pr_cur *>sage_realloc(curves, sizeof(pr_cur)*cur_len) curves[ret][0] = i curves[ret][1] = j ret += 1u if which == 1u: self.numFirst = ret else: self.numSecond = ret return curves cdef do_clift(self): self.cur_lift[0][0] = self.a11 self.cur_lift[0][1] = self.a12 self.cur_lift[0][2] = self.a21 self.cur_lift[0][3] = self.a22 self.cur_lift[0][4] = self.a31 self.cur_lift[0][5] = self.a32 self.cur_lift[0][6] = self.a41 self.cur_lift[0][7] = self.a42 self.cur_lift[0][8] = self.a61 self.cur_lift[0][9] = self.a62 self.cur_lift += 1 cdef do_a4lift(self): cdef int64_t old1,old2 old1 = self.a61 old2 = self.a62 self.do_clift()#0 self.a61-=self.p self.do_clift()#1 self.a62-=self.p self.do_clift()#3 self.a61+=self.p self.do_clift()#2 self.a61 = old1 self.a62 = old2 cdef do_t(self): cdef uint64_t s1 = self.q if self.a11 else 0u cdef uint64_t s2 = self.q if self.a12 else 0u cdef uint64_t t1 = self.q if self.a31 else 0u cdef uint64_t t2 = self.q if self.a32 else 0u cdef uint64_t r1,r2 cdef uint64_t xsqr1 = s1*s1+s2*s2 cdef uint64_t xsqr2 = s2*((s1<<1u)+s2) r1 = self.v*(self.p+xsqr1+self.a21)%self.p r2 = self.v*(xsqr2+self.a22)%self.p cdef uint64_t rpow1 = (r1*r1+r2*r2)%self.p cdef uint64_t rpow2 = r2*((r1<<1u)+r2)%self.p cdef uint64_t temp1 = s2*t2%self.p cdef uint64_t temp2 = (s1*t2+s2*t1+temp1)%self.p temp1 += s1*t1 self.a41 = (self.p+self.A1+3*rpow1-(temp1<<1u)%self.p)%self.p self.a42 = (self.p+self.A2+3*rpow2-(temp2<<1u)%self.p)%self.p temp1 = rpow2*r2 rpow1, rpow2 = temp1+rpow1*r1, temp1+rpow1*r2+rpow2*r1 temp1 = self.A2*r2 temp2 = self.A1*r2+self.A2*r1+temp1 temp1 += self.A1*r1 xsqr1 = (t1*t1+t2*t2)%self.p xsqr2 = t2*((t1<<1u)+t2)%self.p self.a61 = (self.p+self.B1+temp1+rpow1-xsqr1)%self.p self.a62 = (self.p+self.B2+temp2+rpow2-xsqr2)%self.p self.do_a4lift()#0 self.a41-=self.p self.do_a4lift()#1 self.a42-=self.p self.do_a4lift()#3 self.a41+=self.p self.do_a4lift()#2 cdef do_r(self): self.a31 = 0u self.a32 = 0u self.do_t() self.a31 = 1u self.do_t() self.a32 = 1u self.do_t() self.a31 = 0u self.do_t() cdef do_s(self): self.a21 = 0 self.a22 = 0 self.do_r() self.a21 = 1 self.do_r() self.a22 = 1 self.do_r() self.a22 = -1 self.do_r() self.a21 = 0 self.do_r() self.a22 = 1 self.do_r() self.a21 = -1 self.do_r() self.a22 = -1 self.do_r() self.a22 = 0 self.do_r() def __call__(self, cond, ap_list, *args): if len(args) != 2u*self.numPrimes: raise ValueError cdef uint32_t p = self.p cdef uint32_t i cdef uint_cur **big_list_newp = <uint_cur **>sage_malloc(sizeof(uint_cur *)*self.numPrimes) for i in range(self.numPrimes): self.firstp = self.curves_in_GF(1, args[2u*i],i) self.secondp = self.curves_in_GF(2, args[2u*i+1u],i) self.lens[i] = self.numFirst*self.numSecond big_list_newp[i] = self.split_prime_combine(i) sage_free(self.firstp) sage_free(self.secondp) cdef uint_cur *bigger_list_new, *temp self.L = 1u for i in range(self.numPrimes): temp = self.integral_combine(bigger_list_new, big_list_newp[i],i) if i: sage_free(bigger_list_new) sage_free(big_list_newp[i]) bigger_list_new = temp cdef uint32_t L = self.L cdef uint32_t current = 1u cdef uint_cur *e = bigger_list_new gal_cond = gal_con(K(cond)) cN = Integer(K(cond).norm().abs()) cdef uint32_t N = cN tempfact = cN.factor() cdef int_cur *endLift = self.var_lifts+2304 cdef uint32_t Nfactlen = len(tempfact) cdef uint32_t *Nfact, *Nfactend a = K.gen(0) Nfact = <uint32_t *>sage_malloc(sizeof(uint32_t)*Nfactlen) for i in range(Nfactlen): Nfact[i] = tempfact[i][0]**tempfact[i][1] Nfactend = Nfact+Nfactlen mpz_init(self.s) mpz_init(self.t) mpz_init(self.u) for i in range(L): current += 1u self.A1 = e[0][0] self.A2 = e[0][1] self.B1 = e[0][2] self.B2 = e[0][3] self.cur_lift = self.var_lifts self.a11 = 0u self.a12 = 0u self.do_s() self.a11 = 1u self.do_s() self.a12 = 1u self.do_s() self.a11 = 0u self.do_s() self.cur_lift = self.var_lifts while self.cur_lift != endLift: if self.check_disc(self.cur_lift[0],N,Nfact,Nfactend): E = EllipticCurve(K,[self.cur_lift[0][2*j]+a*self.cur_lift[0][2*j+1] for j in range(5)]) c = E.conductor() if c == cond: if self.check_aps_against(E,ap_list): mpz_clear(self.t) mpz_clear(self.s) mpz_clear(self.u) sage_free(Nfact) sage_free(big_list_newp) sage_free(bigger_list_new) return E elif c == gal_cond: print E E = EllipticCurve(K,[gal_con(K(t)) for t in E.ainvs()]) if self.check_aps_against(E,ap_list): mpz_clear(self.t) mpz_clear(self.s) mpz_clear(self.u) sage_free(Nfact) sage_free(big_list_newp) sage_free(bigger_list_new) return E self.cur_lift+=1 e += 1 mpz_clear(self.t) mpz_clear(self.s) mpz_clear(self.u) sage_free(Nfact) sage_free(big_list_newp) sage_free(bigger_list_new) return None cdef bint check_disc(self, int_cur ainv, uint32_t N, uint32_t *fact, uint32_t *factEnd): cdef int32_t b21, b22, b41, b42, b61, b62, b81, b82 cdef int64_t temp b21 = ainv[0]+ainv[1] b22 = ainv[1]*(ainv[0]+b21)+(ainv[3]<<2) b21 += ainv[2]<<2 b81 = b22*ainv[9] b82 = b21*ainv[9]+b22*ainv[8]+b81 b81 += b21*ainv[8] b41 = ainv[1]*ainv[5] b42 = b41+ainv[0]*ainv[5]+ainv[1]*ainv[4] b41 += ainv[0]*ainv[4] temp = ainv[7]*b42 b82 -= temp + ainv[6]*b42 + ainv[7]*b41 b81 -= temp + ainv[6]*b41 b41 += ainv[6]<<1 b42 += ainv[7]<<1 b61 = ainv[4]+ainv[5] b62 = ainv[5]*(ainv[4]+b61) temp = ainv[3]*b62 b82 += temp + ainv[2]*b62 + ainv[3]*b61 b81 += temp + ainv[2]*b61 b61 += ainv[8]<<2 b62 += ainv[9]<<2 temp = ainv[7]*ainv[7] b81 -= temp + ainv[6]*ainv[6] b82 -= temp + ((ainv[6]*ainv[7])<<1) cdef int64_t tmp1,tmp2 tmp1 = b22*b42 tmp2 = tmp1+b22*b41+b21*b42-3*b62 tmp1 += b21*b41-3*b61 temp = b62*tmp2 tmp2 = 9*(temp+b62*tmp1+b61*tmp2) tmp1 = 9*(temp+tmp1*b61) mpz_set_si(self.s,tmp1) mpz_set_si(self.t,tmp2) tmp1 = b41*b41 tmp2 = b42*b42 temp = (b42+3*b41)*tmp2 tmp2 = temp+(tmp2+3*tmp1)*b42<<3 tmp1 = temp+tmp1*b41<<3 if tmp1 >= 0: mpz_sub_ui(self.s,self.s,tmp1) else: mpz_add_ui(self.s,self.s,-tmp1) if tmp2 >= 0: mpz_sub_ui(self.t,self.t,tmp2) else: mpz_add_ui(self.t,self.t,-tmp2) temp = b21*b21+b22*b22 b22 = ((b21<<1)+b22)*b22 b21 = temp tmp1 = b81 tmp2 = b82 temp = b22*tmp2 tmp2 = temp+tmp1*b22+tmp2*b21 tmp1 = temp+b21*tmp1 if tmp1 >= 0: mpz_sub_ui(self.s,self.s,tmp1) else: mpz_add_ui(self.s,self.s,-tmp1) if tmp2 >= 0: mpz_sub_ui(self.t,self.t,tmp2) else: mpz_add_ui(self.t,self.t,-tmp2) if not mpz_sgn(self.s) and not mpz_sgn(self.t): return 0 mpz_sub(self.u,self.s,self.t) mpz_mul(self.u,self.u,self.t) mpz_addmul(self.u,self.s,self.s) if mpz_mod_ui(self.t,self.u,N): return 0 while fact != factEnd: while not mpz_mod_ui(self.t,self.u,fact[0]): mpz_divexact_ui(self.u,self.u,fact[0]) fact += 1 if mpz_sgn(self.u) < 0: mpz_neg(self.u,self.u) return mpz_root(self.u, self.u, 12ul) def ap(self,E,p): return E.change_ring(p.residue_field()).trace_of_frobenius() def check_aps_against(self, E, aps): for i,p in enumerate([11,19,29,31,41,59,61,71,79,89]): if not isinstance(aps[2*i], str): r1,r2 = K.primes_above(p) if self.ap(E,r1) != aps[2*i] or self.ap(E,r2) != aps[2*i+1]: return False return True 

## In this cell, we define a = (1+sqrt(5))/2, import necessary functions, and define the aplist function ## from psage.ellcurve.minmodel.sqrt5 import canonical_model from psage.modform.hilbert.sqrt5.sqrt5 import F as K a = K.gen(0) def ap(E,p): return E.change_ring(p.residue_field()).trace_of_frobenius() def aplist(E): aplist=[] for p in [2,3,5,7,11,19,29,31,41,59,61,71,79,89]: for i in K.primes_above(p): try: aplist.append(ap(E,i)) except ArithmeticError: aplist.append('?') return aplist 

## In this cell, we first construct instances of the SievedSearch class (do[p]) for each split prime p up to ## 89. Then we define a function curve_finder which loops through these split primes and calls do[p]'s __call__ ## method for ap1 and ap2 (where p1 and p2 are primes above the split prime p), or for ap1, ap2, aq1, and aq2 ## (if we want to look in the ring O_K/(pq)). If a curve is found, the function returns its canonical model. from psage.ellcurve.minmodel.sqrt5 import canonical_model do = {} primes = [11,19,29,31,41,59,61,71,79,89] for i,p in enumerate(primes): do[p] = SievedSearch(p) for p2 in primes[i+1:]: do[p*p2] = SievedSearch(p,p2) def curve_finder(cond, aps,imj=False): for i in range(10): pr = primes[i] ap1 = aps[2*i] if not isinstance(ap1, str): ap2 = aps[2*i+1] p1,p2 = K.primes_above(pr) if p1.integral_basis()[1][0] > p2.integral_basis()[1][0]: ap1,ap2=ap2,ap1 if imj == False: print 'Current Search:', pr E = do[pr](cond,aps,ap1,ap2) if E is not None: return canonical_model(E.global_minimal_model()) else: for j in range(i+1,10): pr2 = primes[j] ap3 = aps[2*j] if not isinstance(ap3, str): ap4 = aps[2*j+1] print 'Current search:', pr*pr2 p1,p2 = K.primes_above(pr2) if p1.integral_basis()[1][0] > p2.integral_basis()[1][0]: ap3,ap4=ap4,ap3 E = do[pr*pr2](cond,aps,ap1,ap2,ap3,ap4) if E is not None: return canonical_model(E.global_minimal_model()) return None 

%time ## Here is an example of the SievedSearch class in action ## #f1 = open('/Users/bleveque/Sage/unmatched_728.txt') missing_curves = ["1475 -35*a + 25 [1,1,'?',5,-1,-1,8,7,-9,9,-8,1,6,11,'?','?',-11,-2,9,6,3,-10,10,0]"] for s in missing_curves: S = s.split('[') aps = eval('['+S[1]) aps = aps[4:] print S[0] t = S[0].split(' ') cond = eval(''.join(t[1:])) E = curve_finder(cond,aps,imj=True) if E is None: print 'none' else: print E.conductor(), E.a_invariants() 

1475 -35*a + 25 
Current search: 209
Current search: 319
Current search: 341
Current search: 451

1475 -35*a + 25 
Current search: 209
Current search: 319
Current search: 341
Current search: 451

x = var('x') K.<a> = NumberField(x^2-x-1);K 

Number Field in a with defining polynomial x^2 - x - 1

Number Field in a with defining polynomial x^2 - x - 1

## In this cell we define some helpful functions. ## def grabber(line): return eval('['+line.split('[')[1].split(']')[0]+']') def grabber_two(line): return eval('['+line.split(']')[1].split('[')[1].split(']')[0]+']') def ap(E,p): return E.change_ring(K.residue_field(p)).trace_of_frobenius() def aplist(E): aps = [] for p in [2,3,5,7,11,19,29,31,41,59,61,71,79,89]: for i in K.primes_above(p): try: aps.append(ap(E,i)) except ArithmeticError: aps.append('?') return aps def galois_conj(beta): return beta[0]+beta[1]*(1-a) def check_aps_against(E, aps): Eaps = aplist(E) for i,p in enumerate([11,19,29,31,41,59,61,71,79,89]): if not isinstance(aps[2*i+4], str): r1,r2 = K.primes_above(p) if not isinstance(Eaps[2*i+4],str): if Eaps[2*i+4] != aps[2*i+4]: return False if not isinstance(Eaps[2*i+5],str): if Eaps[2*i+5] != aps[2*i+5]: return False return True #def check_aps_against(E, aps): # for i,p in enumerate([11,19,29,31,41,59,61,71,79,89]): # if not isinstance(aps[2*i+4], str): # r1,r2 = K.primes_above(p) # if ap(E,r1) != aps[2*i+4] or ap(E,r2) != aps[2*i+5]: # return False # return True 

E = EllipticCurve(K,[1,a-1,a+1,-2*a+3,-1]) 

%time check_aps_against(E,['?', -4, -1, 1, -2, 5, -8, -1, '?', '?',-3, 4, 0, -7, 4, 4, -8, 6, 16, 9, 10, -11,-8,-1]) 

True
CPU time: 0.22 s,  Wall time: 0.22 s

True
CPU time: 0.22 s,  Wall time: 0.22 s

## In this cell, we find all curves in our list of known curves (f) which have norm conductor dividing the ## norm conductors of our missing curves. These are candidates for quadratic twists. f = open('/home/bleveque/all_curves_so_far_2000.txt') norm_list = [(1100,['?',4,'?',-6,-4,'?',6,-4,6,10,2,-2,6,2,6,-10,-14,-8,0,-8,4,-8,6,-14])] di={} for N_u,aps in norm_list: di[N_u] = [] for line in f.readlines(): N = eval(line.split()[0]) for N_u,aps in norm_list: if N_u%N == 0: if N<1000: di[N_u].append((N,grabber(line))) di[N_u].append((N,grabber_two(line))) else: di[N_u].append((N,grabber(line))) 

print di[1100] 

[(55, [1, -a + 1, 1, -a, 0]), (55, [1, a, 1, a - 1, 0]), (100, [1, 1, 1,
-3, 1]), (100, [1, 1, 1, -3, 1]), (100, [a + 1, a, a + 1, -25*a - 25,
53*a + 46]), (100, [a + 1, a, a + 1, -25*a - 25, 53*a + 46]), (220, [a +
1, a - 1, a, -a - 3, -2*a + 3]), (220, [a, a + 1, 0, 2*a - 3, -a + 2]),
(220, [1, -1, a, -4*a - 1, 3*a + 3]), (220, [1, -1, a + 1, 3*a - 5, -4*a
+ 6]), (220, [a, -a + 1, 0, -2*a, 0]), (220, [a + 1, 0, 0, 2*a - 2, 0]),
(275, [a, a, a, 332*a - 536, -3529*a + 5713]), (275, [a + 1, a + 1, a,
-330*a - 203, 2994*a + 1853]), (1100, [a + 1, 1, a + 1, a - 1, 1]),
(1100, [a, a, 1, -2, -a + 1]), (1100, [a + 1, -a + 1, 0, 6*a - 10, 12*a
- 22]), (1100, [a, -a - 1, a, 11*a + 3, -16*a - 14]), (1100, [1, -a + 1,
a, 84*a - 138, -400*a + 653]), (1100, [a + 1, -a + 1, 1, 45*a - 49,
-89*a + 152]), (1100, [1, -a, 1, 316*a - 534, 3501*a - 5624])]

[(55, [1, -a + 1, 1, -a, 0]), (55, [1, a, 1, a - 1, 0]), (100, [1, 1, 1, -3, 1]), (100, [1, 1, 1, -3, 1]), (100, [a + 1, a, a + 1, -25*a - 25, 53*a + 46]), (100, [a + 1, a, a + 1, -25*a - 25, 53*a + 46]), (220, [a + 1, a - 1, a, -a - 3, -2*a + 3]), (220, [a, a + 1, 0, 2*a - 3, -a + 2]), (220, [1, -1, a, -4*a - 1, 3*a + 3]), (220, [1, -1, a + 1, 3*a - 5, -4*a + 6]), (220, [a, -a + 1, 0, -2*a, 0]), (220, [a + 1, 0, 0, 2*a - 2, 0]), (275, [a, a, a, 332*a - 536, -3529*a + 5713]), (275, [a + 1, a + 1, a, -330*a - 203, 2994*a + 1853]), (1100, [a + 1, 1, a + 1, a - 1, 1]), (1100, [a, a, 1, -2, -a + 1]), (1100, [a + 1, -a + 1, 0, 6*a - 10, 12*a - 22]), (1100, [a, -a - 1, a, 11*a + 3, -16*a - 14]), (1100, [1, -a + 1, a, 84*a - 138, -400*a + 653]), (1100, [a + 1, -a + 1, 1, 45*a - 49, -89*a + 152]), (1100, [1, -a, 1, 316*a - 534, 3501*a - 5624])]

## Here, we loop through the curves in each entry of the dictionary we created above, and for each one we twist ## by i+j*a for i and j such that i,j<=B and i>=L or j>=L. If the conductor norm and a_p-list of the twist (or ## its conjugate) are equal to the conductor norm and a_p-list of the curve we are searching for, then we print ## the twist. B = 3 L = 0 for N_u,aps in norm_list: print N_u for ct in di[N_u]: for i in [-B..B]: for j in [-B..B]: if i>=L or j>=L: if i!=0 or j!=0: d = i+j*a if N_u % ct[0]*(d.norm())^2 == 0: E = EllipticCurve(K,ct[1]) E_d = E.quadratic_twist(d) N = E_d.conductor().norm() Gainv = [galois_conj(A) for A in E_d.a_invariants()] GE_d = EllipticCurve(K,Gainv) if N == N_u: if check_aps_against(E_d,aps) or check_aps_against(GE_d,aps): print N, E_d.a_invariants(), aplist(E_d) 

1100
1100 (0, -9*a + 2, 0, 200*a - 360, 2560*a - 4080) ['?', 4, '?', -6, -4,
'?', 6, -4, 6, 10, 2, -2, 6, 2, 6, -10, -14, -8, 0, -8, 4, -8, 6, -14]

1100
1100 (0, -9*a + 2, 0, 200*a - 360, 2560*a - 4080) ['?', 4, '?', -6, -4, '?', 6, -4, 6, 10, 2, -2, 6, 2, 6, -10, -14, -8, 0, -8, 4, -8, 6, -14]

%cython ############################################## # The following is code to find curves using # torsion families. It uses cython code to # calculate the discriminant of a curve from # its a-invariants. ############################################## from libc.stdint cimport uint32_t, int32_t, uint64_t, int64_t from sage.schemes.elliptic_curves.constructor import EllipticCurve from psage.modform.hilbert.sqrt5.sqrt5 import F as K from sage.rings.integer cimport Integer cdef int64_t quick_disc_norm(int32_t A, int32_t B, int32_t C, int32_t D, int32_t E, int32_t F, int32_t G, int32_t H, int32_t I, int32_t J): cdef int32_t b21, b22, b41, b42, b61, b62, b81, b82, disc1, disc2, tp1,tp2,tp3,tp4,tp5,tp6 cdef int64_t temp, ret cdef mpz_t s,t,u b21 = A**2+B**2+4*C b22 = B**2+2*A*B+4*D tp1 = A*E+B*F tp2 = A*F+B*(E+F) b41 = tp1+2*G b42 = tp2+2*H tp3 = E**2+F**2 tp4 = 2*E*F+F**2 b61 = tp3+4*I b62 = tp4+4*J b81 = b21*I+b22*J-G*tp1-H*tp1+C*tp3+D*tp4-G**2-H**2 b82 = b21*J+b22*(I+J)-H*(tp1+tp2)-G*tp2+D*(tp3+tp4)+C*tp4-2*G*H-H**2 tp1 = b21**2+b22**2 tp2 = 2*b21*b22+b22**2 tp3 = b41**2+b42**2 tp4 = 2*b41*b42+b42**2 tp5 = b21*b41+b22*b42 tp6 = b21*b42+b22*b41+b22*b42 disc1 = -(tp1)*b81-(tp2)*b82-8*(tp3)*b41-8*(tp4)*b42-27*(b61**2+b62**2)+9*(tp5)*b61+9*(tp6)*b62 disc2 = -(tp1)*b82-(tp2)*(b81+b82)-8*(tp3)*b42-8*(tp4)*(b41+b42)-27*(2*b61*b62+b62**2)+9*(tp5)*b62+9*(tp6)*(b61+b62) ret = disc1**2+disc1*disc2-disc2**2 return ret ################################################ # Below we define functions to compute # curves with given torsion structure using # the parametrizations given in Kubert's # paper. Note: some of these parametrizations # need four values while others only need # two. This can be cleaned up a lot by defining # a function that all of them call to check # the discriminant, print a curve if necessary, # etc... Some of the methods below can also # be sped up by working component-wise and not # introducing a = K.gens(0) until as late as # possible. ################################################ #def check_disc_cond( def family2_new_curve_check(norm_list, uint32_t B, uint32_t L = 0): cdef int32_t u,v,y,z cdef int64_t N_d for u in range(-B,B+1): for v in range(-B,B+1): for y in range(-B,B+1): for z in range(-B,B+1): if abs(u)>=L or abs(v)>=L or abs(y)>=L or abs(z)>=L: N_d = quick_disc_norm(0,0,u,v,0,0,y,z,0,0) if N_d != 0: for N in norm_list: if N_d%N: continue else: a = K.gen(0) E = EllipticCurve(K,[0,u+v*a,0,y+z*a,0]) cond = E.conductor() if cond.norm() in norm_list: print cond.norm(),'---', cond,'---',E.a_invariants() break def family3_new_curve_check(norm_list, uint32_t B, uint32_t L = 0): cdef int32_t u,v,y,z cdef int64_t N_d for u in range(-B,B+1): for v in range(-B,B+1): for y in range(-B,B+1): for z in range(-B,B+1): if abs(u)>=L or abs(v)>=L or abs(y)>=L or abs(z)>=L: N_d = quick_disc_norm(u,v,0,0,y,z,0,0,0,0) if N_d != 0: for N in norm_list: if N_d%N: continue else: a = K.gen(0) E = EllipticCurve(K,[u+v*a,0,y+z*a,0,0]) cond = E.conductor() if cond.norm() in norm_list: print cond.norm(),'---', cond,'---',E.a_invariants() break def family7_new_curve_check(norm_list, uint32_t B, uint32_t L = 0): cdef int32_t u,v cdef int64_t N_d for u in range(-B,B+1): for v in range(-B,B+1): if abs(u)>=L or abs(v)>=L: tp1 = u**2+v**2 tp2 = 2*u*v+v**2 c1 = tp1 - u c2 = tp2 - v b1 = tp1*u + tp2*v -tp1 b2 = tp1*v+tp2*(u+v)-tp2 N_d = quick_disc_norm(1-c1,-c2,-b1,-b2,-b1,-b2,0,0,0,0) #print N_d if N_d != 0: for N in norm_list: if N_d%N: continue else: a = K.gen(0) E = EllipticCurve(K,[1-(c1+c2*a),-b1-b2*a,-b1-b2*a,0,0]) cond = E.conductor() if cond.norm() in norm_list: print cond.norm(),'---', cond,'---',E.a_invariants() break def family5_new_curve_check(norm_list, uint32_t B, uint32_t L = 0): cdef int32_t u,v cdef int64_t N_d for u in range(-B,B+1): for v in range(-B,B+1): if abs(u)>=L or abs(v)>=L: N_d = quick_disc_norm(1-u,-v,-u,-v,-u,-v,0,0,0,0) if N_d != 0: for N in norm_list: if N_d%N: continue else: a = K.gen(0) E = EllipticCurve(K,[1-(u+v*a),-u-v*a,-u-v*a,0,0]) cond = E.conductor() if cond.norm() in norm_list: print cond.norm(),'---', cond,'---',E.a_invariants() break def family2x2_new_curve_check(norm_list, uint32_t B, uint32_t L = 0): cdef int32_t u,v,y,z cdef int64_t N_d for u in range(-B,B+1): for v in range(-B,B+1): for y in range(-B,B+1): for z in range(-B,B+1): if abs(u)>=L or abs(v)>=L or abs(y)>=L or abs(z)>=L: N_d = quick_disc_norm(0,0,u+y,v+z,0,0,u*y+v*z,u*z+v*y+v*z,0,0) if N_d != 0: for N in norm_list: if N_d%N: continue else: a = K.gen(0) E = EllipticCurve(K,[0,(u+y)+(v+z)*a,0,(u*y+v*z)+(u*z+v*(y+z))*a,0]) cond = E.conductor() if cond.norm() in norm_list: print cond.norm(),'---', cond,'---',E.a_invariants() break def family6_new_curve_check(norm_list, uint32_t B, uint32_t L = 0): cdef int32_t u,v cdef int64_t N_d for u in range(-B,B+1): for v in range(-B,B+1): if abs(u)>=L or abs(v)>=L: b1 = u**2+v**2+u b2 = 2*u*v+v**2+v N_d = quick_disc_norm(1-u,-v,-b1,-b2,-b1,-b2,0,0,0,0) if N_d != 0: for N in norm_list: if N_d%N: continue else: a = K.gen(0) E = EllipticCurve(K,[1-(u+v*a),-b1-b2*a,-b1-b2*a,0,0]) cond = E.conductor() if cond.norm() in norm_list: print cond.norm(),'---', cond,'---',E.a_invariants() break def family10_new_curve_check(norm_list, uint32_t B, uint32_t L = 0): cdef int32_t u,v cdef int64_t N_d for u in range(-B,B+1): for v in range(-B,B+1): if abs(u)>=L or abs(v)>=L: a=K.gen(0) f = u+v*a if f!=(f-1)**2: d = f**2/(f-(f-1)**2) c = f*d-f b = c*d N_d = quick_disc_norm(1-c[0],-c[1],-b[0],-b[1],-b[0],-b[1],0,0,0,0) if N_d != 0: for N in norm_list: if N_d%N: continue else: E = EllipticCurve(K,[1-c,-b,-b,0,0]) cond = E.conductor() if cond.norm() in norm_list: print cond.norm(),'---', cond,'---',E.a_invariants() break def family4_new_curve_check(norm_list, uint32_t B, uint32_t L = 0): cdef int32_t u,v cdef int64_t N_d for u in range(-B,B+1): for v in range(-B,B+1): if abs(u)>=L or abs(v)>=L: N_d = quick_disc_norm(1,0,-u,-v,-u,-v,0,0,0,0) if N_d != 0: for N in norm_list: if N_d%N: continue else: a = K.gen(0) E = EllipticCurve(K,[1,-u-v*a,-u-v*a,0,0]) cond = E.conductor() if cond.norm() in norm_list: print cond.norm(),'---', cond,'---',E.a_invariants() break def family8x2_new_curve_check(norm_list, uint32_t B, uint32_t L = 0): cdef int32_t u,v cdef int64_t N_d for u in range(-B,B+1): for v in range(-B,B+1): if abs(u)>=L or abs(v)>=L: a = K.gen(0) alpha = u+v*a d = alpha*(8*alpha+2)/(8*alpha^2-1) b = (2*d-1)*(d-1) c = (2*d-1)*(d-1)/d N_d = quick_disc_norm(c[1-c[0],-c[1],-b[0],-b[1],-b[0],-b[1],0,0,0,0) if N_d != 0: for N in norm_list: if N_d%N: continue else: a = K.gen(0) E = EllipticCurve(K,[1-c,-b,-b,0,0]) cond = E.conductor() if cond.norm() in norm_list: print cond.norm(),'---', cond,'---',E.a_invariants() break def family8_new_curve_check(norm_list, uint32_t B, uint32_t L = 0): cdef int32_t u,v cdef int64_t N_d for u in range(-B,B+1): for v in range(-B,B+1): if abs(u)>=L or abs(v)>=L: a = K.gen(0) d = u+v*a b = (2*d-1)*(d-1) c = (2*d-1)*(d-1)/d N_d = quick_disc_norm(c[1-c[0],-c[1],-b[0],-b[1],-b[0],-b[1],0,0,0,0) if N_d != 0: for N in norm_list: if N_d%N: continue else: a = K.gen(0) E = EllipticCurve(K,[1-c,-b,-b,0,0]) cond = E.conductor() if cond.norm() in norm_list: print cond.norm(),'---', cond,'---',E.a_invariants() break def family6x2_new_curve_check(norm_list, uint32_t B, uint32_t L = 0): cdef int32_t u,v cdef int64_t N_d for u in range(-B,B+1): for v in range(-B,B+1): if abs(u)>=L or abs(v)>=L: a = K.gen(0) alpha = u+v*a c = (10-2*alpha)/(alpha^2-9) b = c+c^2 N_d = quick_disc_norm(c[1-c[0],-c[1],-b[0],-b[1],-b[0],-b[1],0,0,0,0) if N_d != 0: for N in norm_list: if N_d%N: continue else: a = K.gen(0) E = EllipticCurve(K,[1-c,-b,-b,0,0]) cond = E.conductor() if cond.norm() in norm_list: print cond.norm(),'---', cond,'---',E.a_invariants() break def family12_new_curve_check(norm_list, uint32_t B, uint32_t L = 0): cdef int32_t u,v cdef int64_t N_d for u in range(-B,B+1): for v in range(-B,B+1): if abs(u)>=L or abs(v)>=L: a = K.gen(0) t = u+v*a m = (3*t-3*t^2-1)/(t-1) f = m/(1-t) d = m+t c = f*d-f b = c*d N_d = quick_disc_norm(c[1-c[0],-c[1],-b[0],-b[1],-b[0],-b[1],0,0,0,0) if N_d != 0: for N in norm_list: if N_d%N: continue else: a = K.gen(0) E = EllipticCurve(K,[1-c,-b,-b,0,0]) cond = E.conductor() if cond.norm() in norm_list: print cond.norm(),'---', cond,'---',E.a_invariants() break 

norm_list=[1476] 

E = EllipticCurve(K,[1,-a,0,-a,-7*a-3]) E.torsion_subgroup() 

Torsion Subgroup isomorphic to Z/2 associated to the Elliptic Curve
defined by y^2 + x*y = x^3 + (-a)*x^2 + (-a)*x + (-7*a-3) over Number
Field in a with defining polynomial x^2 - x - 1

Torsion Subgroup isomorphic to Z/2 associated to the Elliptic Curve defined by y^2 + x*y = x^3 + (-a)*x^2 + (-a)*x + (-7*a-3) over Number Field in a with defining polynomial x^2 - x - 1

%time family4_new_curve_check(norm_list,60) 

CPU time: 4.78 s,  Wall time: 4.78 s

CPU time: 4.78 s,  Wall time: 4.78 s

%time family2x2_new_curve_check(norm_list,20) 

%time family8_new_curve_check(norm_list,60) 

%time family5_new_curve_check(norm_list,60) 

%time family6_new_curve_check(norm_list,60) 

%time family6x2_new_curve_check(norm_list,60) 

%time family8x2_new_curve_check(norm_list,60) 

%time family12_new_curve_check(norm_list,60)