Are there natural choices of $\sqrt{-1}$ in $\mathbb Z/p\mathbb Z$ for a prime $p\equiv 1\pmod 4$ The definition of "natural" is of course somewhat personal and involves perhaps aesthetic aspects and/or interesting properties. I illustrate it with two examples:
The integer $((p-1)/2)!$ is such a square-root of $-1$ (modulo $p$ for $p$ a prime congruent to $1$ modulo $4$). 
Another choice which is perhaps natural is given by $2a/b\pmod p$ where $p=4a^2+b^2$ with $a,b\in\mathbb N$.
Are there other nice formulae?
An example of a not very natural choice (in my eyes) is $a^{(p-1)/4}\pmod p$ where $a$ is the smallest natural integer generating the multiplicative 
group of invertible elements modulo $p$. However, I accept this answer gladly, if there is some generator $a$ of $(\mathbb Z/p\mathbb Z)^*$ given by a "natural" formula. 
(The motivation for this question was Does (the ideal class of) the different of a number field have a canonical square root? )
 A: Let me mention that Gauss proved that if $p=4n+1$ is a prime, then $p=a^2+b^2$ where $a$, $b$ are the least absolute value remainders of $\binom{2n-1}{n-1}$ and $(2n)!\binom{2n-1}{n-1}$  respectively mod $p$. Although the formulas are simple, the proof is not trivial. 
A: As you probably know, you are also asking for a particular Gaussian prime ideal $(\pi )$ with norm $p$. Which is to say that such a choice of square root of $-1$ is the same as giving a homomorphism from the Gaussian integers to the field with $p$ elements, which is the same as giving its kernel. This accounts for the formulae in $a$ and $b$. The theory of complex multiplication gives a little bit more, namely by identifying the Frobenius automorphism which is subject to some congruence conditions depending on choice of curve with complex multiplication by the Gaussian integers. The formulae with binomial coefficients come out of formulae for the Hasse invariant. 
Edit: To clarify somewhat, given an elliptic curve $C$ over the rational numbers with complex multiplication by the Gaussian integers, it makes good sense to consider its reduction $C$ mod $p$ (away from primes of bad reduction). But it does not make naive sense to ask for the endomorphism ring to reduce mod $p$. In fact this looks like the same issue: the complex square root $i$ of $-1$ acts on C, and if we have a candidate square root $i$ of $-1$ mod $p$, we can then see how to reduce any endomorphism. The natural endomorphism mod $p$ is the Frobenius endomorphism. What I was trying to say amounts to two or three things, about looking at this question:


*

*there is a family of such curves C, not just one;

*making the change of field to C over $Q(i)$ makes the choice required to be one of the two prime ideals above $p$;

*if you lift the Frobenius back to the complex endomorphisms, you get the ways of expressing $i$ mod $p$ via the $a$ and $b$.

