A bunch of cubic Bézier functions

Bezier.cs

using System.Collections.Generic;
using UnityEngine;

public static class Bezier
{
    public static Vector3 Cubic(Vector3 a, Vector3 b, Vector3 c, Vector3 d, float t) {
        float r = 1f - t, r2 = r * r, t2 = t * t;
        return (r2 * r) * a + (3f * r2 * t) * b + (3f * r * t2) * c + (t2 * t) * d;
    }

    public static Vector3 CubicDerivative(Vector3 a, Vector3 b, Vector3 c, Vector3 d, float t) {
        float r = 1f - t, r2 = r * r, t2 = t * t;
        return (-3f * r2) * a + (9f * t2 - 12f * t + 3f) * b + (6f * t - 9f * t2) * c + (3f * t2) * d;
    }

    // Returns closest t and squared distance
    public static Vector2 CubicProject(Vector3 a, Vector3 b, Vector3 c, Vector3 d, Vector3 p, int iterations = 10) {
        float min = 0f, max = 1f;

        var distSqToMin = (p - a).sqrMagnitude;
        var distSqToMax = (p - d).sqrMagnitude;

        while (iterations-- > 0) {
            if (distSqToMin < distSqToMax) {
                max += 0.5f * (min - max);
                distSqToMax = (p - Cubic(a, b, c, d, max)).sqrMagnitude;
            }
            else {
                min += 0.5f * (max - min);
                distSqToMin = (p - Cubic(a, b, c, d, min)).sqrMagnitude;
            }
        }
        return new Vector2(0.5f * (min + max), 0.5f * distSqToMin + 0.5f * distSqToMax);
    }

    // 24 Gauss-Legendre quadrature abcissae and weights, mapped to 0..1
    static readonly Vector2[] GaussLegendre24 = {
        new Vector2(0.4679715535686972f, 0.06396909767337608f),
        new Vector2(0.5320284464313028f, 0.06396909767337608f),
        new Vector2(0.40444056626319186f, 0.06291872817341415f),
        new Vector2(0.5955594337368082f, 0.06291872817341415f),
        new Vector2(0.3424786601519183f, 0.060835236463901696f),
        new Vector2(0.6575213398480817f, 0.060835236463901696f),
        new Vector2(0.28310324618697746f, 0.0577528340268628f),
        new Vector2(0.7168967538130225f, 0.0577528340268628f),
        new Vector2(0.22728926430558022f, 0.05372213505798282f),
        new Vector2(0.7727107356944198f, 0.05372213505798282f),
        new Vector2(0.17595317403151223f, 0.04880932605205694f),
        new Vector2(0.8240468259684878f, 0.04880932605205694f),
        new Vector2(0.12993790421072282f, 0.04309508076597664f),
        new Vector2(0.8700620957892772f, 0.04309508076597664f),
        new Vector2(0.08999900701304853f, 0.03667324070554015f),
        new Vector2(0.9100009929869515f, 0.03667324070554015f),
        new Vector2(0.056792236497799464f, 0.02964929245771839f),
        new Vector2(0.9432077635022005f, 0.02964929245771839f),
        new Vector2(0.0308627239986336f, 0.022138719408709904f),
        new Vector2(0.9691372760013663f, 0.022138719408709904f),
        new Vector2(0.012635722014345263f, 0.014265694314466832f),
        new Vector2(0.9873642779856547f, 0.014265694314466832f),
        new Vector2(0.0024063900014893447f, 0.0061706148999936f),
        new Vector2(0.9975936099985107f, 0.0061706148999936f),
    };

    public static float CubicArcLength(Vector3 a, Vector3 b, Vector3 c, Vector3 d, float t = 1f) {
        var sum = 0f;

        foreach (var strip in GaussLegendre24)
            sum += strip.y * t * CubicDerivative(a, b, c, d, strip.x * t).magnitude;

        return sum;
    }

    public static float CubicArcLengthFast(Vector3 a, Vector3 b, Vector3 c, Vector3 d, float t = 1f) {
        return
            0.1739274225687269f * t * CubicDerivative(a, b, c, d, 0.06943184420297371f * t).magnitude +
            0.32607257743127305f * t * CubicDerivative(a, b, c, d, 0.33000947820757187f * t).magnitude +
            0.32607257743127305f * t * CubicDerivative(a, b, c, d, 0.6699905217924281f * t).magnitude +
            0.1739274225687269f * t * CubicDerivative(a, b, c, d, 0.9305681557970262f * t).magnitude;
    }

    public static float CubicParameterAtArcLength(Vector3 a, Vector3 b, Vector3 c, Vector3 d, float s, int iterations = 10) {
        float t = s / CubicArcLength(a, b, c, d), lower = 0f, upper = 1f;

        while (iterations-- > 0) {
            var error = CubicArcLength(a, b, c, d, t) - s;

            if (Mathf.Abs(error) <= Mathf.Epsilon)
                return t;

            var tCandidate = t - error / CubicDerivative(a, b, c, d, t).magnitude;
            if (error > 0f) {
                upper = t;
                t = tCandidate > lower ? tCandidate : 0.5f * (lower + upper);
            }
            else {
                lower = t;
                t = tCandidate < upper ? tCandidate : 0.5f * (lower + upper);
            }
        }
        return t;
    }

    public static void CubicHull(Vector3 a, Vector3 b, Vector3 c, Vector3 d, float t, List<Vector3> hull) {
        var ab = Vector3.LerpUnclamped(a, b, t);
        var bc = Vector3.LerpUnclamped(b, c, t);
        var cd = Vector3.LerpUnclamped(c, d, t);
        var abbc = Vector3.LerpUnclamped(ab, bc, t);
        var bccd = Vector3.LerpUnclamped(bc, cd, t);
        var abbcbccd = Vector3.LerpUnclamped(abbc, bccd, t);

        hull.Add(a);
        hull.Add(b);
        hull.Add(c);
        hull.Add(d);
        hull.Add(ab);
        hull.Add(bc);
        hull.Add(cd);
        hull.Add(abbc);
        hull.Add(bccd);
        hull.Add(abbcbccd);
    }

    public static void CubicSplit(Vector3 a, Vector3 b, Vector3 c, Vector3 d, float t, List<Vector3> points) {
        var ab = Vector3.LerpUnclamped(a, b, t);
        var bc = Vector3.LerpUnclamped(b, c, t);
        var cd = Vector3.LerpUnclamped(c, d, t);
        var abbc = Vector3.LerpUnclamped(ab, bc, t);
        var bccd = Vector3.LerpUnclamped(bc, cd, t);
        var abbcbccd = Vector3.LerpUnclamped(abbc, bccd, t);

        points.Add(a);
        points.Add(ab);
        points.Add(abbc);
        points.Add(abbcbccd);
        points.Add(abbcbccd);
        points.Add(bccd);
        points.Add(cd);
        points.Add(d);
    }
}