Да, давайте я прокомментирую эти картинки. Это мы рисовали гауссовы суммы

\sum_{n=0}^{p-1} exp( 2πi n^2/p )

И они оказывались равны при p=4k+1 — корню из p, а при p=4k+3 — тому же корню из p, но умноженному на i.
А почему картинки частичных сумм при больших p так выглядят? А вот почему: когда p очень большое, то пока n маленькое, n^2/p при изменении n на 1 меняется мало — всего лишь на ~2n/p. Так что кривая идёт «в одну сторону», потихоньку начиная заворачиваться. Чем быстрее n — тем быстрее, что мы, собственно, и видим. Ну и если сделать замену x = n / \sqrt{p}, то в масштабе «n порядка корня из p» получится практически интегральная сумма Римана для интеграла от exp(2πi x^2), только умноженная на разницу между соседними x — как раз на \sqrt{p}. То есть — практически интеграл от гауссовой плотности, только с мнимой, а не положительной, дисперсией.

Когда n уходит за пределы этого масштаба, сумма начинает дёргаться во все стороны, в итоге стоя на месте. И отсюда получается часть суммы вида

\sqrt{p} * (1+i)/2.

Но. В некоторый момент сдвиги опять начинают идти в одну сторону. И происходит это, что логично, при n около n_0=(p+-1)/2. Потому что там угол между соседними сдвигами, примерно 2π*2n/p, как раз почти обнуляется.

И там будет примерно такая же сумма — только умноженная на
exp(2πi n_0^2/p).

Если p=4k+1, то n_0=2k, соответственно,
4k^2/p = 4k^2/(4k+1) = k - k/(4k+1),
и
exp(2πi n_0^2/p) ~ exp ( - π/2) = -i.

Так что общая сумма примерно должна быть равна

\sqrt{p} * (1+i)/2 * (1-i) = \sqrt{p}.

А вот если p=4k-1, n_0=2k, то

4k^2/p = 4k^2/(4k-1) = k + k/(4k+1),
и
exp(2πi n_0^2/p) ~ exp ( + π/2) = +i.

То есть сдвиг между двумя «натоптанными кругами», где гуляют частичные суммы, для p=4k-1 будет в противоположную сторону — и общая сумма тогда примерно будет равна

\sqrt{p} * (1+i)/2 * (1+i) = i* \sqrt{p}.

И вот множитель i и вылез. А вот то, что равенства точные, а не приближённые, это так увидеть нельзя. Но зато получается объяснить, что же мы видим на картинках.