nt.number theory - What's the "best" proof of quadratic reciprocity?

The question asked for the nicest proof for a first undergraduate course. Has anyone who offered a proposal used their favorite choice in a course? (Obviously the suggestions referring to $K$-theory or Hilbert symbols weren't suggested in that spirit.)
I've taught an undergrad number theory class several times and initially I gave the Gauss sum proof. But I realized afterwards that to the students this truly comes out of nowhere (it seems too magical), so I hunted around for other proofs, preferably some which build on more basic ideas that I could present earlier in the course. The Eisenstein (sine-function) proof doesn't fit that requirement, and Zolotarev seems too far-out if the students have not had group theory (which most have not). So what else is available?

There is a proof due to V. Lebesgue (not H. Lebesgue!) that is based on counting points on hyperspheres mod $p$. It can be found in Ireland-Rosen's book. For an odd prime $p$ and positive integer $n$, let

$$N_n(p) = #{(x_1,dots,x_n) in ({mathbf Z}/(p))^n : x_1^2 + cdots + x_n^2 = 1}.$$

This is the number of mod $p$ points on the sphere in $n$-space mod $p$. Earlier in the course I have the students discover numerically that every number mod $p$ is a sum of two squares. That is,

$$#{(x,y) in ({mathbf Z}/(p))^2 : x^2 + y^2 = a}$$

is positive for every $a$ in ${mathbf Z}/(p)$. This could be shown by the pigeonhole principle, since $x^2$ and $a - y^2$ each take $(p+1)/2$ values mod $p$ and thus have an overlap. But more precisely, if you look at examples, you quickly discover that this 2-variable count is independent of $a$ when $a$ is nonzero (from a more adv. point of view, the independence is because the norm map on unit groups $(({mathbf Z}/(p))[t]/(t^2+1))^times to ({mathbf Z}/(p))^times$ is a homomorphism so its fibers have the same size, but that is a crazy explanation in an elem. number theory course). Enough data suggest what that uniform value is for any nonzero $a bmod p$, and then we prove that in class. With this 2-variable count we return to the hypersphere count and get a simple recursive formula connecting $N_n(p)$ to $N_{n-2}(p)$. If you let $n = q$ be an odd prime, the recursive formula involves $p^{(q-1)/2}$ plus $N_{q-2}(p)$ times a multiple of $q$, so $N_q(p) bmod q$ involves $(frac{p}{q})$. [Note: Although the application will use $N_q(p)$, you must think about $N_n(p)$ for general odd $n$ first since the recursion from $n$ back to $n-2$ makes no sense in general when the number of variables is only an odd prime: $n-2$ usually isn't prime when $n$ is.]

At the same time, the set being counted by $N_n(p)$ is invariant under cyclic shifts of the coordinates. On the very first problem set I have the students discover numerically that the number of cyclic shifts of an $n$-tuple is always a divisor of $n$. So when we let $n = q$ be an odd prime, $N_n(p) = N_q(p)$ is the number of constant $q$-tuples on the unit sphere mod $p$ plus a multiple of $q$. A constant $q$-tuple is basically counting whether or not $q bmod p$ is a square. So $N_q(p) bmod q$ is related to $(frac{q}{p})$.

In the two approaches to counting $N_q(p) bmod q$, one involves $(frac{p}{q})$ and the other involves $(frac{q}{p})$. This implies the QR relation mod $q$, which is actual equality since $1$ is not $-1 bmod q$.

One nice thing about this approach is that it can also be used to prove the supplementary law for $(frac{2}{p})$, by counting

$$#{(x,y) in ({mathbf Z}/(p))^2 : x^2 + y^2 equiv 1 bmod p}$$

in two ways. First there is the exact formula for the count (not just mod $p$) which I mentioned before. Second, most solutions in this count come in packets of size 8 (permute coordinates and change signs to get 8 solutions out of one solution). The exceptions which don't fall into packets of size 8 (when $x$ is $pm y$ mod $p$ or $x$ or $y$ is $0$ mod $p$) depend on whether or not 2 mod $p$ is a square (does $2x^2 equiv 1 bmod p$ have a solution?), and comparing these two formulas mod 8 implies the usual rule for $(frac{2}{p})$.

Since I am able to get the students to work on ideas that are used in the proof much earlier on during the semester (one doesn't need quadratic residues to numerically look at solutions of $x^2 + y^2 equiv a bmod p$, for instance), this proof nicely ties together things they have seen throughout the course. So that's why this is my vote for the best proof to give in an undergrad course.