diff --git a/amicable.sage b/amicable.sage
index 9b3cd71..a9847bd 100644
--- a/amicable.sage
+++ b/amicable.sage
@@ -10,7 +10,7 @@ from itertools import combinations
 
 # p and q should each be ~ L bits.
 
-DEFAULT_TWOADICITY = 21
+DEFAULT_TWOADICITY = 32
 DEFAULT_STRETCH = 0
 
 COEFFICIENT_RANGE = (5,)
@@ -19,14 +19,10 @@ COEFFICIENT_RANGE = (5,)
 ACCEPTABLE_PRIMES = (5,)
 #ACCEPTABLE_PRIMES = Primes()
 
-# <https://eprint.iacr.org/2011/465.pdf>
-# It is well known that if g is neither a square nor a cube in Fp, then all
-# possible group orders an elliptic curve E : y^2 = x^3 + b can have over Fp
-# occur as the order of one of the 6 twists with b \in {1, g, g^2, g^3, g^4, g^5}.
+TWIST_SECURITY = 120
+REQUIRE_PRIMITIVE = True
+REQUIRE_HALFZERO = True
 
-# <https://math.stackexchange.com/questions/127251/when-is-a-not-a-cube-mod-p>:
-#   If p = 2 (mod 3) then all elements are cubes.
-#   If p = 1 (mod 3) then a is a cube iff a^((p-1)/3) = 1.
 
 # <https://cryptojedi.org/papers/pfcpo.pdf> section 2:
 # [...] the order of a curve satisfying the norm equation 3V^2 = 4p - t^2 has one
@@ -39,12 +35,11 @@ ACCEPTABLE_PRIMES = (5,)
 #         = 3(V-1)^2 + 6(V-1) + (t-1)^2 + 2(t-1) + 4
 #       p = 3((V-1)/2)^2 + 3(V-1)/2 + ((t-1)/2)^2 + (t-1)/2 + 1
 #
-# So p-1 will be a multiple of 2^twoadicity, and so will (p+1-t)-1 = (p-1)-(t-1).
+# So p-1 will be a multiple of 2^twoadicity, and so will q-1 for q in
+# { p + 1 - t, p + 1 + (t-3V)/2 }.
 #
-# We'd also like both p and q to be 1 (mod 3), so that we have efficient endomorphisms
-# on both curves. We explicitly check p = 1 (mod 3), and then if t is chosen to be
-# 1 (mod 3) then p+1-t will be 1 (mod 3) (but we must still check q since it
-# is not necessarily that order).
+# We'd also like both p and q to be 1 (mod 6), so that we have efficient endomorphisms
+# on both curves.
 
 def low_hamming_order(L, twoadicity, wid, processes):
     Vlen = (L-1)//2 + 1
@@ -69,16 +64,17 @@ def low_hamming_order(L, twoadicity, wid, processes):
                     if p % 6 == 1 and is_pseudoprime(p):
                         yield (p, t, V)
 
