The problem can be abstractly specified as this:

• Let $p=2^{255} - 19$ be a prime and let $s_0, a, b \in \mathbb{F}_p$ be unknown constants. Let $f(x) = \sum_{k=0}^{14} s_k x^k$ be a degree-14 polynomial over $\mathbb{F}_p$ determined by the following relation:

  • For $i \geq 0$, $s_{i+1} = a s_i + b$. You can choose exactly 14 input values $0 < x_1, \,\dots,\, x_{14} < p$ and obtain their $y$-values $(x_i, y_i) = (x_i, f(x_i)) \in \mathbb{F}_p^2$.

  Using this information, determine $s_0$.

We can encode our setup in SageMath, writing $n = \deg(f) = 14$.

p = 2^255 - 19
n = 14
F = GF(p)
P.<s0, a, b> = PolynomialRing(F) # multivariate polynomial ring for the unknowns
Px.<x> = PolynomialRing(P) # F[s0,a,b][x]

My first instinct was to reach for Lagrange interpolation, which is the normal way of reconstructing a polynomial from sampled points; but for a degree-14 polynomial Lagrange interpolation requires 15 points, and we only have 14. Hmm.

Well, it seems pretty clear at least that it doesn’t really matter what we choose for the $x$-values, they won’t give us any extra information, so we’ll just query the simple case of 1 through 14:

xs = [F(v) for v in [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]]
ys = [
    4004031469935226797678041300513305517447093488292866886390744385810723152992,
    # et cetera
]
ys = [F(v) for v in ys]

If we perform “partial” Lagrange interpolation with 14 points instead of 15, we should get an affine space that represents all possible 15-degree polynomials under the provided constraints, and then maybe we could use the coefficient constraints.

xs.append(0)
ys.append(s0)
basis = []
for j in range(n+1):
    basis.append(
        prod((x-xs[m])*(xs[j]-xs[m]).inverse() for m in range(n+1) if m != j)
    )
undet_poly = sum(ys[i]*basis[i] for i in range(n+1))

Now our undet_poly (which we will denote $g$) is some expression $g \in \mathbb{F}_p [s_0][x]$. We can read off each coefficient:

coefs = [undet_poly.coefficient(i) for i in range(n+1)]
assert coefs == [
    s0,
    14944267905132978895707490166157624493267307919363814442860231562415439532290*s0 + 42734869096716633894829933192697273711142800957529494279467966422605761677192,
    37886807707385071832844031387479480962153490977325705939809345100249480239084*s0 + 46271650052041060360191158582510170254941583025029928660681043770310086080314,
    # et cetera
]

Notice how every entry in the coefficient array is some linear equation on $s_0$. This makes sense; Lagrange interpolation works via first constructing a unique basis and then summing up that basis to produce the final result, so we expect to get something like a linear or affine subspace, and in this case it’s an affine subspace.

If we denote $g$ as $g(x) = \sum_{i=0}^{14} c_i x^i$ where each $c_i$ corresponds to coefs[i], we can represent each $c_i$ as a linear constraint from our array: we write $c_i = m_i s_0 + v_i$ for some known $m_i,v_i \in \mathbb{F}_p$.

ms = [c_i.coefficient({s0: 1}) for c_i in coefs]
vs = [c_i.coefficient({s0: 0}) for c_i in coefs]

This looks like a natural place to apply the constraint $s_{i+1} = a s_i + b$, which translates to $c_{i+1} = a c_i + b$ here. To do this, we first write a function to compute $s_i$ in terms of $s_0$:

# finds s_k in terms of s_0
# ac, bc are the linear coefficients, which default to a,b. we make them params here just in case
def lcgiter(k, s0, ac=a, bc=b):
    if k == 0:
        return s0
    else:
        return lcgiter(k-1, ac*s0+bc, ac, bc)
assert lcgiter(2, s0, ac=4, bc=5) == 4*(4*s0 + 5) + 5

Then $s_k = \operatorname{lcgiter}(k,s_0)$ is a polynomial over $\mathbb{F}_p [a,b,s_0]$, and we can phrase our constraints as $m _i s_0 + v_i = \operatorname{lcgiter}(i,s_0)$. In SageMath, that’s phrased as $m_i s_0 + v_i - \operatorname{lcgiter}(i,s_0) = 0$:

constraints = [ms[i]*s0 + vs[i] - lcgiter(i, s0) for i in range(n+1)]

