I think the method you are looking for is called `division_points`

in Sage.

See the documentation at
http://doc.sagemath.org/html/en/reference/curves/sage/schemes/elliptic_curves/ell_point.html#sage.schemes.elliptic_curves.ell_point.EllipticCurvePoint_field.division_points
which you can find in a web search engine with the query [elliptic curve torsion point sagemath].

Your example, with a guess for what Qp might be.

```
sage: Qp = pAdicField(17)
sage: E = EllipticCurve(Qp, [0, 0, 0, -3267, 45630])
```

Starting with the point of x-coordinate zero as P --- works fine.

```
sage: P = E(0); P
(0 : 1 + O(17^20) : 0)
sage: P.division_points(4)
[(15 + O(17^20) : O(17^10) : 1 + O(17^20)), (0 : 1 + O(17^20) : 0)]
```

Starting with your choice of P, I get an error --- but maybe I have the wrong Qp.

```
sage: P = E([-21, 324]); P
(13 + 15*17 + 16*17^2 + 16*17^3 + 16*17^4 + 16*17^5 + 16*17^6 + 16*17^7 + 16*17^8 + 16*17^9 + 16*17^10 + 16*17^11 + 16*17^12 + 16*17^13 + 16*17^14 + 16*17^15 + 16*17^16 + 16*17^17 + 16*17^18 + 16*17^19 + O(17^20) : 1 + 2*17 + 17^2 + O(17^20) : 1 + O(17^20))
sage: P.division_points(4)
Traceback (most recent call last)
...
PrecisionError: p-adic factorization not well-defined since the discriminant is zero up to the requestion p-adic precision
```

Different choices of Qp give different results.

Here is 5-adic:

```
sage: Qp = pAdicField(5)
sage: E = EllipticCurve(Qp, [0, 0, 0, -3267, 45630])
sage: P = E([-21, 324]); P
(4 + 4*5^2 + 4*5^3 + 4*5^4 + 4*5^5 + 4*5^6 + 4*5^7 + 4*5^8 + 4*5^9 + 4*5^10 + 4*5^11 + 4*5^12 + 4*5^13 + 4*5^14 + 4*5^15 + 4*5^16 + 4*5^17 + 4*5^18 + 4*5^19 + O(5^20) : 4 + 4*5 + 2*5^2 + 2*5^3 + O(5^20) : 1 + O(5^20))
sage: P.division_points(4)
[(4 + 3*5 + 4*5^2 + 5^3 + 5^4 + 5^5 + 3*5^6 + 3*5^7 + 3*5^8 + 4*5^9 + 4*5^10 + 4*5^11 + 5^12 + 3*5^13 + 3*5^14 + 5^15 + 3*5^16 + 3*5^18 + O(5^20) : 1 + 4*5 + 4*5^2 + 4*5^4 + 2*5^5 + 4*5^7 + 2*5^8 + 3*5^9 + 5^10 + 4*5^12 + 4*5^13 + 3*5^14 + 3*5^15 + 4*5^17 + 2*5^18 + 4*5^19 + O(5^20) : 1 + O(5^20)),
(2 + 3*5 + 2*5^2 + 5^3 + 5^4 + 2*5^5 + 2*5^6 + 3*5^7 + 5^8 + 5^9 + 5^11 + 4*5^12 + 5^13 + 3*5^14 + 2*5^16 + 5^17 + 5^19 + O(5^20) : 2 + 5^4 + 2*5^5 + 3*5^6 + 4*5 ...
```

(more)
We are missing a definition of

`Qp`

to be able to play with your code.I think the terminology is as follows.

Points Q such that 4 Q = P are called 4-division points of P on E.

In the special case of P = E(0), they are called 4-torsion points of E.

@slelievre Thank you. That definitely make sense.