How to find the number and coordinates of self-intersections points for a polygon? I have a  self-intersecting polygon defined by $n$ points on the plane:
A <- t(matrix(c(
 0,   0,
 3,   0,
 3,  -1,
 2,  -1,
 2,   2,
 3,   2,
 3,   1,
 0,   1,
 0,   2,
 1,   2,
 1,  -1,
 0,  -1,
 0,   0), nrow=2));

The number of points is less $20$, coordinates are integer.
You can see four points of self-intersections in $(1,1)$, $(2,1)$, $(2,0)$, $(1,0)$ on the figure:
plot(A, col='red', type= 'l', xlim=c(min(A[,1]),max(A[,1])),
      ylim=c(min(A[,2]),max(A[,2])), xlab='x', ylab='y'); 
points(A, col='black', pch = 22); 
grid() 


Question I am looking for an algorithm to define the number and coordinates of self-intersections points.
How to find the number and coordinates of self-intersections points?
 A: An approach that will work for a general self-intersecting polygon in 3-dimensional space is as follows:
For each adjacent vertex pair ${i,j}$, loop over every other set of adjacent vertex pairs ${m,n} : m\ne i\ne j,\hspace {0.5 cm} n\ne i\ne j$ and invoke the following procedure:
The position vector of any point along the line  joining vertices $i$ and $j$ is
${\bf \vec x} = {\bf \vec x_i} + \xi L_{ij} \bf \vec q \tag{EQ. 1}$ 
and likewise that of any point along the line  joining vertices $m$ and $n$ is
${\bf \vec x} = {\bf \vec x_m} + \eta L_{mn}\bf \vec r \tag{EQ. 2}$
where
$$
L_{ij}=\|{\bf \vec x_j -\bf \vec x_i}\|\\
L_{mn}=\|\bf \vec x_n -\bf \vec x_m\|$$
and
$${\bf \vec q}=\frac{1}{L_{ij}}({\bf \vec x_j} -{\bf \vec x_i})\\
{\bf \vec r}=\frac{1}{L_{mn}}({\bf \vec x_n -\bf \vec x_m})$$
At the intersection, we have
${\bf \vec p} + {\bf \vec q}L_{ij}\xi -{\bf \vec r} L_{mn}\eta  =0 \tag{EQ. 3}$ 
where ${\bf \vec p}={\bf \vec x_i} -{\bf \vec x_m}$. Now you can find a unit vector $\bf\vec s$ orthogonal to $\bf\vec r$ as a linear combination of the unit vectors  $\bf\vec q$ and $\bf\vec r$ as follows:
$\bf\vec s = \alpha \bf\vec q + \beta  \bf\vec r\tag{EQ. 4}$
From the condition $\bf\vec s\cdot\bf\vec r=0$ we have
$\alpha \bf\vec q\cdot\bf\vec r +\beta=0\tag{EQ. 5}$
and from the stipulation that $\bf\vec s$ be a unit vector, we have
$\alpha^2 \bf\vec q\cdot\bf\vec q +\beta^2\bf\vec r\cdot\bf\vec r+2\alpha\beta\bf\vec q\cdot\bf\vec r=1\tag{EQ. 6}$
From EQS. 5 and 6, $\beta$ can be eliminated giving an equation for $\alpha$ whereupon only the positive root is relevant. Following this, $\beta$ is given by EQ. 5, and $\bf\vec s$ from EQ. 4
Now, taking the inner product of EQ. 3 with $\bf\vec s$ gives
${\bf \vec s \cdot\bf \vec p} + L_{ij}(\bf \vec s \cdot\bf \vec q)\xi =0\tag{EQ. 7}$
from which you can evaluate $\xi$.  If $\xi<0$ or $\xi>1$, then you abort the procedure (the intersection is outside the lines joining the vertices $i$ and $j$), and proceed to the next pair. 
Otherwise, calculate $\bf\vec x$ from EQ. 1 and then evaluate $\eta$ from
$\eta= \frac{1}{L_{mn}}(\bf\vec x-\bf\vec x_m)\cdot\bf\vec r
\tag{EQ. 8}$
If $0\le\eta\le1$, then $\bf\vec x$ is an intersection.
This is a brute-force approach of $O(n^2)$. A more efficient alternative may be the Bentley-Ottman Algorithm
