What is the probability that the sum of squares of n randomly chosen numbers from $Z_p$ is a quadratic residue mod p?
That is, let $a_1$,..$a_n$ be chosen at random. Then how often is $\Sigma_i a^2_i$ a quadratic residue?
What is the probability that the sum of squares of n randomly chosen numbers from $Z_p$ is a quadratic residue mod p?
That is, let $a_1$,..$a_n$ be chosen at random. Then how often is $\Sigma_i a^2_i$ a quadratic residue?
This probability can be calculated exactly, and indeed it approaches $1/2$ rather quickly — more precisely, for each $p$ it approaches the fraction $(p-1)/(2p)$ of quadratic residues $\bmod p$. This can be proved by elementary means, but perhaps the nicest way to think about it is that if you choose $n$ numbers $a_i$ independently and sum $a_i^2 \bmod p$, the resulting distribution is the $n$-th convolution power of the distribution of a random single square — so its discrete Fourier transform is the $n$-th power of the D.F.T., call it $\gamma$, of the distribution of $a^2 \bmod p$. For this purpose $\gamma$ is normalized so $\gamma(0)=1$. Then for $k \neq 0$ we have $\gamma(k) = (k/p) \gamma(1)$ [where $(\cdot/p)$ is the Legendre symbol], and $$ p \gamma(1) = \sum_{a \bmod p} \exp(2\pi i a^2/p), $$ which is a Gauss sum and is thus a square root of $\pm p$. It follows that $|\gamma(k)| = p^{-1/2}$, from which we soon see that each value of the convolution approaches $1/p$ at the exponential rate $p^{-n/2}$, and the probability you asked for approaches $(p-1)/(2p)$ at the same rate.
As noted above, this result, and indeed the exact probability, can be obtained by elementary means, yielding a (known but not well-known) alternative proof of Quadratic Reciprocity(!). But that's probably too far afield for the present purpose.
The probability depends on the parity of $n$ and the residue of $p$ modulo $4$: it can be calculated in a straightforward way using Gauss sums.
Let $n$ be $2k$ or $2k+1$, and let $p\equiv r\pmod{4}$ where $r=\pm 1$. Then, assuming I made no mistake, the probability equals $$ \frac{p+1}{2p}+\frac{p-1}{2p}(rp)^{-k}. $$
Note that in my calculation I regarded zero as a quadratic residue. If we exclude zero then the final answer will look slightly different, with a main term $\frac{p-1}{2p}$ as Noam Elkies said.
.5
Let me atone for giving too few details by giving too many.
Let $$S=\sum_{a_1=0}^{p-1}\dots\sum_{a_n=0}^{p-1}\sum_{t=0}^{p-1}\sum_{m=0}^{p-1}e^{2\pi im(a_1^2+\cdots+a_n^2-t^2)/p} $$ The innermost sum is $p$ if $a_1^2+\cdots+a_n^2-t^2\equiv0\pmod p$ and zero otherwise, so $S$ counts $2p$ whenever $a_1^2+\cdots+a_n^2$ is a (nonzero) quadratic residue, $p$ whenever it's zero. On the other hand, $$ S=\sum_{m=0}^{p-1}\sum_{a_1=0}^{p-1}\dots\sum_{a_n=0}^{p-1}\sum_{t=0}^{p-1}e^{2\pi im(a_1^2+\cdots+a_n^2-t^2)/p} $$ so $$ S=p^{n+1}+\sum_{m=1}^{p-1}\sum_{a_1=0}^{p-1}\dots\sum_{a_n=0}^{p-1}\sum_{t=0}^{p-1}e^{2\pi im(a_1^2+\cdots+a_n^2-t^2)/p} $$ so $$ S=p^{n+1}+\sum_{m=1}^{p-1}\left(\left(\sum_{a_1=0}^{p-1}e^{2\pi ima_1^2/p}\right)\cdots\left(\sum_{a_n=0}^{p-1}e^{2\pi ima_n^2/p}\right)\left(\sum_{t=0}^{p-1}e^{2\pi imt^2/p}\right)\right) $$ Each of those inner sums is a Gauss sum and known to equal $\sqrt p$ in modulus (more detail: the sum is ${m\overwithdelims()p}\sqrt{{-1\overwithdelims()p}p}$), so $|S-p^{n+1}|\le(p-1)p^{(n+1)/2}$. For $n\gt1$, the main term beats the error term, and you get a good estimate.
Here is a slightly different argument: Let $Q$ be a non degenerate quadratic form over $\mathbb{F}_q$ of rang $n$ and determinant $d$. Let $A(n,d)=|\{x\in \mathbb{F}_q^n:Q(x)=0\}|$. The claim is that $A(n,d)=q^{n-1}+O(q^{n/2})$. For $n>2$ we can write $Q(X)=Q_0(X_1,\ldots,X_{n-2}) +X_{n-1}X_n$, where $Q_0$ is a form of rank $n-2$ in the variables $X_1,\ldots,X_{n-2}$. This decomposition shows instantly that that $A(n,d)=(2 q-1) A(-d,n-2) +(q-1) (q^{n-2}-A(-d,n-2))$. Proceeding by induction we get the estimate $A(n,d)=q^{n-1}+O(q^{n/2})$. (The error term can be computed exactly using Gauss sums). Applying this to the forms $X_1^2+\cdots +X_n^2-X_{n+1}^2$ and $X_1^2+\cdots+ X_n^2$ we get that the desired probability is $(A(n+1,-1)-A(n,1))/(2 q^n)= (q-1)/(2q) +O(q^{-n/2})$.