Once I needed to calculate a certain quadratic form on random vectors: $$ I = \sum_{i,j=1}^N x_i x_j \;, $$ where $x_i = \pm 1$ are random variables, the probabilities of their values are $P[x_i=-1] = P[x_i=+1] = 1/2$. All $x_i$ are statistically independent. Hence $I$ is the integer-valued random variable.
I needed to know a probability of event $I = N$. So for the start I just sampled the vectors $(x_1, \dots x_N)$ and found that the values of the $I$ give a sequence
0, 4, 16, 36, 64, 100, 144, 196, 256, 324, 400, 484, 576, 676, 784, 900, 1024, 1156, 1296, 1600, ...
(the higher $N$ the longer the sequence, of course). As seen, the condition $I = N$ may or may not be satisfied.
Does anyone have an idea about this sequence?
Here, the Python code:
import numpy as np
S = 10000
N = 10
xx = 2 * np.random.randint(2, size=(S, N)) - 1
I = np.array([ np.einsum('ij,i,j', np.ones((N, N)), xx[s], xx[s])
for s in range(S) ])
print(np.unique(I))
Looks like $(2n)^2$. (additional characters)