Now we have a system of polynomial equations in $\mathbb{F}_p [a,b,s_0]$ that we want to solve. Er, how do we do that, again?

I got stuck a little bit at this point. In general, solving polynomial equations over finite fields appears to be quite difficult. To quote (Dell et al., 2024):

Theorem 1.3. If $\text{SETH}$ holds [i.e. the NP-complete problem $k$-SAT cannot be solved in subexponential time $O(2^{o(n)})$ for all $k$, and further as $k \to \infty$ we have the time complexities approaching $O(2^n)$], then for all prime powers $q$ and all rationals $δ > 0$, there exists $d ∈ \mathbb{N}$ such that there is no $O(q^{(1−δ)n})$-time algorithm for $\text{PES}^d_q$ [deciding whether or not any system of polynomial equations of degree at most $d$ over a finite field of size $q$ with $n$ variables has a root or not.]

So this is NP-hard, at least. Not great. Do we have anything better? The multivariate version of Coppersmith’s method won’t work; $s_0,a,b$ are too large.

I struggled with this for half an hour, and then I came up with a stupid idea.

• We have efficient algorithms to solve polynomial equations over $\mathbb{Q}$. What if we “ported” those solutions back into $\mathbb{F}_p$? That is, we found a function $\psi$ that took solutions to $\mathbb{Q}$-polynomials and returned equivalent constants in $\mathbb{F}_p$?

As soon as I thought of it I was 75% convinced it wouldn’t work. The worlds of $\mathbb{Q}$ and $\mathbb{F}_p$ are pretty different: one classic example is that the polynomial $x^4 + 1$ is irreducible over $\mathbb{Q}$ but is reducible in $\mathbb{F}_q$ for every $q$. But I didn’t have any better ideas, so… why not?

The first question is: how do we even “port” solutions over to $\mathbb{F}_p$? It’s a little more difficult than it sounds, because systems of polynomial equations over $\mathbb{Q}$ have solutions in what is called the algebraic numbers $\mathbb{A}$, which isn’t just rationals but involves things like square roots, $i$, and suchlike.¹ I came up with the following construction:

• Given $u \in \mathbb{Z}$, the equivalent in $\mathbb{F}_p$ is passing $u$ through the canonical homomorphism $\mathbb{Z} \to \mathbb{F}_p$, so "$\psi(u) = u$".
• Given $u+v$ or $u v$ for some $u,v \in \mathbb{A}$, then the equivalent in $\mathbb{F}_p$ is $\psi(u+v) = \psi(u) + \psi(v)$ and $\psi(u v) = \psi(u) \psi(v)$ respectively. (This means it’s almost a ring homomorphism.)
• Given $u/v$ for some $u,v \in \mathbb{A}$, then the equivalent in $\mathbb{F}_p$ is $\psi(u/v) = \psi(u) \psi(v)^{-1}$.
• Given $\sqrt[k]u$ for some $u \in \mathbb{A}$, the equivalent in $\mathbb{F}_p$ is $\psi(\sqrt[k]u) = \alpha$ where $\alpha^k = \psi(u) \in \mathbb{F}_p$. (You may wonder how we actually compute $\alpha$. The answer is “call SageMath’s inbuilt algorithm for doing exactly this”.)

The obvious problem with this is cases where e.g. we attempt to take the inverse of $\psi(v) = 0$ or we take the $k$-th root of something that doesn’t have a $k$-th root. But I didn’t have any better ideas, so…

def psi(form):
    try:
        return F(form)
    except:
        pass
    if hasattr(form, 'match'):
        w0, w1 = SR.wild(0), SR.wild(1)
        if t := form.match(w0/w1):
            return psi(t[w0])*(psi(t[w1])^(-1))
        if t := form.match(sqrt(w0)):
            # it turns out k-th roots in generality were not needed
            # for reasons we'll see later
            return psi(t[w0]).sqrt()
        if t := form.match(w0*w1):
            return psi(t[w0])*psi(t[w1])
        if t := form.match(w0+w1):
            return psi(t[w0])+psi(t[w1])
    if (not isinstance(form, Expression)) or form.is_constant():
        return psi(form.real()) + psi(form.imag())*F(-1).sqrt()
    assert False

