I don't know what is the fastest way to solve this, but I think mine is more complicated than it can be. I've keeped a vector for each remainder i mod k and i*i mod k, to quickly find the numbers which give remainder r. Then you can add edge by edge in graph and use union-find algo to check for cycles. Maybe there is a faster math solution?

Yep, graph is very sparse. Almost always (without several trivial cases) you can consider only range [0..2*k]. To avoid any kinds of checking simply get n = max(10, 2*k). For each number in [1..n] try to find cycle with previous ones. Use e.g. disjoint sets for simplicity. Note, that solution always exists!

Is there any mathematical proof for this?

if k=1 =>ans is 3
if k=2 =>ans is 5
if k>2 => you can take numbers 1, k-1, 2*k-1, and it will be a cycle