I haven't understood Pollard's rho. I decided to do away with the polynomial requirement and see what happens. I started with $N = 33.$ I selected at random the numbers $7, 8, 1, 3, 26, 31 \mod\ 33$ as my $x$ sequence. I kept track of the $y$ sequence $\mod 3$ which was $1, 2, 1, ...$ Bingo: we found a collision with $y_1$ and $y_3.$ Let's see if we have $|x_1 - x_3|$ revealing a factor of $33$: $|7 - 1| = 6$ and $\gcd(6, 33) = 3.$
Maybe that was just luck. The fact I chose $3$ as a factor gives me too great of a chance perhaps.
I then picked $13\times19 = 247.$ The random numbers of the $x$ sequence were $57, 128, 231, 24 ...$ The $y$ sequence $\mod\ 13$ was $5, 11, 10, 11 ...$ Found a collision at $y_2, y_4.$ Does $|x_2 - x_4|$ reveal a factor? $|128 - 24| = 104$ and $\gcd(104, 247) = 13.$
Still lucky perhaps.
I then took $N = 57\times257.$ The $x$ sequence was $$10848, 9291, 9501, 4055, 12300, 13061, 12492, 13860.$$ The $y$ sequence, which is $x \mod\ 57$ is $$18, 0, 39, 8, 45, 8, 9, 9, ...$$
We have a collision between $y_4$ and $y_6.$ Does $|x_4 - x_6| = |4055 - 13061| = 9006$ reveal a factor? $\gcd(9006, 14649) = 57.$
Is this still luck?
The polynomial must make a difference and I'm totally ignoring it and still finding factors. Can you guide me into some light here? Thank you.
The Wikipedia article is pretty good at explaining that the function of the polynomial is to be a pseudorandom source of $x_n$s. Since collision relies on the birthday paradox, any other suitable random sequence should serve. The advantage of polynomials is that we can regenerate the same sequence anywhere we like. The Floyd and Brent cycle-finding algorithms allow us to find cycles in the values of the polynomials $\mod N$ keeping track of only two $x_i$ at a time. Finding collisions by your method requires keeping vastly more $x_i$ ($O(\sqrt{N})$, in general).
Note that the runtime of Pollard's rho is frequently listed as $\sqrt{N}$, but that is because one does not know a priori $p$, the least prime factor of $N$, so writing the more accurate runtime, $\sqrt{p}$ is noninformative. Your experiments also utilize information about $p$, so do not represent the situation for general factoring. Reducing modulo $N$ would more accurately capture the situation in general factoring, and this can mask collisions modulo $p$. (For instance, $2$ and $5$ are collisions modulo $3$, but not modulo $33$.)
For $N = 33$, $p = 3$ and $\sqrt{p} = 1.7\dots$. We expect a collision in a multiple of $2$ tries.
For $N = 247$, $p = 13$ and $\sqrt{p} = 3.6\dots$. We expect a collision in a multiple of $4$ tries.
For $N = 14649$, $p = 57$ and $\sqrt{57} = 7.5\dots$. We expect a collision in a multiple of $8$ tries.
I write "in a multiple of" because being $O(\sqrt{p})$ only means that we know the shape of growth. A constant multiple of $\sqrt{p}$ is also $O(\sqrt{p})$. The multiple could be anything, bigger than $1$, less than $1$, we don't know without vastly more understanding of the relevant number theory. Also, all of these examples are "tiny". The "big-O" notation really only describes the asymptotic (very large) behaviour of the function. For these miniscule examples, the behavioue could be very far away from the asymptotic behaviour.