finding 4-torsion point on elliptic curve

edit

I think the method you are looking for is called division_points in Sage.

See the SageMath documentation for division points of elliptic curves 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^8 + 5^9 + 2*5^10 + 3*5^11 + 2*5^12 + 3*5^14 + 4*5^15 + 3*5^16 + 4*5^17 + 4*5^18 + 3*5^19 + O(5^20) : 1 + O(5^20))]

Here is 7-adic:

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