vello_common/src/flatten_simd.rs

// Copyright 2025 the Vello Authors
// SPDX-License-Identifier: Apache-2.0 OR MIT

//! This is a temporary module that contains a SIMD version of flattening of cubic curves, as
//! well as some code that was copied from kurbo, which is needed to reimplement the
//! full `flatten` method.

#[cfg(not(feature = "std"))]
use crate::kurbo::common::FloatFuncs as _;
use crate::kurbo::{CubicBez, Line, ParamCurve, ParamCurveNearest, PathEl, Point, QuadBez};
use crate::{
    flatten::{SQRT_TOL, TOL, TOL_2},
    geometry::RectU16,
    kurbo::Affine,
    tile::Tile,
};
use alloc::vec::Vec;
use bytemuck::{Pod, Zeroable};
use fearless_simd::*;

/// The element of a path made of lines.
///
/// Each subpath must start with a `MoveTo`. Closing of subpaths is not supported, and subpaths are
/// not closed implicitly when a new subpath (with `MoveTo`) is started. It is expected that closed
/// subpaths are watertight in the sense that the last `LineTo` matches exactly with the first
/// `MoveTo`.
///
/// This intentionally allows for non-watertight subpaths, as, e.g., lines that are fully outside
/// of the viewport do not need to be drawn.
///
/// See [`PathEl`] for a more general-purpose path element type.
pub(crate) enum LinePathEl {
    MoveTo(Point),
    LineTo(Point),
}

// Unlike kurbo, which takes a closure with a callback for outputting the lines, we use a trait
// instead. The reason is that this way the callback can be inlined, which is not possible with
// a closure and turned out to have a noticeable overhead.
pub(crate) trait Callback {
    fn callback(&mut self, el: LinePathEl);
}

