The nice thing about Blum-Blum-Shub or Blum-Micali is that they come with a proof of security. Even then, they tend to be impractical, due to performance and side channels.
This one is missing the most important part, the proof. Indeed, a sibling comment notes that empirical results look pretty flawed.
First, I don't think the error term is contributing much to the solution. It almost never affects the high bit. In addition, it isn't fed back into updating the secret vectors, so I think an analysis can pretend it doesn't exist.
The nonlinear step where each value is squared looks questionable to me. You will only produce quadratic residues (https://en.wikipedia.org/wiki/Quadratic_residue) when you square the numbers. This roughly halves the number of possibilities.
So what this really boils down to is this:
You have a matrix G and a vector s and a prime p. You repeatedly compute s' = Gs (Hadamard) Gs mod p. Each time you run this step you are projecting into a space with dimensionality (p/2)^N from a space p^N. My guess is that most operations will get trapped into short cycles.
Using you example values, after 10 iterations it gets to [9, 16, 13, 8]. This then repeats with a cycle length of 20. Given 4 values with p = 17 you could get up to 83520 values before repeating.
In some random tests, 6 values with p=97 enters a cycle at iteration 3802 even though there were 832,972,004,929 values.
6 values with p=271 enters a cycle at iteration 166,684 even though there were 396,109,944,105,121 values.
Comments
The nice thing about Blum-Blum-Shub or Blum-Micali is that they come with a proof of security. Even then, they tend to be impractical, due to performance and side channels.
This one is missing the most important part, the proof. Indeed, a sibling comment notes that empirical results look pretty flawed.
TLDR - this RNG is completely and totally broken.
First, I don't think the error term is contributing much to the solution. It almost never affects the high bit. In addition, it isn't fed back into updating the secret vectors, so I think an analysis can pretend it doesn't exist.
The nonlinear step where each value is squared looks questionable to me. You will only produce quadratic residues (https://en.wikipedia.org/wiki/Quadratic_residue) when you square the numbers. This roughly halves the number of possibilities.
So what this really boils down to is this:
You have a matrix G and a vector s and a prime p. You repeatedly compute s' = Gs (Hadamard) Gs mod p. Each time you run this step you are projecting into a space with dimensionality (p/2)^N from a space p^N. My guess is that most operations will get trapped into short cycles.
Using you example values, after 10 iterations it gets to [9, 16, 13, 8]. This then repeats with a cycle length of 20. Given 4 values with p = 17 you could get up to 83520 values before repeating.
In some random tests, 6 values with p=97 enters a cycle at iteration 3802 even though there were 832,972,004,929 values.
6 values with p=271 enters a cycle at iteration 166,684 even though there were 396,109,944,105,121 values.