Order of Edwards curve and its twist In Mike Hamburg's Ed448-Goldilocks, a new elliptic curve (eprint 2015, WECCS 2015) it is studied untwisted Edwards curves in the prime field $\mathbb F_p$
$$E_d:\,y^2+x^2\,=\,1+d\,x^2\,y^2$$
with large prime $p\equiv3\pmod 4$ and the Legendre symbol $\displaystyle\left(\frac d p\right)=-1$.
The matching &quot;twist&quot; is
$$E'_d:\,y^2-x^2\,=\,1-d\,x^2\,y^2$$
Constant $d$ is chosen with minimal $|d|$ such that the curve's order $|E_d|$ is $4\cdot q$ with $q$ prime, the twist's order is $|E'_d|=4\cdot r$ with $r$ prime, and $q&lt;p/4$.
The paper uses prime $p=4^{224}-2^{224}-1$ and gives $d=-39081$,
$q=2^{446}-\mathtt{8335dc163bb124b65129c96fde933d8d723a70aadc873d6d54a7bb0d_h}$,
$r=2^{446}+\mathtt{0335dc163bb124b65129c96fde933d8d723a70aadc873d6d54a7bb0d_h}$
It holds $|E_d|+|E'_d|=2\cdot p+2$. Update: initially my experiments¹ differed by two but Mike Hamburg kindly pointed my mistake: I did not count the two points at infinity for $E'_d$.
The question (now) boils down to: Why $|E_d|+|E'_d|=2\cdot p+2$ ? And how do we find $d$ given $p$ ?
If the later is by mere enumeration of $d\gets-j\cdot\displaystyle\left(\frac j p\right)$ for incremental $j&gt;0$ and checking $q\gets|E_d|/4$ and $r\gets|E'_d|/4$ are prime, how are these computed²?

¹ With $p\gets4^i-2^i-1$ for $i\in\{4,5\}$, I get
$$\begin{array}{r|rrr}
i&amp;p&amp;d&amp;q&amp;r\\
\hline
4&amp;239&amp;19&amp;59&amp;61\\
5&amp;991&amp;-45&amp;233&amp;263
\end{array}$$
² This may be asking for Schoof–Elkies–Atkin adapted to Edwards curves. Pointer to an implementation also welcome.
Regarding the [B] and [C] parts of the question per the comments:
I'm not sure how exactly did Mike Hamburg find the curve, but from what I know it's usually easier to find the order of the matching Montgomery curve. Recall that Montgomery curves have the form $By^2 = x^3 + Ax^2 + x$. If $B$ is 1, then it fits into the generalized Weierstrass form, and most SEA algorithm implementations work with any curve in the generalized Weierstrass form. (If it's not 1 then you can easily map into a curve with $B = 1$, the same way that short Weierstrass curves can be mapped into $a = -3$)
So basically:

Search for a Montgomery curve matching the criteria;
Then convert it into Edwards form.

One optimization is to instruct SEA to quickly discard curves whose order it knows beforehand that have a small factor (other than 4 or 8), see the tors parameter of the ellsea PARI/GP function, for example.
The paper &quot;A note on high-security general-purpose elliptic curves&quot; has a Magma implementation of the process (though IIRC it uses a slightly different approach). RFC 7748 has a Sage script that also searches for a Montgomery curve (though it will probably be much slower, since it doesn't seem to support that optimization).