flag_printer picoCTF 2024 Solution

Description

1. Step 1Recognise the shape: polynomial interpolation

   Observation

   encoded.txt holds about 1.77 million x-y integer pairs over a finite field. Those are sample points of an unknown polynomial, not raw data, so the way back is interpolation.

   The challenge ships a transcript of (x, f(x)) pairs over GF(7514777789). A unique polynomial of degree n-1 passes through n distinct points; finding its coefficients is the entire challenge. Once you have them, the bytes of those coefficients spell out the flag.

   Learn more

   Why a polynomial. Plain text is a sequence of bytes; a polynomial is a sequence of coefficients. The challenge author treated the flag-printer's memory image as polynomial coefficients a_0, a_1, ..., a_{n-1} over GF(p), then sampled f(x_i) = sum a_k * x_i^k mod p at n distinct x_i values and shipped the table. Recovering a_0..a_{n-1} is exactly polynomial interpolation.

   Why GF(p). Working in a finite field keeps the coefficients bounded (everything stays under p, which fits in 33 bits). Over the rationals, an interpolating polynomial of degree 1.77 million would have coefficients with billions of digits each, blowing up storage and arithmetic. p = 7514777789 is just over 2^32, large enough to encode any byte sequence as an integer modulo p with no ambiguity, small enough for fast modular arithmetic.

   Why the result is a BMP. When you reconstruct a_0..a_{n-2} and dump them as raw bytes (lowest 8 bits of each, or the first byte of the little-endian encoding), the byte stream starts with 42 4D(BMP magic) and renders as an image of the flag. The challenge name "flag printer" is the giveaway: the polynomial is the flag-printer's firmware, and decoding it prints the flag visually.

2. Step 2Why naive Lagrange is hopeless

   Observation

   With 1.77 million points and a Lagrange formula that scales quadratically, the direct approach never finishes. This needs divide and conquer.

   Use the divide-and-conquer interpolation algorithm (Modern Computer Algebra ch. 10): split the points in half, compute vanishing polynomials Z1, Z2 for each half, recurse, and combine as f1*Z2 + f2*Z1. Total complexity O(n log^2 n).

   Learn more

   The recursion. Given (X, Y) with |X| = n, partition into halves (X1, Y1) and (X2, Y2). Define vanishing polynomials:

   Z1(x) = product over xk in X1 of (x - xk)
   Z2(x) = product over xk in X2 of (x - xk)

   Z1 is zero at every x in X1 and nonzero on X2; symmetrically for Z2. Then the interpolating polynomial is:

   f(x) = f1(x) * Z2(x) + f2(x) * Z1(x)
   
   where f1 interpolates (X1, Y1 / Z2(X1))   <- adjusted Y values
         f2 interpolates (X2, Y2 / Z1(X2))

   On X1 the second term vanishes (because Z1 is zero there) and the first term contributes f1(x) * Z2(x). Pre-dividing Y1 by Z2(X1) exactly cancels the Z2 factor, giving f(x_i) = Y_i. Symmetric on X2.

   Why it ends up O(n log^2 n). The provided code builds each Z with a sequential product loop (O(n^2) polynomial multiplications), which is fast enough in practice with SageMath's compiled arithmetic and 2-way parallelism. A true subproduct tree would reduce this to O(n log^2 n), which is what FLINT implements internally. Evaluating one polynomial at many points (multipoint evaluation) is also O(n log^2 n). Recursion depth is log n. Total for a fully optimised implementation: O(n log^2 n), roughly n * 21^2 = ~440 * n field ops at this scale.

   Library choices. SageMath's R.lagrange_polynomial uses the naive algorithm and is unusable here. FLINT's fmpz_mod_poly_interpolate_fmpz_vec ships the divide-and-conquer version and finishes in around 25 seconds; pure SageMath with a hand-rolled recursion takes ~1 minute. Plain Python or NumPy will not finish.