/// See the docs for the kurbo implementation of flattening:
/// <https://docs.rs/kurbo/latest/kurbo/fn.flatten.html>
///
/// This version works using a similar approach but using f32x4/f32x8 SIMD instead.
#[inline(always)]
pub(crate) fn flatten<S: Simd>(
    simd: S,
    path: impl IntoIterator<Item = PathEl>,
    // Note: we explicitly pass the `Affine` here instead of using `map` on
    // the iterator because it seems like the `map` function isn't always inlined
    // properly, even when annotating the closure with `#[inline(always)]`.
    // See https://github.com/linebender/vello/pull/1600.
    affine: Affine,
    callback: &mut impl Callback,
    flatten_ctx: &mut FlattenCtx,
    cull_bbox: RectU16,
) {
    flatten_ctx.flattened_cubics.clear();

    // For the culling performed here to be correct, the top y coordinate of the cull bbox must be
    // aligned to strip row boundaries. Consider the alternative: for example, a strip row starting
    // at y=8, a cull bbox starting at y=10, and (nearly) horizontal geometry at y=9 that does not
    // cross into the cull bbox and does not extend above y=8 (and note this is all in a y-down
    // coordinate space).
    //
    // The culling performed here would remove that geometry. However, as it does not extend above
    // the strip row, it does not add coarse winding to the strip, and does not produce a sparse
    // fill. Yet, if this is the top part of, say, a rectangle that is partially visible, the
    // geometry is necessary for tiling to emit intermediate tiles that are necessary for the
    // bottom part of the strip to get filled.
    //
    // Therefore, we align `top` to strip row boundaries, such that all intermediate tiles are
    // produced.
    //
    // This is not necessary for the bottom. Consider a path where a segment extends just below the
    // cull bbox, but remains in the same strip row. The path is closed, so if anything is to be
    // rendered at all, there will be other geometry above that edge of the cull bbox. If that
    // geometry extends above the row, there will be coarse winding for a sparse fill. If not, it
    // there will be geometry to generate the intermediate tiles.
    let left = cull_bbox.x0 as f64;
    let top = ((cull_bbox.y0 / Tile::HEIGHT) * Tile::HEIGHT) as f64;
    let right = cull_bbox.x1 as f64;
    let bottom = cull_bbox.y1 as f64;

    let mut path = path.into_iter();
    let Some(first_el) = path.next() else {
        return;
    };
    let first_el = affine * first_el;
    let PathEl::MoveTo(start_pt) = first_el else {
        debug_assert!(
            matches!(first_el, PathEl::MoveTo(_)),
            "Non-empty paths must begin with `PathEl::MoveTo`, got {first_el:?}"
        );
        return;
    };

    let mut start_pt = start_pt;
    let mut last_pt = start_pt;
    callback.callback(LinePathEl::MoveTo(start_pt));

    for el in path {
        match affine * el {
            PathEl::MoveTo(p) => {
                if last_pt != start_pt {
                    callback.callback(LinePathEl::LineTo(start_pt));
                }
                last_pt = p;
                start_pt = p;
                callback.callback(LinePathEl::MoveTo(p));
            }
            PathEl::LineTo(p) => {
                last_pt = p;
                callback.callback(LinePathEl::LineTo(p));
            }
            PathEl::QuadTo(p1, p2) => {
                let p0 = last_pt;
                let line = Line::new(p0, p2);
                // If the quadratic Bézier is fully to the right, top, or bottom of the culling
                // bbox, it does not impact pixel coverage or winding. We can ignore it. The
                // following checks that conservatively by checking whether the bounding box of the
                // Bézier's control points is fully outside the culling bbox.
                if [p0, p1, p2].into_iter().all(|p| p.x > right)
                    || [p0, p1, p2].into_iter().all(|p| p.y < top)
                    || [p0, p1, p2].into_iter().all(|p| p.y > bottom)
                {
                    callback.callback(LinePathEl::MoveTo(p2));
                }
                // The following checks two things. First, if the quadratic Bézier is fully to the
                // left of the culling bbox, it may affect pixel coverage and winding, but its
                // exact shape does not matter. It can be emitted as a line segment [p0, p2].
                //
                // Second, an upper bound on the shortest distance of any point on the quadratic
                // Bézier curve to the line segment [p0, p2] is 1/2 of the control-point-to-line-segment
                // distance.
                //
                // The derivation is similar to that for the cubic Bézier (see below). In
                // short:
                //
                // q(t) = B0(t) p0 + B1(t) p1 + B2(t) p2
                // dist(q(t), [p0, p1]) <= B1(t) dist(p1, [p0, p1])
                //                       = 2 (1-t)t dist(p1, [p0, p1]).
                //
                // The maximum occurs at t=1/2, hence
                // max(dist(q(t), [p0, p1] <= 1/2 dist(p1, [p0, p1])).
                //
                // The following takes the square to elide the square root of the Euclidean
                // distance.
                else if [p0, p1, p2].into_iter().all(|p| p.x < left)
                    || line.nearest(p1, 0.).distance_sq <= 4. * TOL_2
                {
                    callback.callback(LinePathEl::LineTo(p2));
                } else {
                    let q = QuadBez::new(p0, p1, p2);
                    let params = q.estimate_subdiv(SQRT_TOL);
                    let n = ((0.5 / SQRT_TOL * params.val).ceil() as usize).max(1);
                    let step = 1.0 / (n as f64);
                    for i in 1..n {
                        let u = (i as f64) * step;
                        let t = q.determine_subdiv_t(&params, u);
                        let p = q.eval(t);
                        callback.callback(LinePathEl::LineTo(p));
                    }
                    callback.callback(LinePathEl::LineTo(p2));
                }
                last_pt = p2;
            }
            PathEl::CurveTo(p1, p2, p3) => {
                let p0 = last_pt;
                let line = Line::new(p0, p3);
                // If the cubic Bézier is fully to the right, top, or bottom of the culling bbox,
                // it does not impact pixel coverage or winding. We can ignore it. The following
                // checks that conservatively by checking whether the bounding box of the Bézier's
                // control points is fully outside the culling bbox.
                if [p0, p1, p2, p3].into_iter().all(|p| p.x > right)
                    || [p0, p1, p2, p3].into_iter().all(|p| p.y < top)
                    || [p0, p1, p2, p3].into_iter().all(|p| p.y > bottom)
                {
                    callback.callback(LinePathEl::MoveTo(p3));
                }
                // The following checks two things. First, if the cubic Bézier is fully to the left
                // of the culling bbox, it may affect pixel coverage and winding, but its exact
                // shape does not matter. It can be emitted as a line segment [p0, p3].
                //
                // Second, an upper bound on the shortest distance of any point on the cubic Bézier
                // curve to the line segment [p0, p3] is 3/4 of the maximum of the
                // control-point-to-line-segment distances.
                //
                // With Bernstein weights Bi(t), we have
                // c(t) = B0(t) p0 + B1(t) p1 + B2(t) p2 + B3(t) p3
                // with t from 0 to 1 (inclusive).
                //
                // Through convexivity of the Euclidean distance function and the line segment,
                // we have
                // dist(c(t), [p0, p3]) <= B1(t) dist(p1, [p0, p3]) + B2(t) dist(p2, [p0, p3])
                //                      <= (B1(t) + B2(t)) max(dist(p1, [p0, p3]), dist(p2, [p0, p3]))
                //                       = 3 ((1-t)t^2 + (1-t)^2t) max(dist(p1, [p0, p3]), dist(p2, [p0, p3])).
                //
                // The inner polynomial has its maximum of 1/4 at t=1/2, hence
                // max(dist(c(t), [p0, p3])) <= 3/4 max(dist(p1, [p0, p3]), dist(p2, [p0, p3])).
                //
                // The following takes the square to elide the square root of the Euclidean
                // distance.
                else if [p0, p1, p2, p3].into_iter().all(|p| p.x < left)
                    || f64::max(
                        line.nearest(p1, 0.).distance_sq,
                        line.nearest(p2, 0.).distance_sq,
                    ) <= 16. / 9. * TOL_2
                {
                    callback.callback(LinePathEl::LineTo(p3));
                } else {
                    let c = CubicBez::new(p0, p1, p2, p3);
                    let max = flatten_cubic_simd(simd, c, flatten_ctx);

                    for p in &flatten_ctx.flattened_cubics[1..max] {
                        callback.callback(LinePathEl::LineTo(Point::new(p.x as f64, p.y as f64)));
                    }
                }
                last_pt = p3;
            }
            PathEl::ClosePath => {
                if last_pt != start_pt {
                    callback.callback(LinePathEl::LineTo(start_pt));

                    // Kurbo says: "If `quad_to` [or another drawing op] is called immediately
                    // after `close_path` then the current subpath starts at the initial point of
                    // the previous subpath."
                    //
                    // Hence, we set `last_pt` back to the just-closed subpath's `start_pt`.
                    last_pt = start_pt;
                }
            }
        }
    }

    if last_pt != start_pt {
        callback.callback(LinePathEl::LineTo(start_pt));
    }
}

