> **is there an algorithm to generate a random point on the curve?**

By having a good random source like `/dev/urandom` you can generate a random point $R$ on a curve at least by two methods;

1. **By using random scalar and scalar multiplication**

 1. Choose a generator point $P$ on the curve
 2. Get random integer between $0 < k < \text{Order of the Group}$
 3. Calculate $R =[k]P$ by scalar multiplication, prefable by [double-and-add][1] algorithm.

2. **Using random integer with sample rejection on the curve equation.**

 1. Choose a random integer $x \in \mathbb{Z_{25}}$.
 2. Check the equation $y^2 = x^3 + x + 1$ has a solution, i.e. it is a [quadratic residue][2].
 
 If there is no solution return to 1. step ( reject the sample)

 Else we have two solutions, $y$ and $-y$
 3. Now use the random number to select $\bar{y}$ as $y$ or $-y$, or toss a coin.
 4. Form the random point on the curve as $R=(x,\bar{y})$

3. **If you are using SageMath then it has a function**

 `random_element()` and the explanation is [given as][3] 
 
 _Return a random point on this elliptic curve, uniformly chosen among all rational points._

 SageMath `random_element()` function uses the 2. method with an addition that it can also return the point at infinity, $\mathcal{O}$

> **In my implementation points on this field are represented by polynomials, if that is relevant.**

The $\mathbb{F}_{25}$ is an Extension Field and it is usual to represent the elements of the field by the polynomial representation which is very advantageous on operations. As a result, the coordinates of points have polynomial representation.

**From comments of OP** 
> I used Sage to calculate the order of each point: (sagecell.sagemath.org/…) and it can be seen that even though the cardinality of the group is 32 the biggest order is 16. Is it just a matter of choosing an EC with prime cardinality?

The SageMath construction is not correct; 

`E = EllipticCurve(GF(25, 'x'),1, 1)` provides this curve $y^2 = x^3 + x + 2$ not $y^2 = x^3 + x + 1$

One can construct as below;

```
 E = EllipticCurve(GF(25, 'x'), [1, 1])
 print(E)
 print(E.abelian_group())
 print ( E.cardinality())
 for i in E.points():
 print (i.order())
```
The construction `EllipticCurve(Finite Field,[a,b])` is for the standard Weierstrass from $y^2 = x^3+ a x + b$ where $a,b \in F$ in which the curve defined. You only need to provide $a,b$, and $F$. The set of points commonly represented as $E_{a,b}(F)$ 

The `abelian_group()` functionality prints the abelian group structure of the group of points on this elliptic curve. The result is

$$\mathbb{Z_3} \oplus \mathbb{Z_9}$$

One should not call this function for a large group. There is a theory behind this;

**Theorem:** Let $E$ ben elliptic curve group over the finite field $\mathbb{F_q}$. Then $$E(\mathbb{F_q}) \simeq \mathbb{Z_p} \text{ or } \mathbb{Z_{n_1}} \oplus \mathbb{Z_{n_2}}$$ for some integer $n \geq 1$, or for some integers $n_1,n_2 \geq 1$ with $n_1|n_2$.

 [1]: https://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication#Double-and-add
 [2]: https://en.wikipedia.org/wiki/Quadratic_residue
 [3]: https://doc.sagemath.org/html/en/reference/curves/sage/schemes/elliptic_curves/ell_finite_field.html