ag.algebraic geometry - How do you explicitly compute the p-torsion points on a general elliptic curve in Weierstrass form?

Here is an attempt to answer my own question, using the "division
polynomials" of Kevin's and Jared's answers. It is probably the maximally naive idea, and I do not claim it
works, though it's not clear to me that it can't.
I've community wiki-ed this answer, as it's probably junk anyway ...

Fix $A,Bin mathbb{Z}$, obtaining an elliptic curve $E$ over $S=Spec
mathbb{Z}[Delta^{-1}]$, and a prime $pgeq 5$ not dividing the
discriminant $Delta=Delta(A,B)$. The
division polynomial $psi_p(t)$ is supposed to be a polynomial of
degree $d=(p^2-1)/2$ over $mathbb{Z}$, whose roots are the values
$t(P)=x(P)/z(P)$ as $P$ ranges over the points of exact order $p$ in
$E$.

Turn $psi_p$ into a homogeneous polynomial $g$ of degree $d$ in
$mathbb{Z}[x,y,z]$, so that $g(x,y,1)=psi_p(x)$. The polynomial
deterimes a curve $C=(g)$ in $P^2/S$, and thus a closed subscheme
$D=Ecap C$ of $E$. Over $mathbb{Z}[Delta^{-1},p^{-1}]$, $D$ should
consist of
the points of exact order $p$ on $E$ (with multiplicity $1$), together
with the basepoint of $E$ with multiplicity $d$.

Claim. $D$ is an effective Cartier
divisor on $E/S$, of degree $3d$.

Proof. I don't know. I need to prove that $Dto
S$ is flat, the main concern being the behavior over
$mathbb{Z}_{(p)}$. I don't even know if this is really plausible in general.

Let's pretend we somehow know $D$ is an effective Cartier divisor on $E$
relative to the base $S$. There is another relative Cartier
divisor, namely
$$ D' = E[p] quad + quad (d-1)[0]$$
where $E[p]$ is the $p$-torsion subgroup scheme of $E$, and $[0]$ is the
degree one divisor of the basepoint of $E$.
It seems clear that away from characteristic $p$, the divisors $D$
and $D'$ are equal. Equality of effective divisors on a smooth curve
is a closed condition, so they should be equal over all of
$S$.

The divisor I want is thus $D''=D-d[0]=D'-d[0]$ (which is still
effective). Then it's really easy to find a homogeneous polynomial
$h$ of degree $2d/3$ which defines $D''$; because of the form of the
Weierstrass equation, you can produce it from $g$ by hand, and if my claim is true you can produce it globally, i.e., with coefficients in $mathbb{Z}[A,B]$.

I've carried out this out in the case $p=5$, and it appears to "work".
That is, I get an answer which appears sane for general $A$ and $B$, and
which for some explicit cases I've tried of $A,Bin mathbb{Z}$ appears to give me a flat $Dto S$.
For instance, if $E/mathbb{Z}[6^{-1}]$ is the curve with $(A,B)=(0,1)$ (which reduces to a supersingular curve at $p=5$), I find
$$h= 729z^{8}-1350x^4z^4+360x^6z^2+5x^8.$$