# TokyoWesterns - `sqrt` (216pt / 45 solves)

We are given the output of the following script:

```python
from Crypto.Util.number import bytes_to_long, isPrime
from secret import flag, p

def encrypt(m, k, p):
return pow(m, 1 << k, p)

assert flag.startswith("TWCTF{")
assert len(flag) == 42
assert isPrime(p)

k = 64
pt = bytes_to_long(flag.encode())
ct = encrypt(pt, k, p)

with open("output.txt", "w") as f:
f.write(str(ct) + "\n")
f.write(str(p) + "\n")
```

**TL;DR:** We get $c \equiv m^{2^{64}} \pmod p$.

Normally, we could just compute $\phi(p) = p - 1$ and invert the exponent
$e = 2^{64}$ modulo $\phi(p)$. However, we discover that the prime chosen
by the authors has
$$ 2^{30} \mid p-1. $$

Thus, the best we can do is
$$ed' \equiv 2^{30} \pmod{p-1},$$
$$c' \equiv c^{d'} \equiv m^{2^{30}} \pmod p.$$

We can then try to calculate the square root of $c'$ modulo $p$ -- that's
easy when $p$ is prime. However, just like in $\mathbb R$, there are two valid
results -- $\pm \sqrt{c'}$. Because $p \equiv 1 \pmod 4$, *both* results
will be quadratic residues, and thus we need to choose between the two
results after each of the 30 sqrt operations.

$2^{30}$ is not prohibitively big, so we decide to brute-force it.
One optimization we can make here is that
$$\sqrt{-\sqrt{c'}} \equiv \sqrt{-1} \cdot \sqrt{\sqrt{c'}}.$$
We can generalize this with higher roots of $-1$. Thus, a solution looks
like this:

```sage
ct = 5602276430032875007249509644314357293319755912603737631044802989314683039473469151600643674831915676677562504743413434940280819915470852112137937963496770923674944514657123370759858913638782767380945111493317828235741160391407042689991007589804877919105123960837253705596164618906554015382923343311865102111160
p = 6722156186149423473586056936189163112345526308304739592548269432948561498704906497631759731744824085311511299618196491816929603296108414569727189748975204102209646335725406551943711581704258725226874414399572244863268492324353927787818836752142254189928999592648333789131233670456465647924867060170327150559233

K = Zmod(p)
ct = K(ct)

print(hex(p-1))

k = 1<<64
m = 1<<30
q = (p-1)//m
assert q*m == p - 1

e = ZZ(pow(k, -1, q))
#print(e*k*m%(p-1))

ct64 = ct
ct30 = ct^(e*m)

cta = ct30.sqrt()
adj = K(-1)
for _ in range(29):
cta = cta.sqrt()
adj = adj.sqrt()

print('cta =', cta)
print('adj =', adj)
print('p =', p)
```

On a beefier machine, we then bruteforce in pure Python:
```python
import sys
cta = 2870488233024552312211912744984364723362092073401853563020153002240999628038594079978068371406100016947499112660863712384845343971877800926813383568986079674887394443071221968018307459337467225923659696329457180694735414816160530949316647077716726636742679728616967912214567159925834523302439092575850221797093
adj = 3039464130592063363268172491709664258669230391606827037044490175611151346859585269501224540364274702937290928429337888612793876451080615986165398334636621206940191889565327484346806401927835678635869985173316792564575533161011240534381880069725881090619803496841776907428150698981211903041608207198343473228015
p = 6722156186149423473586056936189163112345526308304739592548269432948561498704906497631759731744824085311511299618196491816929603296108414569727189748975204102209646335725406551943711581704258725226874414399572244863268492324353927787818836752142254189928999592648333789131233670456465647924867060170327150559233

threshold = 2**340
cores = 32
percore = 2**30 // cores
n = int(sys.argv[1])
cta = cta * pow(adj, n*percore, p)
for _ in range(percore):
if cta < threshold:
print(cta)
break
cta = cta * adj % p
# 46118657867175077974333582503480021752138026928648315026752964583384841511234055522567479767514232701
```
This obtains the result in about 5 minutes.