Why is the remainder uniformly distributed when 1,2,3,... are divided by an irrational number? Let remainder $r$ be defined as
$$
r = n - pq
$$
where $n \in \mathbb{N}$ is the dividend , $q \in \mathbb{R}$ is the divisor, and $p = \mathrm{floor}(n/q)$.
I calculated the remainders by dividing by $\pi$ for $1,2,\dots,100000$, then I got the reminders seem to follow a uniform distribution on $[0,\pi)$.

When the divisor is rational, it's obvious that the remainder is limited to some certain values.
When it comes to a irrational divisor, it makes intuitive sense to me that the remainders can have any values on $[0,q)$, and follow a uniform distribution. However, I have no idea how to prove this.
EDIT: Added histogram.
 A: Dividing your formula by $q=\pi$ gives
$$
n \cdot \frac{1}{\pi}
= p + \frac{r}{\pi}
,
$$
so $p$ is the integer part,
and $\frac{r}{\pi}$ is the fractional part, of $n \cdot \frac{1}{\pi}$.
As $n$ runs through the integers, this fractional part $\frac{r}{\pi}$ is uniformly distributed on $[0,1)$
by Weyl's equidistribution theorem,
since $\frac{1}{\pi}$ is irrational. Hence $r$ is uniformly distributed on $[0,\pi)$.
A: As said in the answer of @HansLundmark this is Weyl's equidistribution theorem. It is an application of the ergodic theorem combined with the theorem that an irrational rotation of the circle is ergodic, and here's some details of how to reduce it to those two theorems.
Consider the unit circle $S^1$ in the complex plane. Given $q \in \mathbb{R}$, consider the function which rotates the circle $S^1$ through the angle $\frac{2 \pi}{q}$:
$$R(z) = e^{2 \pi i / q} z
$$
Define the iterates of this function by induction to be
$$R^n(z) = R(R^{n-1}(z))
$$
Your remainder function $r(n) \in [0,q)$ is the unique function such that
$$R^n(1) = e^{2 \pi i \, r(n)/q}
$$
This sequence $R^n(1)$ is the "orbit" of $1=1+0i$ under the action of the function $R(z)$.
So your question comes down to asking: Why is the subset $\{R(x), R^2(x),\ldots,R^n(x)\}$ equally distributed in the circle as $n \to +\infty$?
Now let's bring in the ergodic theory. 
Let $\mu$ denote the Borel measure on $S^1$ which assigns to an interval of angle $\alpha$ a measure of $\alpha/2\pi$. The theorem referred to above says that since $\frac{1}{q}$ is irrational, the transformation $R$ is ergodic with respect to $\mu$, and that means: for all $\mu$-measurable subsets $A \subset S^1$, if $R(A)=A$ then $\mu(A)=0$ or $1$.
Next we apply the Ergodic Theorem which, using that $R$ is ergodic with respect to $\mu$, says: For every $\mu$-integrable function $f : S^1 \to \mathbb{R}$ and for $\mu$-almost every $x \in S^1$ we have
$$\lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^n f(R^i(x)) = \int_{S^1} f \, d\mu
$$
So finally we can use this to see why an orbit is equally distributed. Take an angular interval $A \subset S^1$ of angle $\alpha$, and let $\chi_A$ be the characteristic function of $A$. The expression $\frac{1}{n} \sum_{i=1}^n \chi_A(R^i(x))$ (under the limit sign on the left hand side) is simply the proportion of the points in the set $\{R(x), R^2(x),\ldots,R^n(x)\}$ which lie in $A$. What "uniformly distributed" should mean is that limit of this proportion as $n \to \infty$ (the left hand side) should be equal to $\alpha/2\pi=\mu(A)=\int_{A} d\mu = \int_{S^1} \chi_A$ (the right hand side). And they are indeed equal, that is exactly what the Ergodic Theorem says.
A: Ok, I'm sure I am making an error but I can't spot it.  I understand (because I just looked it up) that this is typically proven by ergodic methods but it still seems elementary to me.  I expect I am making a blunder and would appreciate having someone show me where.
We wish to show that if two intervals within $[0,q]$ have the same length then the remainder is equally likely to be in either of them.
Suppose the desired interval has length $\Delta$.  Since $q$ is irrational we know that, given any positive $\Delta$ we can find an integer $m$ with $(mq)$ less than $\Delta$ (where $(x)$ denotes the remainder $x\;mod (q)$.  More generally, if $N$ is a (large) integer we can find an integer $m(N)$ such that $(m(N)q)<\frac{\Delta}{N}$.  We now partition the natural numbers by congruence classes mod (m(N)).  Take any interval $[a,b]$ of length $\Delta$. We count the number of congruence classes $m_i\;\;mod(m(N))$ such that $(m_i q)\in[a,b]$ and get $\left[\frac{a}{(m(N)q)}\right]-\left[\frac{b}{(m(N)q)}\right]$.  But this is always within $2$ of $\frac{a-b}{(m(N)q)}$ and as $N\rightarrow\infty$ we see independence of the placement of $a$ and $b$.  
