$y^2 = x^3 + ax + b\bmod p$ is the Short Weierstrass equation. The theory behind it is here
Using Bezout’s Theorem, it can be shown that every irreducible cubic has a flex (a point where the tangent intersects the curve with multiplicity three) or a singular point (a point where there is no tangent because both partial derivatives are zero). [Reducible cubics consist of a line and a conic, which are easy to study.]
An irreducible cubic with a flex can be affinely transformed into a Weierstrass equation:
$$Y^2 + a_1 XY + a_3 Y = X^3 + a_2 X^2 + a_4 X + a_6$$
And this can be converted into a short Weierstrass equation if $p$ is larger than $3$. The transformation can be achieved by
Change $Y\longrightarrow y-(a_1x+a_3)/2$ and $X \longrightarrow x$
$$y^2 = x^3 - (1/4 a_1^2 +a_2 )x^2 + (1/2 a_1 a_3)x + (a_4+ 1/4 a_3^2 + a_6)$$
with renaming
\begin{align}
A & = - (1/4 a_1^2 +a_2 )\\
B & = 1/2 a_1 a_3\\
C & = a_4+ 1/4 a_3^2 + a_6
\end{align}
now the new equation has of the form $$y^2=x^3+Ax^2+Bx+C.$$
And now change $x\longrightarrow x-A/3$, $$y^2 = x^3 - \left(1/3A^2x + Bx \right)+ \left( 2/27 A^3 - 1/3AB + C \right)$$
with renaming;
\begin{align}
a &= 1/3A^2 + B\\
b & = 2/27 A^3 - 1/3AB + C\\
\end{align}
now the new equation has of the form $$y^2=x^3+ax+b.$$
The reason for $p \neq 2,3$ can be seen from the divisions. Also one should note that the two changes of variables are invertible, therefore their composition is too. This concludes that the rational points on both equations are in bijective correspondence.
For $p=2$ and $p=3$, there are other short forms. In Cryptography we work in either binary extension fields ($\mathbb{F}_{2^m}$), which are no more secure, large prime fields ($\mathbb{F}_{p}$), and their extensions ($\mathbb{F}_{p^m}$).
If the discriminant is zero ($\Delta = 0$), then the curve is singular and the DLP is easy. Geometrically, this means that the graph has cusps, self-intersections (nodes), or isolated points.
Cusp; when $\Delta=0$ and $a=0$ ( $y^2 = x^3$);
Intersection or Node; when $\Delta=0$ and $a\neq 0$ ( $y^2 = x^3 +3x +1$);
Acnode (isolated point); $y^2 + x^3 + x^2$, at $(0,0)$
Note these plots are not in finite fields where we cannot see such figures easily. A curve over a finite field is just points when we draw on the plane as below, and we cannot see the above singular points; however, we can still calculate them.
The above is the point of the curve over $\operatorname{GF}(307)$ with the curve equation $ y^2 = x^3 + x^2 + x$ and this curve's group of points is isomorphic to $ \mathbb Z_{160} \oplus \mathbb Z_2$
And, Wolfram has a nice page about those singular points.
This is not completely true. Some curves have prime order (i.e., $\#E(\mathbb F_p)$ is a prime), and they are called prime curves. All of the elements other than the identity can be a generator. This is the direct result of Lagrange's theorem on the group theory, i.e., the order of the subgroup generated by an element must divide the order of the group. To see what kind of groups can be formed by the points of a curve over a finite field see the theorem on the bottom.
Some curves don't have prime order, like the Curve25519 and Curve448. Curve25519 has a co-factor $h=8$ i.e. $h = \frac{\#E(\mathbb F_p)}{n}$ where $n$ is the order of the generator point $G$. In this case, we chose the largest prime to form the subgroup and select a generator for it.
There is also twist security related to the order that one must consider. Although the NIST P-224 has co-factor 1 ( i.e. it has prime order), it doesn't have the twisted security.
Theorem: Let $E$ be an 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$.
If the order of the curve is prime then it must be $E(\mathbb{F_q}) \simeq \mathbb{Z_p}$.
SageMath code for the transformation. First part;
var('x y z a b c')
f = y^2 - x^3- a*x^2-b*x-c
print(f)
g = f.subs(x == x-a/3)
print(g)
g.expand()
The second part;
var('X Y Z x y z a1 a2 a3 a4 a5 a6')
f = Y^2 + a1 * X*Y + a3*Y - X^3 - a2*X^2 - a4*X - a6
f
g = f.subs(Y == Y - (a1*X +a3)/2)
print(g)
g.expand()