// The below methods are copied from kurbo and needed to implement flattening of normal quad curves.

/// An approximation to $\int (1 + 4x^2) ^ -0.25 dx$
///
/// This is used for flattening curves.
fn approx_parabola_integral(x: f64) -> f64 {
    const D: f64 = 0.67;
    x / (1.0 - D + (D.powi(4) + 0.25 * x * x).sqrt().sqrt())
}

/// An approximation to the inverse parabola integral.
fn approx_parabola_inv_integral(x: f64) -> f64 {
    const B: f64 = 0.39;
    x * (1.0 - B + (B * B + 0.25 * x * x).sqrt())
}

impl FlattenParamsExt for QuadBez {
    #[inline(always)]
    fn estimate_subdiv(&self, sqrt_tol: f64) -> FlattenParams {
        // Determine transformation to $y = x^2$ parabola.
        let d01 = self.p1 - self.p0;
        let d12 = self.p2 - self.p1;
        let dd = d01 - d12;
        let cross = (self.p2 - self.p0).cross(dd);
        let x0 = d01.dot(dd) * cross.recip();
        let x2 = d12.dot(dd) * cross.recip();
        let scale = (cross / (dd.hypot() * (x2 - x0))).abs();

        // Compute number of subdivisions needed.
        let a0 = approx_parabola_integral(x0);
        let a2 = approx_parabola_integral(x2);
        let val = if scale.is_finite() {
            let da = (a2 - a0).abs();
            let sqrt_scale = scale.sqrt();
            if x0.signum() == x2.signum() {
                da * sqrt_scale
            } else {
                // Handle cusp case (segment contains curvature maximum)
                let xmin = sqrt_tol / sqrt_scale;
                sqrt_tol * da / approx_parabola_integral(xmin)
            }
        } else {
            0.0
        };
        let u0 = approx_parabola_inv_integral(a0);
        let u2 = approx_parabola_inv_integral(a2);
        let uscale = (u2 - u0).recip();
        FlattenParams {
            a0,
            a2,
            u0,
            uscale,
            val,
        }
    }

