Find cubic Bézier control points given four points What I need is to generate an SVG file while having a series of (x,y) ready.
P0(x0,y0)
P1(x1,y1)
P2(x2,y2)
P3(x3,y3)
P4(x4,y4)
P5(x5,y5)
...

I need to make a Bézier curve from them so I need to calculate control points (mid-points) for them. This image explains what I am exactly looking for:

I have points P0, P1, P2 and P3 ready. What I need is to calculate control points C1 and C2. The curve does not pass them. But it bends toward them.
I need a formula which gives me C1 and C2 in clear direct form:
C1= fomula1 (P0,P1,P2,P3)
C2= fomula2 (P0,P1,P2,P3)

I was thinking about least square method and some other methods but I had no idea how to implement them.
References: Animation, SVG
 A: Your problem, as stated, does not have a unique solution. Suppose that point $P_j$ is at location $(3j, 0)$, for each integer $j$, so that they're equi-spaced on the $x$-axis. Now let $y$ be any real number. Then by adding control points at locations
$$(6i+1, y)\\ (6i+2, y)\\ (6i + 4, -y)\\ (6i+5, -y)$$
for each integer $i$, 
you get two "control points" between any two of your original points. For instance, near the origin, for $y = 2$, we have points
$$
(-3, 0) \leftarrow (\mathrm{one~ of~ the}~ P_i)\\
(-2, -2)\\
(-1, -2)\\
(0, 0) \leftarrow (\mathrm{one~ of~ the}~ P_i)\\
(1, 2)\\
(2, 2)\\
(3, 0) \leftarrow (\mathrm{one~ of~ the}~ P_i)
$$
These determine two Bezier segments that glue up nicely at the origin, with a slope of $2$ at the origin. 
You may say "But it's obvious that the control points in this case should be on the $x$-axis!" and I say "but your problem statement doesn't require it." Indeed, I chose this example because it was easy to write, but given any set of $P_i$, I can again find an infinitely family of ways to place the intermediate control points so as to join the $P_i$ with Bezier segments. 
I'm going to suggest that you consider looking at Catmull-Rom splines, which are piecewise cubics passing through a sequence of points like your $P_i$. Each segment of a CR-spline can be expressed as a Bezier curve, because the Bezier basis functions span the space of cubic curves. One detailed reference on this is Computer Graphics: Principles and Practice, 3rd edition, of which I am a coauthor, but there are plenty of other references as well. 
Here are somewhat brief details on CR spline construction from a sequence of points $M_0, M_1, \ldots$. I'm going to describe how to find the control points for the part of the curve between $M_1$ and $M_2$, so as to avoid any negative indices. The four control points will be $P_0, P_1, P_2, P_3$. Two of these are easy:
\begin{align}
P_0 &= M_1 \\
P_3 &= M_2
\end{align}
so that the Bezier curve starts and ends at $M_1$ and $M_2$, respectively. 
The other two are only slightly trickier. We compute
\begin{align}
v_1 &= \frac{1}{2} (M_2 - M_0)\\
v_2 &= \frac{1}{2} (M_3 - M_1)
\end{align}
which are the velocity vectors at $M_1$ and $M_2$. We then have 
\begin{align}
P_1 &= P_0 + \frac{1}{3} v_1 = M_1 + \frac{1}{6} (M_2 - M_0)\\
P_2 &= P_4 - \frac{1}{3} v_2 = M_2 - \frac{1}{6} (M_3 - M_1)
\end{align}
Applying these rules in the example I gave earlier, with 
$$
M_i = (3i, 0)
$$
we have
\begin{align}
M_0 &= (0, 0)\\
M_1 &= (3, 0)\\
M_2 &= (6, 0)\\
M_3 &= (9, 0)
\end{align}
so that 
\begin{align}
P_0 &= (3, 0)\\
P_3 &= (6, 0)\\
v_1 &= \frac{1}{2}((6, 0) - (0, 0)) = (3, 0)\\
v_2 &= \frac{1}{2}((9, 0) - (3, 0)) = (3, 0)\\
P_1 &= P_0 + \frac{1}{3} v_1 = (3, 0) + (1, 0) = (4, 0)\\
P_2 &= P_3 - \frac{1}{3} v_2 = (6, 0) - (1, 0) = (5, 0)
\end{align}
as expected. 
Hint for the start and end points: Assuming you have a sequence of points 
$$
M_0, M_1, \ldots, M_{n-1}
$$
you can let
\begin{align}
M_{-1} &= M_0 - (M_1 - M_0) = 2M_0 - M_1 \\
M_n &= M_{n-1} + (M_{n-1} - M_{n-2}) = 2M_{n-1} - M_{n-2}
\end{align}
to extend your list just enough that the CR scheme above provides interpolation all the way from $M_0$ to $M_{n-1}$. 
