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 $(p1)/(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 $(p1)/(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 wellknown) 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{p1}{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{p1}{2p}$ as Noam Elkies said. 


.5 Let me atone for giving too few details by giving too many. Let $$S=\sum_{a_1=0}^{p1}\dots\sum_{a_n=0}^{p1}\sum_{t=0}^{p1}\sum_{m=0}^{p1}e^{2\pi im(a_1^2+\cdots+a_n^2t^2)/p} $$ The innermost sum is $p$ if $a_1^2+\cdots+a_n^2t^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}^{p1}\sum_{a_1=0}^{p1}\dots\sum_{a_n=0}^{p1}\sum_{t=0}^{p1}e^{2\pi im(a_1^2+\cdots+a_n^2t^2)/p} $$ so $$ S=p^{n+1}+\sum_{m=1}^{p1}\sum_{a_1=0}^{p1}\dots\sum_{a_n=0}^{p1}\sum_{t=0}^{p1}e^{2\pi im(a_1^2+\cdots+a_n^2t^2)/p} $$ so $$ S=p^{n+1}+\sum_{m=1}^{p1}\left(\left(\sum_{a_1=0}^{p1}e^{2\pi ima_1^2/p}\right)\cdots\left(\sum_{a_n=0}^{p1}e^{2\pi ima_n^2/p}\right)\left(\sum_{t=0}^{p1}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 $Sp^{n+1}\le(p1)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^{n1}+O(q^{n/2})$. For $n>2$ we can write $Q(X)=Q_0(X_1,\ldots,X_{n2}) +X_{n1}X_n$, where $Q_0$ is a form of rank $n2$ in the variables $X_1,\ldots,X_{n2}$. This decomposition shows instantly that that $A(n,d)=(2 q1) A(d,n2) +(q1) (q^{n2}A(d,n2))$. Proceeding by induction we get the estimate $A(n,d)=q^{n1}+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^2X_{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)= (q1)/(2q) +O(q^{n/2})$. 