    #[inline(always)]
    fn determine_subdiv_t(&self, params: &FlattenParams, x: f64) -> f64 {
        let a = params.a0 + (params.a2 - params.a0) * x;
        let u = approx_parabola_inv_integral(a);
        (u - params.u0) * params.uscale
    }
}

trait FlattenParamsExt {
    fn estimate_subdiv(&self, sqrt_tol: f64) -> FlattenParams;
    fn determine_subdiv_t(&self, params: &FlattenParams, x: f64) -> f64;
}

// Everything below is a SIMD implementation of flattening of cubic curves.
// It's a combination of https://gist.github.com/raphlinus/5f4e9feb85fd79bafc72da744571ec0e
// and https://gist.github.com/raphlinus/44e114fef2fd33b889383a60ced0129b.

// TODO(laurenz): Perhaps we should get rid of this in the future and work directly with f32,
// as it's the only reason we have to pull in proc_macros via the `derive` feature of bytemuck.
#[derive(Clone, Copy, Debug, Default, Zeroable, Pod)]
#[repr(C)]
struct Point32 {
    x: f32,
    y: f32,
}

struct FlattenParams {
    a0: f64,
    a2: f64,
    u0: f64,
    uscale: f64,
    /// The number of `subdivisions * 2 * sqrt_tol`.
    val: f64,
}

/// This limit was chosen based on the pre-existing GitHub gist.
/// This limit should not be hit in normal operation, but _might_ be hit for very large
/// transforms.
const MAX_QUADS: usize = 16;

/// The context needed for flattening curves.
#[derive(Default, Debug)]
pub struct FlattenCtx {
    // The +4 is to encourage alignment; might be better to be explicit
    even_pts: [Point32; MAX_QUADS + 4],
    odd_pts: [Point32; MAX_QUADS],
    a0: [f32; MAX_QUADS],
    da: [f32; MAX_QUADS],
    u0: [f32; MAX_QUADS],
    uscale: [f32; MAX_QUADS],
    val: [f32; MAX_QUADS],
    n_quads: usize,
    /// Reusable buffer for flattened cubic points.
    flattened_cubics: Vec<Point32>,
}

#[inline(always)]
fn is_finite_simd<S: Simd>(x: f32x4<S>) -> mask32x4<S> {
    let simd = x.simd;

    let x_abs = x.abs();
    let reinterpreted = u32x4::from_bytes(x_abs.to_bytes());
    simd.simd_lt_u32x4(reinterpreted, u32x4::splat(simd, 0x7f80_0000))
}

/// An approximation to $\int (1 + 4x^2) ^ -0.25 dx$
///
/// This is used for flattening curves.
///
/// SIMD version of [`approx_parabola_integral`].
#[inline(always)]
fn approx_parabola_integral_simd<S: Simd, F: SimdFloat<S, Element = f32>>(x: F) -> F {
    let simd = x.witness();

    const D: f32 = 0.67;
    const D_POWI_4: f32 = 0.201_511_2;

    let temp = F::splat(x.witness(), 0.25)
        .mul_add(x * x, F::splat(simd, D_POWI_4))
        .sqrt()
        .sqrt();
    let denom = temp + (1. - D);
    x / denom
}

/// An approximation to the inverse parabola integral.
///
/// SIMD version of [`approx_parabola_inv_integral`].
#[inline(always)]
fn approx_parabola_inv_integral_simd<S: Simd>(x: f32x8<S>) -> f32x8<S> {
    let simd = x.simd;

    const B: f32 = 0.39;
    const ONE_MINUS_B: f32 = 1.0 - B;

    let temp = f32x8::splat(simd, 0.25).mul_add(x * x, B * B).sqrt();
    let factor = f32x8::splat(simd, ONE_MINUS_B) + temp;
    x * factor
}

#[inline(always)]
fn pt_splat_simd<S: Simd>(simd: S, pt: Point32) -> f32x8<S> {
    let p_f64: f64 = bytemuck::cast(pt);
    f64x4::splat(simd, p_f64).bitcast()
}

