Can you please explain a bit more how to make this speed up?

We were able to get accepted by simply optimizing the O(n^3 log n) approach.

What is the proof for "image stabilizes" fact? Specially, when we are working in mod prime?

There is additional trick here. We can't use baby step giant step as a blackbox. We need to use fact that we know that our result will be small, so that even though it is modulo q, it works in (because p is upper bound on result we are searching for in this moment) instead of usual (and we know p can be upper bounded by 10⁹ because cases when are easy (this is where this comes from, because we can bound and then we solve BSGS in ))

Btw, I don't think that Tonelli-Shanks is pretty simple xD

Well, I think our is quite ok.

It's quite easy to get a range that is non-decreasing. Let's do a big loop, inside which do following. In another Get next two such ranges, and merge them. If one left, to nothing with it. If only one left total we are already won.

The only problem we had, is we wrote solution for <=, not <, and got problems with stability.

You can do a mergesort. There are only 18 levels of recursion which can be emulated by 18 (identical except for variables) nested loops.

Let's consider the image on a 2D plane. Now assuming that upper right corner is on (2^k, 2^k) and lower left corner is on (0, 0). We can do recursion to assign coordinate to each square. It would be like this:

go(x1, y1, x2, y2):

int x;

cin >> x;

if x == 1:

squares.add({ x1, y1, x2, y2 });

if x == 2:

go(x1, (y1 + y2) / 2, (x1 + x2) / 2, y2);

go((x1 + x2) / 2, (y1 + y2) / 2, x2, y2);

go(...);

go(...);

Now we can iterate over the squares and check whether it connected with some previous squares or not. So we have a graph and there we need to find the number of connected components (we can use DSU instead of storing whole graph in memory).

But the problem is that 2^k is really large to store this number in a long long variable. We can create struct to store coordinates is a binary representation (storing only powers).

https://ideone.com/7ZgrFU

We solved it differently: For every square store 8 pointers (4 for squares inside), and 4 for adjacent squares that are bigger. We can maintain this data just by doing dfs/bfs + using dsu for simplicity. It is O(n).

/***
2 0 2 0 1 1 0 2 0 1 0 1 1
***/
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef long double ld;
int num = 0;
bool ok[20000];
struct DSU
{
	int sz[20000];
	int pa[20000];
	DSU()
	{
		for (int i = 0; i < 20000; i++)
		{
			sz[i] = 1;
			pa[i] = i;
		}
	}
	int root(int a)
	{
		if (a == pa[a])
			return a;
		return pa[a] = root(pa[a]);
	}
	void merge(int a, int b)
	{
		if (a == -1 || b == -1)
			return;
		a = root(a);
		b = root(b);
		if (a == b)
			return;
		if (sz[a] > sz[b])
			swap(a, b);
		sz[b] += sz[a];
		pa[a] = b;
	}
};
DSU dsu;
struct Square
{
	int id;
	int x;
	Square *a = NULL, *b = NULL, *c = NULL, *d = NULL;
	Square *x_ = NULL, *y_ = NULL, *z_ = NULL, *t_ = NULL;
	Square()
	{
		id = num++;
		cin >> x;
		if (x == 2)
		{
			a = new Square();
			b = new Square();
			c = new Square();
			d = new Square();
			b->x_ = a;
			c->x_ = d;

			c->y_ = b;
			d->y_ = a;

			a->z_ = b;
			d->z_ = c;

			a->t_ = d;
			b->t_ = c;
		}
		if (x == 1)
		{
			ok[id] = true;
		}
	}
	void fix()
	{
		if (x == 2)
		{
			{
				d->t_ = t_;
				if (t_ != NULL && t_->x == 2)
				{
					d->t_ = t_->a;
				}
				d->x_ = x_;
				if (x_ != NULL && x_->x == 2)
				{
					d->x_ = x_->c;
				}
			}
			{
				c->t_ = t_;
				if (t_ != NULL && t_->x == 2)
				{
					c->t_ = t_->b;
				}
				c->z_ = z_;
				if (z_ != NULL && z_->x == 2)
				{
					c->z_ = z_->d;
				}
			}
			{
				b->z_ = z_;
				if (z_ != NULL && z_->x == 2)
				{
					b->z_ = z_->a;
				}
				b->y_ = y_;
				if (y_ != NULL && y_->x == 2)
				{
					b->y_ = y_->c;
				}
			}
			{
				a->y_ = y_;
				if (y_ != NULL && y_->x == 2)
				{
					a->y_ = y_->d;
				}
				a->x_ = x_;
				if (x_ != NULL && x_->x == 2)
				{
					a->x_ = x_->b;
				}
			}
			a->fix();
			b->fix();
			c->fix();
			d->fix();
		}
		if (x == 1)
		{
			ok[id] = true;
		}
	}
	void answer()
	{
		if (x == 2)
		{
			a->answer();
			b->answer();
			c->answer();
			d->answer();
		}
		if (x == 1)
		{
			if (x_)
			{
				if (x_->x == 1)
					dsu.merge(id, x_->id);
			}
			if (y_)
			{
				if (y_->x == 1)
					dsu.merge(id, y_->id);
			}
			if (z_)
			{
				if (z_->x == 1)
					dsu.merge(id, z_->id);
			}
			if (t_)
			{
				if (t_->x == 1)
					dsu.merge(id, t_->id);
			}
		}
	}
};
int main()
{
	ios_base::sync_with_stdio(false), cin.tie(NULL), cout.tie(NULL);
	Square A = Square();
	A.fix();
	A.answer();
	int ans = 0;
	for (int i = 0; i < 20000; i++)
		if (ok[i] && dsu.root(i) == i)
			ans++;
	cout << ans << "\n";
}

Backtracking was one possibility, you could also for every triplet of pieces (with order) and all their rotations, reflections and displacements try to create a corner, compute the border of the corner and check whether there is a matching corner from the other three pieces.