p = 8444461749428370424248824938781546531375899335154063827935233455917409239041
F = GF(p)
z = F(2841681278031794617739547238867782961338435681360110683443920362658525667816)
# First, solve p - 1 = 2**n * m where n >= 1 and m odd
n = 47
m = (p - 1) / 2**n
# Next, compute g = z^m where z is a quadratic non-residue in Fp
g = z ** m
# Finally, define l_0... l_{k-1} > 0 such that l_0 + .. + l_{k-1} = n-1
k = 6
l = [7, 7, 8, 8, 8, 8]
w = 8
# g powers lookup table: Indexed by power of two, then nu
powers_needed_for_t_i = [0, 8, 16, 24, 32]
powers_needed_for_t_k_over_2 = [40]
# This is Table 1 and 2 from Sarkar 2020 (combined)
gtab = {}
for power_of_two in powers_needed_for_t_i + powers_needed_for_t_k_over_2:
gtab_i = []
for nu in range(0, 2**w):
exponent = nu * F(2) ** power_of_two
gtab_i.append(g**exponent)
gtab[power_of_two] = gtab_i
# s lookup table: Indexed by g^{\nu * 2^{n-l}}
# This is Table 3 from Sarkar 2020
s_lookup_table = {}
for nu in range(0, 2**w):
s_lookup_table[g**(-1 * nu * 2**(n-w))] = nu
def sqrt_ratio(num, den):
s = den**(2**n - 1)
t = s**2 * den
w = (num * t)**((m-1)/2) * s
v = w * den;
uv = w * num;
assert v == (num / den)**((m - 1) / 2)
assert uv == (num / den)**((m + 1) / 2)
return table_based_findSqRoot_sarkar(uv, v)
def table_based_findSqRoot_sarkar(uv, v):
"""Implements decaf377 Sarkar 2020 based method"""
x = uv * v # Such that x = u**m
x5 = x
x4 = x5 ** (2**l[5])
x3 = x4 ** (2**l[4])
x2 = x3 ** (2**l[3])
x1 = x2 ** (2**l[2])
x0 = x1 ** (2**l[1])
# i = 0
# Here we want to look up s_0 = q_0 * g**{n-7}, but our table has
# powers of g**{n-8}, so we're actually looking up q_prime = 2*q_0.
q_0_prime = s_lookup_table[x0] # Since x0 = alpha0
t = q_0_prime
# i = 1
alpha_1 = x1 * gtab[32][t]
#breakpoint()
q_1_prime = s_lookup_table[alpha_1]
# The left shift values come from the final expression for t:
# t_6 = q_0_prime + q_1_prime * 2^7 + ... + q_5 * 2^39
t += q_1_prime << 7
# i = 2
alpha_2 = x2 * gtab[24][t & 0xFF] * gtab[32][t >> 8]
assert (t >> 8) & 0xFF == q_1_prime // 2
q_2 = s_lookup_table[alpha_2]
t += q_2 << 15
# i = 3
alpha_3 = x3 * gtab[16][t & 0xFF] * gtab[24][(t >> 8) & 0xFF] * gtab[32][(t >> 16) & 0xFF]
assert (t >> 16) & 0xFF == q_2 // 2
q_3 = s_lookup_table[alpha_3]
t += q_3 << 23
# i = 4
alpha_4 = x4 * gtab[8][t & 0xFF] * gtab[16][(t >> 8) & 0xFF] * gtab[24][(t >> 16) & 0xFF] * gtab[32][(t >> 24) & 0xFF]
assert (t >> 24) & 0xFF == q_3 // 2
q_4 = s_lookup_table[alpha_4]
t += q_4 << 31
# i = 5
alpha_5 = x5 * gtab[0][t & 0xFF] * gtab[8][(t >> 8) & 0xFF] * gtab[16][(t >> 16) & 0xFF] * gtab[24][(t >> 24) & 0xFF] * gtab[32][(t >> 32) & 0xFF]
assert (t >> 32) & 0xFF == q_4 // 2
q_5 = s_lookup_table[alpha_5]
t += q_5 << 39
assert t == q_0_prime + q_1_prime * 2**7 + q_2 * 2**15 + q_3 * 2**23 + q_4 * 2**31 + q_5 * 2**39 # Lemma 4 assertion
t = (t + 1) >> 1;
# Take 8 bits at a time, e.g. (t & 0xFF) is taking the last 8 bits of t to yield a value from 0-255, which are the allowed values in each g lookup table
gamma = gtab[0][t & 0xFF] * gtab[8][(t >> 8) & 0xFF] * gtab[16][(t >> 16) & 0xFF] * gtab[24][(t >> 24) & 0xFF] * gtab[32][(t >> 32) & 0xFF] * gtab[40][(t >> 40)]
y = uv * gamma
if q_0_prime % 2 != 0:
y *= z**((1-m) / 2)
return y
def test_sqrt(n):
for _ in range(n):
num = F.random_element()
den = F(1)
u = num / den
res = sqrt_ratio(num, den)
if u.is_square():
assert res**2 == u
else:
assert res**2 == u * z
#test_sqrt(10000)
test_sqrt(100)