#[inline(always)]
fn eval_cubics_simd<S: Simd>(simd: S, c: &CubicBez, n: usize, result: &mut FlattenCtx) {
    result.n_quads = n;
    let dt = 0.5 / n as f32;

    // TODO(laurenz): Perhaps we can SIMDify this better
    let p0p1 = f32x4::from_slice(
        simd,
        &[c.p0.x as f32, c.p0.y as f32, c.p1.x as f32, c.p1.y as f32],
    );
    let p2p3 = f32x4::from_slice(
        simd,
        &[c.p2.x as f32, c.p2.y as f32, c.p3.x as f32, c.p3.y as f32],
    );

    let split_single = |input: f32x4<S>| {
        let t1: f64x2<S> = input.bitcast();
        let p0 = simd.zip_low_f64x2(t1, t1);
        let p1 = simd.zip_high_f64x2(t1, t1);

        let p0: f32x4<S> = p0.bitcast();
        let p1: f32x4<S> = p1.bitcast();

        (f32x8::block_splat(p0), f32x8::block_splat(p1))
    };

    let (p0_128, p1_128) = split_single(p0p1);
    let (p2_128, p3_128) = split_single(p2p3);

    // Use Horner's method to evaluate p(t) as ((a*t + b)*t + c)*t + d. This allows hoisting
    // coefficients out of the loop, and then evaluating as three sequential FMAs. Estrin's method
    // would expose more ILP, but the increase in work on such small polynomials (cubic here and
    // quad below) makes it a net slowdown.
    let coeff_a = (p1_128 - p2_128).mul_add(3.0, p3_128 - p0_128);
    let coeff_b = p1_128.mul_add(-2.0, p0_128 + p2_128) * 3.0;
    let coeff_c = (p1_128 - p0_128) * 3.0;
    let coeff_d = p0_128;

    let iota = f32x8::from_slice(simd, &[0.0, 0.0, 2.0, 2.0, 1.0, 1.0, 3.0, 3.0]);
    let step = iota * dt;
    let mut t = step;
    let t_inc = f32x8::splat(simd, 4.0 * dt);

    let even_pts: &mut [f32] = bytemuck::cast_slice_mut(&mut result.even_pts);
    let odd_pts: &mut [f32] = bytemuck::cast_slice_mut(&mut result.odd_pts);

    for i in 0..n.div_ceil(2) {
        let evaluated = (coeff_a.mul_add(t, coeff_b))
            .mul_add(t, coeff_c)
            .mul_add(t, coeff_d);

        let (low, high) = simd.split_f32x8(evaluated);

        low.store_slice(&mut even_pts[i * 4..][..4]);
        high.store_slice(&mut odd_pts[i * 4..][..4]);

        t += t_inc;
    }

    p3_128.store_slice(&mut even_pts[n * 2..][..8]);
}

