Commit a7926c46 authored by MASSON Simon's avatar MASSON Simon

add final expo formulas file

parent 0a360a19
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), ]),
# For i==1, 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), ]),
]
# 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 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("")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment