Code snippetst from blender source to interpolate it's own bezier types

BKE_curve_forward_diff_bezier.c.inc

/* forward differencing method for bezier curve */
void BKE_curve_forward_diff_bezier(float q0, float q1, float q2, float q3, float *p, int it, int stride)
{
	float rt0, rt1, rt2, rt3, f;
	int a;

	f = (float)it;
	rt0 = q0;
	rt1 = 3.0f * (q1 - q0) / f;
	f *= f;
	rt2 = 3.0f * (q0 - 2.0f * q1 + q2) / f;
	f *= it;
	rt3 = (q3 - q0 + 3.0f * (q1 - q2)) / f;

	q0 = rt0;
	q1 = rt1 + rt2 + rt3;
	q2 = 2 * rt2 + 6 * rt3;
	q3 = 6 * rt3;

	for (a = 0; a <= it; a++) {
		*p = q0;
		p = (float *)(((char *)p) + stride);
		q0 += q1;
		q1 += q2;
		q2 += q3;
	}
}

blender_bpy_mathutils_interpolate_bezier.c.inc

PyDoc_STRVAR(M_Geometry_interpolate_bezier_doc,
".. function:: interpolate_bezier(knot1, handle1, handle2, knot2, resolution)\n"
"\n"
"   Interpolate a bezier spline segment.\n"
"\n"
"   :arg knot1: First bezier spline point.\n"
"   :type knot1: :class:`mathutils.Vector`\n"
"   :arg handle1: First bezier spline handle.\n"
"   :type handle1: :class:`mathutils.Vector`\n"
"   :arg handle2: Second bezier spline handle.\n"
"   :type handle2: :class:`mathutils.Vector`\n"
"   :arg knot2: Second bezier spline point.\n"
"   :type knot2: :class:`mathutils.Vector`\n"
"   :arg resolution: Number of points to return.\n"
"   :type resolution: int\n"
"   :return: The interpolated points\n"
"   :rtype: list of :class:`mathutils.Vector`'s\n"
);
static PyObject *M_Geometry_interpolate_bezier(PyObject *UNUSED(self), PyObject *args)
{
	VectorObject *vec_k1, *vec_h1, *vec_k2, *vec_h2;
	int resolu;
	int dims;
	int i;
	float *coord_array, *fp;
	PyObject *list;

	float k1[4] = {0.0, 0.0, 0.0, 0.0};
	float h1[4] = {0.0, 0.0, 0.0, 0.0};
	float k2[4] = {0.0, 0.0, 0.0, 0.0};
	float h2[4] = {0.0, 0.0, 0.0, 0.0};


	if (!PyArg_ParseTuple(args, "O!O!O!O!i:interpolate_bezier",
	                      &vector_Type, &vec_k1,
	                      &vector_Type, &vec_h1,
	                      &vector_Type, &vec_h2,
	                      &vector_Type, &vec_k2, &resolu))
	{
		return NULL;
	}

	if (resolu <= 1) {
		PyErr_SetString(PyExc_ValueError,
		                "resolution must be 2 or over");
		return NULL;
	}

	if (BaseMath_ReadCallback(vec_k1) == -1 ||
	    BaseMath_ReadCallback(vec_h1) == -1 ||
	    BaseMath_ReadCallback(vec_k2) == -1 ||
	    BaseMath_ReadCallback(vec_h2) == -1)
	{
		return NULL;
	}

	dims = max_iiii(vec_k1->size, vec_h1->size, vec_h2->size, vec_k2->size);

/// --- we are curious the thing that happens below: 

	for (i = 0; i < vec_k1->size; i++) k1[i] = vec_k1->vec[i];
	for (i = 0; i < vec_h1->size; i++) h1[i] = vec_h1->vec[i];
	for (i = 0; i < vec_k2->size; i++) k2[i] = vec_k2->vec[i];
	for (i = 0; i < vec_h2->size; i++) h2[i] = vec_h2->vec[i];

	coord_array = MEM_callocN(dims * (resolu) * sizeof(float), "interpolate_bezier");
	for (i = 0; i < dims; i++) {
		// this function is defined in BKE_curve.h
		BKE_curve_forward_diff_bezier(k1[i], h1[i], h2[i], k2[i], coord_array + i, resolu - 1, sizeof(float) * dims);
	}

/// --- 

	list = PyList_New(resolu);
	fp = coord_array;
	for (i = 0; i < resolu; i++, fp = fp + dims) {
		PyList_SET_ITEM(list, i, Vector_CreatePyObject(fp, dims, Py_NEW, NULL));
	}
	MEM_freeN(coord_array);
	return list;
}