#[inline(always)]
fn estimate_subdiv_simd<S: Simd>(simd: S, sqrt_tol: f32, ctx: &mut FlattenCtx) {
    let n = ctx.n_quads;

    let even_pts: &mut [f32] = bytemuck::cast_slice_mut(&mut ctx.even_pts);
    let odd_pts: &mut [f32] = bytemuck::cast_slice_mut(&mut ctx.odd_pts);

    for i in 0..n.div_ceil(4) {
        let p0 = f32x8::from_slice(simd, &even_pts[i * 8..][..8]);
        let p_onehalf = f32x8::from_slice(simd, &odd_pts[i * 8..][..8]);
        let p2 = f32x8::from_slice(simd, &even_pts[(i * 8 + 2)..][..8]);
        let x = p0 * -0.5;
        let x1 = p_onehalf.mul_add(2.0, x);
        let p1 = p2.mul_add(-0.5, x1);

        p1.store_slice(&mut odd_pts[(i * 8)..][..8]);

        let d01 = p1 - p0;
        let d12 = p2 - p1;
        let d01x = simd.unzip_low_f32x8(d01, d01);
        let d01y = simd.unzip_high_f32x8(d01, d01);
        let d12x = simd.unzip_low_f32x8(d12, d12);
        let d12y = simd.unzip_high_f32x8(d12, d12);
        let ddx = d01x - d12x;
        let ddy = d01y - d12y;
        let d02x = d01x + d12x;
        let d02y = d01y + d12y;
        // (d02x * ddy) - (d02y * ddx)
        let cross = ddx.mul_add(-d02y, d02x * ddy);

        let x0_x2_a = {
            let (d01x_low, _) = simd.split_f32x8(d01x);
            let (d12x_low, _) = simd.split_f32x8(d12x);

            simd.combine_f32x4(d12x_low, d01x_low) * ddx
        };
        let temp1 = {
            let (d12y_low, _) = simd.split_f32x8(d12y);
            let (d01y_low, _) = simd.split_f32x8(d01y);

            simd.combine_f32x4(d12y_low, d01y_low)
        };
        let x0_x2_num = temp1.mul_add(ddy, x0_x2_a);
        let x0_x2 = x0_x2_num / cross;
        let (ddx_low, _) = simd.split_f32x8(ddx);
        let (ddy_low, _) = simd.split_f32x8(ddy);
        let dd_hypot = ddy_low.mul_add(ddy_low, ddx_low * ddx_low).sqrt();
        let (x0, x2) = simd.split_f32x8(x0_x2);
        let scale_denom = dd_hypot * (x2 - x0);
        let (temp2, _) = simd.split_f32x8(cross);
        let scale = (temp2 / scale_denom).abs();
        let a0_a2 = approx_parabola_integral_simd(x0_x2);
        let (a0, a2) = simd.split_f32x8(a0_a2);
        let da = a2 - a0;
        let da_abs = da.abs();
        let sqrt_scale = scale.sqrt();
        let temp3 = simd.or_i32x4(x0.bitcast(), x2.bitcast());
        let mask = simd.simd_ge_i32x4(temp3, i32x4::splat(simd, 0));
        let noncusp = da_abs * sqrt_scale;
        // TODO: should we skip this if neither is a cusp? Maybe not worth branch prediction cost
        let xmin = sqrt_tol / sqrt_scale;
        let approxint = approx_parabola_integral_simd(xmin);
        let cusp = (sqrt_tol * da_abs) / approxint;
        let val_raw = simd.select_f32x4(mask, noncusp, cusp);
        let finite_mask = is_finite_simd(val_raw);
        let val = simd.select_f32x4(finite_mask, val_raw, f32x4::splat(simd, 0.0));
        let u0_u2 = approx_parabola_inv_integral_simd(a0_a2);
        let (u0, u2) = simd.split_f32x8(u0_u2);
        let uscale_a = u2 - u0;
        let uscale = 1.0 / uscale_a;

        a0.store_slice(&mut ctx.a0[i * 4..][..4]);
        da.store_slice(&mut ctx.da[i * 4..][..4]);
        u0.store_slice(&mut ctx.u0[i * 4..][..4]);
        uscale.store_slice(&mut ctx.uscale[i * 4..][..4]);
        val.store_slice(&mut ctx.val[i * 4..][..4]);
    }
}

#[inline(always)]
fn output_lines_simd<S: Simd>(
    simd: S,
    ctx: &mut FlattenCtx,
    i: usize,
    x0: f32,
    dx: f32,
    n: usize,
    start_idx: usize,
) {
    let p0 = pt_splat_simd(simd, ctx.even_pts[i]);
    let p1 = pt_splat_simd(simd, ctx.odd_pts[i]);
    let p2 = pt_splat_simd(simd, ctx.even_pts[i + 1]);

    const IOTA2: [f32; 8] = [0., 0., 1., 1., 2., 2., 3., 3.];
    let iota2 = f32x8::from_slice(simd, IOTA2.as_ref());
    let x = iota2.mul_add(dx, f32x8::splat(simd, x0));
    let da = f32x8::splat(simd, ctx.da[i]);
    let mut a = da.mul_add(x, f32x8::splat(simd, ctx.a0[i]));
    let a_inc = 4.0 * dx * da;
    let uscale = f32x8::splat(simd, ctx.uscale[i]);
    let u0 = f32x8::splat(simd, ctx.u0[i]);

    // See Horner comment above
    let coeff_a = p1.mul_add(-2.0, p0) + p2;
    let coeff_b = (p1 - p0) * 2.0;
    let coeff_c = p0;

    let out: &mut [f32] = bytemuck::cast_slice_mut(&mut ctx.flattened_cubics[start_idx..]);

    for j in 0..n.div_ceil(4) {
        let u = approx_parabola_inv_integral_simd(a);
        let t = (u - u0) * uscale;
        let p = coeff_a.mul_add(t, coeff_b).mul_add(t, coeff_c);
        p.store_slice(&mut out[j * 8..][..8]);
        a += a_inc;
    }
}