+
 def near_powerof2_order(L, twoadicity, wid, processes):
     trailing_zeros = twoadicity+1
     Vbase = isqrt((1<<(L+1))//3) >> trailing_zeros
-    for Voffset in symmetric_range(10000, base=wid, step=processes):
+    for Voffset in symmetric_range(100000, base=wid, step=processes):
         V = ((Vbase + Voffset) << trailing_zeros) + 1
         assert(((V-1)/2) % (1 << twoadicity) == 0)
         tmp = (1<<(L+1)) - 3*V^2
         if tmp < 0: continue
         tbase = isqrt(tmp) >> trailing_zeros
-        for toffset in symmetric_range(10000):
+        for toffset in symmetric_range(100000):
             t = ((tbase + toffset) << trailing_zeros) + 1
             assert(((t-1)/2) % (1<<twoadicity) == 0)
             if t % 6 != 1:
@@ -87,16 +83,12 @@ def near_powerof2_order(L, twoadicity, wid, processes):
             assert(p4 % 4 == 0)
             p = p4//4
             assert(p % (1<<twoadicity) == 1)
+            if REQUIRE_HALFZERO and p>>(L//2) != 1<<(L - 1 - L//2):
+                continue
+
             if p > 1<<(L-1) and p % 6 == 1 and is_pseudoprime(p):
                 yield (p, t, V)
 
-def find_nonsquare_noncube(p):
-    for g_int in xrange(2, 1000):
-        g = Mod(g_int, p)
-        if g^((p-1)//3) != 1 and g^((p-1)//2) != 1:
-            return g
-    return None
-
 def symmetric_range(n, base=0, step=1):
     for i in xrange(base, n, step):
         yield -i
@@ -109,34 +101,100 @@ def find_nice_curves(strategy, L, twoadicity, stretch, wid, processes):
 
         for (q, qdesc) in ((p + 1 - t,          "p + 1 - t"),
                            (p + 1 + (t-3*V)//2, "p + 1 + (t-3*V)/2")):
-            if q > 1<<(L-1) and q % 6 == 1 and q % (1<<twoadicity) == 1 and is_prime(q) and is_prime(p):
-                bp = find_coefficient(p, q)
+            if REQUIRE_HALFZERO and q>>(L//2) != 1<<(L - 1 - L//2):
+                continue
+
+            if q not in (p, p+1, p-1) and q > 1<<(L-1) and q % 6 == 1 and q % (1<<twoadicity) == 1 and is_prime(q) and is_prime(p):
+                (Ep, bp) = find_curve(p, q)
                 if bp == None: continue
-                bq = find_coefficient(q, p)
+                (Eq, bq) = find_curve(q, p)
                 if bq == None: continue
-                gp = find_nonsquare_noncube(p)
-                if gp == None: continue
-                gq = find_nonsquare_noncube(q)
-                if gq == None: continue
 
-                aq = gq**((q-1)//3)
-                assert(aq**3 == Mod(1, q))
-                ap = gp**((p-1)//3)
-                assert(ap**3 == Mod(1, p))
-                yield (p, q, bp, bq, ap, aq, qdesc)
+                sys.stdout.write('*')
+                sys.stdout.flush()
 
-def find_coefficient(p, q):
+                primp = (Mod(bp, p).multiplicative_order() == p-1)
+                if REQUIRE_PRIMITIVE and not primp: continue
+                primq = (Mod(bq, q).multiplicative_order() == q-1)
+                if REQUIRE_PRIMITIVE and not primq: continue
+
+                twsecp = twist_security(p, q)
+                if twsecp < TWIST_SECURITY: continue
+                twsecq = twist_security(q, p)
+                if twsecq < TWIST_SECURITY: continue
+
+                secp = curve_security(order=q)
+                secq = curve_security(order=p)
+
+                zetap = GF(p).zeta(3)
+                zetap = min(zetap, zetap^2)
+                assert(zetap**3 == Mod(1, p))
+
+                zetaq = GF(q).zeta(3)
+                P = Ep.gens()[0]
+                zP = endo(Ep, zetap, P)
+                if zP != int(zetaq)*P:
+                    zetaq = zetaq^2
+                    assert(zP == int(zetaq)*P)
+                assert(zetaq**3 == Mod(1, q))
+
+                Q = Eq.gens()[0]
+                assert(endo(Eq, zetaq, Q) == int(zetap)*Q)
+
+                embeddivp = embedding_divisor(p, q)
+                embeddivq = embedding_divisor(q, p)
+                twembeddivp = twist_embedding_divisor(p, q)
+                twembeddivq = twist_embedding_divisor(q, p)
+
+                yield (p, q, bp, bq, zetap, zetaq, qdesc, primp, primq, secp, secq, twsecp, twsecq,
+                       embeddivp, embeddivq, twembeddivp, twembeddivq)
+
+def endo(E, zeta, P):
+   (xp, yp) = P.xy()
+   return E(zeta*xp, yp)
+
+def find_curve(p, q):
     for b in COEFFICIENT_RANGE:
         E = EllipticCurve(GF(p), [0, b])
         if E.count_points() == q:
-            return b
-    return None
+            return (E, b)
+    return (None, None)
 
 def find_lowest_prime(p):
     for r in ACCEPTABLE_PRIMES:
         if gcd(p-1, r) == 1:
             return r
 
+pi_12 = (pi/12).numerical_approx()
+
+def curve_security(order):
+    sys.stdout.write('!')
+    sys.stdout.flush()
+    r = factor(order)[-1][0]
+    return log(pi_12 * r, 4)
+
+def twist_security(p, q):
+    return curve_security(2*(p+1) - q)
+
+def embedding_divisor(p, q):
+    sys.stdout.write('#')
+    sys.stdout.flush()
+    assert(gcd(p, q) == 1)
+    Z_q = Integers(q)
+    u = Z_q(p)
+    d = q-1
+    V = factor(d)
+    for (v, k) in V:
+        while d % v == 0:
+            if u^(d/v) != 1: break
+            d /= v
+
+    return (q-1)/d
+
+def twist_embedding_divisor(p, q):
+    return embedding_divisor(p, 2*(p+1) - q)
+
+
 def format_weight(x, detail=True):
     X = format(abs(x), 'b')
     if detail:
@@ -186,20 +244,30 @@ def worker(*args):
         print_exc()
 
 def real_worker(*args):
-    for (p, q, bp, bq, ap, aq, qdesc) in find_nice_curves(*args):
+    for (p, q, bp, bq, zetap, zetaq, qdesc, primp, primq, secp, secq, twsecp, twsecq,
+         embeddivp, embeddivq, twembeddivp, twembeddivq) in find_nice_curves(*args):
         output  = "\n"
         output += "p   = %s\n" % format_weight(p)
         output += "q   = %s\n" % format_weight(q)
         output += "    = %s\n" % qdesc
-        output += "α_p = %s (mod p)\n" % format_weight(int(ap), detail=False)
-        output += "α_q = %s (mod q)\n" % format_weight(int(aq), detail=False)
+        output += "ζ_p = %s (mod p)\n" % format_weight(int(zetap), detail=False)
+        output += "ζ_q = %s (mod q)\n" % format_weight(int(zetaq), detail=False)
 
-        output += "Ep/Fp : y^2 = x^3 + %d (%ssquare)\n" % (bp, "" if Mod(bp, p).is_square() else "non")
-        output += "Eq/Fq : y^2 = x^3 + %d (%ssquare)\n" % (bq, "" if Mod(bq, q).is_square() else "non")
+        output += "Ep/Fp : y^2 = x^3 + %d\n" % (bp,)
+        output += "Eq/Fq : y^2 = x^3 + %d\n" % (bq,)
 
         output += "gcd(p-1, %d) = 1\n" % find_lowest_prime(p)
         output += "gcd(q-1, %d) = 1\n" % find_lowest_prime(q)
 
+        output += "%d is %ssquare and %sprimitive in Fp\n" % (bp, "" if Mod(bp, p).is_square() else "non", "" if primp else "non")
+        output += "%d is %ssquare and %sprimitive in Fq\n" % (bq, "" if Mod(bp, q).is_square() else "non", "" if primq else "non")
+
+        output += "Ep security = %.1f, embedding degree = (q-1)/%d\n" % (secp, embeddivp)
+        output += "Eq security = %.1f, embedding degree = (p-1)/%d\n" % (secq, embeddivq)
+
+        output += "Ep twist security = %.1f, embedding degree = (2p + 1 - q)/%d\n" % (twsecp, twembeddivp)
+        output += "Eq twist security = %.1f, embedding degree = (2q + 1 - p)/%d\n" % (twsecq, twembeddivq)
+
         print(output)  # one syscall to minimize tearing
 
 main()