I like this and will certainly review it positively soon. I have a few questions first:

1. A few more comments might be useful. What you do for many primes l is find a point P of order l defined over and extension (if necessary), compute phi(P) and phi^2(P) and use discrete logs to find a linear relation between them, which tells you the trace mod l.
2. In the two lines where you use objgen() you don't need the ring so could just use the function polygen().
3. I like the CM example! But if the endo in the example was really 1+2*i then its trace would be 2 not 4, so perhaps the endo is actually 2+i. I did not do the calculation.
4. In the code block which finds P of order l you use a splitting field except over number fields where you just take an extension. Why not always take a simple extension?
5. It might be useful to have a function which takes a curve E and a prime l and returns a point P of order l, defined over the base field if such a P exists and otherwise over a suitable extension, found as in your code. That function could go in somewhere like ell_generic or ell_field, and you could call it here. I'm sure there will be aother times when this would be useful.

Thanks for your work on this!

Some more comments, possibly more interesting:

1. Over fields of characteristic zero the trace is just phi.scaling_factor().trace(), I suggest using that instead. In characteristic p this does give the trace mod p, which gives at least partial information and even full information if some inequality holds (so that the Hasse interval has length < p).
2. I think that in the loop over l you should skip l which divide the degree. Otherwise the point P which is l-torsion could be in ker(phi), then Q=0 and the discrete log will fail.
3. Also in char. 0 for curves without CM phi must be [n] where d=n^2 and t=2n, but that is not good enough as it does not give you the sign of n. Given point 1 this can be ignored.

Correction to previous point: if s=phi.scaling_factor() and d=phi.degree() then trace=s+d/s. That's in char.0 -- as before this gives the trace mod p in char p, provided d mod p is not 0.

My comments above were by way of review

Thank you for the very useful feedback, I think I've addressed it all in the current branch. I also have a question: Do you happen to know if the optimization with the scaling factor has appeared in the literature? I (re)discovered this trick a little while ago and assumed it was new, but evidently this was either known or obvious to you.

Do you mean the s+d/s formula? No, I don't know a reference, it just seemed fairly obvious to me. Perhaps I thought of it while writing notes for a course a few years ago.

Does this need a new review? If so I will.

It's the same as before, only rebased onto the latest beta. I think you hadn't yet reviewed the changes I made following your first batch of comments, though.

OK thanks, I thought that maybe I had something yet to do. I'll do it soon, hopefully this week while at LuCaNT (pity you are not here. @yyyyx4 !)

On Linux 32-bit:

**********************************************************************
File "src/sage/schemes/elliptic_curves/ell_field.py", line 2062, in sage.schemes.elliptic_curves.ell_field.point_of_order
Failed example:
    P = point_of_order(E, 3); P
Expected:
    (x : -Y : 1)
Got:
    (x : Y : 1)
**********************************************************************

That's funny. I thought we had dealt with this issue of nondeterministic lifting of X coordinates.

Indeed. The inconsistency stems from the following behaviour, which has nothing to do with any of the changes here (or earlier):

sage: K.<Y> = NumberField(x^3 - x^2 + 7*x + 7)
sage: Y <= -Y
False
sage: -Y <= Y
False

Therefore, in the failing doctest, the .lift_x() in point_of_order() has no way to identify the "smaller" point among (x,Y) and (x,-Y), so it returns a non-deterministic choice after all.

I'm not sure how to resolve this: It looks like this behavior is intended for number fields, since there's usually no mathematically meaningful ordering on their elements (in the sense of ordered fields). On the other hand, finite fields, which are equally un‑ordered, simply apply lexicographic comparison anyway. I suppose what this highlights is that Sage should — but currently doesn't — distinguish between mathematically meaningful orderings and simple comparison functions for practical purposes (sorting, deterministic choice). Fixing all that seems far beyond the scope of this ticket, though...

I think we have to either output just the x-coordinate (OK for a test, not so good for the documentation), or add a random tag,

Documentation preview for this PR (built with commit c039809; changes) is ready! 🎉

Thanks, I added a random tag. Since nothing else changed since you looked at this, I'll set it back to "positive review".