What is the distribution of primes modulo $n$? Let $n\geq 2$ and let $k$ be "considerably larger" than $n$ (like some large multiple of $n$). Then for each $i$ such that $0<i<n$ and $\gcd(i,n)=1$ let's define
$$c_i=\left|\{p_j\;|\; p_j\equiv i \mod n,\;\mbox{where $p_j$ is the $j$-th prime, $1\leq j\leq k$}\}\right|$$
so $c_i$ represents how many of the first $k$ primes are congruent to $i$ modulo $n$.
What can we say about $c_i$s, that is about the distribution of the first $k$ primes modulo $n$?
I thought the distribution will be seemingly random, and that is mostly true - for example $c_i$s are always very close together. But there are observable non-random patterns. For example for $n=3$, for various $k$s I've tried  (up to $10^6$) I always got $c_2>c_1$. If this were just some kind of random discrepancy due to distribution of small primes, it would eventually vanish for large $k$s, which I don't observe.
 A: Dirichlet's theorem on primes in arithmetic progression tells us that the proportion of primes will be the same, for values of $i$ that are coprime to $n$ (and 0 otherwise).
However, Chebyshev's bias tells us that numerically there are more primes with a quadratic non-residue, than a residue, when you're counting primes up to $N$.
A: I once tried to answer the particular case of primes $4n\pm 1$ to myself. 
(I think it will strongly help to read the answers by Raymond, Raymond  and Greg. At the end of the last answer there's also a link to the chat, where we continued the discussion.)
Here is how far I got with an explicit formula for the number of primes of the form $4n+3$ below $x$, $\pi^*(x;4,3)$, expressed in terms of (sums of) sums of Riemann's $R$ functions over roots of Riemann's $\zeta$ resp. Dirichlet $\beta$ function:
\begin{align*}
\Pi^*(x;4,3)
&= \pi^*(x;4,3) + \tfrac12 \sum_{\substack{b\pmod 4 \\ b^2\equiv
3\pmod 4}} \pi^*(x^{1/2};4,b) + \tfrac13 \sum_{\substack{c\pmod q \\
c^3\equiv 3\pmod 4}} \pi^*(x^{1/3};4,c) + \cdots \\
\end{align*}
Then I try to complete things by adding several up
\begin{align*}
\Pi^*(x;4,3) &= \tfrac11\pi^*(x;4,3) + \tfrac13 \pi^*(x^{1/3};4,3) + \cdots \\
\tfrac12\Pi^*(x^{1/2};4,3) &= \tfrac12\pi^*(x^{1/2};4,3) + \tfrac16
\pi^*(x^{1/6};4,3) + \cdots \\
\tfrac14\Pi^*(x^{1/4};4,3) &= \tfrac14\pi^*(x^{1/4};4,3) +
\tfrac1{12} \pi^*(x^{1/12};4,3) + \cdots \\
&\vdots&\\
\hline\\
\tag{1}\sum_{k=0}^\infty
2^{-k}\Pi^*(x^{2^{-k}};4,3)&=\sum_{m=1}^\infty \tfrac1m
\pi^*(x^{1/m};4,3)
\end{align*}
Using Möbuis inversion I'll get
\begin{align*}
\pi^*(x;4,3)&=\sum_{m=1}^\infty \tfrac{\mu(m)}m\sum_{k=0}^\infty
2^{-k}\Pi^*(x^{2^{-k}/m};4,3)\\
\tag{2}&=\sum_{k=0}^\infty 2^{-k}\sum_{m=0}^\infty
\tfrac{\mu(m)}m\Pi^*(x^{2^{-k}/m};4,3)
\end{align*}
Now I use
\begin{align*}
\Pi^*(x^{2^{-k}};4,3)&=\frac1{\phi(4)} \sum_{\chi\pmod 4}
\overline{\chi(3)}\Pi^*(x^{2^{-k}},\chi)\\
\tag{3}&=\frac12 \left( \Pi^*(x^{2^{-k}},\chi_1)-
\Pi^*(x^{2^{-k}},\chi_2) \right)
\end{align*}
and then
\begin{align*}
\tag{$4_1$}\Pi^*(x^{2^{-k}},\chi_k)&=\operatorname{li}(x^{1/2^{k}})-\sum_{\rho_\zeta}
\operatorname{li}(x^{\rho_\zeta/2^k})\text{ if $k=1$}\\
\tag{$4_2$}&=\phantom{\operatorname{li}(x^{1/2^{k}})}-\sum_{\rho_\beta}
\operatorname{li}(x^{\rho_\beta/2^k})\text{ if $k=2$}\\
\end{align*}
which gives
\begin{align*}
\tag{3'}\Pi^*(x^{2^{-k}};4,3)&=\frac12 \left(
\operatorname{li}(x^{1/2^{k}})-\sum_{\rho_\zeta}
\operatorname{li}(x^{\rho_\zeta/2^k}) +\sum_{\rho_\beta}
\operatorname{li}(x^{\rho_\beta/2^k}) \right)
\end{align*}
so finally
\begin{align*}
\pi^*(x;4,3)&=\sum_{k=0}^\infty 2^{-k}\sum_{m=0}^\infty
\tfrac{\mu(m)}m\frac12 \left(
\operatorname{li}(x^{1/2^{k}})-\sum_{\rho_\zeta}
\operatorname{li}(x^{\rho_\zeta/2^k}) +\sum_{\rho_\beta}
\operatorname{li}(x^{\rho_\beta/2^k}) \right)\\
\tag{5}&=\sum_{k=0}^\infty 2^{-k-1}\left(
\operatorname{R}(x^{1/2^{k}})-\sum_{\rho_\zeta}
\operatorname{R}(x^{\rho_\zeta/2^k}) +\sum_{\rho_\beta}
\operatorname{R}(x^{\rho_\beta/2^k}) \right)
\end{align*} 
I would be very, very glad to read your opinion...
