# This is slightly unfortunate, as reloads of attached files will overwrite symbols defined in the running session. # TODO: Can we go finer grain ? from sage.all_cmdline import * # import sage library from sage.rings.factorint import factor_trial_division import sage.rings.integer from enumerate_sparse_T import * import itertools from collections import defaultdict import re import pprint import random from hashlib import md5 import time import pickle #Zimmermann trick max_prime = 10**6 primes_list = prime_range(max_prime) prod_primes = prod(primes_list) seed=0 # utilities def md5int(T): return Integer(md5(str(T)).hexdigest(), 16) def get_prod_pi(k): """ return the product of all primes that are congruent to 1 mod k, or that divide k. """ prod_pi = prod(k.prime_factors()) for p in primes_list: if p%k == 1: prod_pi *= p return prod_pi def phex(t): if abs(t) < 10: return str(t) else: return "%s0x%s" % ("-" if t < 0 else "", abs(t).hex()) class sequence_with_stateful_iterators(object): def __init__(self, S): self.S = S def __iter__(self): return self.iterator(self) class iterator(object): def __init__(self, range_parent): self.S = range_parent.S self.i = 0 self.previous_state = self.getstate() def __iter__(self): return self def __next__(self): return self.next() def next(self): self.previous_state = self.getstate() if self.i >= len(self.S): raise StopIteration i = self.i self.i += 1 return self.S[i] def getstate(self): return self.i def setstate(self, state): self.previous_state = state self.i = state class range_with_stateful_iterators(object): def __init__(self, x0, x1): self.x0 = x0 self.x1 = x1 def __iter__(self): return self.iterator(self) class iterator(object): def __init__(self, range_parent): self.x0 = range_parent.x0 self.x1 = range_parent.x1 self.x = self.x0 self.previous_state = self.getstate() def __iter__(self): return self def __next__(self): return self.next() def next(self): self.previous_state = self.getstate() if self.x >= self.x1: raise StopIteration x = self.x self.x += 1 return x def getstate(self): return self.x def setstate(self, state): self.previous_state = state self.x = state class random_chooser_for_integer_of_low_weight(object): def __init__(self, source, weight, x0, x1, random_state=None): self.f = source.successor self.random_state = random_state self.weight = weight self.x0 = x0 self.x1 = x1 self.previous_state = self.getstate() def __iter__(self): return self def __next__(self): return self.next() def next(self): self.previous_state = self.getstate() while True: if self.random_state: r=self.random_state.randrange(self.x0,self.x1) else: r=randrange(self.x0,self.x1) r=self.f(Integer(r), self.weight) if r < self.x1: return r def getstate(self): if self.random_state: return self.random_state.getstate() else: return None def setstate(self, state): self.previous_state = state if self.random_state: self.random_state.setstate(state) else: # maybe warn, or even die ? pass class cap_range(object): def __init__(self, r, cap): self.r = r self.cap = cap def __iter__(self): return self.iterator(self) def __len__(self): return self.cap class iterator(object): # probably want to have that inherit from the original type, right ? # or maybe not. After all, self.i is part of the state, so we've got # to do something with it. def __init__(self, parent): self.ri = iter(parent.r) self.cap = parent.cap self.i = 0 self.previous_state = self.getstate() def __next__(self): return self.next() def next(self): self.previous_state = self.getstate() if self.i >= self.cap: raise StopIteration self.i += 1 return next(self.ri) def __iter__(self): return self def getstate(self): return (self.i, self.ri.getstate()) def setstate(self, state): i,s = state self.i = i self.ri.setstate(s) def range_with_strategy(x0, x1, split=None, strategy=None, random_state=None): """ given a strategy string and a pair of integers describing an interval [x0, x1), return a range object whose iterators follows the requested strategy. The syntax of the strategy argument is a comma-separated string of tokens which may be: - random: abide by the other parameters, except that we do random picks. - hamming=X: restrict to Hamming weight equal to X (int, positive) - 2-naf=X: restrict to 2-NAF weight equal to X (int, positive) - x0=X: restrict to x>=x0 (intersects with the provided argument) - x1=X: restrict to x= self.parent.cap: raise StopIteration self.current += 1 x0 = self.parent.x0 x1 = self.parent.x1 if self.parent.random_state: return self.parent.random_state.randrange(x0, x1) else: return randrange(x0, x1) def getstate(self): if self.parent.random_state: return self.parent.random_state.getstate() else: return None def setstate(self, state): self.previous_state = state if self.parent.random_state: self.parent.random_state.setstate(state) else: # maybe warn, or even die ? pass def __init__(self, x0, x1, cap, random_state = None): self.cap = cap self.random_state = random_state self.x0 = x0 self.x1 = x1 def __iter__(self): return self.rsc_iterator(self) def __len__(self): return self.cap source = None rand_c = False up_to = False weight=0 cap = None if strategy is None: strategy="" # We want h in the interval [x0, x1). # First check for explicit instructions on x0 and x1 for ss in strategy.split(','): m=re.match("x0=(0x[0-9a-f]+|\d+)",ss) if m: x0=max(x0,int(m.group(1),0)) continue m=re.match("x1=(0x[0-9a-f]+|\d+)",ss) if m: x1=min(x1,int(m.group(1),0)) continue m=re.match("max=(0x[0-9a-f]+|\d+)",ss) if m: x1=min(x1,int(m.group(1),0)) x0=max(x0,1-int(m.group(1),0)) continue if split: i,n = split dx = x1 - x0 x1 = x0 + ((i + 1) * dx) // n x0 = x0 + (i * dx) // n for ss in strategy.split(','): if ss == "random": rand_c = True continue m=re.match("cap=(\d+)",ss) if m: cap=int(m.group(1)) continue m=re.match("(hamming|2-naf)(<)?=(\d+)",ss) if m: assert source is None up_to=m.group(2)=='<' weight=int(m.group(3)) if m.group(1) == 'hamming': source = enumerator_factory_hamming elif m.group(1) == '2-naf': source = enumerator_factory_2naf continue if source is None: if rand_c: return random_state_chooser_in_range(x0, x1, cap, random_state) else: # want a stateful iterator in this range. # return xsrange(x0,x1) r = range_with_stateful_iterators(x0,x1) elif not rand_c: assert source is not None r = enumerator_factory(source, weight, x0, x1, up_to=up_to) else: # Does not distinguish between <= weight and = weight. r = random_chooser_for_integer_of_low_weight(source, weight, x0, x1, random_state=random_state) if cap is not None: r = cap_range(r, cap) return r def attempt_to_factor(n, known=[], max_trialdiv=10**6, max_B1=10**4): """ tries without too much effort to factor n. Return a list of primes with multiplicities, and a list of composites without multiplicities INPUT: n: integer to (attempt to) factor known: list of already known factors max_trialdiv (integer): bound for trial division max_B1 (integer): bound for ECM OUTPUT: primes: a list of primes with multiplicities (tuples) composites': a list of composites without multiplicities """ primes=[] composites=[] for p in known: assert p.is_prime(proof=False) pp = gcd(n, p) if pp==1: continue n = n // pp primes.append((pp,1)) # we have one single composite n at this point td = n.factor(limit=max_trialdiv) for p,e in td: if p.is_prime(proof=False): primes.append((p,e)) else: assert e==1 composites.append(p) # theoretically ECM should work even with integers that are not # squarefree. It's only partly true though, see e.g.: # ECM(100).factor(1009^2*prod([random_prime(4000,lbound=1024) for i in range(3)])) # (maybe try a few times, but at least with sage 8.5 it won't take # long until you get an infinite loop). B1=400 nc=20 while B1 < max_B1 and composites: oc=composites composites=[] for n in oc: L=defaultdict(int) # print "Trying ECM with B1=%d, c=%d (max_B1=%d)" % (B1, nc, max_B1) global seed seed = md5int(seed) e=ECM(B1=B1, sigma=seed) # print "echo %d | %s" % (n, " ".join(e._cmd)) try: fn=e.one_curve(n, B1=B1, c=nc) except IndexError: # bug in sage ECM interface. We probably found the input # number N anyway, so let's move on. # Example: # sage: ECM().one_curve(3808350204649, B1=400000, sigma=42) fn=[n] except KeyboardInterrupt: raise except: print "ECM failed on %d" % n raise for p in fn: L[p]+=1 assert prod(L)==n for p,e in L.items(): if p.is_prime(proof=False): primes.append((p,e)) elif p>1: assert e==1 composites.append(p) if oc == composites: B1 = int(B1 * sqrt(5)) nc = int(nc * sqrt(2)) return primes,composites def chain_alternate_iterators(gp, gm, with_zero=False): gp = iter(gp) gm = iter(gm) while True: try: yield next(gp) yield next(gm) except StopIteration: break while True: try: yield next(gp) except StopIteration: break while True: try: yield next(gm) except StopIteration: break # class for holding the final results. # # run "make check" to make sure that all tests below fulfill their # promises. class CocksPinchVariantResult(object): """ sage: C=CocksPinchVariantResult(6,3,34359607296,5,ht=0x101,hy=-2,max_B1=1000) sage: C.E2(factor=True)["text_factorization"] '2^2 * 3 * 19 * 73 * 163 * 33637 * p48 * r' sage: C=CocksPinchVariantResult(6,3,0x600100002,5,ht=0x428,hy=-0x639,allowed_cofactor=420,max_B1=600) sage: C.is_small_subgroup_secure() True sage: C.is_twist_small_subgroup_secure() True sage: C.is_G2_small_subgroup_secure() True sage: C.is_twist_G2_small_subgroup_secure() True sage: C=CocksPinchVariantResult(5,10000000019,0xeffff80100000000,1,ht=0x40200,hy=0xfe0000000000001) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure()) (True, True) sage: C=CocksPinchVariantResult(8,4,0x27d80,7,ht=-0x451,hy=-0x481) sage: C.set_test_info(allowed_size_cofactor=10) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(8,4,0x29072,7,ht=0x9bf,hy=-0x10e) sage: C.set_test_info(allowed_size_cofactor=10) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(8,4,0x29f24,7,ht=-0x289,hy=0x53f) sage: C.set_test_info(allowed_size_cofactor=10) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(8,4,0x2a1c8,3,ht=0x53f,hy=-0x437) sage: C.set_test_info(allowed_size_cofactor=10) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(8,4,0x27d80,7,ht=-0x451,hy=-0x481) sage: C.set_test_info(allowed_size_cofactor=10) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(8,4,0x2617e,5,ht=-0xd93,hy=0x305) sage: C.set_test_info(allowed_size_cofactor=10) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(8,4,0x28f86,3,ht=0x8cf,hy=0x2e0) sage: C.set_test_info(allowed_size_cofactor=10) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(6,3,0x600081000,1,ht=-0x191,hy=0x7e2) sage: C.set_test_info(allowed_size_cofactor=10) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(5,10000000147,0xe000000000008000,1,ht=3,hy=0x11e36418c7c8b454,max_B1=600) sage: C.set_test_info(allowed_size_cofactor=10) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(6,3,0xfffffffffffffff00000000000000000,1,ht=0x43fff,hy=0xffffffffff800007fffe,allowed_size_cofactor=10,max_B1=600) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(6,3,0xff800000000000200000000000000000,1,ht=-1,hy=0xffffff823ffffe008000,allowed_cofactor=420,allowed_size_cofactor=10,max_B1=600) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(6,3,0xffe00008000000000000000000000000,1,ht=-1,hy=0xffbfffe3f80200000000,allowed_cofactor=420,allowed_size_cofactor=10,max_B1=600) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(6,3,0xffe00008000000000000000000000000,1,ht=-1,hy=0xfffffd0010001ffc0000,allowed_cofactor=420,allowed_size_cofactor=10,max_B1=600) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(6,3,0xefffffffffffffe00000000000000000,1,ht=-1,hy=0xffbbffffffffffffc020,allowed_cofactor=420,allowed_size_cofactor=10,max_B1=600) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(7,20,0x5ec7fc01ff8,4,ht=-3,hy=1,allowed_size_cofactor=10,max_B1=600) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, False, False) sage: C=CocksPinchVariantResult(8,4,0xffffffffffffffc0,1,ht=-0x1821f,hy=-0x1fdc,allowed_cofactor=1232,allowed_size_cofactor=10,max_B1=600) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(8,4,0xffffffffeff7c200,5,ht=5,hy=-0xd700,allowed_cofactor=420,allowed_size_cofactor=10,max_B1=600) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) sage: C=CocksPinchVariantResult(8,4,0xffdffffc7ffffc00,3,ht=5,hy=0xc5f4,allowed_cofactor=420,allowed_size_cofactor=10,max_B1=600) sage: (C.is_small_subgroup_secure(), C.is_twist_small_subgroup_secure(), C.is_G2_small_subgroup_secure(), C.is_twist_G2_small_subgroup_secure()) (True, True, True, True) """ def __init__(self,k,D,T,i,ht=Integer(0),hy=Integer(0),max_poly_coeff=0,pre=False,l=1, allowed_cofactor = 1, allowed_size_cofactor = 5, max_trialdiv=10**6, max_B1=10**4, new_semantics=False): kl = k * l fD = -fundamental_discriminant(-D) self.k = k self.kl = kl self.fD = Integer(fD) self.T = Integer(T) self.i = i self.r = cyclotomic_polynomial(kl)(T) self.max_poly_coeff = max_poly_coeff self.allowed_size_cofactor = allowed_size_cofactor self.allowed_cofactor = allowed_cofactor self.max_trialdiv = max_trialdiv self.max_B1 = max_B1 r = self.r K = FiniteField(r, proof = False) self.t0 = Integers()(K(T)**i + 1) self.y0 = Integers()(K(self.t0-2)/sqrt(K(-fD))) # We like to have the smallest representative, in particular # because this matches what we do formally when we take the # representative that is given by the polynomial expression # having the smallest degree. if abs(r-self.t0) < abs(self.t0): self.t0 -= r if abs(r-self.y0) < abs(self.y0): self.y0 -= r self.y0 = abs(self.y0) # Determination of the lifted (t,y) from the solution mod r # (t0,y0) and the cofactors (ht,hy). Note that t0 and y0 are # chosen as centered integer representatives of the values mod r # (that is, both are in the interval [-(r-1)/2,(r-1)/2]). # t is (t0 + cofactor * r), where cofactor is: # 2*ht if -fD % 4 == 0 and t0 is even # 2*ht+1 if -fD % 4 == 0 and t0 is odd # ht if -fD % 4 == 1 # y is (y0 + cofactor * r), where cofactor is: # hy if -fD % 4 == 0 # 2*hy if -fD % 4 == and (t-y0) is even # 2*hy+1 if -fD % 4 == and (t-y0) is odd self.ht = Integer(ht) self.hy = Integer(hy) self.t = self.t0 + ht * r self.y = self.y0 + hy * r if pre: # compatibility mode with previous semantics. htformat = "%s" hyformat = "%s" if -fD % 4 == 0: htformat = "2 * %s" self.t += ht * r if self.t0 % 2 == 1: htformat = "(2 * %s + 1)" self.t += r elif -fD % 4 == 1: hyformat = "2 * %s" self.y += hy * r if (self.t - self.y0) % 2 == 1: hyformat = "(2 * %s + 1)" self.y += r old_call = "ht=%s,hy=%s,pre=True" % (phex(self.ht), phex(self.hy)) new_ht_expr = htformat % self.ht new_hy_expr = hyformat % self.hy new_ht = eval(preparse(new_ht_expr)) new_hy = eval(preparse(new_hy_expr)) new_call = "ht=%s=%s,hy=%s=%s" % (new_ht_expr,phex(new_ht),new_hy_expr,phex(new_hy)) # raise ValueError("Replace call with %s by %s" % (old_call, new_call)) # print("Replace call with %s by %s" % (old_call, new_call)) self.ht = new_ht self.hy = new_hy t = self.t y = self.y if (t**2 + self.fD*y**2) % 4 != 0: raise ValueError("1/2*((t0+%d*r)+(y0+%d*r)*sqrt(-%d)) not an algebraic integer: " %(ht,hy,self.fD)) self.p = (t**2 + self.fD*y**2)//4 p = self.p # Also compute all twists. self._twists=[] self._twists.append({ 'name': "", 'card': p+1-t }) self._twists.append({ 'name': "t", 'card': p+1+t }) # do we have quartic twists ? Only possible if the endomorphism ring # contains 4-th roots of unity, hence -fD==-4 if self.fD == 4: assert (2*y)**2 == 4*p-t**2 self._twists.append({ 'name': "q2", 'card': p+1+2*y }) self._twists.append({ 'name': "q3", 'card': p+1-2*y }) # As for sextic twist, this is only for -fD==-3 elif self.fD == 3: assert 3*y**2 == 4*p-t**2 tr3a = (t-3*y) // 2 tr3b = (t+3*y) // 2 self._twists.append({ 'name': "s2", 'card': p+1+tr3a }) self._twists.append({ 'name': "s3", 'card': p+1+tr3b }) self._twists.append({ 'name': "s4", 'card': p+1-tr3a }) self._twists.append({ 'name': "s5", 'card': p+1-tr3b }) self._E2 = None self._twist_E2 = None if pre: print "C=%s" % repr(self) def u(self): return self.T def T(self): return self.T def _prepare_E2(self): if self._E2 is not None: return k = self.k p = self.p fD = self.fD r = self.r t = self.t y = self.y ZP=ZZ['x'] x=ZP.gen() d = 2 if fD == 4: d = 4 elif fD == 3: d = 6 d0 = gcd(k, d) k0 = k // d0 # Then G2 is actually defined as an order r subgroup of a degree # d0 twist E2 of the same curve, over GF(p^k0) # First compute that curve over GF(p^k0). tk0 = (p**k0+1) - (x**2-t*x+p).resultant(x**k0-1) K = NumberField(x**2+fD,'z') z = K.gen() pi = (t+y*z)/2 pibar = (t-y*z)/2 assert tk0 == pi**k0 + pibar**k0 assert tk0 == (pi**k0).trace() yk0 = ZZ((pi**k0 - pibar**k0)/z) assert 4*p**k0 == tk0**2+fD*yk0**2 # We expect that exactly one of the non-trivial degree d0 twists # will have cardinal divisible by r. In fact, we can determine # this more directly, but more or less for the same cost anyway. # p^k0 will be a d0-th root of unity mod r. So one of T and -T if # d0==d==4, one of T and -T^2 if d0==d==6, and one of T^2 and # -T if d0==3 and d==6. One still has to determine p^k0 mod r # anyway, and identify the automorphism that multiplies by T with # the multiplication by some torsion unit. which_twist = [] if d0 == 1: which_twist.append((0, (p**k0+1) - tk0)) assert (p**k0+1 - tk0) % r**2 == 0 elif d0 == 2: # multiply pi by 1 which_twist.append((1, (p**k0+1) + tk0)) elif d0 == 4: assert fD == 4 # multiply pi^k0 by i -- as sqrt(-fD)=2i, a 2 pops out i = z/2 assert i**2 + 1 == 0 assert -2*yk0 == (pi**k0 * i).trace() which_twist.append((2, (p**k0+1) + 2 * yk0)) which_twist.append((3, (p**k0+1) - 2 * yk0)) elif d0 == 3 or d0 == 6: assert fD == 3 # multiply pi^k0 by j+1 j = (-1+z)/2 assert j**2 + j + 1 == 0 ak0 = (tk0-3*yk0) // 2 assert ak0 == (pi**k0 * (j+1)).trace() bk0 = (tk0+3*yk0) // 2 assert bk0 == (pi**k0 * -j).trace() if d0 == 6: which_twist.append((4, (p**k0+1) - ak0)) which_twist.append((5, (p**k0+1) - bk0)) else: # Only multiplication by j and j^2 -- the others are # sextic twists. which_twist.append((2, (p**k0+1) + ak0)) which_twist.append((3, (p**k0+1) + bk0)) # Let E2 be the curve of which G2 is a subgroup. # # Let Fpk0 = F_{p^{k0}} and Fpk=F_{p^k} # # E has a group of r-torsion # Ek0 is E extended to Fpk0 # Ek is E extended to Fpk which is a degree d0 extension of Fpk0 # Ek0 has d0 twists. # #Ek is the product of the # of these d0 twists # Ek has another subgroup of order r. # This subgroup of order r exists in one of the twists of Ek0 # (conceivably E0 itself if d0==1, in which case k0==0 and our # extension to Fpk0 did all the work). # # Let ' denotes the (possibly trivial) twisting action such that # Ek0' has this extra subgroup of order r. E2 is this Ek0' # # Because d0 = k/k0 = gcd(d,k), no divisor x of k0 is such that # d0*x s a divisor of d. Hence if ' is not trivial and k0>1, # there is no twist of E (or of an intermediate curve between E # and Ek0) that Ek0' is an extension of. On the contrary, if # k0==1, then E2 is no other than E', and we know its cardinal # already. # # It follows that #E2 is divisible by r*#E if ' is trivial, and # by r only otherwise. # # Now if ~ denotes the quadratic twist, E2~ = Ek0'~. Now '~ might # be trivial, if and only if d0==2. In this case, #E2~ is # divisible by #E. In all other cases, #E2~ has no known factor. good_twist = [ x for x in which_twist if x[1] % r == 0 ] assert len(good_twist) == 1 action, card = good_twist[0] self._E2 = { 'name': "2", 'card': card, 'd': d0, 'action': action } tcard = 2 * (p**k0+1) - good_twist[0][1] self._twist_E2 = { 'name': "2t", 'card': tcard, 'd': d0 } def set_test_info(self, allowed_cofactor = 1, allowed_size_cofactor = 5, max_trialdiv=10**6, max_B1=10**4): self.allowed_size_cofactor = allowed_size_cofactor self.allowed_cofactor = allowed_cofactor self.max_trialdiv = max_trialdiv self.max_B1 = max_B1 def E2(self, factor=False): self._prepare_E2() tt = self._E2 return self._factor_and_return(tt, factor=factor) def twist_E2(self, factor=False): self._prepare_E2() tt = self._twist_E2 return self._factor_and_return(tt, factor=factor) def _factor_and_return(self, tt, factor=False): if not factor: return tt if tt.has_key("factorization"): if tt["max_B1"] >= self.max_B1 or not tt["factorization"][1]: return tt tt["max_B1"] = self.max_B1 cof = tt["card"] known_primes=[] known_composites=[] if tt is self._twists[0]: cof = cof // self.r elif self._E2['d'] == self.k and tt is self._twists[self._E2['action']]: cof = cof // self.r elif tt is self._E2: if tt['d'] == self.k: # Then E2 is a twist of E, let's delegate our work there! uu = self.twist(tt['action'], factor=True) tt["factorization"] = uu["factorization"] tt["text_factorization"] = uu["text_factorization"] return tt if tt['d'] == 1: # Here we're a direct extension. We know that #E will # divide. cof = cof // self._twists[0]["card"] cof = cof // self.r elif tt is self._twist_E2 and tt['d'] == 1 and self.k % 2 == 1: # then -pi^k is (-pi)^k, so the twist of the extension is the # extension of the twist ! cof = cof // self._twists[1]["card"] elif tt is self._twist_E2 and tt['d'] == 2: cof = cof // self._twists[0]["card"] tt["factorization"] = attempt_to_factor(cof, max_trialdiv=self.max_trialdiv, max_B1=self.max_B1) primes,composites = tt["factorization"] ps = [ ("%d" % pf if pf.nbits() < 30 else "p%d" % pf.nbits()) + (("^" + str(e)) if e>1 else "") for pf,e in primes ] cs = [ ("%d" % pf if pf.nbits() < 30 else "c%d" % pf.nbits()) for pf in composites] tt["text_factorization"] = " * ".join(ps+cs) if tt is self._twists[0]: ps.append((r,1)) tt["text_factorization"] += " * r" elif self._E2['d'] == self.k and tt is self._twists[self._E2['action']]: ps.append((r,1)) tt["text_factorization"] += " * r" elif tt is self._E2: if tt['d'] == 1: uu = self.twist(0, factor=True) ups, ucs = uu["factorization"] ps += ups cs += ucs tt["text_factorization"] += " * (%s)" % uu["text_factorization"] ps.append((r,1)) tt["text_factorization"] += " * r" elif tt is self._twist_E2 and tt['d'] == 1 and self.k % 2 == 1: uu = self.twist(1, factor=True) ups, ucs = uu["factorization"] ps += ups cs += ucs tt["text_factorization"] += " * (%s)" % uu["text_factorization"] elif tt is self._twist_E2 and tt['d'] == 2: uu = self.twist(0, factor=True) ups, ucs = uu["factorization"] ps += ups cs += ucs tt["text_factorization"] += " * (%s)" % uu["text_factorization"] return tt def twist(self, i, factor=False): """ returns the i-th twist (i<2 in the normal case, but can go to i<4 or i<6 depending on the discriminant). INPUT: factor (boolean): whether we want to get the (partially) factored group order. The result is cached and will be reused for later calls. Defaults to False. OUTPUT: a dictionary with the following members: name: a text identifier for which twist we're talking of (one "", "t", "q1", "q3", "s1", "s2", "s4", "s5") card: the group order of the twist factorization: primes and composites in the factorization (see attempt_to_factor) text_factorization: a human-readable version of the latter. max_B1: the latest B1 that was used. """ self._prepare_E2() # I don't think this does a copy -- modifications in the called # functions should be by reference, I believe. tt = self._twists[i] return self._factor_and_return(tt, factor=factor) # "twist-secure" is generally understood as having a cardinal that is # (smallish times a prime). In our application, we'll hardly ever get this, # at least on the original curve. # # "twist-small-subgroup-secure" is a looser check. We still allow a limited # size contribution of small primes, but then we require all remaining # cofactors to be primes larger than the bit length of r # # By default we only check the curve and its quadratic twist, but can do so # for other twists as well. def is_twist_small_subgroup_secure(self, all_twists = False): """ Checks whether the group order is free of small divisors beyond the ones of size allowed_size_cofactor. (Here, small divisor means smaller than the bit length of r.) The number of twists that are checked is normally 2, unless the positional parameter all_twists is set, in which case all defined twists are checked. """ for i in range(len(self._twists) if all_twists else 2): if not self.is_small_subgroup_secure(twist_index=i): return False return True def is_small_subgroup_secure(self, twist_index = 0): """ Checks whether the group order is free of small divisors beyond the ones in allowed_cofactor. (Here, small divisor means smaller than the bit length of r.) """ tt = self.twist(twist_index, factor=True) primes,composites = tt["factorization"] return self._is_safe(primes, composites) def is_G2_small_subgroup_secure(self, twist_index = 0): tt = self.E2(factor=True) primes,composites = tt["factorization"] return self._is_safe(primes, composites) def is_twist_G2_small_subgroup_secure(self, twist_index = 0): tt = self.twist_E2(factor=True) primes,composites = tt["factorization"] return self._is_safe(primes, composites) def _is_safe(self, primes, composites, tolerance=10): if [ p for p in composites if p.nbits() < 2 * self.r.nbits() - tolerance ]: # a composite that is less than twice the bit length of r # surely leads to a prime factor that is smaller than r. return False tp = {p:e for p,e in primes} for p,e in factor(self.allowed_cofactor): if p in tp: tp[p] = max(0,tp[p]-e) allowed_size_cofactor = self.allowed_size_cofactor for p in tp: if tp[p] * ZZ(p).nbits() < allowed_size_cofactor : allowed_size_cofactor -= tp[p] * ZZ(p).nbits() tp[p] = 0 tp = {p:e for p,e in tp.items() if e > 0 and Integer(p).nbits() < self.r.nbits() - tolerance} if tp: return False if composites: # We still have composites, and we can't rule out the # possibility that there are small factors. We couldn't find # them, though. What we want isn't completely clear. return False return True def __repr__(self): """ return an unambiguous representation of the object """ main="CocksPinchVariantResult(%d,%d,%s,%d" % (self.k,self.fD,phex(self.T),self.i) if self.ht != 0: main += ",ht=" + phex(self.ht) if self.hy != 0: main += ",hy=" + phex(self.hy) if self.kl != self.k: main += ",l=%d" % (self.kl // self.k) if self.allowed_cofactor != 1: main += ",allowed_cofactor=%d" % self.allowed_cofactor if self.allowed_size_cofactor != 5: main += ",allowed_size_cofactor=%d" % self.allowed_size_cofactor if self.max_trialdiv != 10**6: main += ",max_trialdiv=%d" % self.max_trialdiv if self.max_B1 != 10**4: main += ",max_B1=%d" % self.max_B1 main += ")" return main def __str__(self): """ return a human-readable representation of the object. This is meant to be copy-pastable into sage for further inspection. """ k = self.k kl = self.kl fD = self.fD T = self.T i = self.i r = self.r y = self.y t = self.t p = self.p ht = self.ht hy = self.hy y0 = self.y0 t0 = self.t0 # Temporarily put more effort into ECM, just for printing. saved_max_B1 = self.max_B1 self.max_B1 = 600 dt0 = t0 - ((T**i+1) % r) y0base = ZZ((t0-2)/sqrt(Integers(r)(-fD))) if r - y0base < y0base: y0base = r - y0base dy0 = y0 - y0base assert dt0 in [0,-r] assert dy0 in [0,-r] s=[ "", "# Cocks-Pinch pairing-friendly curve of embedding degree %d:" % k, "C=" + repr(self), "fD = %d" % fD, "k = %d" % k, "p = 0x%s # %d bits" % (p.hex(), p.nbits()), "rho = %.2f" % float(log(p,r)), "T = 0x%s" % T.hex(), "r = cyclotomic_polynomial(%d)(T)" % kl, "r = 0x%s # %d bits" % (r.hex(), r.nbits()), "i = %d" % i, "t0 = (T**i+1) % r" + (" - r" if dt0==-r else ""), "t0 = " + phex(t0), "ht = " + phex(ht) + " # %d bits, HW=%d, HW_{2-NAF}=%d" % (ht.nbits(), len(bit_positions(ht)), len(bit_positions_2naf(ht))), "t = t0 + ht * r", "t = " + phex(t), "y0 = ZZ((t0-2)/sqrt(Integers(r)(-fD)))" + (" - r" if dy0==-r else ""), "y0 = " + phex(y0), "hy = " + phex(hy) + " # %d bits, HW=%d, HW_{2-NAF}=%d" % (hy.nbits(), len(bit_positions(hy)), len(bit_positions_2naf(hy))), "y = y0 + hy * r", "y = " + phex(y), "is_small_subgroup_secure = %s" % self.is_small_subgroup_secure(), "is_twist_small_subgroup_secure = %s" % self.is_twist_small_subgroup_secure(), "is_G2_small_subgroup_secure = %s" % self.is_G2_small_subgroup_secure(), "is_twist_G2_small_subgroup_secure = %s" % self.is_twist_G2_small_subgroup_secure(), "is_twist_small_subgroup_secure_all = %s" % self.is_twist_small_subgroup_secure(all_twists = True), ] for i in range(len(self._twists)): tt = self.twist(i) s.append("card_E%s = %d" % (tt["name"], tt["card"])) tt = self.E2() s.append("card_E%s = %d" % (tt["name"], tt["card"])) tt = self.twist_E2() s.append("card_E%s = %d" % (tt["name"], tt["card"])) for i in range(len(self._twists)): tt = self.twist(i, factor = True) s.append("# card_E%s = %s" % (tt["name"], tt["text_factorization"])) tt = self.E2(factor = True) s.append("# card_E%s = %s" % (tt["name"], tt["text_factorization"])) tt = self.twist_E2(factor = True) s.append("# card_E%s = %s" % (tt["name"], tt["text_factorization"])) self.max_B1 = saved_max_B1 return "\n".join(s) # search proper class CocksPinchVariantSearch(object): def __init__(self, k, Drange, lambdar, lambdap, verbose=False, output_file="", T_choice="", hty_choice="", max_poly_coeff=0, l=1, required_cofactor=1, allowed_cofactor = 1, allowed_size_cofactor=5, allowed_automatic_cofactor=None, # defaults to allowed_cofactor check_small_subgroup_secure=0, restrict_i=None, parallel=None, seed=None, statefile=None, ): """ INPUT: - k (int, positive): embedding degree - Drange (iterable): discriminants to check (use e.g. [D0, D1] or xsrange(D0, D1)) - lambdar (int, positive): bitzise of prime order subgroup r, typically 256 bits. lambdar-1 < log_2 r < lambdar - lambdap (int, positive): max accepted bitzise of p - T_choice (string): strategy to select T this is a comma-separated list of tokens which may be: - random: abide by the other parameters, except that we do random picks. - hamming=X: restrict to Hamming weight equal to X (int, positive) - 2-naf=X: restrict to 2-NAF weight equal to X (int, positive) - max_poly_coeff (int): maximum value of alpha such that x^k - alpha is irreducible - verbose (boolean): verbose mode - output_file (string): filename where to write the parameters - hty_choice (string): strategy to select hy in y = y0 + hy*r and ht in t = t0 + ht*r; same syntax as T_choice - l (int, positive): consider r = Phi_{k*l}(T) instead of r = Phi_k(T) (if D=4 then l=4/gcd(4,k), if D=3 then l=3/gcd(3,k)) - required_cofactor (int, positive): #E(Fp) = p+1-tr = r*h and required_cofactor | h. Typically required_cofactor = 2,4,2,4,3 to allow Montgomery, Edwards, twisted Edwards, Jacobi Quartic or Hessian representation. - allowed_size_cofactor (int, positive): #E(Fp) = p+1-tr = r*h and the curve is subgroup secure with respect to allowed_size_cofactor, that is, h < 2**allowed_size_cofactor. - restrict_i: examine only these exponents - parallel (string,int,int): (s,i,n) for 0<=i>self.out_file, "# START [T:%s ; hty:%s] %s" % (self.T_choice, self.hty_choice, time.asctime()) self.out_file.flush() class range_class_on_T(object): def __init__(self, search_parent): self.search_parent = search_parent self.k = self.search_parent.k self.kl = self.k * self.search_parent.l self.phi = cyclotomic_polynomial(self.kl) self.prod_pi = get_prod_pi(Integer(self.kl)) # We want: 2^(lambdar-1) <= Phi(T) < 2^lambdar # Phi has no real roots. As a function over the reals, it is # strictly increasing for values >= 1. Therefore we may # simply compute bounds by solving the equation. There has to # be only one root above 1, each time. lambdar = self.search_parent.lambdar RR = RealField(lambdar+10) self.Tmin = ceil((self.phi - 2**(lambdar-1)).roots(RR)[-1][0]) self.Tmax = ceil((self.phi - 2**(lambdar)).roots(RR)[-1][0]) print "Explore T range %d .. %d [%s]" % (self.Tmin, self.Tmax, self.search_parent.T_choice) self.raw_Trange = range_with_strategy( self.Tmin, self.Tmax, strategy=self.search_parent.T_choice, random_state=self.search_parent.random_state) def __iter__(self): return self.iterator(self) class iterator(object): def __init__(self, range_parent): self.range_parent = range_parent self.search_parent = range_parent.search_parent self.it = iter(range_parent.raw_Trange) self.previous_state = self.getstate() def getstate(self): return self.it.getstate() def setstate(self, state): self.previous_state = state return self.it.setstate(state) def __next__(self): return self.next() def next(self): self.previous_state = self.getstate() Tp = self.search_parent.parallel["T"] fail = self.search_parent.fail phi = self.range_parent.phi lambdar = self.search_parent.lambdar while True: T = next(self.it) if Tp and md5int(T) % Tp[1] != Tp[0]: continue if phi(T).nbits() != lambdar: fail['phi(T) bits'] += 1 continue #if self.verbose: # print ("# log_2 T = {}, T = {}".format(T.nbits(), hex(T))) r = phi(T) if r.nbits() > lambdar: break if gcd(r, self.range_parent.prod_pi) != 1: fail['r trialdiv'] += 1 continue if not r.is_pseudoprime() : fail['r milrab'] += 1 continue break # T^k is 1 mod r, so that k divides phi(r)==r-1 assert r % self.search_parent.k == 1 self.search_parent.r = r self.search_parent.K = FiniteField(r, proof = False) # Not clear it makes sense to set this value. Probably not. self.search_parent.T = T return T class range_class_on_D(object): def __init__(self, search_parent): self.search_parent = search_parent self.Drange = sequence_with_stateful_iterators(search_parent.Drange) def __iter__(self): return self.iterator(self) class iterator(object): def __init__(self, range_parent): self.range_parent = range_parent self.search_parent = range_parent.search_parent self.it = iter(self.range_parent.Drange) self.previous_state = self.getstate() self.rf = [ p for p,v in self.search_parent.required_cofactor.factor() if p > 2 ] def __iter__(self): return self def getstate(self): return self.it.getstate() def setstate(self, state): self.previous_state = state return self.it.setstate(state) def __next__(self): return self.next() def next(self): self.previous_state = self.getstate() while True: D = next(self.it) fD = -fundamental_discriminant(-D) if self.search_parent.allowed_automatic_cofactor % 3 != 0: if fD % 3 == 2: # In this case, (t^2+fD*y^2)/4 == t^2-y^2 hence # (t:y) is either (0:1) or (1:0). # If (t:y) = (1:0), one of E and Et is # automatically divisible by 3. # If (t:y) = (0:1), we have p==2 mod 3, and then # p+1\pm t is 0 mod 3 so that both E and Et are # divisible by 3. self.search_parent.fail['unwanted automatic factor ell=3'] += 1 continue if self.search_parent.k % 2 == 1: ell=2*self.search_parent.k+1 if ell.is_prime(): if self.search_parent.allowed_automatic_cofactor % ell != 0: if legendre_symbol(-D, ell) == 1: self.search_parent.fail['unwanted automatic factor ell=2k+1'] += 1 continue # The code in the article requires that we iterate # over fundamental discriminants. As we skip # square-free integers above, we can retain # essentially the same semantics as before. if not self.search_parent.K(-fD).is_square(): self.search_parent.fail['-D not square'] += 1 continue # All odd prime divisors of the cardinal of E must be # such that -D is a square. if prod([1+legendre_symbol(-D,p) for p in self.rf]) == 0: self.search_parent.fail['required_cofactor imposssible'] += 1 continue if not is_squarefree(D): self.search_parent.fail['skip D'] += 1 continue break self.search_parent.fD = fD # Not clear it makes sense to set this value. Probably not. self.search_parent.D = D return D class range_class_on_i(object): def __init__(self, search_parent): self.search_parent = search_parent k = self.search_parent.k l = self.search_parent.l kl = k * l list_i = list({(i*l) % kl for i in range(1,kl) if gcd(i,kl) == 1}) if self.search_parent.restrict_i is not None: list_i = [i for i in list_i if i in self.search_parent.restrict_i] self.list_i = sequence_with_stateful_iterators(list_i) def __iter__(self): return self.iterator(self) class iterator(object): def __init__(self, range_parent): self.range_parent = range_parent self.search_parent = range_parent.search_parent self.it = iter(self.range_parent.list_i) self.previous_state = self.getstate() def __iter__(self): return self def getstate(self): return self.it.getstate() def setstate(self, state): self.previous_state = state return self.it.setstate(state) def __next__(self): return self.next() def next(self): self.previous_state = self.getstate() T = self.search_parent.T fD = self.search_parent.fD r = self.search_parent.r K = self.search_parent.K ZZ = Integers() i = next(self.it) t0 = K(T)**i + 1 y0 = K(t0-2)/sqrt(K(-fD)) # Lift arbitrarily. Anyway we'll iterate over multiple # possible representatives. # The normalisation choice that we do in final_expo_k68 # (at least) is that we use the least positive integer # representative of y0=\pm(t0-2)*inv_sqrt_D # # (as for t0, we have no sign indetermination, so we # simply choose the representative of smallest absolute # value, and that may mean a negative integer) t0 = ZZ(t0) y0 = ZZ(y0) if abs(r-t0) < abs(t0): t0 -= r if abs(r-y0) < abs(y0): y0 -= r y0 = abs(y0) # We want to constrain the bit length of t^2+fD*y^2{{{ # with t = t0 + ht * r and y = y0 + hy * r # to be exactly in the range: # 2^(logp+1)<=(t^2+D*y^2)<2^(logp+2) # that is, we're actually enumerating over # t/2+y/2*sqrt(-D), except that we want that to be an # algebraic integer. Therefore, this entails: # if D==0 mod 4 then t must be even # if D==1 mod 4 then t-y must be even.}}} PP = Integer(2)**(self.search_parent.lambdap+2) # {{{ t1 is an adjusted t0 t1 = t0 rt = r if -fD % 4 == 0: rt = 2*r if t1 % 2 == 1: t1 += r # }}} # t1 and rt are such that # t0 + ht * r is actually constructed as t1 + pre_ht * # rt, which allows ht to be one of 2*pre_ht, or # 2*pre_ht+1, or pre_ht. # # we'll do the same for hy, by the way. # {{{ adjust tmin..tmax # constraints on t only, so that we have a maybe # non-empty range for y # # # 2^(logp+1) <= t^2+D*y^2 < 2^(logp+2) # t^2 < 2^(logp+2) # -sqrt(2^(logp+2)) < t1+pre_ht*rt < sqrt(2^(logp+2)) # (-sqrt(2^(logp+2))-t1)/rt < ht < (sqrt(2^(logp+2))-t1)/rt self.search_parent.pre_htmin = Integer(ceil((-sqrt(PP)-t1)/rt)) self.search_parent.pre_htmax = Integer(ceil( (sqrt(PP)-t1)/rt)) # }}} self.search_parent.rt = rt # self.search_parent.ry = ry self.search_parent.PP = PP self.search_parent.t0 = t0 self.search_parent.t1 = t1 self.search_parent.y0 = y0 # y1 will make sense only when t is chosen. # Not clear it makes sense to set this value. Probably not. self.search_parent.i = i return i class range_class_on_hty_choices(object): def __init__(self, search_parent): """ parse the hty_choice argument, and deduce a sequence of choice strategies that we will process in order """ self.search_parent = search_parent hty_tokens = self.search_parent.hty_choice.split(',') ht_only_tokens = [] hy_only_tokens = [] hty_common_tokens = [] expand_tokens = None for tok in hty_tokens: mm = re.match("(hamming|2-naf)<=(\d+)", tok) if mm: expand_tokens=[] for w in range(int(mm.group(2))+1): tc="%s=%d" % (mm.group(1),w) hc="%s<=%d" % (mm.group(1),int(mm.group(2))-w) expand_tokens.append((tc,hc)) continue mm = re.match("hy:(.*)", tok) if mm: hy_only_tokens.append(mm.group(1)) continue mm = re.match("ht:(.*)", tok) if mm: ht_only_tokens.append(mm.group(1)) continue hty_common_tokens.append(tok) ht_only_tokens += hty_common_tokens hy_only_tokens += hty_common_tokens if not expand_tokens: tc = ",".join(ht_only_tokens) th = ",".join(hy_only_tokens) all_hty_choices = [ (",".join(ht_only_tokens), ",".join(hy_only_tokens)) ] else: all_hty_choices = [] for et,eh in expand_tokens: tc = ",".join([et] + ht_only_tokens) th = ",".join([eh] + hy_only_tokens) all_hty_choices.append((tc, th)) print "Will enumerate cofactors ht and hy based on the following schedule:" for ht_choice, hy_choice in all_hty_choices: pt = ht_choice if not pt: pt = "(normal range order)" py = hy_choice if not py: py = "(normal range order)" if self.search_parent.parallel['ht']: i,n = self.search_parent.parallel['ht'] pt += ", split %d/%d" % (i,n) if self.search_parent.parallel['hy']: i,n = self.search_parent.parallel['hy'] py += ", split %d/%d" % (i,n) print "\t%s and %s" % (pt, py) self.all_hty_choices = sequence_with_stateful_iterators(all_hty_choices) def __iter__(self): return self.iterator(self) class iterator(object): def __init__(self, range_parent): self.range_parent = range_parent self.search_parent = range_parent.search_parent self.it = iter(self.range_parent.all_hty_choices) self.previous_state = self.getstate() def __iter__(self): return self def getstate(self): return self.it.getstate() def setstate(self, state): self.previous_state = state return self.it.setstate(state) def __next__(self): return self.next() def next(self): self.previous_state = self.getstate() # XXX XXX At this point, we know: # self.search_parent.rt // self.search_parent.r (2 or 1) # (self.search_parent.t1 - self.search_parent.t0) // self.search_parent.r (0 or 1) # these two values could lead us to modify the "max" # parameter, and maybe (more broadly) tinker with how we are going to # pass arguments to range_with_strategy pre_ht_choice, pre_hy_choice = next(self.it) # Not clear it makes sense to set this value. Probably not. self.search_parent.pre_ht_choice = pre_ht_choice # Not clear it makes sense to set this value. Probably not. self.search_parent.pre_hy_choice = pre_hy_choice return (pre_ht_choice, pre_hy_choice) class cartesian_product_with_stateful_iterators(object): """ must use ranges that have stateful iterators as well. This is somewhat special as we do not assume that ranges are totally independent, and in fact we mean them not to. Thus we must have generated something from the first iterator in order to call next() on the second one. (the range _objects_ are independent, just the range _iterators_, and precisely their next() functions, depend on one another). This has some implications on how we do next() below. """ def __init__(self, *ranges): self.ranges = ranges def __iter__(self): return self.iterator(self) class iterator(object): def __init__(self, range_parent): self.range_parent = range_parent self.it = [ ] self.current = [ ] def __iter__(self): return self def getstate(self): return [ i.getstate() for i in self.it ] @property def previous_state(self): return [ i.previous_state for i in self.it ] def setstate(self, s): self.current = [] self.it = [ iter(r) for r in self.range_parent.ranges[:len(s)]] for i,ss in enumerate(s): print "Resetting iterator %d with state %s" % (i, ss) self.it[i].setstate(ss) # If previous_state was correctly implemented on all # sub-iterators, then we will not get a StopIteration # # We do not fetch values from the iterators right # now, for two reasons: # 1 - it is important to not fetch from the last # iterator, or we'll end up skipping one tuple from # the iteration. # 2 - it is best to ensure the invariant that # self.current has length either 0 or n. # It's simplest if we hand over the generation to # next() alone. def __next__(self): return self.next() def next(self): n = len(self.range_parent.ranges) if len(self.current) == n: _ = self.current.pop() else: assert not self.current while len(self.current) < n: k = len(self.current) try: if len(self.it) == k: self.it.append(iter(self.range_parent.ranges[k])) self.current.append(next(self.it[k])) except StopIteration: # Then the iterator on the k-th component is # consumed... _ = self.it.pop() # and we're done with our (k-1)-st value, so that # must fetch the next value. If k==0, then this # is the end. if not self.current: raise _ = self.current.pop() return tuple(self.current) def run(self): it = iter(self.range_all) if self.statefile: try: with open(self.statefile, "r") as f: state = pickle.load(f) it.setstate(state) print "Restored state from %s: %s" % (self.statefile, state) except IOError: pass except EOFError: pass for foo in it: if self.statefile: oldstate = it.previous_state with open(self.statefile, "w") as f: pickle.dump(oldstate, f) if self.verbose: print "Saved state to %s: %s" % (self.statefile, oldstate) try: self.subloop() except KeyboardInterrupt: print "Interrupted" return # ss = it.previous_state # it.setstate(ss) # bar = next(it) # assert foo == bar print "DONE [T:%s ; hty:%s] %s" % (self.T_choice, self.hty_choice, time.asctime()) print "generated %d good curves among %d candidate curves in %f s" % (self.ng, self.nc, time.time() - self.time0) pprint.pprint(dict(self.fail)) if self.out_file: print >>self.out_file, "# DONE [T:%s ; hty:%s] %s" % (self.T_choice, self.hty_choice, time.asctime()) print >>self.out_file, "# generated %d good curves among %d candidate curves in %f s" % (self.ng, self.nc, time.time() - self.time0) print >>self.out_file, re.sub("^","# ",pprint.pformat(dict(self.fail)), re.MULTILINE) + "\n" self.out_file.flush() self.out_file.close() def check_automatic_factors(self, E): g = self.gcd_per_curve if g['ncurves'] >= 16: return True n0 = E.twist(0)['card'] // self.r n1 = E.twist(1)['card'] n2 = E.E2()['card'] // self.r if E.E2()['action'] == 0: n2 = n2 // self.r n3 = E.twist_E2()['card'] g['E'] = gcd(g['E'], n0) g['Et'] = gcd(g['Et'], n1) g['E2'] = gcd(g['E2'], n2) g['E2t'] = gcd(g['E2t'], n3) g['all'] = gcd(g['all'], prod([x // gcd(x, self.allowed_automatic_cofactor) for i,x in enumerate([n0,n1,n2,n3])])) g['checked'] = gcd(g['checked'], prod([x // gcd(x, self.allowed_automatic_cofactor) for i,x in enumerate([n0,n1,n2,n3]) if self.check_small_subgroup_secure & (2**i) != 0])) g['ncurves'] += 1 if g['ncurves'] < 16: # True means: keep checking. return True print "Gcd of the curve orders after 16 curves generated [allowed_automatic_cofactor=%d, allowed_cofactor=%d, allowed_size_cofactor=%d, required_cofactor=%d]:" % (self.allowed_automatic_cofactor, self.allowed_cofactor, self.allowed_size_cofactor, self.required_cofactor) pprint.pprint(dict(g)) # Time for a final check. # We don't check the 'all' value, as it doesn't seem to be a nuisance # to have it. Note that guarding against 'all' being divisible by some # factors needs some care, esp. for ell=2*k+1 # # Note however that if we insist on finding curves that avoid small # factors for all of E, E2, Et, E2t, then we must pay attention and do # this work. # # --> use the 'all' marker anyway, because we want the output to # be filtered to produce curves that are valid for all cases. for kk in ['E', 'E2', 'Et', 'E2t', 'all']: if gcd(g[kk], self.allowed_automatic_cofactor) != g[kk]: print "#########################################" print "## Aborting this search, as we have non-trivial algebraic factors on the curves that are larger than allowed_automatic_cofactor (%d):" % self.allowed_automatic_cofactor pprint.pprint(dict(g)) print "#########################################" return False return True def subloop(self): """ This explores only on t and y (actually on ht and hy) This is the place where we could consider balancing the criteria on ht and hy: maybe cap the sum of the Hamming weights, or the sum of the bit sizes. Note that when sqrt(-fD) has no nice algebraic expression in GF(r), then we may choose ht and hy essentially freely (only subject to size constraints). This implies that we're likely to spend a *LOT* of time in this loop, which lessens the win of having made the upper layers of the iteration interruptible. Sure we have a state, but if it's for taking only one big run of the subloop here, we might as well not have any. """ k = self.k T = self.T D = self.D fD = self.fD i = self.i r = self.r K = self.K pre_htmin = self.pre_htmin pre_htmax = self.pre_htmax rt = self.rt # ry = self.ry PP = self.PP t0 = self.t0 t1 = self.t1 y0 = self.y0 pre_ht_choice = self.pre_ht_choice pre_hy_choice = self.pre_hy_choice fail = self.fail print "T=%s i=%d # Explore pre_ht range is (raw) %d .. %d [%s]" % (phex(T), i, pre_htmin, pre_htmax, pre_ht_choice) pre_htgen = itertools.chain( chain_alternate_iterators( range_with_strategy( Integer(0), self.pre_htmax+1, strategy=pre_ht_choice, split=self.parallel['ht'], random_state=self.random_state), itertools.imap(lambda x:-x, range_with_strategy( Integer(1), 1-self.pre_htmin, strategy=pre_ht_choice, split=self.parallel['ht'], random_state=self.random_state) ) ) ) old_range_text = "" q_ht = rt // r r_ht = (t1 - t0) // r for pre_ht in pre_htgen: t = t1 + pre_ht * rt ht = q_ht * pre_ht + r_ht assert t1 + pre_ht * rt == t0 + ht * r if gcd(t-1, self.required_cofactor) != 1: # This might score many many times. fail['required_cofactor impossible'] += 1 continue self.gcd_per_curve = defaultdict(int) # {{{ y1 is an adjusted version of y0 y1 = y0 ry = r if -fD % 4 == 1: ry = 2*r if (t - y0) % 2 == 1: y1 += r # }}} # y1 and ry are such that # y0 + hy * r is actually constructed as y1 + pre_hy * # ry, which allows hy to be one of 2*pre_hy, or # 2*pre_hy+1, or pre_hy. ## {{{ See whether we have reasons that prevent p ## from being prime. ZP=ZZ['x'] x=ZP.gen() p_poly = ZP((t**2 + fD*(y1+x*ry)**2) / 4) if p_poly.content() != 1: fail['p-content'] += 1 continue # p has degree 2. It might be x^2+x mod 2, in which # case it's not good either. no need to do this for # 3, since 3 is larger than the degree if list(p_poly.change_ring(GF(2))) == [0,1,1]: fail['p-algebraic-content'] += 1 continue if not p_poly.is_irreducible(): # should never happen. fail['p not irreducible'] += 1 continue # }}} # {{{ Compute necessary congruence classes for hy # cmodk_set = set((p_poly-1).change_ring(Integers(k)).roots(multiplicities=False)) # As of sage 8.5, integer root finding modulo # composites does not work, as shown by the following # example: # sage: Integers(6)['x']([0,3]).roots(multiplicities=False) # [] cmodk_set = [ ZZ(kk) for kk in Integers(k) if p_poly(kk) == 1 ] assert len(cmodk_set) == len(set(cmodk_set)) if not cmodk_set: fail['p mod k impossible'] += 1 continue kq = k for q in [2,3,5,7]: if k % q == 0: # Then since p is going to be 1 mod k, it # will also be 1 mod q, so there's no point # in constraining the congruence classes mod # k more than we're already doing. continue # otherwise check whether we have roots forbidden = p_poly.roots(GF(q), multiplicities=False) if not forbidden: continue qkeep = [ZZ(v) for v in (set(GF(q))-set(forbidden))] cmodk_set = [ CRT([xk,xq],[kq,q]) for xk in cmodk_set for xq in qkeep ] kq = kq * q assert len(cmodk_set) == len(set(cmodk_set)) # }}} # {{{ adjust hymin and hymax. deduce the generator that we will iterate on # we'll consider y = y1 + hy * ry # we want # 2^(logp-1) <= (t^2+D*y^2)/4 < 2^logp # PP/2 <= (t^2+D*y^2) < PP # PP/2 - t^2 <= D*y^2 < PP - t^2 # sqrt((PP/2 - t^2)/D) - y1 / ry <= pre_hy < sqrt((PP - t^2)/D) - y1 / ry # OR (if y is negative): # -sqrt((PP/2 - t^2)/D) >= y1 + pre_hy * ry > -sqrt(PP - t^2)/D) # (sqrt((PP/2 - t^2)/D) + y1) / ry <= -pre_hy < (sqrt(PP - t^2)/D) + y1) / ry if PP < t**2: continue pre_hymax = 1+floor(((sqrt((PP - t**2)/fD) - y1)/ry)) mpre_hymax = 1+floor(((sqrt((PP - t**2)/fD) + y1)/ry)) if PP/2 < t**2: pre_hygen = range_with_strategy(-mpre_hymax, pre_hymax, strategy=pre_hy_choice, split=self.parallel['hy'], random_state=self.random_state) range_text = "%d .. %d" % (-mpre_hymax, pre_hymax) pre_hysize = pre_hymax+mpre_hymax else: pre_hymin = ceil(((sqrt((PP/2 - t**2)/fD) - y1)/ry)) mpre_hymin = ceil(((sqrt((PP/2 - t**2)/fD) + y1)/ry)) range_text = "%d .. %d + %d .. %d" % (-mpre_hymax, -mpre_hymin, pre_hymin, pre_hymax) pre_hysize = pre_hymax-pre_hymin + mpre_hymax-mpre_hymin # We want to do the negative low-weight before # the positive large-weight, for sure... mm = re.match("(hamming|2-naf)<=(\d+)", pre_hy_choice) if not mm: # Simple things will do ii0 = range_with_strategy( pre_hymin, pre_hymax, strategy=pre_hy_choice, split=self.parallel['hy'], random_state=self.random_state) ii1 = range_with_strategy( mpre_hymin, mpre_hymax, strategy=pre_hy_choice, split=self.parallel['hy'], random_state=self.random_state) ii1 = itertools.imap(lambda x:-x, ii1) pre_hygen = itertools.chain(ii0, ii1) else: # Then we must do something special pre_hygen = iter([]) for w in range(int(mm.group(2))+1): hc="%s=%d" % (mm.group(1),w) ii0 = range_with_strategy( pre_hymin, pre_hymax, strategy=hc, split=self.parallel['hy'], random_state=self.random_state) ii1 = range_with_strategy( mpre_hymin, mpre_hymax, strategy=hc, split=self.parallel['hy'], random_state=self.random_state) ii1 = itertools.imap(lambda x:-x, ii1) pre_hygen = itertools.chain(pre_hygen, ii0, ii1) # }}} if pre_hysize == 0: fail['empty hy range'] += 1 continue if range_text == old_range_text: range_text = "" else: old_range_text = range_text range_text = " # Explore hy range %s [%s]" % (range_text, pre_hy_choice) print "\tT=0x%s, -D=%d, i=%d, ht=%d%s" % (T.hex(), -fD, i, ht, range_text) # print "Looking at %d congruence classes among %d" % (len(cmodk_set), kq) # arrange to print only rarely # if self.verbose and sum(fail.values()) % 32 == 0: if self.verbose: pprint.pprint(dict(fail)) q_hy = ry // r r_hy = (y1 - y0) // r for pre_hy in pre_hygen: if pre_hy % kq not in cmodk_set: # don't even record this. continue hy = q_hy * pre_hy + r_hy y = y1 + ry * pre_hy assert y == y0 + r * hy assert (t**2 + fD*y**2) % 4 == 0 p = (t**2 + fD*y**2)//4 #if self.verbose: # print("# p: {} bits, r: {} bits".format(p.nbits(), r.nbits())) if p.nbits() < (self.lambdap-20): fail['p too small'] += 1 continue if p % k != 1: fail['p mod k'] += 1 continue if ((p+1-t) % self.required_cofactor) != 0: fail['no required cofac'] += 1 continue if p.nbits() > (self.lambdap+1) : fail['p too large'] += 1 continue if gcd(prod_primes,p) > 1: # This should never happen for primes for # which we've explicitly computed the # forbidden congruence classes. # for q,v in gcd(210,p).factor(): # assert kq % q != 0 # fail["%d^%d"%(q,v)] += 1 fail['p trialdiv'] += 1 continue if not p.is_pseudoprime(): fail['p milrab'] += 1 continue if self.max_poly_coeff > 0 : boo = False for alpha in range(1, self.max_poly_coeff): if (x**k - alpha).is_irreducible() : boo = True break if (x**k + alpha).is_irreducible() : boo = True break if not(boo) : fail['large poly'] += 1 E = CocksPinchVariantResult(k,fD,T,i,ht=ht,hy=hy,l=self.l) if not self.check_automatic_factors(E): # There are various effects that condition the # factors that we encounter in the different # cardinals. # # - if an odd prime ell divides #E2 and not #E, then # one of the eigenvalues of the Frobenius is a k-th # root of unity mod ell, whence ell is 1 mod k. (and # then mod 2k). # # - if the following conditions are met: # a) k is odd # b) ell=2*k+1 is prime (and, by (a), -1 qnr mod ell) # c) -D is a square mod ell # then pi=y*sqrt(-D)/2 maps to a square in one of the # residue fields mod ell. Then pi^k is a 2k-th power # mod ell, and therefore is equal to 1, so that ell # divides #E2 (maybe pi maps to 1 in one of these # fields, in which case ell divide #E -- we're not # saying that elll divides #E2 another time in that # case). # # - if t is 0 mod a prime ell, then ell|p+1-t and # ell|p+1+t are equivalent. # # changing t seems sufficient to counter these # effects. But we might as wel reset # self.gcd_per_curve less often (e.g. right at the # beginning of this function), and return right away # to move on to the next i at least, if not the next # D. break self.nc += 1 if self.nc % 65536 == 0: print "generated %d candidate curves in %f s" % (self.nc, time.time() - self.time0) pprint.pprint(dict(fail)) sys.stdout.flush() #if self.verbose: # print repr(E) E.set_test_info(allowed_cofactor = self.allowed_cofactor, allowed_size_cofactor=self.allowed_size_cofactor, max_B1=0) try: csgs = self.check_small_subgroup_secure if (csgs & 1) and not E.is_small_subgroup_secure(): fail['E not sgs'] += 1 continue if (csgs & 2) and not E.is_twist_small_subgroup_secure(): fail['E not twist secure'] += 1 continue if (csgs & 4) and not E.is_G2_small_subgroup_secure(): fail['E not G2 secure'] += 1 continue if (csgs & 8) and not E.is_twist_G2_small_subgroup_secure(): fail['E not twist G2 secure'] += 1 continue # E.set_test_info(max_B1=600) # if not E.is_small_subgroup_secure(): # fail['E not sgs'] += 1 # continue # Is there really a point in continuing with # ECM at this point ? After all, the results # we'll be chiefly interested in will have # been printed at this point anyway. All we # can do at this point is try in vain to # catch a factor which we have basically very # very little chance to catch. # # if not E.is_small_subgroup_secure(): # fail['E not sgs'] += 1 # continue self.ng += 1 print E if self.out_file: print >>self.out_file, E self.out_file.flush() except KeyboardInterrupt: raise except: print "Problem with %s" % repr(E) raise if self.out_file: # print >>self.out_file # print >>self.out_file, "# DONE T=0x%s D=%d i=%d hty_choice=(%s,%s) %s" % (T.hex(), D, i, ht_choice, hy_choice, time.asctime()) self.out_file.flush() def CocksPinchVariant(*args, **kw): """ See CocksPinchVariantSearch """ S = CocksPinchVariantSearch(*args, **kw) S.run()