#[inline(always)]
fn flatten_cubic_simd<S: Simd>(simd: S, c: CubicBez, ctx: &mut FlattenCtx) -> usize {
    let n_quads = estimate_num_quads(c, TOL as f32);
    eval_cubics_simd(simd, &c, n_quads, ctx);
    let tol = (TOL as f32) * (1.0 - TO_QUAD_TOL);
    let sqrt_tol = tol.sqrt();
    estimate_subdiv_simd(simd, sqrt_tol, ctx);
    let sum: f32 = ctx.val[..n_quads].iter().sum();
    let n = ((0.5 * sum / sqrt_tol).ceil() as usize).max(1);
    let target_len = n + 4;
    if target_len > ctx.flattened_cubics.len() {
        ctx.flattened_cubics.resize(target_len, Point32::default());
    }

    let step = sum / (n as f32);
    let step_recip = 1.0 / step;
    let mut val_sum = 0.0;
    let mut last_n = 0;
    let mut x0base = 0.0;

    for i in 0..n_quads {
        let val = ctx.val[i];
        val_sum += val;
        let this_n = val_sum * step_recip;
        let this_n_next = 1.0 + this_n.floor();
        let dn = this_n_next as usize - last_n;
        if dn > 0 {
            let dx = step / val;
            let x0 = x0base * dx;
            output_lines_simd(simd, ctx, i, x0, dx, dn, last_n);
        }
        x0base = this_n_next - this_n;
        last_n = this_n_next as usize;
    }

    ctx.flattened_cubics[n] = ctx.even_pts[n_quads];

    n + 1
}

#[inline(always)]
fn estimate_num_quads(c: CubicBez, accuracy: f32) -> usize {
    let q_accuracy = (accuracy * TO_QUAD_TOL) as f64;
    let max_hypot2 = 432.0 * q_accuracy * q_accuracy;
    let p1x2 = c.p1.to_vec2() * 3.0 - c.p0.to_vec2();
    let p2x2 = c.p2.to_vec2() * 3.0 - c.p3.to_vec2();
    let err = (p2x2 - p1x2).hypot2();
    let err_div = err / max_hypot2;

    estimate(err_div)
}

const TO_QUAD_TOL: f32 = 0.1;

#[inline(always)]
fn estimate(err_div: f64) -> usize {
    // The original version of this method was:
    // let n_quads = (err_div.powf(1. / 6.0).ceil() as usize).max(1);
    // n_quads.min(MAX_QUADS)
    //
    // Note how we always round up and clamp to the range [1, max_quads]. Since we don't
    // care about the actual fractional value resulting from the powf call we can simply
    // compute this using a precomputed lookup table evaluating 1^6, 2^6, 3^6, etc. and simply
    // comparing if the value is less than or equal to each threshold.

    const LUT: [f64; MAX_QUADS] = [
        1.0, 64.0, 729.0, 4096.0, 15625.0, 46656.0, 117649.0, 262144.0, 531441.0, 1000000.0,
        1771561.0, 2985984.0, 4826809.0, 7529536.0, 11390625.0, 16777216.0,
    ];

    #[expect(clippy::needless_range_loop, reason = "better clarity")]
    for i in 0..MAX_QUADS {
        if err_div <= LUT[i] {
            return i + 1;
        }
    }

    MAX_QUADS
}

#[cfg(test)]
mod tests {
    use crate::flatten_simd::{MAX_QUADS, estimate};

    fn old_estimate(err_div: f64) -> usize {
        let n_quads = (err_div.powf(1. / 6.0).ceil() as usize).max(1);
        n_quads.min(MAX_QUADS)
    }

    // Test is disabled by default since it takes 10-20 seconds to run, even in release mode.
    #[test]
    #[ignore]
    fn accuracy() {
        for i in 0..u32::MAX {
            let num = f32::from_bits(i);

            if num.is_finite() {
                assert_eq!(old_estimate(num as f64), estimate(num as f64), "{num}");
            }
        }
    }
}