assert psi(3) == F(3)
assert psi(3/4) == F(3)*(F(4)^(-1))
assert psi(3*i) == F(-9).sqrt()
assert psi(sqrt(3))^2 == F(3)
assert psi(3/4)*F(4) == F(3)
assert psi(3*i/sqrt(2))^2 * 2 + 9 == F(0)

We then solve the system of equations² over $\mathbb Q$…

s0var, avar, bvar = var('s_0 a b')
constraints_qq = [(ms[i].change_ring(ZZ)*s0var + vs[i].change_ring(ZZ) - lcgiter(i, s0var, ac=avar, bc=bvar)) for i in range(4)]
solns = solve(constraints_qq, [s0var,avar,bvar])
# solns[i][0].rhs() contains the i-th solution for s[i]

…and hope it translates to $\mathbb{F}_p$.

s0_guess = psi(solns[0][0].rhs())
print('guessed s0 is', s0_guess)

We then attempt this against the server, and… huh? It actually got the correct value!

It was unclear to me how this worked, but I was in the middle of a CTF, so I just submitted this and moved on. But I came back to this much later for curiosity’s sake, and figured out a couple of very interesting things that makes this work.

Let’s try this formally.

The formal description

This entire thing above relies on our magical $\psi$. Can we obtain something more intuitive and/or formalizable? Let’s try a restriction $\phi: \mathbb{Q} \to \mathbb{F}_p$ of our $\psi$ that only works on the rational numbers. Specifically, we’ll pass integers through the canonical homomorphism: for any $a/b \in \mathbb{Q}$, we have

$$\phi(a/b) := a b^{-1} \in \mathbb{F}_p.$$

This is almost a ring homomorphism.³ If you calculate it out, it turns out it preserves addition, multiplication, and sends $\phi(0) = 0$ and $\phi(1) = 1$. So it’s a ring homomorphism, right?

Not quite, it’s not actually well-defined! If $p$ divides $b$ (which I will write as $p\,|\,b$), then $b^{-1}$ is not defined. That’s a bit sad, isn’t it? We’ve got this function that almost perfectly looks like a ring homomorphism except for a few strange cases. Can we make it into one?

Definition. The $p$-localization of the integers, denoted $\mathbb{Z}_{\langle p \rangle}$, is the ring

$$\mathbb{Z}_{\langle p \rangle} = \{a/b : a, b \in \mathbb{Z} \land p \cancel| b\}.$$

(This is a subring of $\mathbb{Q}$.) Now our $\phi: \mathbb{Z}_{\langle p \rangle} \to \mathbb{F}_p$ is in fact a well-defined ring homomorphism, and if we want we can find $\ker(\phi)$ in order to exhibit an isomorphism

$$\mathbb{Z}_{\langle p \rangle}/\ker(\phi) \cong \mathbb{F}_p.$$

That’s kind of cool. We can also extend $\phi: \mathbb{Z}_{\langle p \rangle} \to \mathbb{F}_p$ to a homomorphism between the polynomial rings $\phi: \mathbb{Z}_{\langle p \rangle}[x] \to \mathbb{F}_p[x]$ (by simply defining $\phi(x) = x$.)

Our first glimmer at the idea that there’s something more going on here comes from the following result:

Theorem. (“Casting from $\mathbb{Z}_{\langle p \rangle}$ to $\mathbb{F}_p$”) Let $f \in \mathbb{Z}_{\langle p \rangle}[x]$. If there exists a root $r \in \mathbb{Z}_{\langle p \rangle}$ of $f$, so $f(r) = 0$, then $\phi(f)(\phi(r)) = 0$. [That is, $r$ represents a root of $\phi(f)$ too.]

The proof for this is just symbolic manipulation: if $f = \sum a_k x^k$, then

$$\begin{align*} \phi(f) (\phi(r)) &= \phi\left(\left. \sum_k a_k x^k\right)\right |_{x=\phi(r)} \\ &= \left. \left(\sum_k \phi(a_k) x^k \right)\right|_{x=\phi(r)} \\ &= \left(\sum_k \phi(a_k) \phi(r)^k\right) \\ &= \sum_k \phi(a_k r^k) \\ &= \phi\left(\sum_k a_k r^k\right) \\ &= \phi(f(r)) = \phi(0) = 0. \end{align*}$$

We can of course generalize this to multivariate polynomials. This is no different from roots in $\mathbb{Z}$ implying roots mod $m$, or in general ring homomorphisms giving a different way to calculate roots. But this doesn’t explain why our solution worked. Let’s look at the calculated solutions in $\mathbb{A}$ to our system of polynomial equations:

