work

4 files changed, 393 insertions, 265 deletions

diff --git a/README.md b/README.md
index c0e920b..d93012d 100644
--- a/README.md
+++ b/README.md
@@ -1 +1,3 @@
 # forbidden salamanders
+
+The cryptography in this library is not audited and should not be used in production.
diff --git a/aesgcmanalysis.py b/aesgcmanalysis.py
index db919e1..0b81522 100644
--- a/aesgcmanalysis.py
+++ b/aesgcmanalysis.py
@@ -1,10 +1,12 @@
-import random, struct, hmac
+import random, struct, hmac, itertools
 from Crypto.Cipher import AES
 import numpy as np
 
 ## Computation in GF(2^128)/(x^128 + x^7 + x^2 + x^1 + 1)
+## Elements are represented as an integer n where (n & (1 << i)) is the coefficient for x^i.
 
 class X:
+    """A helper class for specifying elements of GF(2^n)."""
     def __pow__(self, y):
         return 1 << y
 x = X()
@@ -54,7 +56,7 @@ def gf128_mul(a, b):
     return p
 
 def gf128_divmod(a, b):
-    """Compute q, r such that a = qb + r.
+    """Compute q, r such that a = qb + r; gf128_deg(r) < gf128_deg(b).
     >>> assert gf128_divmod(x**20 + x**5 + 1, x**13 + x**2) == (x**7, x**9 + x**5 + 1)
     >>> assert gf128_divmod(x**90, x**80) == (x**10, 0)
     """
@@ -127,87 +129,131 @@ def gf128_exp(a, n):
         rec = gf128_mul(rec, a)
     return rec
 
+def gf128_sqrt(x):
+    """Compute y such that y^2 = x: the inverse Frobenius automorphism.
+    Reference: https://math.stackexchange.com/questions/943417/square-root-for-galois-fields-gf2m
+    >>> assert gf128_sqrt(1) == 1
+    >>> assert gf128_sqrt(x**2 + 1) == x**1 + 1
+    >>> assert gf128_sqrt(x**3 + 1) == x**126 + x**123 + x**120 + x**117 + x**114 + x**111 + x**108 + x**105 + x**102 + x**99 + x**96 + x**93 + x**90 + x**87 + x**84 + x**81 + x**78 + x**75 + x**72 + x**69 + x**66 + x**63 + x**62 + x**60 + x**59 + x**57 + x**56 + x**54 + x**53 + x**51 + x**50 + x**48 + x**47 + x**45 + x**44 + x**42 + x**41 + x**39 + x**38 + x**36 + x**35 + x**33 + x**32 + x**30 + x**29 + x**27 + x**26 + x**24 + x**23 + x**21 + x**20 + x**18 + x**17 + x**15 + x**14 + x**12 + x**11 + x**9 + x**8 + x**6 + x**3 + 1
+    """
+    return gf128_exp(x, 2**127)
+
 ## Computation in GF(2^128)[X]/(x^128 + x^7 + x^2 + x^1 + 1)
+## Elements are represented as arrays where the ith element is the coefficient for x^i.
+
+def gf128poly_lead(x):
+    """Compute the leading coefficient of x.
+    >>> assert gf128poly_lead([x**3, x**1, 1]) == 1
+    >>> assert gf128poly_lead([x**3, x**1, x**2, x**4]) == x**4
+    >>> assert gf128poly_lead([1]) == 1
+    """
+    assert len(x) > 0
+    return x[-1]
 
-# Get rid of trailing zeros?
 def collapse(x):
     if len(x) == 0:
-        return x
-    if x[-1] == 0:
+        return []
+    if gf128poly_lead(x) == 0:
         return collapse(x[:-1])
-# CAN WE REMOVE MOD?
-    return [gfmod(y, m) for y in x]
+    return [gf128_mod(y, gcm_modulus) for y in x]
 
-def gf128polyring_deg(x):
+def gf128poly_deg(x):
+    """Compute the degree of x.
+    >>> assert gf128poly_deg([x**3, x**1, 1]) == 2
+    >>> assert gf128poly_deg([x**2, x**3]) == 1
+    >>> assert gf128poly_deg([x**15]) == 0
+    >>> assert gf128poly_deg([1]) == 0
+    >>> assert gf128poly_deg([]) == -1
+    """
     return len(x) - 1
-# assert gf128deg([x**1, x**1, x**1]) == 2
-def gf128polyring_add(a, b):
-    return collapse([gfadd(x, y) for x, y in zip_longest(a, b, fillvalue=0)])
-# assert gf128add([x**1, x**1, x**1, x**1], [x**2, x**3, x**4]) == [x**1+x**2, x**1+x**3, x**1+x**4, x**1]
-def gf128polyring_sub(a, b):
-    return gf128add(a, b)
-
-def gf128polyring_cmul(c, a):
-    return collapse([gfmul(x, c) for x in a])
-
-def gf128polyring_mul(a, b):
-    a = a.copy()
-    b = b.copy()
-    
+
+def gf128poly_add(a, b):
+    """Compute a + b.
+    >>> assert gf128poly_add([x**3, x**1, 1], [x**2 + 1, x**3]) == [x**3 + x**2 + 1, x**1 + x**3, 1]
+    >>> assert gf128poly_add([x**2, x**3], [x**5, x**12, x**2, x**4]) == [x**5 + x**2, x**12 + x**3, x**2, x**4]
+    >>> assert gf128poly_add([x**15], []) == [x**15]
+    >>> assert gf128poly_add([x**15], [x**15]) == []
+    """
+    return collapse([gf128_add(x, y) for x, y in itertools.zip_longest(a, b, fillvalue=0)])
+
+def gf128poly_mul(a, b):
+    """Compute ab.
+    >>> assert gf128poly_mul([x**3, x**1, 1], [x**2 + 1, x**3]) == [x**5 + x**3, x**6 + x**3 + x**1, x**4 + x**2 + 1, x**3]
+    >>> assert gf128poly_mul([x**2, x**3], [x**5, x**12, x**2, x**4]) == [x**7, x**14 + x**8, x**15 + x**4, x**6 + x**5, x**7]
+    >>> assert gf128poly_mul([x**15], []) == []
+    >>> assert gf128poly_mul([x**15], [1]) == [x**15]
+    >>> assert gf128poly_mul([0, x**15], [0, x**15]) == [0, 0, x**30]
+    >>> assert gf128poly_mul([x**70], [1, x**70, 1]) == [x**70, x**19 + x**14 + x**13 + x**12, x**70]
+    """
     p = []
-    while gf128deg(a) >= 0:
-        p = gf128add(p, gf128cmul(a[0], b))
+    while gf128poly_deg(a) >= 0:
+        p = gf128poly_add(p, [gf128_mul(x, a[0]) for x in b])
         a = a[1:]
         b = [0] + b
     return collapse(p)
 
-def gf128polyring_divmod(a, b):
-    q, r = [], []
-    while gf128deg(a) >= gf128deg(b):
-        highA = a[-1]
-        highB = b[-1]
-        highBinv = 1 if highB == 1 else gfmodinv(highB, m)
-        #print(highA, highBinv)
-        qq = gfmodmul(highA, highBinv, m)
-        #qq, rr = gfdivmod(highA, highB)
-        #assert rr == 0
-        de = gf128deg(a) - gf128deg(b)
+def gf128poly_divmod(a, b):
+    """Compute q, r such that a = qb + r; gf128poly_deg(r) < gf128poly_deg(b).
+    >>> assert gf128poly_divmod([x**5 + x**3, x**6 + x**3 + x**1, x**4 + x**2 + 1, x**3], [x**2 + 1, x**3]) == ([x**3, x**1, 1], [])
+    >>> assert gf128poly_divmod([x**7, x**14 + x**8, x**15 + x**4, x**6 + x**5, x**7], [x**2, x**3]) == ([x**5, x**12, x**2, x**4], [])
+    >>> assert gf128poly_divmod([x**15], [1]) == ([x**15], [])
+    >>> assert gf128poly_divmod([x**15], [x**15]) == ([1], [])
+    >>> assert gf128poly_divmod([0, 0, x**30], [0, x**15]) == ([0, x**15], [])
+    >>> assert gf128poly_divmod([x**70, x**19 + x**14 + x**13 + x**12, x**70], [1, x**70, 1]) == ([x**70], [])
+    >>> assert gf128poly_divmod([x**18, x**18, x**3], [x**3, x**3]) == ([x**15 + 1, 1], [x**3])
+    """
+    assert b != [], 'divide by zero'
+    q = []
+    while gf128poly_deg(a) >= gf128poly_deg(b):
+        highA = gf128poly_lead(a)
+        highB = gf128poly_lead(b)
+        highBinv = 1 if highB == 1 else gf128_inv(highB)
+        qq = gf128_mul(highA, highBinv)
+        de = gf128poly_deg(a) - gf128poly_deg(b)
         new_term = [0]*(de+1)
         new_term[de] = qq
-        q = gf128add(q, new_term)
-        fact = gf128mul(b, new_term)
-        a = gf128sub(a, fact)
-        
+        q = gf128poly_add(q, new_term)
+        fact = gf128poly_mul(b, new_term)
+        a = gf128poly_add(a, fact)
     return collapse(q), collapse(a)
 
-def InverseFrobeniusAutomorphismGF128(x):
-    return gfmodexp(x, 2**127, m) # https://math.stackexchange.com/questions/943417/square-root-for-galois-fields-gf2m
-
-def gf128polyring_sqrt(p):
-    p = p.copy()
+def gf128poly_sqrt(p):
+    """Compute q such that q^2 = p for p such that p' = 0.
+    >>> assert gf128poly_sqrt([1, 0, 1, 0, 1, 0, 1]) == [1, 1, 1, 1]
+    >>> assert gf128poly_sqrt([1, 0, 1, 0, 0, 0, 1]) == [1, 1, 0, 1]
+    >>> assert gf128poly_sqrt([1, 0, x**4]) == [1, x**2]
+    >>> assert gf128poly_sqrt([1, 0, x**3]) == [1, x**126 + x**123 + x**120 + x**117 + x**114 + x**111 + x**108 + x**105 + x**102 + x**99 + x**96 + x**93 + x**90 + x**87 + x**84 + x**81 + x**78 + x**75 + x**72 + x**69 + x**66 + x**63 + x**62 + x**60 + x**59 + x**57 + x**56 + x**54 + x**53 + x**51 + x**50 + x**48 + x**47 + x**45 + x**44 + x**42 + x**41 + x**39 + x**38 + x**36 + x**35 + x**33 + x**32 + x**30 + x**29 + x**27 + x**26 + x**24 + x**23 + x**21 + x**20 + x**18 + x**17 + x**15 + x**14 + x**12 + x**11 + x**9 + x**8 + x**6 + x**3]
+    """
+    assert gf128poly_formal_derivative(p) == []
     q = []
     for x in p:
         if x == 0 or x == 1:
             q.append(x)
         else:
-            inv = InverseFrobeniusAutomorphismGF128(x)
-            q.append(inv)
-    #print(q)
-    p = q
-    for i in range(2,len(p),2):
-        p[i//2] = gfadd(p[i//2], p[i])
-        p[i] = 0
-    return collapse(p)
+            q.append(gf128_sqrt(x))
+    for i in range(2, len(q), 2):
+        q[i//2] = q[i]
+        q[i] = 0
+    return collapse(q)
 
-def gf128polyring_monic(p):
-    p = p.copy()
-    finalCoeff = p[-1]
-    if finalCoeff == 1:
+def gf128poly_monic(p):
+    """Compute p/p_lead for p such that gf128poly_deg(p) >= 0.
+    >>> assert gf128poly_monic([x**54, x**2]) == [x**52, 1]
+    >>> assert gf128poly_monic([x**2, x**3]) == [x**127 + x**6 + x**1 + 1, 1]
+    """
+    assert len(p) > 0
+    lead = gf128poly_lead(p)
+    if lead == 1:
         return p
-    inv = gfmodinv(finalCoeff, m)
-    return [gfmodmul(x, inv, m) for x in p]
-
-def gf128polyring_formal_derivative(p):
+    lead_inv = gf128_inv(lead)
+    return [gf128_mul(x, lead_inv) for x in p]
+
+def gf128poly_formal_derivative(p):
+    """Compute p'.
+    >>> assert gf128poly_formal_derivative([]) == []
+    >>> assert gf128poly_formal_derivative([x**25]) == []
+    >>> assert gf128poly_formal_derivative([x**25, x**3 + x**1, x**10, x**7]) == [x**3 + x**1, 0, x**7]
+    """
     q = []
     for i, x in enumerate(p):
         if i == 0:
@@ -219,128 +265,158 @@ def gf128polyring_formal_derivative(p):
         q.append(e)
     return collapse(q)
 
-def gf128polyring_sff(f):
+def gf128poly_modexp(a, n, m):
+    """Compute a^n (mod m).
+    >>> assert gf128poly_modexp([x**2, x**3, x**4], 0, [x**3, x**2, 1]) == [1]
+    >>> assert gf128poly_modexp([x**2, x**3, x**4], 1, [x**3, x**2, 1, x**5]) == [x**2, x**3, x**4]
+    >>> assert gf128poly_modexp([x**2, x**3, x**4], 5, [x**3, x**2, 1]) == [x**39 + x**38 + x**37 + x**36 + x**35 + x**33 + x**32 + x**30 + x**27 + x**26 + x**25 + x**24 + x**21 + x**20 + x**15 + x**10, x**38 + x**36 + x**35 + x**33 + x**32 + x**31 + x**26 + x**24 + x**23 + x**22 + x**21 + x**20 + x**14 + x**11]
+    """
+    def modmul(a, b, m):
+        q, r = gf128poly_divmod(gf128poly_mul(a, b), m)
+        return r
+
+    if n == 0:
+        return [1]
+    rec = gf128poly_modexp(a, n//2, m)
+    rec = modmul(rec, rec, m)
+    if n % 2 == 1:
+        rec = modmul(rec, a, m)
+    return rec
+
+def gf128poly_monic_gcd(a, b):
+    """Compute g such that g | a, g | b, gf128poly_lead(g) = 1;
+	for any h such that h | a and h | b, gf128poly_deg(g) >= gf128poly_deg(h).
+    >>> assert gf128poly_monic_gcd([x**3, x**2, x**4], [x**3, x**2, x**4]) == [x**127 + x**6 + x**1 + 1, x**127 + x**126 + x**6 + x**5 + x**1, 1]
+    >>> assert gf128poly_monic_gcd([1], [1]) == [1]
+    >>> assert gf128poly_monic_gcd([x**5 + x**2], []) == [x**5 + x**2]
+    >>> assert gf128poly_monic_gcd([x**3, 0, 0, x**6], [x**3, x**8 + x**4, x**9 + x**7, x**8]) == [x**127 + x**6 + x**1 + 1, 1]
+    """
+    assert gf128poly_deg(a) >= 0 or gf128poly_deg(b) >= 0
+    if gf128poly_deg(a) < 0:
+        return b
+    if gf128poly_deg(b) < 0:
+        return a
+    q, r = gf128poly_divmod(b, a)
+    ret = gf128poly_monic_gcd(r, a)
+    return gf128poly_monic(ret)
+
+def gf128poly_square_free_factorization(f):
+    """Compute the square-free factorization of a monic polynomial f.
+    Reference: https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields#Square-free_factorization
+    >>> a = [x**1, x**2, 1]
+    >>> b = [x**2, x**3+1, 1]
+    >>> c = [x**3, x**4, x**5, 1]
+    >>> asq = gf128poly_mul(a, a)
+    >>> acb = gf128poly_mul(asq, a)
+    >>> bsq = gf128poly_mul(b, b)
+    >>> acbbsq = gf128poly_mul(acb, bsq)
+    >>> p = gf128poly_mul(acbbsq, c)
+    >>> assert gf128poly_square_free_factorization([1]) == [([1], 1)]
+    >>> assert gf128poly_square_free_factorization(p) == [([x**3, x**4, x**5, 1], 1), ([x**1, x**2, 1], 3), ([x**2, x**3 + 1, 1], 2)]
+    """
     R = []
     assert f[-1] == 1, 'monic'
-    fprime = gf128formal_derivative(f)
-    #print('fprime', fprime)
-    c = gf128gcd(f, fprime)
-    #print('c', c)
-    w, r = gf128divmod(f, c)
-    #print('w', w)
+    if f == [1]:
+        return [([1], 1)]
+    fprime = gf128poly_formal_derivative(f)
+    c = gf128poly_monic_gcd(f, fprime)
+    w, r = gf128poly_divmod(f, c)
     assert r == []
     i = 1
     while w != [1]:
-        y = gf128gcd(w, c)
-        #print('y', i, y)
-        fac, r = gf128divmod(w, y)
-        #print('fac', i, fac)
+        y = gf128poly_monic_gcd(w, c)
+        fac, r = gf128poly_divmod(w, y)
         assert r == []
         R.append((fac, i))
         w = y
-        c, r = gf128divmod(c, y)
-        #print('c', i, c, r)
+        c, r = gf128poly_divmod(c, y)
         assert r == []
         i = i + 1
-    if c != [1]: 
-        #print('sqrt of', ds(c))
-        c = gf128sqrt(c)
-        c = gf128monic(c)
-        #print('rec', ds(c))
-        Rrec = gf128sff(c)
+    if c != [1]:
+        c = gf128poly_sqrt(c)
+        c = gf128poly_monic(c)
+        Rrec = gf128poly_square_free_factorization(c)
         R += [(fac, i*2) for (fac, i) in Rrec]
-    return R
-    # Frobenius automorphism???
-
-def gf128polyring_modmul(a, b, m):
-    c = gf128mul(a, b)
-    q, r = gf128divmod(c, m)
-    return r
-def gf128polyring_modexp(a, n, m):
-    if n == 0:
-        return [1]
-    rec = gf128modexp(a, n//2, m)
-    rec2 = gf128modmul(rec, rec, m)
-    if n % 2 == 1:
-        rec2 = gf128modmul(rec2, a, m)
-    return rec2
-def ddf(f):
+    return [(fac, mult) for fac, mult in R if fac != [1]]
+
+def gf128poly_distinct_degree_factorization(f):
+    """Compute the distinct-degree factorization of a square-free monic polynomial f.
+    Reference: https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields#Distinct-degree_factorization
+    >>> a = [x**1, 1]
+    >>> b = [x**2, 1]
+    >>> ab = gf128poly_mul(a, b)
+    >>> assert gf128poly_distinct_degree_factorization(ab) == {((x**3, x**2 + x**1, 1), 1)}
+    """
     i = 1
     S = set()
     f_ = f.copy()
-    
-    while gf128deg(f_) >= 2 * i:
-        qq = gf128modexp([0, 1], 2**(128*i), f_)
-        qq = gf128add(qq, [0, 1])
-        g = gf128gcd(f_, qq)
-        #print('gcd', ds(f_), '||', qq, '||', ds(g))
+    while gf128poly_deg(f_) >= 2*i:
+        qq = gf128poly_modexp([0, 1], 2**(128*i), f_)
+        qq = gf128poly_add(qq, [0, 1])
+        g = gf128poly_monic_gcd(f_, qq)
         if g != [1]:
             S.add((tuple(g), i))
-            q, r = gf128divmod(f_, g)
+            q, r = gf128poly_divmod(f_, g)
             assert r == []
             f_ = q
         i += 1
     if f_ != [1]:
-        #print('last')
-        S.add((tuple(f_), gf128deg(f_)))
+        S.add((tuple(f_), gf128poly_deg(f_)))
     if len(S) == 0:
         return {(tuple(f), 1)}
     else:
         return S
 
-exp = (2**128-1)//3
-def gfrand():
-    return random.randint(0, 2**128-1)
-def gf128polyring_rand(degree):
-    qq= [gfrand() for _ in range(degree+1)]
-    return qq
-def edf(f, d):
-    n = gf128deg(f)
+def gf128poly_equal_degree_factorization(f, d):
+    """Compute the equal-degree factorization of a square-free monic polynomial f whose factors are all of degree d.
+    Uses the Cantor-Zassenhaus algorithm.
+    Reference: https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields#Equal-degree_factorization
+    >>> a = [x**1, 1]
+    >>> b = [x**2, 1]
+    >>> ab = gf128poly_mul(a, b)
+    >>> assert gf128poly_equal_degree_factorization(ab, 1) == {(x**2, 1), (x**1, 1)}
+    """
+    def gf128_rand():
+        return random.randint(0, 2**128-1)
+
+    def gf128poly_rand(degree):
+        return [gf128_rand() for _ in range(degree+1)]
+
+    exp = (2**128-1)//3
+    n = gf128poly_deg(f)
     r = n//d
     S = {tuple(f)}
     while len(S) < r:
-        h = gf128rand(n-1)
-        g = gf128gcd(h, f)
+        h = gf128poly_rand(n - 1)
+        g = gf128poly_monic_gcd(h, f)
         if g == [1]:
-            g = gf128add(gf128modexp(h, exp, f), [1])
+            g = gf128poly_add(gf128poly_modexp(h, exp, f), [1])
         for u in list(S):
             u = list(u)
-            if gf128deg(u) == d:
+            if gf128poly_deg(u) == d:
                 continue
-            gugcd = gf128gcd(g, u)
+            gugcd = gf128poly_monic_gcd(g, u)
             if gugcd != [1] and gugcd != u:
                 S.remove(tuple(u))
                 S.add(tuple(gugcd))
-                qq, rr = gf128divmod(u, gugcd)
+                qq, rr = gf128poly_divmod(u, gugcd)
                 assert rr == []
                 S.add(tuple(qq))
-        #print("Iteration", S)
     return S
 
-def factorize(f):
-    #print('Computing factors over GF(2^128)/[y](x^128+x^7+x^2+x+1) of', ds(f))
+def gf128poly_factorize(f):
+    """Compute the factors of a polynomial f. Does not return multiplicity.
+    """
     factors = set()
-    f = gf128monic(f)
-    #print('Reduced to monic polynomial', ds(f))
-    fs = gf128sff(f)
+    f = gf128poly_monic(f)
+    fs = gf128poly_square_free_factorization(f)
     for p, _ in fs:
-        #print('Computed square-free factor', ds(p))
-        qs = ddf(p)
+        qs = gf128poly_distinct_degree_factorization(p)
         for q, d in qs:
-            #print('Computed distinct-degree factor', ds(q), 'of degree', d)
-            rs = edf(list(q), d)
-            #for r in rs:
-            #   print('Computed equal-degree factor', ds(r))
+            rs = gf128poly_equal_degree_factorization(list(q), d)
             factors |= rs
     return factors
 
-def monic_roots(factors):
-    roots = []
-    for factor in factors:
-        if gf128deg(factor) == 1:
-            roots.append(factor[0])
-    return roots
-
 ## AES-GCM
 
 def reverse_mask(b):
@@ -371,9 +447,17 @@ def gf128_to_bytes(n):
     assert len(s) == 16
     return s
 
+def pad(xs):
+    return xs + b'\x00' * ((16 - len(xs)) % 16)
+
+def build_blocks(aad, c):
+    padded_aad = pad(aad)
+    padded_c = pad(c)
+    length = struct.pack('>QQ', len(aad)*8, len(c)*8)
+    return padded_aad + padded_c + length
+
 def gmac(h, s, aad, c):
-    assert len(n) == 12
-    blocks = build_blockseq(ad, c)
+    blocks = build_blocks(aad, c)
     g = 0
     for i in range(len(blocks)//16):
         block = bytes_to_gf128(blocks, i*16, (i + 1)*16)
@@ -394,7 +478,7 @@ def gctr(k, nonce, x):
     ctr = (int.from_bytes(nonce, 'big') << 32) + 2
     i = 0
     while len(y) < len(x):
-        c += xor(x[i*16:(i + 1)*16], ecb_encrypt(k, ctr.to_bytes(16, 'big')))
+        y += xor(x[i*16:(i + 1)*16], ecb_encrypt(k, ctr.to_bytes(16, 'big')))
         ctr += 1
         i += 1
     return y
@@ -402,7 +486,7 @@ def gctr(k, nonce, x):
 def authentication_key(k):
     return bytes_to_gf128(ecb_encrypt(k, b'\x00'*16))
 
-def blind(k, n):
+def blind(k, nonce):
     return bytes_to_gf128(ecb_encrypt(k, nonce + b'\x00\x00\x00\x01'))
 
 def gcm_encrypt(k, nonce, aad, m, mac_bytes=16):
@@ -413,70 +497,68 @@ def gcm_encrypt(k, nonce, aad, m, mac_bytes=16):
 def gcm_decrypt(k, nonce, aad, c, mac, mac_bytes=16):
     m = gctr(k, nonce, c)
     expected_mac = gmac(authentication_key(k), blind(k, nonce), aad, c)[:mac_bytes]
-    assert hmac.compare_digest(mac, exp_mac), f'Expected MAC {exp_mac.hex()}, got {mac.hex()}'
+    assert hmac.compare_digest(mac, expected_mac)
     return m
 
 ## Forbidden Attack
 
-def compute_polydiff(ad1, ct1, mac1, ad2, ct2, mac2):
-    bs1 = build_blockseq(ad1, ct1)
-    bs2 = build_blockseq(ad2, ct2)
+def compute_forbidden_polynomial(aad1, aad2, c1, c2, mac1, mac2):
+    bs1 = build_blocks(aad1, c1)
+    bs2 = build_blocks(aad2, c2)
     if len(bs1) < len(bs2):
         bs1 = b'\x00'*(len(bs2)-len(bs1)) + bs1
     else:
         bs2 = b'\x00'*(len(bs1)-len(bs2)) + bs2
     assert len(bs1) == len(bs2)
-    #print(bs1.hex())
-    #print(bs2.hex())
     f = []
     N = len(bs1)//16
     for i in range(N):
-        b1 = binToGF128(bs1[i*16:(i+1)*16])
-        b2 = binToGF128(bs2[i*16:(i+1)*16])
-        c = gfadd(b1, b2)
-        f.append(c)
-    mac1gf = binToGF128(mac1)
-    mac2gf = binToGF128(mac2)
-    macSgf = gfadd(mac1gf, mac2gf)
-    f.append(macSgf)
+        b1 = bytes_to_gf128(bs1[i*16:(i+1)*16])
+        b2 = bytes_to_gf128(bs2[i*16:(i+1)*16])
+        f.append(gf128_add(b1, b2))
+    f.append(gf128_add(bytes_to_gf128(mac1), bytes_to_gf128(mac2)))
     return list(reversed(f))
 
-def forge_gcm_tag(h, origAd, origCt, newAd, newCt, tag, n):
-    # First recover s for the nonce
-    tag_ = gmac(h, 0, origAd, origCt, n)
-    s = gfadd(binToGF128(tag), binToGF128(tag_))
-    tag = gmac(h, s, newAd, newCt, n)
-    return tag
-def gcm_repeatednonce_attack(factors, n, ad1, ct1, mac1, ad2, ct2, mac2):
-    #f=compute_polydiff(ad1, ct1, mac1, ad2, ct2, mac2)
-    #factors = factorize(f)
-    roots = monic_roots(factors)
-    h = None
-    tags = []
-    for root in roots:
-        forged_tag = forge_gcm_tag(root, ad1, ct1, ad2, ct2, mac1, n)
-        tags.append(forged_tag)
-    return tags
-# tag_guesses = gcm_repeatednonce_attack(factors, n, ad1, ct1, mac1, ad2, ct2, mac2)
-
-
-def d(n):
-    if n == 0:
-        return '0'
-    s = []
-    k = 0
-    while n > 0:
-        if n & 1 == 0:
-            k += 1
-            n >>= 1
-            continue
-        if k == 0:
-            s.append('1')
-        elif k == 1:
-            s.append('x')
-        else:
-            s.append('x^' + str(k))
-        n >>= 1
-        k += 1
-    return ' + '.join(reversed(s))
-
+def nonce_reuse_recover_secrets(nonce, aad1, aad2, c1, c2, mac1, mac2):
+    f = compute_forbidden_polynomial(aad1, aad2, c1, c2, mac1, mac2)
+    factors = gf128poly_factorize(f)
+    secrets = []
+    for factor in factors:
+        if gf128poly_deg(factor) == 1:
+            h = factor[0]
+            zero_tag = gmac(h, 0, aad1, c1)
+            s = gf128_add(bytes_to_gf128(mac1), bytes_to_gf128(zero_tag))
+            secrets.append((h, s))
+    return secrets
+
+def demo():
+    k = b"tlonorbistertius"
+    nonce = b"jorgelborges"
+    m1 = b"The universe (which others call the Library)"
+    aad1 = b"The Anatomy of Melancholy"
+    m2 = b"From any of the hexagons one can see, interminably"
+    aad2 = b"Letizia Alvarez de Toledo"
+
+    c1, mac1 = gcm_encrypt(k, nonce, aad1, m1)
+    c2, mac2 = gcm_encrypt(k, nonce, aad2, m2)
+
+    # Recover the authentication key and blind from public information
+    possible_secrets = nonce_reuse_recover_secrets(nonce, aad1, aad2, c1, c2, mac1, mac2)
+
+    # Forge the ciphertext
+    m_forged = b"As was natural, this inordinate hope"
+    assert len(m_forged) <= len(m1)
+    c_forged = xor(c1, xor(m1, m_forged))
+    aad_forged = b"You who read me, are You sure of understanding my language?"
+
+    # Check possible candidates for authentication key
+    succeeded = False
+    for h, s in possible_secrets:
+        mac_forged = gmac(h, s, aad_forged, c_forged)
+        try:
+            assert gcm_decrypt(k, nonce, aad_forged, c_forged, mac_forged) == m_forged
+            succeeded = True
+            print(c_forged.hex(), mac_forged.hex())
+        except AssertionError:
+            pass
+    assert succeeded
diff --git a/index.html b/index.html
index 7b4572b..1fd52b1 100644
--- a/index.html
+++ b/index.html
@@ -35,10 +35,10 @@
 			Choose one of the following missions.
 		</p>
         <p>
-            <strong><a href="#">Nonce reuse</a>.</strong> Due to rising entropy
-            prices, Roseacrucis has started to reuse nonces. You must perform the
-            Forbidden Attack in order to recover the authentication key and
-            forge arbitrary ciphertext.
+            <strong><a href="/forbidden-salamanders/nonce-reuse">Nonce
+            reuse</a>.</strong> Due to rising entropy prices, Roseacrucis has
+            started to reuse nonces. You must perform the Forbidden Attack in order to
+            recover the authentication key and forge arbitrary ciphertext.
         </p>
         <p>
             <strong><a href="#">Nonce truncation</a>.</strong> The sorcerer
@@ -82,7 +82,7 @@
 				interpreted as the set of polynomials with coefficients in \(\mathbb{F}_2\)
 				of degree less than \(128\). Multiplication
 				is performed modulo \(m\). This field is of characteristic 2;
-				e.g., \((\alpha^5 + \alpha+1)+(\alpha^5 + \alpha+1) = 0\).
+				e.g., \((\alpha^5 + 1)+(\alpha^5 + 1) = 0\).
 			</p>
 			<p>
 				We interpret 16-byte blocks as elements in \(\mathbb{K}\)
@@ -95,40 +95,40 @@
 				the block.
 			</p>
 			<p>
-				12-byte nonces are interpreted as 96-bit integers in big-endian byte order.
-			</p>
+                12-byte nonces are interpreted as 96-bit integers in big-endian
+                byte order. Let \(\operatorname{Byte} = [0, 2^8-1]\) and
+                \(x_i\) refer to the \(i\)th 16-byte chunk of the bytestring
+                \(x\).
+            </p>
             <p>
-                Let \(\operatorname{Byte} = [0, 2^8-1]\).
+                \(\operatorname{encode_{big}}(x, n)\) encodes an integer \(x\) into \(n\) bytes in big-endian
+                byte order. \(\operatorname{pad_n}(x, p)\) pads the length of
+                the bytestring \(x\) to the nearest multiple of \(n\) with the
+                byte \(p\). \(\operatorname{AES}(k, x)\) refers to
+                the <a href="https://en.wikipedia.org/wiki/Advanced_Encryption_Standard">128-bit AES block cipher</a>.
             </p>
 			<br>
 			<div class="algorithm">
                 <p>\(\operatorname{GMAC}(h\in \mathbb{K}, s\in \mathbb{K}, aad\in \operatorname{Byte}^{y}, c\in \operatorname{Byte}^{z})\)</p>
 				<ol class="algorithm-code">
-                    <li>\( aad' = \operatorname{chunk}_{16}(aad, \operatorname{pad}=\mathtt{0x00}) \)</li>
-                    <li>\( c' = \operatorname{chunk}_{16}(c, \operatorname{pad}=\mathtt{0x00}) \)</li>
-					<li>\( len = \operatorname{encode_{big}}(128\vert aad' \vert, 8) \mathbin\Vert \operatorname{encode_{big}}(128\vert c'\vert, 8) \)</li>
-					<li>\( blocks = aad' \mathbin\Vert c' \mathbin\Vert (len) \mathbin\Vert (s) \)</li>
-					<li>\( \operatorname{return} \sum\limits_{i=1}^{\vert blocks\vert} blocks_{\vert blocks \vert-i} h^{i-1}\)</li>
+					<li>\( len = \operatorname{encode_{big}}(8y, 8) \mathbin\Vert \operatorname{encode_{big}}(8z, 8) \)</li>
+					<li>\( blocks = \operatorname{pad}_{16}(aad, 0) \mathbin\Vert \operatorname{pad}_{16}(c, 0) \mathbin\Vert len \mathbin\Vert s \)</li>
+                    <li>\( N = \frac{\vert blocks \vert}{16} \)</li>
+                    <li>\( \operatorname{return} \sum\limits_{i=1}^{N} blocks_{N-i} h^{i-1}\)</li>
 				</ol>
 			</div>
 			<br>
 			<br>
 			<div class="algorithm">
-                <p>\(\operatorname{GCM}(k\in \operatorname{Byte}^{16}, n\in \operatorname{Byte}^{12}, aad\in \operatorname{Byte}^{y}, m\in \operatorname{Byte}^{z})\)</p>
+                <p>\(\operatorname{AES-GCM}(k\in \operatorname{Byte}^{16}, n\in \operatorname{Byte}^{12}, aad\in \operatorname{Byte}^{y}, m\in \operatorname{Byte}^{z})\)</p>
 				<ol class="algorithm-code">
-					<li> \( r = \mathop{\Vert}\limits_{n'=2^{32}n+2}^{2^{32}n+2^{32}-1} \operatorname{AES-ECB}(k, n') \)</li>
+					<li> \( r = \mathop{\Vert}\limits_{n'=2^{32}n+2}^{2^{32}n+2^{32}-1} \operatorname{AES}(k, n') \)</li>
 					<li> \( c = r \oplus m \) </li>
-					<li> \( h = \operatorname{AES-ECB}(k, 0) \) </li>
-					<li> \( s = \operatorname{AES-ECB}(k, 2^{32}n + 1) \) </li>
+					<li> \( h = \operatorname{AES}(k, 0) \) </li>
+					<li> \( s = \operatorname{AES}(k, 2^{32}n + 1) \) </li>
 					<li> \( \operatorname{return} c, \operatorname{GMAC}(h, s, aad, c) \)</li>
 				</ol>
 			</div>
-            <p>
-                The authentication key \( h \) is independent of the
-                nonce \( n \). The constant term \( s \) acts as a blind to
-                hide the confidential block data in the MAC. Finally, note
-                that the polynomial computation reverses the order of the blocks.
-            </p>
 		</details>
     <!-- <script id="MathJax-script" async src="/forbidden-salamanders/static/mathjax.js"></script> -->
 	<!-- <script type="text/x-mathjax-config"> -->
diff --git a/nonce-reuse.html b/nonce-reuse.html
index 0795386..21c794d 100644
--- a/nonce-reuse.html
+++ b/nonce-reuse.html
@@ -48,7 +48,10 @@
             nonces are reused, the keystream will be identical, allowing us to
             recover plaintext via
             <a href="https://samwho.dev/blog/toying-with-cryptography-crib-dragging/">
-            crib dragging</a>, which makes this attack particularly effective.
+            crib dragging</a>, which makes this attack particularly effective:
+			\[
+				c' = c \oplus m \oplus m'.
+			\]
         </p>
         <p>
             However, we still need to compute a new MAC over the forged ciphertext.
@@ -65,35 +68,32 @@
             \[
                 mac + mac' = (s + s') + (len + len')h + (c_1 + c'_1)h^2 + (c_0+c'_0)h^3,
             \]
-            Since \(s = s'\) and \(x+x=0 \in \mathbb{F}_{2^{128}}\), we are
+            Since \(s = s'\) and \(x+x=0\) in \(\mathbb{F}_{2^{128}}\), we are
             left with the polynomial equation
             \[
                 0 = (mac + mac') + (len + len')h + (c_1 + c'_1)h^2 + (c_0+c'_0)h^3
             \]
             where all variables are known other than \(h\). Thus, recovering \(h\)
-            is a simple matter of finding the roots by <a href="https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields">factoring the
-            polynomial</a> using a computer algebra system such as SageMath.
+            is a matter of finding the roots by <a href="https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields">factoring the
+            polynomial</a>.
         </p>
         <p>
             We plug \(h\) back into the first equation to recover \(s\),
             and finally, we can forge the MAC for arbitary ciphertext under the
-            same nonce.
+            same nonce. Note that there may be multiple possible monomial roots;
+            in this case, one can check each possibility online.
         </p>
         </details>
         <br>
         <form>
             <div>
-            <input type="button" value="Autofill example">
-            </div>
-
-            <div>
             <label for="key">Key (16 bytes in hex)</label>
-            <input id="key" type="text">
+            <input id="key" type="text" value="59454c4c4f575f5355424d4152494e45">
             </div>
 
             <div>
             <label for="nonce">Nonce (12 bytes in hex)</label>
-            <input id="nonce" type="text">
+            <input id="nonce" type="text" value="4a4f5247454c424f52474553">
             </div>
 
             <div>
@@ -102,27 +102,17 @@
             </div>
 
             <div>
-            <label for="aad1">First additional authenticated data (in ASCII)</label>
-            <input id="aad1" type="text">
-            </div>
-
-            <div>
             <label for="m2">Second message (in ASCII)</label>
             <input id="m2" type="text">
             </div>
 
             <div>
-            <label for="aad2">Second additional authenticated data (in ASCII)</label>
-            <input id="aad2" type="text">
+            <label for="mf">Forged message; shorter than the first message (in ASCII)</label>
+            <input id="mf" type="text">
             </div>
 
             <div>
-            <label for="c">Forged ciphertext (in hex)</label>
-            <input id="c" type="text">
-            </div>
-
-            <div>
-            <input type="button" value="Compute authentication key and forge MAC">
+            <input type="button" value="Recover authentication key and forge MAC">
             </div>
         </form>
         <div>
@@ -130,46 +120,100 @@
         <input id="h" type="text" disabled value="deadbeef">
         </div>
         <div>
+        <label for="mf">Forged ciphertext, assuming first message is known</label>
+        <input id="mf" type="text">
+        </div>
+        <div>
         <label for="mac">Forged MAC</label>
         <input id="mac" type="text" disabled value="deadbeef">
         </div>
         <br>
+		<details>
+			<summary>
+                Show me the code.
+			</summary>
         <pre>
-from <a href="#">aesgcmanalysis</a> import xor, gcm_encrypt, nonce_reuse_recover_secrets, forge_gmac
-k = b"YELLOW_SUBMARINE"
-nonce = b"JORGELBORGES"
-m1 = b"""The universe (which others call the Library) is composed of an indefinite and
-perhaps infinite number of hexagonal galleries, with vast air shafts between,
-surrounded by very low railings."""
+from <a href="/git/forbidden-salamanders">aesgcmanalysis</a> import xor, gcm_encrypt, gcm_decrypt, nonce_reuse_recover_secrets
+
+k = b"tlonorbistertius"
+nonce = b"jorgelborges"
+m1 = b"The universe (which others call the Library)"
 aad1 = b"The Anatomy of Melancholy"
-m2 = b"From any of the hexagons one can see, interminably, the upper and lower floors."
+m2 = b"From any of the hexagons one can see, interminably"
 aad2 = b"Letizia Alvarez de Toledo"
-c1, mac1 = gcm_encrypt(k, nonce, m1, aad1)
-c2, mac2 = gcm_encrypt(k, nonce, m2, aad2)
+
+c1, mac1 = gcm_encrypt(k, nonce, aad1, m1)
+c2, mac2 = gcm_encrypt(k, nonce, aad2, m2)
 
 # Recover the authentication key and blind from public information
-h, s = nonce_reuse_recover_secrets(nonce, aad1, aad2, c1, c2, mac1, mac2)
+possible_secrets = nonce_reuse_recover_secrets(nonce, aad1, aad2, c1, c2, mac1, mac2)
 
 # Forge the ciphertext
-m_forged = b"As was natural, this inordinate hope was followed by an excessive depression."
-c_forged = xor(c2, xor(m2, m_forged))
+m_forged = b"As was natural, this inordinate hope"
+assert len(m_forged) <= len(m1)
+c_forged = xor(c1, xor(m1, m_forged))
 aad_forged = b"You who read me, are You sure of understanding my language?"
-mac_forged = forge_gmac(nonce, h, s, c_forged, aad_forged)
-assert gcm_decrypt(k, nonce, c_forged, aad_forged, mac_forged) == m_forged
-        </pre>
-		<details>
-			<summary>
-                Show me the code.
-			</summary>
-        </details>
+
+# Check possible candidates for authentication key
+succeeded = False
+for h, s in possible_secrets:
+    mac_forged = gmac(h, s, aad_forged, c_forged)
+    try:
+        assert gcm_decrypt(k, nonce, aad_forged, c_forged, mac_forged) == m_forged
+        succeeded = True
+        print(c_forged.hex(), mac_forged.hex())
+    except AssertionError:
+        pass
+assert succeeded</pre></details>
 		<details>
 			<summary>
                 Show me the math.
 			</summary>
-            see homepage
-            The forged ciphertext can be computed as a xor from original ciphertext,
-            which should be particularly easy as nonce reuse reveals the XOR and then
-            crib dragging.
+            <p>
+                A description of the construction of GMAC can be found at the <a
+                href="/forbidden-salamanders">mission homepage</a>.
+            </p>
+            <p>
+                Once the polynomial difference is computed, one can use SageMath
+                to compute the factors:
+            </p>
+            <pre>
+K = GF(2**128, name='x', modulus=x^128+x^7+x^2+x+1)
+x = K.gen()
+S = PolynomialRing(K, 'y')
+y = S.gen()
+p = (x^127 + x^125 + x^121 + x^118 + x^117 + x^116 + x^113 + x^111 + x^108 + x^105 + x^102 + x^99 + x^97 + x^95 + x^89 + x^87 +
+ x^85 + x^84 + x^81 + x^78 + x^73 + x^71 + x^69 + x^67 + x^65 + x^63 + x^62 + x^59 + x^57 + x^55 + x^54 + x^53 + x^52 + x^51 +
+x^49 + x^47 + x^46 + x^45 + x^43 + x^41 + x^39 + x^38 + x^37 + x^36 + x^33 + x^29 + x^28 + x^25 + x^21 + x^20 + x^17 + x^15 + x
+^13 + x^9 + x^7 + x^4 + x^3 + x)*y^5 + (x^127 + x^125 + x^121 + x^118 + x^117 + x^116 + x^113 + x^111 + x^108 + x^105 + x^102 +
+ x^99 + x^97 + x^95 + x^89 + x^87 + x^85 + x^84 + x^81 + x^78 + x^73 + x^71 + x^69 + x^67 + x^65 + x^63 + x^62 + x^59 + x^57 +
+x^55 + x^54 + x^53 + x^52 + x^51 + x^49 + x^47 + x^46 + x^45 + x^43 + x^41 + x^39 + x^38 + x^37 + x^36 + x^33 + x^29 + x^28 + x
+^25 + x^21 + x^20 + x^17 + x^15 + x^13 + x^9 + x^7 + x^4 + x^3 + x)*y^4 + (x^127 + x^125 + x^124 + x^123 + x^122 + x^120 + x^11
+9 + x^118 + x^115 + x^112 + x^111 + x^110 + x^109 + x^108 + x^106 + x^101 + x^100 + x^99 + x^98 + x^96 + x^93 + x^88 + x^87 + x
+^84 + x^83 + x^82 + x^78 + x^77 + x^74 + x^71 + x^70 + x^69 + x^68 + x^65 + x^62 + x^61 + x^60 + x^57 + x^55 + x^53 + x^50 + x^
+49 + x^47 + x^46 + x^44 + x^43 + x^41 + x^40 + x^39 + x^37 + x^36 + x^35 + x^34 + x^33 + x^31 + x^30 + x^25 + x^19 + x^16 + x^1
+4 + x^13 + x^12 + x^11 + x^10 + x^5 + x^2 + x + 1)*y^3 + (x^124 + x^123 + x^121 + x^120 + x^119 + x^118 + x^117 + x^116 + x^115
+ + x^113 + x^112 + x^110 + x^107 + x^106 + x^104 + x^103 + x^101 + x^99 + x^98 + x^95 + x^94 + x^83 + x^82 + x^81 + x^80 + x^79
+ + x^77 + x^76 + x^75 + x^73 + x^72 + x^67 + x^64 + x^63 + x^62 + x^61 + x^59 + x^57 + x^56 + x^54 + x^53 + x^51 + x^49 + x^48
++ x^46 + x^45 + x^44 + x^42 + x^36 + x^35 + x^34 + x^33 + x^31 + x^28 + x^22 + x^21 + x^17 + x^12 + x^11 + x^10 + x^8 + x^5 + x
+^4 + 1)*y^2 + (x^120 + x^119 + x^56 + x^55)*y + (x^127 + x^126 + x^125 + x^124 + x^123 + x^118 + x^117 + x^116 + x^115 + x^114
++ x^113 + x^105 + x^103 + x^98 + x^96 + x^94 + x^91 + x^90 + x^89 + x^87 + x^86 + x^85 + x^83 + x^81 + x^80 + x^78 + x^77 + x^7
+0 + x^66 + x^63 + x^62 + x^59 + x^58 + x^54 + x^53 + x^52 + x^50 + x^47 + x^45 + x^44 + x^42 + x^40 + x^39 + x^37 + x^35 + x^34
+ + x^30 + x^29 + x^27 + x^26 + x^25 + x^23 + x^18 + x^16 + x^14 + x^13 + x^12 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3
++ 1)
+for factor, _ in p.factor():
+    if factor.degree() == 1:
+        print('Authentication key:', factor - y)</pre>
+		<p>
+			However, the library powering this demonstration implements <a href="https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields">polynomial factoring over finite fields</a>, which is an edifying exercise.
+		</p>
+		<p>
+			We present advice for those who wish to implement polynomial factorization as well:
+		</p>
+		<ul>
+			<li>The gcd of two polynomials is unique only up to multiplication by a non-zero constant because &ldquo;greater&rdquo; is defined for polynomials in terms of degree. When used in algorithms, gcd refers to the <em>monic</em> gcd, which is unique.</li>
+			<li>The <a href="https://math.stackexchange.com/a/943626/1084004">inverse Frobenius automorphism</a> (i.e., square root) in \(\mathbb{F}_{2^{128}}\) is given by \(\sqrt{x} = x^{2^{127}})\).</li>
+		</ul>
         </details>
 <script>
 MathJax = {