Computing the endomorphism ring of an elliptic curve over a finite field (using SAGE) $ 
\newcommand{\End}{\mathrm{End}}
\newcommand{\Gal}{\mathrm{Gal}}
\newcommand{\kb}{\overline{k}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\F}{\mathbb{F}}
\newcommand{\Q}{\mathbb{Q}}
$
I would like to have an algorithm (possibly very inefficient) that computes the endomorphism ring of a given elliptic curve $E$ over a finite field $k$.
For simplicity, we shall assume that $E$ is ordinary (to avoid maximal orders in quaternion algebras...), so it is enough to compute the conductor of $\End(E)$ in the imaginary quadratic field $K := \Q(\pi)$, where $q = |k|$ and $\pi = \sqrt{a_q^2 - 4q}$. We know that this conductor divides $[O_K : \Z[\pi]]$, the latter being quite easy to compute in SAGE I suppose.
But now, is there a way to check whether, for a given $f \mid [O_K : \Z[\pi]]$, we have $\Z + f O_K = \End(E)$ ? This is where I don't know how to proceed.
I am aware of Kohel's thesis, which involves isogeny graphs, but I'm not sure if one can implement this on SAGE easily.
Ideally, I want to reproduce the table on page 303 here, which lists $\End_{\Bbb F_7}(E)$ for all (isomorphism classes of) elliptic curves over $\Bbb F_7$.
 A: I couldn't resist trying to code up the algorithm hinted at in my comment.
Sage code at:
https://pastebin.com/3ANkqfZp
This is v2 of the code; version 1 used only prime-degree endomorphisms, and thus failed to distinguish between the orders $\mathbf{Z}[\sqrt{-15}]$ and $\mathbf{Z}[\tfrac{1 + \sqrt{-15}}{2}]$ (since the sets of primes arising as norms of elements in these two rings are the same). So it failed on $y^2 = x^3 + x + 8$ over $\mathbf{F}_{19}$. The new implementation uses endomorphisms of all degrees (piecing them together from isogenies of prime degree).
For instance, here we recover the values of $c$ in the table from p.303 in your linked article, using the functions defined in the script linked above:
sage: listofcurves = []
sage: K = GF(7)
sage: for a in K:
....:     for b in K:
....:         if 4*a^3 + 27*b^2 == 0:
....:             continue
....:         E = EllipticCurve([0, 0, 0, a, b])
....:         if not any(E.is_isomorphic(F) for F in listofcurves):
....:             print(E.ainvs(), endomorphism_conductor(E))
....:             listofcurves.append(E)
(0, 0, 0, 0, 1) 1
(0, 0, 0, 0, 2) 1
(0, 0, 0, 0, 3) 1
(0, 0, 0, 0, 4) 1
(0, 0, 0, 0, 5) 1
(0, 0, 0, 0, 6) 1
(0, 0, 0, 1, 0) 2
(0, 0, 0, 1, 1) 1
(0, 0, 0, 1, 3) 1
(0, 0, 0, 1, 4) 1
(0, 0, 0, 1, 6) 1
(0, 0, 0, 3, 0) 1
(0, 0, 0, 3, 1) 2
(0, 0, 0, 3, 2) 3
(0, 0, 0, 3, 3) 1
(0, 0, 0, 3, 4) 1
(0, 0, 0, 3, 5) 3
(0, 0, 0, 3, 6) 2

PS: Important caveat: this computes the endomorphisms of $E$ over the given field $k$ (not over $\overline{k}$). If $E$ is ordinary this makes no difference (all endomorphisms are defined over $k$). In the supersingular case, it makes a difference: if $E$ is supersingular, but $\operatorname{Frob}_E$ isn't in $\mathbf{Z}$, we don't get a quaternion order, but an imaginary quadratic order (the centraliser of the Frobenius in $\operatorname{End}_{\overline{k}}(E)$). If $E$ is supersingular and the Frobenius is an integer, which can only happen if $k$ has even degree over its prime field, then the algorithm fails.
