from final_expo_k57 import * def horner_list(f,T): if f == 0: return [] L=horner_list(f//T,T) L.append(f%T) return L def horner(f,T): if f//T == 0: return f return "(%s)*%s+%s" % (horner(f//T,T),T,f%T) def count_formula_k8(i, c): (T, h_t, h_y,) = c.parent()._first_ngens(3) c_list = horner_list(c, T) # formulas for c in case k=8, using the simplification h_u=(h_t+1)/2 # (i=1,3,5,7 -- note that from i to i+4, it's only a matter of changing # signs). # let r = u^2+y^2 # ((((r-u+1/4)*T+y )*T -y+1/4)*T +u-1)*T +r # ((((r-u+1/4)*T+u-1/2)*T -y+1/4)*T -y )*T +r # ((((r-u+1/4)*T-y )*T -y+1/4)*T -u+1)*T +r # ((((r-u+1/4)*T-u+1/2)*T -y+1/4)*T +y )*T +r # # notice the following cyclic transformations on the coefficients # that do change between the cases i=1 and i=3. Alas, this doesn't # extend to all coefficients # # (u,y) -> (1-y,u-1/2) -> (3/2-u,1/2-y) -> (1/2+y,1-u) -> (u, y) # R2 = QQ['U, V, u, y']; (U, V, u, y,) = R2._first_ngens(4) cl2 = [x(h_t=2*u-1,h_y=y) for x in c_list] e0, e1, e2, e3, e4 = cl2 print "reduced horner: ", cl2 T=(2*U+V)*2 # our target is the following exponent (here for i==7) # target = ((((u^2 + y^2 - u + 1/4)*T+-u + 1/2)*T+-y + 1/4)*T+y)*T+u^2 + y^2 # we will compute the input raised to the following exponents: # c0 d0 c1 d1 c2 d2 c3 d3 c4 d4 # # with: c0 = u**2 + y**2 - u d0 = 1 assert c0 + d0/4 == e0 assert c0 + d0/4 == u**2 + y**2 - u + 1/4 c02 = c0*2 c1 = (c02*2 + d0) * U + c02*V if i==7: c1 = c1 - u d1 = 1 elif i==3: c1 = c1 + u d1 = -1 elif i==1: c1 = c1 + y d1 = 0 elif i==5: c1 = c1 - y d1 = 0 d1 = d1 + d0*V assert c1 + d1/2 == (c0 + d0/4) * T + e1 c12 = c1*2 c2 = (c12 + d1) * 2*U + (c12 + d1) * V - y d2 = 1 assert c2 + d2/4 == (c1 + d1/2) * T + e2 c22 = c2*2 c3 = (c22*2 + d2) * U + c22*V if i == 7: c3 = c3 + y elif i == 3: c3 = c3 - y elif i == 1: c3 = c3 + u - 1 elif i == 5: c3 = c3 - u + 1 d3 = d2*V assert c3 + d3/2 == (c2 + d2/4) * T + e3 c4 = (c3*2 + d3) * (2*U+V) + u**2 + y**2 assert c4 == (c3 + d3/2) * T + e4 counts = counter() # 2y 2u 2m a = counts.get_element() ai = a^-1 ay = a ** y ayi = ay^-1 ay2 = ay ** y au = a ** u aui = au^-1 au2 = au ** u ar = au2 * ay2 r = ar * aui s = a assert r.val == c0 assert s.val == d0 # 1T 2m 1m# r = r ** 2 r = (r ** 2 * s) ** U * r ** V if i == 7: r = r * aui s = a * s ** V ## this multiplication happens *only* if V==1 elif i == 3: r = r * au s = ai * s ** V ## this multiplication happens *only* if V==1 elif i == 1: r = r * ay s = s ** V elif i == 5: r = r * ayi s = s ** V assert r.val == c1 assert s.val == d1 # 1T 2m r = r ** 2 r = r * s r = (r ** 2) ** U * r ** V r = r * ayi s = a assert r.val == c2 assert s.val == d2 # 1T 2m r = r ** 2 r = (r ** 2 * s) ** U * r ** V if i == 7: r = r * ay elif i == 3: r = r * ayi elif i == 1: r = r * au r = r * ai elif i == 5: r = r * aui r = r * a s = s ** V assert r.val == c3 assert s.val == d3 # 1T 3m r = r ** 2 r = r * s ## r = (r ** 2) ** U * r ** V r = r * ar assert r.val == c4 counts.mark_free('I') counts.rename('2', 's') # Use the fact that 1T = 1U + 1V + 1M + 2s nT = counts._counts['U'] counts._counts['U'] -= nT counts._counts['V'] -= nT counts._counts['M'] -= nT counts._counts['s'] -= 2*nT counts._counts['T'] = nT print counts.text_counts() # This prints 11M + 2u + 4T + 2V + 2y ; V is actually just a variable # assignment, and one of the M can be elided if V=0, at least in the # cases i=3 and i=7, so we get: # 2y 2u 4T 10m 1m# # for i=1 and i=5, it seems that we don't have the m# def count_formula_k6(i, tr, c): (T, h_t, h_y,) = c.parent()._first_ngens(3) R2 = QQ['U, V, u, w']; (U, V, u, w,) = R2._first_ngens(4) # Be smart. # For k=6, the expressions of t and y are: # t = T + 1 +h_t*r # y = 1/3*T^2 - 2/3*T +h_y*r # t = T + 1 +h_t*r # y = -1/3*T^2 - 2/3 +h_y*r # t = -T + 2 +h_t*r # y = 1/3*T^2 - 2/3*T + 1 +h_y*r # t = -T + 2 +h_t*r # y = -1/3*T^2 + 1/3 +h_y*r # but if we reduce to the parity bit of t+y only: # parity = T + 1 +h_t + T^2 +h_y = 1 + h_t + h_y # parity = T + 1 +h_t - T^2 +h_y = 1 + h_t + h_y # parity = -T +h_t + T^2 + 1 +h_y = 1 + h_t + h_y # parity = -T +h_t - T^2 + 1 +h_y = 1 + h_t + h_y # so that in all cases, we have either h_t odd and h_y # even, or the converse. for parity_ht in range(2): if parity_ht == 0: # This one expresses the result as a function of: # u = h_t/2 # w = (h_y-z)/2 z = -1 if tr == 0 else 1 # so we're assuming that h_t is even and h_y is odd. new_c=horner_list(3*c(h_t=2*u,h_y=(2*w+z),T=tr+3*U),U) new_c[1] /= 3 new_c[0] /= 9 assert sum([x(u=h_t/2,w=(h_y-z)/2)*(T-tr)^j for j,x in enumerate(reversed(new_c))]) == 3*c elif parity_ht == 1: # This one expresses the result as a function of: # u = (h_t-z)/2 # w = h_y/2 # so we're assuming that h_t is odd and h_y is even. z = -1 new_c=horner_list(3*c(h_t=2*u+z,h_y=2*w,T=tr+3*U),U) new_c[1] /= 3 new_c[0] /= 9 assert sum([x(u=(h_t-z)/2,w=h_y/2)*(T-tr)^j for j,x in enumerate(reversed(new_c))]) == 3*c print "formula for 3*c in case T mod 3 == %d and h_t mod 2 == %d" % (tr, parity_ht) print new_c # u and w are always integers. # order is: # i == 1, T mod 3 == 0, h_t even u = h_t/2 w = (h_y+1)/2 # i == 1, T mod 3 == 0, h_t odd u = (h_t+1)/2 w = h_y/2 # i == 1, T mod 3 == 1, h_t even u = h_t/2 w = (h_y-1)/2 # i == 1, T mod 3 == 1, h_t odd u = (h_t+1)/2 w = h_y/2 # i == 5, T mod 3 == 0, h_t even u = h_t/2 w = (h_y+1)/2 # i == 5, T mod 3 == 0, h_t odd u = (h_t+1)/2 w = h_y/2 # i == 5, T mod 3 == 1, h_t even u = h_t/2 w = (h_y-1)/2 # i == 5, T mod 3 == 1, h_t odd u = (h_t+1)/2 w = h_y/2 # we're going to horner-exponentiate this with exponent T-(T mod 3), # which we know is a multiple of 3. Only thing is that we do not want # to deal with (T-(T mod 3))/3, because its 2-naf weight can be far # apart the weight of T. # [1 + r - 6*w, -r + 3*u + 3*w, r + 3*u - 9*w + 3] # [1 + r - 3*u + 3*w, -r + 6*u - 6*w - 3, r] # [1 + r + 6*w, r + 3*u + 3*w, r + 6*u] # [1 + r - 3*u - 3*w, r - 6*w, r + 3*u - 9*w] # [1 + r - 6*w, -r - 3*u + 3*w, r + 6*u] # [1 + r - 3*u + 3*w, -r - 6*w, r + 3*u + 9*w] # [1 + r + 6*w, r - 3*u + 3*w, r + 3*u + 9*w + 3] # [1 + r - 3*u - 3*w, r - 6*u - 6*w + 3, r] # where r is 3*u^2 + 9*w^2 e0, e1, e2 = new_c a = 1 a3 = 3 a3u = a3 * u a6u = a3u * 2 a3w = a3 * w a6w = a3w * 2 a9w = a6w + a3w a9w2 = a9w * w a3u2 = a3u * u ar = a3u2 + a9w2 c0 = a + ar rplus3 = ar + a3 if parity_ht == 1: c0 = c0 - a3u if parity_ht == 0 and tr == 0: c0 = c0 - a6w elif parity_ht == 1 and tr == 0: c0 = c0 + a3w elif parity_ht == 0 and tr == 1: c0 = c0 + a6w elif parity_ht == 1 and tr == 1: c0 = c0 - a3w assert c0 == e0 if i == 1: if tr == 0: if parity_ht == 0: c1 = -ar + a3u + a3w c2 = rplus3 + a3u - a9w else: c1 = -rplus3 + a6u - a6w c2 = ar else: if parity_ht == 0: c1 = ar + a3u + a3w c2 = ar + a6u else: c1 = ar - a6w c2 = ar + a3u - a9w else: if tr == 0: if parity_ht == 0: c1 = -ar - a3u + a3w c2 = ar + a6u else: c1 = -ar - a6w c2 = ar + a3u + a9w else: if parity_ht == 0: c1 = ar - a3u + a3w c2 = rplus3 + a3u + a9w else: c1 = rplus3 - a6u - a6w c2 = ar assert c1 == e1 assert c2 == e2 counts = counter() a = counts.get_element() a3 = a^2 * a a3u = a3 ^ u a6u = None a3w = a3 ^ w a6w = a3w ^ 2 a9w = a6w * a3w a9w2 = a9w ^ w a3u2 = a3u ^ u ar = a3u2 * a9w2 acc = a * ar if i == 1: if tr == 0: if parity_ht == 0: ar3 = ar * a3 acc = acc * a6w^-1 acc = acc^U acc = acc * ar^-1 * a3u * a3w acc = acc^U acc = acc * ar3 * a3u * a9w^-1 else: ar3 = ar * a3 tmp = a3u^-1 * a3w acc = acc * tmp acc = acc^U acc = acc * (ar3 * tmp^2)^-1 acc = acc^U acc = acc * ar else: if parity_ht == 0: acc = acc * a6w acc = acc^U ar3u = ar * a3u acc = acc * ar3u * a3w acc = acc^U acc = acc * ar3u * a3u else: acc = acc * (a3u * a3w)^-1 acc = acc^U acc = acc * ar * a6w^-1 acc = acc^U acc = acc * ar * a3u * a9w^-1 else: if tr == 0: if parity_ht == 0: acc = acc * a6w^-1 acc = acc^U ar3u = ar * a3u acc = acc * (ar3u)^-1 * a3w acc = acc^U acc = acc * ar3u * a3u else: acc = acc * a3u^-1 acc = acc * a3w acc = acc^U acc = acc * (ar * a6w)^-1 acc = acc^U acc = acc * ar * a3u * a9w else: if parity_ht == 0: ar3 = ar * a3 acc = acc * a6w acc = acc^U acc = acc * ar * a3u^-1 * a3w acc = acc^U acc = acc * ar3 * a3u * a9w else: ar3 = ar * a3 a3u3w = a3u * a3w acc = acc * (a3u3w)^-1 acc = acc^U acc = acc * ar3 * (a3u3w^2)^-1 acc = acc^U acc = acc * ar assert acc.val == (e0 * U + e1) * U + e2 counts.mark_free('I') counts.rename('2', 's') print counts.text_counts() ######################################## # summary of counts. Haven't been able to find lower counts. Note # that it's not uniform, maybe we should take that as a hint that # there's some better way... # 12M + 2s + 2u + 2w + 2U # 10M + 3s + 2u + 2w + 2U # 10M + 2s + 2u + 2w + 2U # 11M + 2s + 2u + 2w + 2U # 10M + 2s + 2u + 2w + 2U # 11M + 2s + 2u + 2w + 2U # 12M + 2s + 2u + 2w + 2U # 10M + 3s + 2u + 2w + 2U def formulas(k): R = QQ['T, h_t, h_y']; (T, h_t, h_y,) = R._first_ngens(3) r = cyclotomic_polynomial(k)(T) K = CyclotomicField(k, names=('w',)); (w,) = K._first_ngens(1) if k==6: D=3 elif k==8: D=4 inv_sqrt_D = (1/sqrt(K(-D))).polynomial() ld = inv_sqrt_D.list()[-1] assert inv_sqrt_D.degree() < euler_phi(k) # choose positive leading coefficient if ld < 0: inv_sqrt_D = -inv_sqrt_D inv_sqrt_D = inv_sqrt_D(T) if k==6: subfamilies=[ # For i==1, t0 == T+1 # Recall that T=2 mod 3 is forbidden since r=Phi_6(T). # The representatives of (t0-2)*inv_sqrt_D in 1/3 * Z are: # (2T^2-3T+1)/3 : outside range # (T^2-2T)/3 = T*(T-2)/3 : good if T = 0,2 mod 3 (hence only 0) # (-T-1)/3 : good if T is 2 mod 3, so *never good* ! # (-T^2-2)/3 : good if T is 1, 2 mod 3 (hence only 1) (1, T+1, [(0,6,(T*(T-2))/3), (1,6,(1-T**2)/3-1), ]), # what might correspond to CocksPinchVariant is: #(1, T+1, [(0,6,(T*(T-2))/3), (1,6,(2+T**2)/3), ]), # For i==5, t0 == 2-T # The representatives of (t0-2)*inv_sqrt_D in 1/3 * Z are: # (T-2*T^2)/3 = T*(1-2*T)/3 : outside range # (1-T^2)/3 : good if T is 1 or 2 mod 3 (hence only 1) # (2-T)/3 : good if T is 2 mod 3, so *never good* ! # (T^2-2*T+3)/3 : good if T is 0 or 2 mod 3 (hence only 0) (5, 2-T, [ (0,3,1+T*(T-2)/3), (1,3,(1-T**2)/3), ]), # what might correspond to CocksPinchVariant is: #(5, 2-T, [ (0,3,1-T*(T-2)/3), (1,3,-(1-T**2)/3), ]), ] # Note that congruence classes on ht and hy will force p to be an # integer, even though it seems to have a 1/4 in the denominator. # E.g. in the case i=1, T=1 mod 3, t0=T+1, y0=-(T^2+2)/3, we have # sage: (4*p(T,h_t,h_y)).change_ring(GF(2)).factor() # (h_t + h_y + 1)^2 * (T^2 + T + 1)^2 # meaning that the denominator cancels only if h_t+h_y+1 is even (and # in fact it is also sufficient). (note that (T^2 + T + 1) is # necessarily even since r is prime). # # -> as a consequence, the formula should not be specific to one # congruence class of T mod 4 elif k==8: subfamilies=[ (1, T+1, [(0,2,(T-1)*T**2/2)]), (3, T**3+1, [(0,2,-(T+1)*T/2)]), (5, -T+1, [(0,2,(-T-1)*T**2/2)]), (7, -T**3+1, [(0,2,(-T+1)*T/2)]), ] else: # just for completeness. This ignores the fact that the # denominator of the expression may force some congruence classes subfamilies=[] for i in [i for i in range(k) if gcd(i,k) == 1]: t0 = (T**i+1) % r y0 = ((t0-2)*inv_sqrt_D) % r subfamilies.append((i, t0, [(0,1,y0)])) for i,t0,y0class in subfamilies: for tr,tq,y0 in y0class: assert (y0 - (t0-2)*inv_sqrt_D) % r == 0 assert (D * y0**2 + (t0-2)**2) % r == 0 t = t0 + h_t*r y = y0 + h_y*r p = (t**2 + D*y**2)/4 print("#" * 40) print("k={} D={} i={} T%{}={}".format(k,D,i,tq,tr)) print("t = {} +h_t*r".format(t0)) print("y = {} +h_y*r".format(y0)) print("coeffs 4*p, T^{} first, constant last:".format(2*euler_phi(k))) for pp,j in [(((4*p).coefficient({T:j})).factor(),j) for j in range(2*euler_phi(k), -1, -1)]: if j > 1: print("({})*T^{} + ".format(pp,j)) elif j == 1: print("({})*T + ".format(pp)) elif j == 0: print("({})".format(pp)) print("") order_red = (p+1-t0) print("p+1-t0 = {}".format(order_red.factor())) assert order_red % r == 0 c = (order_red // r) print "denominator cancellation condition for c: ", (4*c(T,h_t,h_y)).change_ring(GF(2)).factor() print("coeffs c, T^{} 1st, constant last:".format(euler_phi(k))) for cc,j in [((c.coefficient({T:j})).factor(),j) for j in range(euler_phi(k), -1, -1)]: if j > 1: print("({})*T^{} + ".format(cc,j)) elif j == 1: print("({})*T + ".format(cc)) elif j == 0: print("({})".format(cc)) print("horner c:") print(horner(c, T)) if k==8: print("horner c using h_u=(h_t+1)/2:") print(str(horner(c(h_t=2*h_t-1),T)).replace('h_t', 'h_u')) count_formula_k8(i, c) elif k==6: count_formula_k6(i, tr, c) print("")