3. Step 3Implement the recursion, dump coefficients as BMP

   Observation

   The challenge is called flag printer, and every recovered coefficient fits in a byte. The coefficients are image data: write them straight out as a BMP and look at it.

   Run the Sage solver below, pipe coefficients to a BMP file, and open the file in any image viewer. The flag is rendered as image content.

   bash

   sage solve.sage   # ~1 minute on a modern laptop

   Copied!

   bash

   xdg-open output.bmp

   Copied!

   bash

   # or convert to PNG for sharing:

   Copied!

   bash

   convert output.bmp output.png && xdg-open output.png

   Copied!

   Expected output

   picoCTF{...}

   What didn't work first

   Tried: Opening output.bmp while running bytearray(coeffs) on coefficients that include values larger than 255

   A coefficient outside 0 to 255 makes bytearray() raise an overflow error and no file is written. That almost always means the recursion swapped which half gets Z1 and which gets Z2. Y1 divides out Z2, the vanishing product over X2, and Y2 divides out Z1, then the halves recombine as solve(X1, Y1)*Z2 + solve(X2, Y2)*Z1.

   Tried: Using R.lagrange_polynomial() on the full 1.77 million point dataset in SageMath instead of the divide-and-conquer recursion

   SageMath's lagrange_polynomial is the naive quadratic algorithm. At 1.77 million points that is roughly 3 * 10^12 field multiplications, which will not finish on consumer hardware. The script only calls it for subproblems of ten points or fewer, where it is fast, and splits everything above that on the vanishing polynomials Z1 and Z2.

   Learn more

   solve.sage (saved alongside encoded.txt):

   import multiprocessing as mp
   
   p = 7514777789
   X, Y = [], []
   for line in open('encoded.txt').read().strip().split('\n'):
       x, y = line.split(' ')
       X.append(int(x))
       Y.append(int(y))
   
   K = GF(p)
   R = PolynomialRing(K, 'x')
   
   def compZ(X):
       x = R.gen()
       Z = K(1)
       for xk in X:
           Z *= (x - xk)
       return Z
   
   def comp(X, Y, Xother):
       Z = compZ(Xother)
       Y = [y / Z(x) for x, y in zip(X, Y)]
       return Y, Z
   
   def solve(X, Y):
       n = len(Y)
       if n <= 10:
           return R.lagrange_polynomial(list(zip(X, Y)))
       nhalf = n // 2
       X1, Y1 = X[:nhalf], Y[:nhalf]
       X2, Y2 = X[nhalf:], Y[nhalf:]
       if nhalf > 10000:
           with mp.Pool(2) as pool:
               r1 = pool.apply_async(comp, (X1, Y1, X2))
               r2 = pool.apply_async(comp, (X2, Y2, X1))
               Y1, Z2 = r1.get()
               Y2, Z1 = r2.get()
       else:
           Y1, Z2 = comp(X1, Y1, X2)
           Y2, Z1 = comp(X2, Y2, X1)
       return solve(X1, Y1) * Z2 + solve(X2, Y2) * Z1
   
   f = solve(X, Y)
   coeffs = f.coefficients(sparse=False)
   # The interpolated polynomial carries one coefficient more than the image has bytes,
   # so drop the top one. Sanity check: the first two bytes written must be 0x42 0x4D.
   open('output.bmp', 'wb').write(bytearray(coeffs[:-1]))
   print(f'wrote {len(coeffs) - 1} bytes')

   Why the coefficients are already bytes. The encoder maps each byte of the BMP file to a coefficient by zero-extending it to fit in GF(p). So the recovered coefficients live in [0, 256) exactly; bytearray(coeffs) just packs the integer-valued coefficients back to bytes. If your recovered coeffs contain values outside [0, 256), your interpolation is wrong (off-by-one in the recursion is the usual culprit).

   Sanity-test on a slice first. Before running the full n = 1.77M interpolation, run the recursion on the first 1000 pairs and check f(x_i) == y_i for those points. If the recursion is broken, you find out in 0.1 seconds instead of 60.

   See the Python for CTF post for the small-scale Lagrange building blocks and the FLINT manual for the C-level interface if Sage is too slow.

Flag

The flag is rendered as image content inside output.bmp once the polynomial is interpolated; open the BMP in any viewer to read it. The exact characters depend on the encoded.txt your instance ships, but the methodology (divide-and-conquer interpolation with vanishing polynomials, then dump coefficients as bytes) is what unlocks any instance of this challenge.