A: Even $((p-1)/2)!\pmod p$ is a "computationally difficult" square root of $-1\pmod p$, this seems to be the simplest uniform formula. Less universal is $a^{(p-1)/4}\pmod p$ for a quadratic nonresidue $a$ modulo $p$, which is easy $2^{(p-1)/4}\pmod p$ for $p\equiv 5\pmod 8$ as Victor mentions. For $p\equiv1\pmod8$, you simply search for the least quadratic nonresidue $a$ modulo $p$; as far as I remember, there is a bound of the form $a<\operatorname{const}\cdot \log^2p$ assuming GRH. Computing the Legendre--Jacobi symbol $(a/p)$ has polynomial complexity.
I cannot be original in what I indicated above, so let me add some "complicated" (but curious) stuff on a discrete analogue of the formula
$$
(1-x)^{-1/2}=\sum_{n=0}^\infty \binom{2n}{n}\biggl(\frac x4\biggr)^n.
$$
This is the subject of Lemma 3 in arXiv:1004.4337.
Denote $q(x)=q_p(x)=(x^{p-1}-1)/p$ the Fermat quotient of $x\in\mathbb Z_p^*$.
For a prime $p>3$, let $x$ be a rational number such that both
$x$ and $1-x$ are not divisible by $p$ and $1-x$ is a quadratic residue modulo $p$. Take $y$ such
that $y^2\equiv 1-x\pmod p$. Then
$$
\sum_{n=0}^{p-1}\binom{2n}n\biggl(\frac x4\biggr)^n
\equiv 1+\frac{px}{2(1-x)}\bigl(-q(x)+q(y+1)(y+1)-q(y-1)(y-1)\bigr)\pmod{p^2}.
$$
Taking $p\equiv1\pmod 4$ and $x=2$ one obtains a closed form for the left-hand side by means of $y$ which is your square root of $-1$ modulo $p$. Isn't it a curious result?
A: Let me address a more general question. Let $q$ and $p$ be primes. We can ask whether or not an element $x \in \mathbb{F}_p$ has a $q^{th}$ root, and if so, what are such roots (I'll assume throughout that $x \ne 0$). If $p \not \equiv 1$ mod $q$ then EVERY $x$ has a unique $q^{th}$ root given by $x^{a}$ with $a$ given by $a=\frac{{1 - (p - 1)^{q - 1} }}{q}$. If $p \equiv 1$ mod $q$, then $x$ is a $q^{th}$ root if and only if $x^{\frac{{p - 1}}{q}}=1$ in $\mathbb{F}_p$ and if it has a $q^{th}$ root, it has $q$ of them. If in addition, $p \not \equiv 1$ mod $q^2$, one such root can be given as $x^{b}$ where $b = \frac{{1 - \left( {\frac{{p - 1}}{q}} \right)^{q(q-1)}}}{q}$ and the others are all $q^{th}$ roots of unity times this (I admit there is some choices made even to write this formula down). This unfortunately rules out exactly your case, which is the most complicated one, namely when $p \equiv 1$ mod $q^2$ (for you $q=2$). To my knowledge, in this case, the only way to write down a $q^{th}$ root involves PICKING an arbitrary $y \in \mathbb{F}_p$ which is not a $q^{th}$ power and expressing the answer in terms of this (and the formula becomes a bit messy). So, your question seems to boil down to whether or not there is a "natural candidate" for a non-square in $\mathbb{F}_p$. Of course, you could always take the "smallest" such non-square- but there are other natural choices (e.g. the "largest").
In case you are interested in how to use such a non-$q^{th}$ power $y$ to find the $q^{th}$ root of $x$, here it goes-
WARNING, I'm fairly sure there is an easier way. This is just what I figured out for fun 6 years ago, when I was just learning abstract algebra, and is copied out of my notebook. I'm pretty sure it can be simplified, but, this at least works:
Let $n$ be the largest $n$ such that $q^{n}$ divides $p-1$. Let $$T(X) = x^{ - \left( {\frac{{p - 1}}{{q^n }}} \right)^{q^n \left( {q - 1} \right)} }.$$
The order of this is $q^{m}$ for some $m$ less than $n$. Let $z = y^{\frac{{p - 1}}{{q^{m + 1} }}}$. Let $b$ be the smallest positive integer such that $z^{bq}=T(x)$. Let $l = \frac{{1 - \left( {\frac{{p - 1}}{{q^n }}} \right)^{q^n \left( {q - 1} \right)} }}{q}$.
Then the $q^{th}$ roots of $x$ are given by the formula
$$z^{b(p - 2) + kq^m } x^l $$
with $k$ varying from $0$ to $q-1$.
To see how this works observe that $T:\mathbb{F}_p^{*} \to \mathbb{F}_p^{*}$ is a group hom onto the $q^{n}$th roots of unity, so $T(x)=:u$ is a $q^{n}$th root of unity and $x^{l}$ is the $q^{th}$ root of $x \cdot u$.
NOTE: If $x$ is a $q^{n}$th power, then $T(x)=1$ and the formula $x^{l}$ provides a $q^{th}$ root WITHOUT making a choice for $y$. HOWEVER, this never happens when $q=2$ and $x=-1$ because $\frac{{p - 1}}{{2^n }}$ is always odd hence $-1$ to that power is never $1$ so $-1$ is never a $2^{n}$th power. In fact, $T(-1)$ will always be $-1$ to an odd number, hence, $-1$, so, this becomes useless as $x \cdot u= 1$.