assert solns[0][0].rhs() == 39/640059776070390094052306977533233511092604974660525289112584478198039887405245516535313864670250703723826685386742649473790744636126116060802233887780560
    *I*sqrt(429122645378233701296997495540159683865547910093220985770373383068324630497859733686567258077432424497353442676283400008363754666322262870389855882732021296173375620223858447371334461025946899987665926401904006640013983854905620778281091749741863580387784735415962295887183345877535637768905955156060359) - 483752157230349238924830438452111260772096395601048258144831061752835632956860155710293579160016884871928338135207597017194378422059113065051179141020461/640059776070390094052306977533233511092604974660525289112584478198039887405245516535313864670250703723826685386742649473790744636126116060802233887780560

There’s an I (SageMath $i = \sqrt{-1}$) and a sqrt in there! That’s definitely not over $\mathbb{Z}_{\langle p \rangle}$. What happened?

Hmm, maybe we can try to find the most obvious reason this could be true?

Theorem. (“Casting from $\mathbb{Z}_{\langle p \rangle}[\sqrt a]$ to $\mathbb{F}_p$”) Let $a \in \mathbb{Z}_{\langle p \rangle}$ be such that $\phi(a)$ is a quadratic residue in $\mathbb{F}_p$ (i.e. $\exists \beta \in \mathbb{F}_p$ such that $\beta^2 = \phi(a)$). Then there exists a ring homomorphism $\psi: \mathbb{Z}_{\langle p \rangle}[\sqrt a] \to \mathbb{F}_p$.

This is actually true!

Proof: Define $\psi(\sqrt a) = \beta$ in $\mathbb{F}_p$ and $\psi(t) = \phi(t)$ for $t \in \mathbb{Z}_{\langle p \rangle}$. Then, given any $t + s \sqrt a \in \mathbb{Z}_{\langle p \rangle}[\sqrt{a}]$, define

$$\psi(t + s \sqrt a) = \psi(t) + \psi(s) \psi(\sqrt a).$$

It’s not hard to check this is in fact a well-defined ring homomorphism. $\square$

In fact, this kind of construction works if you simply just keep extending: there exists a homomorphism $\mathbb{Z}_{\langle p \rangle} [\sqrt{a_1},\, \dots,\, \sqrt{a_k}] \to \mathbb{F}_p$ when all $a_1,\, ...,\, a_k$ are quadratic residues in $\mathbb{F}_p$.

So, since our system of equations had a root over $\mathbb{Z}_{\langle p \rangle} [\sqrt{-1}, \sqrt{a}]$ for some gargantuan $a$ being a quadratic residue of $\mathbb{F}_p$, it also had a root over $\mathbb{F}_p$, which was exactly what we wanted!

So why did I get lucky enough that this construction actually worked in the real CTF? It doesn’t seem to be luck at all. In fact, if I keep randomly generating polynomials and cracking them my code has a 100% rate of finding a solution (concluded from 1000 samples). This is extremely surprising to me, and I have no idea why this happens. (At a first guess, it could be because a polynomial being reducible over $\mathbb{Z}_{\langle p \rangle}$ implies it’s also reducible over $\mathbb{F}_p$? But this doesn’t explain why we never run into polynomials that don’t exist in $\mathbb{Z}_{\langle p \rangle}$: that is, polynomials which have a coefficient with $p$ present as a factor of the denominator. I would love to know the answer, if anyone can figure it out!)

Postscript

Later, after the CTF had concluded, someone else posted this solution:

Ideal(constraints).groebner_basis() # gives you the exact values you want

I should have probably thought of this. Sigh. Oh well.

Footnotes

1. To be more specific, we have $\mathbb{A} \subsetneq \mathbb{C}$. If $X$ is the set of all values achievable by some combination of addition, multiplication, additive/multiplicative inverses, and $k$-th roots, then we have $X \subsetneq \mathbb{A}$ (because of things like the unsolvability of quintics in radicals.) ↩

2. Why am I restricting to the first 4 constraints (range(4))? Because the unsolvability of the quintic is well-known, and I want exact solutions in radicals. At this point I just hoped I would get something reasonable out of this. ↩

3. Quick fun exercise: can you rule out $\phi$ being a ring homomorphism without doing any computations? There’s a one-sentence argument. ↩