The two ways that you mention, plus looking your curve up in a database (as Matt E. suggests), are the most practical and efficient ways I can think of. Here are two other ways to do this, one practical (but not as efficient) and one which is by no means practical (I think):
Use division polynomials. If $E$ is defined over $\mathbb{Q}$, then $E(\mathbb{Q})$ can only have torsion points of order up to certain bounds. If $P$ is a point of prime order, then the order is $2$, $3$, $5$, or $7$ (by Mazur's theorem). You can find the division polynomials for each of these primes, and see if they have any roots in $\mathbb{Q}$. A root in $\mathbb{Q}$ would give you a rational $x$-coordinate of a torsion point of said order. You can also find the division polynomials of order $p^n$. If a point $P$ has order $p^n$, then the order is at most $8$, $9$, $5$, or $7$, so you only really need to find the $8$th, $9$th, $5$th and $7$th division polynomials, and factor those, to find all the torsion points on $E(\mathbb{Q})$.
Not practical, and conjectural: Use the Birch and Swinnerton-Dyer conjecture. BSD provides a conjectural formula for a certain Taylor coefficient of the $L$-series of the elliptic curve in question. The size of the torsion subgroup appears in the denominator of the coefficient... However, this coefficient is a real number. The real period $\Omega_E$ can always be calculated. If you could calculate the regulator $R_E$ of your elliptic curve (for instance, if you know the rank is $0$, then $R_E=1$), then you could calculate $\Omega_E\cdot R_E$, calculate the coefficient computationally, divide the result by $\Omega_E\cdot R_E$ and obtain a rational number, whose denominator is a divisor of the square of the size of the torsion subgroup of $E$ (notice that there may be some cancellation with the numerator!). This information would then allow you to find the torsion points on $E$ (see previous bullet point).
Example (with division polynomials): Let $E$ be the curve 53b3 in Cremona's database, with Weierstrass equation
$$E: y^2 +xy+y= x^3-x^2-14x+29.$$
Let us first find division polynomials:
($p=2$) The 2nd division polynomial, defining the $x$-coordinates of $2$-torsion points, is given by $4x^3 - 3x^2 - 54x + 117$. This is irreducible, so there are no $2$-torsion points defined over $\mathbb{Q}$.
($p=3$) The 3rd division polynomial is
$$3x^4 - 3x^3 - 81x^2 + 351x - 270,$$
and it factors as $3(x - 1)(x^3 - 27x + 90)$. Thus, there is a point of order $3$ with $x$-coordinate $1$. By plugging into the equation of $E$, we find that $(1,3)$ and $(1,-5)$ are points of order $3$ on $E$. Since we found a point of order $3$, we need to keep going and look for points modulo $9$, $27$ (can't happen over $\mathbb{Q}$), etc...
($9$) The $9$ division polynomial has degree $40$ so we won't reproduce it here. However, it has a bunch of linear factors in its factorization:
$$(x - 9)(x - 3)(x - 1)(x + 3)\cdot (\text{higher order factors}).$$
The points with $x=1$ are of no interest to us, since they are points of order $3$. The other coordinates $x=\pm 3$ and $x=9$ correspond to points of order exactly 9, so we need to check if their $y$-coordinates are in $\mathbb{Q}$. Indeed, they do correspond to points of order $9$, namely
$$(3,1), (-3,7), (9,-29),(9,19),(-3,-5), (3,-5).$$
Notice, however, that $(3,1)$ generates all of them. We move on to order $27$, since we found points of order $9$.
- ($27$) The 27th division polynomial is of degree $364$, but all the linear factors in its factorization already appeared in the 9th division polynomial, so there are no points of order $27$ (again, there are NO points of order $27$ on elliptic curves over $\mathbb{Q}$).
($p=5$) The 5th division polynomial is of degree $12$ and irreducible. Thus, there are no points of order $5$.
($p=7$) The 7th division polynomial is of degree $24$ and irreducible. Thus, there are no points of order $7$.
Hence, $E(\mathbb{Q})_\text{tors} \cong \mathbb{Z}/9\mathbb{Z}$.
Example (with BSD and L-functions): Let $E$ be the curve 53b3 in Cremona's database, with Weierstrass equation
$$E: y^2 +xy+y= x^3-x^2-14x+29.$$
First, we perform a $2$-descent on $E$ to calculate the $2$-Selmer group. It turns out to be trivial. This means two things: the rank is $0$, and Sha is trivial. In particular, the regulator is $R_E=1$. Moreover:
- The real period is $\Omega_E=3.09156554910300755665231500284\cdots$
- The Tamagawa numbers are $c_2=9$ and $c_3=3$ for $p=2$ and $3$ respectively.
- The value of the L-function at $1$ can be calculated to be
$$L(E,1)\approx 1.03052184970100251888410500095\cdots$$
If we believe BSD in this case, then we must have:
$$(\# E(\mathbb{Q})_\text{tors})^2 = \frac{\# \text{Sha} \cdot \Omega_E\cdot R_E\cdot \prod c_p}{L(E,1)} = \frac{1\cdot (3.09156554910300755665231500284\cdots)\cdot 1 \cdot 27}{1.03052184970100251888410500095\cdots} =81.000000000000000000000000000\cdots.$$
Thus, since $(\# E(\mathbb{Q})_\text{tors})^2$ must be an integer, it must be $81$, and $\# E(\mathbb{Q})_\text{tors} = 9$. Since $ E(\mathbb{Q})_\text{tors}$ is an abelian group, there are only two options: $\mathbb{Z}/3\mathbb{Z}\times \mathbb{Z}/3\mathbb{Z}$ or $\mathbb{Z}/9\mathbb{Z}$. However, the former is imposible because if $E[n]$ is defined over a field $F$, then $F$ must contain the $n$th roots of unity. But $\mathbb{Q}$ does not contain $\sqrt{-3}$. Hence, it must be $\mathbb{Z}/9\mathbb{Z}$.
Magma code Here is Magma code for all I did above:
E:=EllipticCurve("54b3");E;
DivisionPolynomial(E,2);
Factorization(DivisionPolynomial(E,2));
DivisionPolynomial(E,3);
Factorization(DivisionPolynomial(E,3));
IsPoint(E,1);
P:=E![1,3,1];
2*P; 3*P;
DivisionPolynomial(E,9);
Factorization(DivisionPolynomial(E,9));
IsPoint(E,3);
IsPoint(E,9);
P:=E![3,1,1];
P,2*P,3*P,4*P,5*P,6*P,7*P,8*P,9*P;
Degree(DivisionPolynomial(E,27));
// Factorization(DivisionPolynomial(E,27));
DivisionPolynomial(E,5);
Factorization(DivisionPolynomial(E,5));
DivisionPolynomial(E,7);
Factorization(DivisionPolynomial(E,7));
TwoSelmerGroup(E);
RealPeriod(E);
TamagawaNumbers(E);
L:=LSeries(E);
Evaluate(L,1);
27*RealPeriod(E)/Evaluate(L,1);
The key fact you need is this: if $E/K$ is an elliptic curve which has complex multiplication over $K$, then the associated $\ell$-adic Galois representations
$$\rho_{E,\ell}:G_K\to\mathrm{GL}_2(\overline{\mathbb Q}_\ell)$$
are reducible for all primes $\ell$. Indeed, $\mathrm{End}(E)\otimes\overline{\mathbb Q}_\ell$ embeds into the endomorphism ring of $\rho_{E,\ell}$, and the former is not a division ring.
In your case, your curve $E: y^2 = x^3-1$ does not have CM over $\mathbb Q$, but it picks up extra endomorphisms over $K=\mathbb Q(\zeta_3)$. It follows that $\rho_{E,\ell}$ is irreducible, but $\rho_{E,\ell}|_{G_K}$ becomes reducible, and representation theory kicks in to show that
$$\rho_{E,\ell}\cong \mathrm{Ind}_{G_K}^{G_\mathbb Q}\chi_\ell\qquad(*)$$
where $\chi_\ell$ is a character -- actually it is the Galois character corresponding to the Hecke grossencharacter which is attached to $E$ by the theory of CM.
Finally, one can show using pure representation theory that $(*)$ is equivalent to the statement
$$\rho_{E,\ell}\cong\rho_{E,\ell}\otimes\theta\qquad(**)$$
where $\theta:G_\mathbb Q\to\overline{\mathbb Q}_\ell^\times$ is the one dimensional representation which has as its kernel $G_K$. Explicitly, $\theta$ is a lift of the unique non-trivial Dirichlet character of conductor $3$.
In particular:
- If $p\equiv 2\pmod 3$, then $\theta(\mathrm{Frob}_p) = -1$ and it follows from $(**)$ that $a_p = -a_p$, so $a_p = 0$.
- If $p\equiv 1\pmod 3$, then it follows from $(*)$ and the usual formula for the trace of an induced representation that $a_p = 2\chi_\ell(\mathrm{Frob}_p)$, which turns out to be the value you'd expect.
Best Answer
A point $P = (x_0,y_0)$ is two torsion if and only if $y_0= 0$, you can see this from the definition of doubling on an elliptic curve, the tangent line goes to infinity only for a vertical tangent line, which happens only for points with $y$-coordinate zero (take the partial derivatives of the defining equation). So to find the two torsion points we simply have to factor $x^3 + ax + b$ over the rationals.
More generally however this sort of analysis of tangent lines to find torsion points gets harder, but there is a theorem by Lutz-Nagell (https://en.wikipedia.org/wiki/Nagell%E2%80%93Lutz_theorem) that classifies torsion points in a very explicit way for elliptic curves over the rationals. This theorem always gives us a finite list of possible points that we can check if they are torsion (if they even define points on the curve) by repeatedly adding one point to itself. Note that as soon as we get non-integer coordinates for a multiple of our point we know it cannot have been torsion so we can stop.