pure Python Bézier curve implementation I came up with this recursive pure-Python implementation of De Casteljau's algorithm for computing points on a Bézier curve:

def bezier_curve(control_points, number_of_curve_points):
    return [
        bezier_point(control_points, t)
        for t in (
            i / (number_of_curve_points - 1) for i in range(number_of_curve_points)
        )
    ]


def bezier_point(control_points, t):
    if len(control_points) == 1:
        result, = control_points
        return result
    control_linestring = zip(control_points[:-1], control_points[1:])
    return bezier_point([(1 - t) * p1 + t * p2 for p1, p2 in control_linestring], t)


Assumptions about control_points

The elements of control_points represent the control points of the Bézier curve. They must be of the same type of mutually compatible types fulfilling the following rules:


points shall all be of the same dimension.
Multiplying a point by a scalar shall result in a point of the same
dimension and with a value according to vector-scalar multiplication
(i.e., multiply each of the point's Cartesian coordinates with the scalar)
Adding two points shall result in a point of the same dimension and
with a value according to vector addition (i.e., component-wise
addition of the points' Cartesian coordinates)


Some examples that work as control_points:


list of turtle.Vec2D
list of complex
list of numpy.array with shape (2,)
numpy.array with shape (n, 2), where n is the number of control points


(tuples instead of lists work, too. Probably any sequencey container will work.)

Why pure Python?

Because I want this to be usable in a QGIS plugin, but SciPy, NumPy, etc. (usually) aren't available to QGIS plugins. As I'm unsure which Python libraries are available in QGIS (and the answer to that seems to be platform-dependent), I'd like to avoid external libraries (those that would have to be installed with pip or one of its alternatives) completely.

Using standard library functions should be fine, so if any part of the implementation could benefit from them, please point that out.

What I'd like to know in this review


Could / should the readability and comprehensibility of this implementation be improved?
Did I go foul of any performance (computation speed, memory usage, etc.) no gos? (It doesn't need to be super-fast, but it shouldn't be unnecessarily slow if I can avoid it.)


Performance for low degrees (e.g., degree 2, i.e. cubic Bézier with three control points per curve) will probably be more relevant than performance at high degrees (many control points per curve)
Performance for large outputs (large number_of_curve_points) might be relevant

About the destructuring assignment result, = control_points to unpack the single point while at the same time making sure it really is exactly one point


Is this idiomatic in Python (i.e. "pythonic")?
Is this readable and comprehensible enough or too obscure?
Is there any good alternative that's an expression, i.e. that could be used directly in the return statement without going through an assignment? (control_points[0] is an expression but doesn't fail when there's more than one element in control_points.)
Unnecessary Generator

You've got an unnecessary generator expression here:

def bezier_curve(control_points, number_of_curve_points):
    return [
        bezier_point(control_points, t)
        for t in (
            i / (number_of_curve_points - 1) for i in range(number_of_curve_points)
        )
    ]


You don't need to generate i / (n-1); you could simply pass that argument to the bezier_point() function:

def bezier_curve(control_points, number_of_curve_points):
    return [ bezier_point(control_points, i / (number_of_curve_points - 1))
             for i in range(number_of_curve_points)
           ]


Slight optimization: instead of computing number_of_curve_points - 1 \$O(N)\$ times (pure Python will not cache the result), precompute it:

def bezier_curve(control_points, number_of_curve_points):
    last_point = number_of_curve_points - 1
    return [ bezier_point(control_points, i / last_point )
             for i in range(number_of_curve_points)
           ]


Tail Recursion

Python does not do Tail Call Optimization, so with M control points, you will recursively enter and exit M calls, for each of the N points along your curve.  That is M*N unnecessary stack frame entry/exits.  You should do the looping yourself:

def bezier_point(control_points, t):
    while len(control_points) &gt; 1:
        control_linestring = zip(control_points[:-1], control_points[1:])
        control_points = [(1 - t) * p1 + t * p2 for p1, p2 in control_linestring]
    return control_points[0]


Since we loop while len(control_points) &gt; 1, it should be guaranteed that control_points will only have one point when the loop exits, so return control_points[0] is safe.  The exception is if the function is called with zero control points, but then control_points[0] will properly fail with an IndexError.