What you need is something called "Boehm's algorithm" (after its originator, Wolfgang Boehm). It has a simple geometric interpretation, and drawing a few pictures should make it clear. There is a pretty good explanation (with pictures) in this document.
The algorithm is based on a process called "knot insertion". You keep inserting knots into the b-spline curve until each knot has multiplicity 3. Then, the b-spline control points of this refined curve give you the Bezier control points of its segments.
So, if you're writing code to do this, one approach is to write a knot insertion function first, and then call it repeatedly.
There is knot insertion code here.
It looks like your set of endpoints and control points can be
any set of points in the plane.
This means that the $order$ of the points
is critical,
so that the generated curve goes through the points
in a specified order.
This is much different than the ordinary interpolation problem,
where the points
of the form $(x_i, y_i)$
are ordered so that
$x_i < x_{i+1}$.
As I read your desire,
if you gave a set of points on a circle
ordered by the angle of the line
from the center to each point,
you would want the result to be
a curve close to the circle.
There are a number of ways this could be done.
I will assume that
you have $n+1$ points
and your points are $(x_i, y_i)_{i=0}^n$.
The first way I would do this
is to separately parameterize the curve
by arc length,
with $d_i = \sqrt{(x_i-x_{i-1})^2+(y_i-y_{i-1})^2}$
for $i=1$ to $n$,
so $d_i$ is the distance from
the $i-1$-th point to the
$i$-th point.
For a linear fit,
for each $i$ from $1$ to $n$,
let $t$ go from
$0$ to $d_i$
and construct separate curves
$X_i(t)$ and $Y_i(t)$
such that
$X_i(0) = x_{i-1}$,
$X_i(d_i) = x_i$,
and
$Y_i(0) = y_{i-1}$,
$Y_i(d_i) = y_i$.
Then piece these together.
For a smoother fit,
do a spline curve
through each of
$(T_i, x_i)$
and
$(T_i, y_i)$
for $i=0$ to $n$,
where
$T_0 = 0$
and
$T_i = T_{i-1}+d_i$.
To get a point for any $t$ from
$0$ to $T_n$,
find the $i$ such that
$T_{i-1} \le t \le T_i$
and then,
using the spline fits
for $x$ and $y$
(instead of the linear fit),
get the $x$ and $y$ values from their fits.
Note that
$T_i$ is the cumulative length
from $(x_0, y_0)$
to $(x_i, y_i)$,
and $T_n$ is the total length of the line segments
joining the consecutive points.
To keep the curves from
not getting too wild,
you might look up "splines under tension".
Until you get more precise,
this is as far as I can go.
Best Answer
It is not clear whether the "control points" you need is for Bezier curve or for B-spline curve. Since many graphic toolkits only take quadratic/cubic Bezier curves, I will assume this is what you want.
For creating multiple Bezier curves interpolating a given set of data points, you can go with Catmull-Rom spline interpolation or natural spline interpolation. Both schemes will produce a cubic Bezier curve in between each two data points but the natural spline interpolation will require solving a linear equation set.
If you really want to do B-spline interpolation, please refer this link for more details. You will have to decide a knot vector in advance, which is typically derived from the parameters for the points so that the linear equation matrix is not ill-conditioned.