find-stretch.sage: Calculate stretch shifts for various block sizes.

[ocb-tv] / find-stretch.sage

diff --git a/find-stretch.sage b/find-stretch.sage
new file mode 100644 (file)
index 0000000..cdd408b
--- /dev/null
+++ b/find-stretch.sage
@@ -0,0 +1,55 @@
+#! /usr/local/bin/sage
+### -*-python-*-
+
+k = GF(2)
+
+def memoize(f):
+  memo = dict()
+  def _(*args):
+    try:
+      return memo[args]
+    except KeyError:
+      r = f(*args)
+      memo[args] = r
+      return r
+  return _
+
+@memoize
+def Abase(w, c):
+  return matrix(k, [[j == i for j in xrange(w)]
+                    for i in xrange(w)] +
+                   [[j == i or j == i + c for j in xrange(w)]
+                    for i in xrange(w)])
+
+@memoize
+def A(w, c, i): return Abase(w, c)[i:i + w, 0:w]
+
+def B(w, c, i, j): return A(w, c, i) + A(w, c, j)
+
+def urank(w, c, i): return A(w, c, i).rank()
+
+def domlen(w, c):
+  for i in xrange(w):
+    if urank(w, c, i) < w: return ZZ(i)
+    for j in xrange(i):
+      if B(w, c, i, j).rank() < w: return ZZ(i)
+  return ZZ(w) #unlikely
+
+def shifts(w):
+  for c1 in xrange(0, 8):
+    for c8 in xrange(w - 8, -8, -8):
+      yield c8 + c1
+
+def stretch_shift(w):
+  want_bits = (w - 1).nbits()
+  best_bits, best_dom, best_c = 0, 0, -1
+  for c in shifts(w):
+    d = domlen(w, c)
+    bits = d.nbits()
+    if bits == want_bits: return c, d
+    elif bits > best_bits: best_bits, best_dom, best_c = bits, d, c
+  return best_c, best_dom
+
+for w in [64, 96, 128, 192, 256]:
+  c, dom = stretch_shift(w)
+  print '%3d: %3d [%d]' % (w, c, dom)