new files

final_expo_k57.py

+from collections import defaultdict
+
+# This companion file can be used to derive and check generic formulas
+# for the final exponentiation for prime k.
+
+# There is no algorithm here. Just inspection of formulas, and
+# optimization work. It does seem that number of necessary inversions
+# varies, but the rest seems to be pretty constant.
+#
+# I have no reason to imagine that I haven't missed some optimization.
+#
+# Results:
+# 
+# sage: load("final-expo-k57.sage")
+# cost for k=5 i=1: 1c + 3T + 7M + 3p
+# cost for k=5 i=2: 1I + 1c + 3T + 7M + 3p
+# cost for k=5 i=3: 2I + 1c + 3T + 7M + 3p
+# cost for k=5 i=4: 1I + 1c + 3T + 7M + 3p
+# cost for k=7 i=1: 1c + 5T + 11M + 5p
+# cost for k=7 i=2: 1I + 1c + 5T + 11M + 5p
+# cost for k=7 i=3: 1I + 1c + 5T + 11M + 5p
+# cost for k=7 i=4: 2I + 1c + 5T + 11M + 5p
+# cost for k=7 i=5: 2I + 1c + 5T + 11M + 5p
+# cost for k=7 i=6: 1I + 1c + 5T + 11M + 5p
+
+# This structure keeps track of all the applied operations.
+class counter():
+    def __init__(self):
+        self._counts=defaultdict(int)
+    class element:
+        def __init__(self, count, val=1):
+            self.val=val
+            self.count_dict = count
+        def __mul__(self, other):
+            w = self.count_dict._MUL(self.val, other.val)
+            return counter.element(self.count_dict, w)
+        def __pow__(self, e):
+            if e == -1:
+                w = self.count_dict._INV(self.val)
+            else:
+                w = self.count_dict._POW(self.val, e)
+            return counter.element(self.count_dict, w)
+        def frob(self, e, j):
+            w = self.count_dict._FROB(self.val, e, j)
+            return counter.element(self.count_dict, w)
+        def __str__(self):
+            return str(self.val)
+        def __repr__(self):
+            return repr(self.val)
+    def get_element(self):
+        return counter.element(count=self)
+    def counts(self):
+        return self._counts
+    def _POW(self, a, e):
+        assert a != 0
+        self._counts[str(e)] += 1
+        return a * e
+    def _FROB(self, a, e, j):
+        assert a != 0
+        self._counts[str(e)] += 1
+        return a * e**j
+    def _MUL(self, a, b):
+        assert a != 0 and b != 0
+        self._counts['M'] += 1
+        return a + b
+    def _INV(self, a):
+        assert a != 0
+        self._counts['I'] += 1
+        return -a
+    def mark_free(self, v):
+        self._counts.pop(str(v), None)
+    def rename(self, v, w):
+        n=self._counts.pop(str(v), None)
+        if n is not None:
+            self._counts[str(w)]=n
+
+    def text_counts(self):
+        return " + ".join(["%s%s" % (v,k) for k,v in self._counts.items() if v])
+
+def print_final_expo_k57():
+    ZZc = PolynomialRing(ZZ, names=('c',)); (c,) = ZZc._first_ngens(1)
+    RTc = PolynomialRing(ZZc, names=('T',)); (T,) = RTc._first_ngens(1)
+    RTcp = PolynomialRing(RTc, names=('p',)); (p,) = RTcp._first_ngens(1)
+
+    for k in [5,7]:
+        phi = cyclotomic_polynomial(k)
+
+        # We're relying on k be prime several times. For example, Frobenius-es
+        # don't simplify in subfields.
+        assert Integer(k).is_prime()
+
+        for i in range(1,k):
+            P_p = T**i + c*phi(T)
+            assert phi(P_p) % phi(T) == 0
+            P_expo = phi(P_p) // phi(T)
+            C=[]
+            E = P_expo
+            while E != 0:
+                E, C0 = E.quo_rem(P_p)
+                C.append(C0)
+            assert len(C) == k-1
+            Cp = RTcp(C)
+            assert P_expo == Cp(P_p)
+
+            counts = counter()
+
+            a=counts.get_element()
+            cc=[0 for x in C]
+            res=a
+
+            if i == 1:
+                cc[-1] = a**c
+                for j in range(1, k-1):
+                    cc[-1-j] = cc[-j]**T * cc[-1]
+                cc[0] = cc[0] * a
+            elif i==k-1:
+                cc[-1] = a*a**c
+                cc[0] = a * (cc[-1]**T)**-1
+                for j in range(1, k-2):
+                    cc[j] = cc[0] * cc[j-1]**T
+            if k == 5:
+                if i == 2:
+                    # [-c*T^3 - T + 1,
+                    #  -c*T^3 + (-c - 1)*T + 1,
+                    #   c*T^2 + c + 1,
+                    #   c]
+                    cT0 = cc[3] = a**c
+                    cT = cT0**T
+                    cT2 = cT**T
+                    ca = cT0 * a
+                    cc[2] = cT2 * ca
+                    b = (cc[2]**T)**-1
+                    cc[1] = b * a
+                    cc[0] = cc[1] * cT
+                elif i == 3:
+                    # c*p^3
+                    # + (c*T^3 + T^2 - T + c)*p^2
+                    # + (c*T^3 + T^2 - T + c + c*T + 1)*p
+                    #  - c*T^2 - T + 1
+                    cT0 = cc[3] = a**c
+                    cT = cT0**T
+                    cTa = cT * a
+                    ai = a**-1
+                    b = cTa**T * ai
+                    cc[0] = b**-1
+                    cc[2] = b**T * cT0
+                    cc[1] = cc[2] * cTa
+            if k==7:
+                if i==2:
+                    # [-c*T^5 - T + 1,
+                    #  -c*T^5 - c*T^3 - T + 1,
+                    #  -c*T^5 - c*T^3 - T + 1 -c*T,
+                    #   c*T^4 + c*T^2 + c + 1,
+                    #   c*T^2 + c,
+                    #   c]
+                    cT0 = cc[5] = a**c
+                    cT = cT0**T
+                    c1 = cT0 * a
+                    cT2 = cT**T
+                    cT2c = cc[4] = cT2 * cT0
+                    cT3T = cT2c ** T
+                    cT4T2 = cT3T ** T
+                    cc[3] = cT4T2 * c1
+                    cc[2] = (cc[3]**T)**-1 * a
+                    cc[1] = cc[2] * cT
+                    cc[0] = cc[2] * cT3T
+                elif i==3:
+                    #         - c*T^4                             - T + 1,
+                    #         - c*T^4                       - c*T - T + 1,
+                    #  -c*T^5 - c*T^4 - T^2         - c*T + 1,
+                    #  -c*T^5 - c*T^4 - T^2 - c*T^2 - c*T + 1,
+                    #   c*T^3 + c + 1,
+                    #   c]
+                    cT0 = cc[5] = a**c
+                    cT = cT0**T
+                    cT2 = cT**T
+                    cT3 = cT2**T
+                    cT3c = cT3 * cT0
+                    cc[4] = cT3c * a
+                    b = cc[4]**-1
+                    bT = b**T
+                    cc[1] = bT*a
+                    cc[0] = cc[1]*cT
+                    cc[3] = cc[1]**T*cc[1]
+                    cc[2] = cc[3] * cT2
+                elif i==4:
+                    # target=[-c*T^3 - T + 1,
+                    #          c*T^5 + c*T^4 + T^3 + c*T^2 + (c - 1)*T + c + 1,
+                    #          c*T^5 + c*T^4 + T^3 +         (c - 1)*T + c,
+                    #          c*T^4 + T^2 + (c - 1)*T + c,
+                    #          c*T^4 + T^2 - T + c,
+                    #          c]
+                    cT0 = cc[5] = a**c
+                    cT = cT0**T
+                    cT2 = cT**T
+                    cT2a = cT2*a
+                    cT3T = cT2a**T
+                    cc[0] = cT3T**-1 * a
+                    cc[4] = (cc[0]**T)**-1 * cT0
+                    cc[3] = cc[4] * cT
+                    cc[2] = cc[4]**T*cc[4]
+                    cc[1] = cc[2]*cT2a
+                elif i==5:
+                    # target=[-c*T^2 - T + 1,
+                    #         -c*T^4 - T^3 + (-c + 1)*T^2 - T + 1,
+                    #          c*T^5 + T^4 + (c - 1)*T^3 + T^2 + (c - 1)*T + c + 1,
+                    #          c*T^5 + T^4 + (c - 1)*T^3 + T^2 - T + c,
+                    #          c*T^5 + T^4 - T^3 + c,
+                    #          c]
+                    cT0 = cc[5] = a**c
+                    cT = cT0**T
+                    cTa = cT*a
+                    cT2T = cTa**T
+                    cc[0] = cT2T**-1 * a
+                    b = cc[0]**T
+                    bT = b**T
+                    cc[1] = bT * cc[0]
+                    cc[3] = (cc[1]**T)**-1*cT0
+                    cc[2] = cc[3] * cTa
+                    cc[4] = cc[3] * b
+
+            if isinstance(cc[0], counter.element):
+                res = cc[0]
+                for j in range(1, k-1):
+                    res = res * cc[j].frob(p, j)
+
+            if res.val == Cp:
+                print "cost for k=%d i=%d: %s" % (k, i, counts.text_counts())
+            else:
+                print "k=%d i=%d: wrong formula" % (k,i)
+                print res.val
+                print Cp
+                print C
+
+# Other nice expressions, but I'm not sure they're any useful.
+#
+# This uses phi(p)/r = phi(T^i)/phi(T) + (phi(p)-phi(T^i))/(p-T^i)
+#
+# While the expression above writes down nicely for i=1, it is much less
+# so for i>4.
+# R.<p,T,c>=QQ[]
+#
+# k=5
+# phi_k=cyclotomic_polynomial(k)
+# r=phi_k(T)
+# for i in range(1,k):
+#     cx=((phi_k(p)-phi_k(T^i))/(p-T^i))
+#     assert cx.denominator()==1
+#     cx=cx.numerator()
+#     [cx.coefficient({p:j})%r for j in range(k-1)]
+# 
+# CC.<T>=QQ[]
+# 
+# [T^3 + T^2 + T + 1, T^2 + T + 1, T + 1, 1]
+# c3=x
+# c2=c3*T+x
+# c1=c2*T+x
+# c0=c1*T+x
+# [c0,c1,c2,c3]
+# TMTMTMFFF
+# 
+# [-T^3, -T^3 - T, T^2 + 1, 1]
+# a1=x*T
+# a2=a1*T
+# c2=a2+x
+# a3=-c2
+# c1=a3*T
+# c0=c1+a1
+# c3=x
+# [c0,c1,c2,c3]
+# TTMITMFFF
+# 
+# [-T^2, T^3 + T + 1, T^3 + 1, 1]
+# a1=x*T 
+# a2=a1*T
+# a3=a2*T
+# c0=-a2
+# c2=x+a3
+# c1=c2+a1
+# c3=x
+# [c0,c1,c2,c3]
+# TTTMMIFFF
+# 
+# [-T, -T^2 - T, -T^3 - T^2 - T, 1]
+# a1=x*T
+# c0=-a1
+# c1=c0*T+c0
+# c2=c1*T+c0
+# [c0,c1,c2,c3]
+# TITMTMFFF
+
+# attach("code/CocksPinchVariant.py")
+# load("redo-bar-duq.sage")