# Solving Wii Phit via Pell's Equations

## Overview

This is a challenge from the [2021 Cyber Apocalpyse CTF](https://www.hackthebox.eu/cyber-apocalypse-ctf-2021). In this challenge, we are given a flag encrypted under a strange RSA modulus (of the form $p^3q$):

N = p**3 * q
e = 0x10001
c = pow(bytes_to_long(FLAG),e,N)

print(f'Flag: {hex(c)}')

# Hint

w = 25965460884749769384351428855708318685345170011800821829011862918688758545847199832834284337871947234627057905530743956554688825819477516944610078633662855
x = p + 1328
y = p + 1329
z = q - 1

assert w*(x*z + y*z - x*y) == 4*x*y*z

We are also given a "hint" about $p$ and $q$, a Diophantine equation $p$ and $q$ satisfy. This is in fact all the information we are given about $p$ and $q$ -- we are not even given $N$ (although we can roughly estimate its size from $c$).

This writeup will discuss how to solve this by reducing this Diophantine equation to a [Pell's equation](https://en.wikipedia.org/wiki/Pell%27s_equation). I will emphasize that this is *not* a particularly good method for solving this specific challenge: the computations are fairly tedious, and at the end it will essentially reduce to guessing the correct form of the solution (see [here](https://blog.cryptohack.org/cyber-apocalypse-2021#wii-phit) for the much easier intended solution). Nonetheless this technique is sometimes applicable in other situations, so hopefully this write-up will still have some use.

## Reducing to a Pell's equation.

We start with the diophantine equation $W(xz + yz - xy) = 4xyz$ (I will capitalize $W$ to make it clear that it is a constant). Moreover, we know that in the solution we seek that $y = x+1$. Substituting this in, we thus get:

$$W(xz + (x+1)z - x(x+1)) = 4x(x+1)z.$$

The first thing to notice is that this is a linear equation in $z$. So let's solve for $z$. We can rearrange this to get:

$$z = \frac{Wx(x+1)}{W(2x+1) - 4x(x+1)}$$

Since $z$ is an integer, we want the RHS of this equation to be an integer. The RHS of this equation is the quotient of two quadratics in $x$. We can make our life slightly easier by performing polynomial division to remove the quadratic term in the numerator:

$$4z = -W + \frac{W^2(2x+1)}{W(2x+1) - 4x(x+1)}$$

We now want $W^2(2x+1)$ to be divisible by $W(2x+1) - 4x(x+1)$. Let's write

$$W^2(2x + 1) = k(W(2x+1) - 4x(x+1))$$

for some integer $k$ (note that this means we can write the equation two lines above as $4z = -W + k$). This is a quadratic equation in $x$. Rearranging, we get:

$$4kx^2 + (2W^2 - 2Wk + 4k)x + (W^2 - Wk) = 0$$

In order for this quadratic to have integer roots, its discriminant must be a perfect square. Letting $a^2$ denote the discriminant of this quadratic equation, we can write:

$$(2W^2 - 2Wk + 4k)^2 - 4(W^2-Wk)(4k) = a^2.$$

This is finally a quadratic in both $k$ and $a$; i.e. it is a [quadratic form](https://en.wikipedia.org/wiki/Quadratic_form). This means that we can rewrite it to make it into a diagonal quadratic form. More specifically, after expanding the LHS, writing it as a quadratic in $k$, and completing the square, we can rewrite the above equation as

$$((2W^2 + 8)k - 2W^3)^2 + 16W^4 = (W^2 + 4)a^2.$$

Let $b = (2W^2 + 8)k - 2W^3$. We've finally arrived at our Pell's equation:

$$b^2 - (W^2 + 4)a^2 = -16W^4.$$

Before we continue to solving it, there is one more simplification we can make here that I did not make when originally solving this challenge, but that will help things somewhat. Note that both $b$ and $-16W^4$ are even, but $W^2 + 4$ is odd (the given $W$ is odd). This means $a$ must be divisible by $2$. Writing $\alpha = a/2$ and $\beta = b/2 = (W^2 + 4)k - W^3$, this simplifies further to

$$\beta^2 - (W^2 + 4)\alpha^2 = -4W^4$$

## Solving a Pell's equation

There are many good online resources that describe how to solve Pell's equations (a quick search turns up [this](https://crypto.stanford.edu/pbc/notes/ep/pell.html), which seems to be a nice reference). Here is a very quick, high-level summary of what's going on:

- The expression $x^2 - Dy^2$ can be thought of as the "norm" of the number $x + y\sqrt{D}$. There is a precise algebraic sense in which this is true, but the relevant fact for us is that if you write $N(x + y\sqrt{D}) = x^2 - Dy^2$, then $N$ "respects multiplication": $N((x_1 + y_1\sqrt{D})(x_2 + y_2\sqrt{D})) = N(x_1 + y_1\sqrt{D})N(x_2 + y_2\sqrt{D}).$
- In particular, the integer solutions to $x^2 - Dy^2 = 1$ form an abelian group (the group of units in $\mathbb{Z}[\sqrt{D}]$). It turns out this group is generated by a single element, called the "fundamental solution". It is possible to find this solution by looking at the continued fraction expansion of $\sqrt{D}$. This isn't guaranteed to be efficient (it depends on the period of the continued fraction expansion of $\sqrt{D}$), but it is guaranteed to be efficient if there exists a small solution $(x, y)$ (that takes polylogarithmic time in $|x|$ and $|y|$).
- If you want to find solutions to $x^2 - Dy^2 = N$, first find some solution $(x_0, y_0)$ to $x^2 - Dy^2 = N$. Then you can multiply $(x_0, y_0)$ with any solution to $x^2 - Dy^2 = 1$ (mapping them to elements of $\mathbb{Z}[\sqrt{D}]$ first) to get a bunch of solutions to $x^2 - Dy^2 = N$.
- Unfortunately, unlike the $N = 1$ case, there isn't a single "fundamental solution" to $x^2 - Dy^2 = N$; there can be many, and in general it depends on the factorization of $N$. Moreover (as far as I know), there isn't an especially good way to find such solutions aside from guessing or enumerating all $x$ and $y$ in some range.
- However, for the special case of $N = -1$, there is a continued fraction approach to find a solution of $x^2 - Dy^2 = -1$ if one exists. This will be useful to us.

In our problem, $D = W^2 + 4$. Let's start by finding fundamental solutions to $\beta^2 - D\alpha^2 = \pm 1$. As mentioned, this can be done via continued fractions methods (see the above link). Here is a quick implementation in Sage:

# finds a solution (x, y) to x^2 - Dy^2 = 1 (if one exists)
def solve_pell (D):
    cf = continued_fraction(sqrt(D))
    for i in range(10):
        denom = cf.denominator(i)
        numer = cf.numerator(i)
        if numer^2 - D*denom^2 == 1:
            return numer, denom
    return None, None

# finds a solution (x, y) to x^2 - Dy^2 = -1 (if one exists)
def solve_pell_neg(D):
    cf = continued_fraction(sqrt(D))
    for i in range(10):
        denom = cf.denominator(i)
        numer = cf.numerator(i)
        if numer^2 - D*denom^2 == -1:
            return numer, denom
    return None, None

D = w^2 + 4
bu, au = solve_pell(D)
bn, an = solve_pell_neg(D)

From these fundamental solutions we can now generate all solutions to $\beta^2 - D\alpha^2 = \pm 1$. Let's now focus on the value of $N$. As mentioned, there's not an especially quick general way to find solutions to $\beta^2 - D\alpha^2 = \pm N$. But in our case, $D$ and $N$ have a lot of structure, and there are a couple solutions we can immediately eyeball.

Specifically, recall that we have the equation $\beta^2 - (W^2 + 4)\alpha^2 = -4W^4$. If we substitute $(\beta, \alpha) = (W, 1)$, we find a solution with norm $-4$; by scaling $\alpha$ and $\beta$ up by $W^2$ (which scales up the norm by $W^4$) and multiplying by any of our norm $1$ solutions, we can find a family of solutions to our equation. If we substitute $(\beta, \alpha) = (2, 1)$, we find a solution with norm $-W^2$; we can likewise scale this up by $2W$ and multiply it with any norm $1$ solution. Finally, we can scale up any of our solutions with norm $-1$ by $2W^2$ to find a solution with norm $-4W^4$ (and get another family of solutions this way). I don't immediately see any other classes of solutions, but it would not be surprising if they exist (there are almost definitely other classes of solutions but with ugly parameterizations).

Note that although all these solutions are valid solutions to $\beta^2 - D\alpha^2 = -4W^4$, they don't necessarily lead to valid values for $x$ and $z$. Importantly, $x$ and $z$ must be integral, positive, and be specific offsets from a prime. We can enumerate solutions from these families and check them. When we do, we quickly find a solution (coming from that last family of solutions).

Interestingly, this solution doesn't use any powers of the fundamental norm 1 solution (although it does use the norm -1 solution, so the continued fraction step was in some sense necessary). That said, given the official writeup, there should be an explicit polynomial form for the norm -1 solution as well (i.e. that doesn't require continued fractions to find).

## Code

from Crypto.Util.number import *

w = 25965460884749769384351428855708318685345170011800821829011862918688758545847199832834284337871947234627057905530743956554688825819477516944610078633662855

# finds a solution (x, y) to x^2 - Dy^2 = 1 (if one exists)
def solve_pell (D):
    cf = continued_fraction(sqrt(D))
    for i in range(10):
        denom = cf.denominator(i)
        numer = cf.numerator(i)
        if numer^2 - D*denom^2 == 1:
            return numer, denom
    return None, None

# finds a solution (x, y) to x^2 - Dy^2 = -1 (if one exists)
def solve_pell_neg(D):
    cf = continued_fraction(sqrt(D))
    for i in range(10):
        denom = cf.denominator(i)
        numer = cf.numerator(i)
        if numer^2 - D*denom^2 == -1:
            return numer, denom
    return None, None



D = w^2 + 4
bu, au = solve_pell(D)
bn, an = solve_pell_neg(D)

FAMILIES = [(w^3, w^2), (4*w, 2*w), (2*w^2*bn, 2*w^2*an)]


def next_pell(b, a):
    return (b*bu + D*a*au), (b*au + a*bu)

def check_sol(beta, alpha):
    b = 2*beta
    a = 2*alpha

    assert (b*b - (w*w + 4)*a*a) == -16*(w**4)

    k = (b + 2*w^3)/(2*w^2 + 8)
    z = (k-w)/4

    B = (2*w^2 - 2*w*k + 4*k)
    x = (-B + a)/(8*k)

    y = x+1
    assert w*(x*z + y*z - x*y) == 4*x*y*z

    if int(x) !=x:
        return None
    
    if int(z) !=z:
        return None

    if int(x) == x and int(z) == z:
        print("Candidate: ", x, z)
        return (x, z)


def find_xz():
    candidates = []
    for b, a in FAMILIES:
        for i in range(5):
            ans = check_sol(b, a)
            if ans is not None:
                candidates.append(ans)
            b, a = next_pell(b, a)
    return candidates
        

print(find_xz())

X = 12982730442374884692175714427854159342672585005900410914505931459344379272923599916417142168935973617313528952765371978277344412909738758472305039316831427 
Z = 4376511920801673769046982367789644084746600661635151104602579081967083768976309788885633491753761209012042953502416064276555378570438196809829053232168930363213412874907199642703512833211163084612185095343783021125928854760406110492494016250237689683218940269389627326164130063600050024089561671913076715913062539181517022517918888557966172136200366952844819842691032260278270002343732370611629511112677638724984890716752906711541254394792207074699964511035503380

def solve(x, z):
    q = z+1
    p = x-1328

    FLAG = 0x12f47f77c4b5a72a0d14a066fedc80ba6064058c900a798f1658de60f13e1d8f21106654c4aac740fd5e2d7cf62f0d3284c2686d2aac261e35576df989185fee449c20efa171ff3d168a04bce84e51af255383a59ed42583e93481cbfb24fddda16e0a767bff622a4753e1a5df248af14c9ad50f842be47ebb930604becfd4af04d21c0b2248a16cdee16a04b4a12ac7e2161cb63e2d86999a1a8ed2a8faeb4f4986c2a3fbd5916effb1d9f3f04e330fdd8179ea6952b14f758d385c4bc9c5ae30f516c17b23c7c6b9dbe40e16e90d8734baeb69fed12149174b22add6b96750e4416ca7addf70bcec9210b967991e487a4542899dde3abf3a91bbbaeffae67831c46c2238e6e5f4d8004543247fae7ff25bbb01a1ab3196d8a9cfd693096aabec46c2095f2a82a408f688bbedddc407b328d4ea5394348285f48afeaafacc333cff3822e791b9940121b73f4e31c93c6b72ba3ede7bba87419b154dc6099ec95f56ed74fb5c55d9d8b3b8c0fc7de99f344beb118ac3d4333eb692710eaa7fd22

    N = p**3 * q
    e = 0x10001
    d = inverse(e, (p-1)*p*p*(q-1))
    pt = pow(FLAG, d, N)
    print(long_to_bytes(pt))

solve(X, Z)