aada 0.1.1

//! 单条三次贝塞尔曲线的离散化：把曲线转成 (x, y) 折线点表，供时间采样查询。
//!
//! 三种采样方式输出**同构**的点表（首尾都精确为端点 P0/P3，x 非降）：
//! - [`CurveSampling::Differenced`]：参数 τ 等间隔差分采样（spine 官方算法，点数 = segs+1）
//! - [`CurveSampling::Adaptive`]：自适应细分，折线到曲线的最大几何误差 ≤ eps
//! - [`CurveSampling::TimeUniform`]：时间轴 x 等间隔，值用牛顿迭代精确求取
//!
//! 三者的消费端一致：构建期 [`prebuild_curve`] 生成点表，运行期 [`sample_curve`] 二分查表线性插值。



/// 二维点：x = 时间坐标，y = 值坐标
pub type Pt2 = [f32; 2];

/// 三次贝塞尔曲线（4 个控制点）
#[derive(Debug, Clone, Copy)]
pub struct Bezier {
    pub p0: Pt2,
    pub p1: Pt2,
    pub p2: Pt2,
    pub p3: Pt2,
}

impl Bezier {
    /// 参数求值（τ ∈ `[0,1]`），幂基三次多项式
    pub fn eval(&self, tau: f32) -> Pt2 {
        let t = tau.clamp(0.0, 1.0);
        let u = 1.0 - t;
        let a = u * u * u;
        let b = 3.0 * u * u * t;
        let c = 3.0 * u * t * t;
        let d = t * t * t;
        [
            a * self.p0[0] + b * self.p1[0] + c * self.p2[0] + d * self.p3[0],
            a * self.p0[1] + b * self.p1[1] + c * self.p2[1] + d * self.p3[1],
        ]
    }

    /// 解 x(τ) = target_x 求 τ（时间单调时唯一），牛顿迭代 + 二分回退
    fn solve_tau(&self, target_x: f32) -> f32 {
        let mut tau =
            ((target_x - self.p0[0]) / (self.p3[0] - self.p0[0])).clamp(0.0, 1.0);
        for _ in 0..8 {
            let x = self.x(tau);
            let dx = self.dx(tau);
            if dx.abs() < 1e-12 {
                break;
            }
            tau = (tau - (x - target_x) / dx).clamp(0.0, 1.0);
            if (x - target_x).abs() < 1e-5 {
                break;
            }
        }
        if (self.x(tau) - target_x).abs() > 1e-4 {
            // 牛顿未收敛（x 分量非单调等）时二分回退
            let (mut lo, mut hi) = (0.0f32, 1.0f32);
            for _ in 0..64 {
                let mid = 0.5 * (lo + hi);
                if self.x(mid) < target_x {
                    lo = mid;
                } else {
                    hi = mid;
                }
            }
            tau = 0.5 * (lo + hi);
        }
        tau
    }

    /// 给定时间 x 精确求值（解 x(τ)=x 后取 y，时间单调区间内为上升分支）
    pub fn value_at_time(&self, x: f32) -> f32 {
        let tau = self.solve_tau(x);
        self.eval(tau)[1]
    }

    /// 反求时间 x 对应的曲线参数 τ，取官方查询语义（上升分支首个交点）：
    /// 线性扫描 τ 找 x(τ) 首次 ≥ x 的段，再二分收敛；全程未达则返回 τ=1。
    /// 回绕曲线在 x=t1 处有两个分支，此函数保证取上升分支的 τ，而非 τ=1 的终点。
    pub fn solve_tau_upper(&self, target_x: f32) -> f32 {
        if self.x(0.0) >= target_x {
            return 0.0;
        }
        const SCAN: usize = 200;
        let mut prev = 0.0f32;
        for k in 1..=SCAN {
            let tau = k as f32 / SCAN as f32;
            if self.x(tau) >= target_x {
                let (mut lo, mut hi) = (prev, tau);
                for _ in 0..80 {
                    let mid = 0.5 * (lo + hi);
                    if self.x(mid) < target_x {
                        lo = mid;
                    } else {
                        hi = mid;
                    }
                }
                return 0.5 * (lo + hi);
            }
            prev = tau;
        }
        1.0
    }

    fn x(&self, tau: f32) -> f32 {
        let u = 1.0 - tau;
        u * u * u * self.p0[0]
            + 3.0 * u * u * tau * self.p1[0]
            + 3.0 * u * tau * tau * self.p2[0]
            + tau * tau * tau * self.p3[0]
    }

    fn dx(&self, tau: f32) -> f32 {
        let u = 1.0 - tau;
        3.0 * u * u * (self.p1[0] - self.p0[0])
            + 6.0 * u * tau * (self.p2[0] - self.p1[0])
            + 3.0 * tau * tau * (self.p3[0] - self.p2[0])
    }
}

/// 曲线离散化方式
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum CurveSampling {
    /// 参数 τ 等间隔差分采样（spine 官方算法）：segs 段，输出 segs+1 个点（含端点）
    Differenced { segs: usize },
    /// 自适应细分：折线到曲线的最大几何误差 ≤ eps
    Adaptive { eps: f32 },
    /// 时间轴均匀采样：x 等间隔 segs 段，输出 segs+1 个点（值精确求取）
    TimeUniform { segs: usize },
}

/// 把一条曲线按指定方式离散成折线点表（首尾为端点，x 非降）
pub fn prebuild_curve(curve: &Bezier, method: CurveSampling) -> Vec<Pt2> {
    match method {
        CurveSampling::Differenced { segs } => differenced(curve, segs.max(1)),
        CurveSampling::Adaptive { eps } => adaptive(curve, eps),
        CurveSampling::TimeUniform { segs } => time_uniform(curve, segs.max(1)),
    }
}

/// spine 官方差分采样：参数 τ 等间隔（步长 h = τ_max/segs），用差分递推（纯加法）生成点。
/// τ_max = 反求 t1 对应的上升分支 τ（回绕曲线不再采到 τ=1 的 p3 终点分支），
/// 最后点 = 曲线在该 τ 处的采样信息（eval(τ_max)）。
fn differenced(curve: &Bezier, segs: usize) -> Vec<Pt2> {
    let tau_max = curve.solve_tau_upper(curve.p3[0]);
    let h = tau_max / segs as f32;
    let h2 = h * h;
    let h3 = h2 * h;
    let init = |p0: f32, p1: f32, p2: f32, p3: f32| {
        let a = p3 - 3.0 * p2 + 3.0 * p1 - p0;
        let b = 3.0 * (p2 - 2.0 * p1 + p0);
        let c = 3.0 * (p1 - p0);
        let d1 = a * h3 + b * h2 + c * h;
        let d2 = 6.0 * a * h3 + 2.0 * b * h2;
        let d3 = 6.0 * a * h3;
        (d1, d2, d3)
    };
    let (mut d1x, mut d2x, d3x) = init(curve.p0[0], curve.p1[0], curve.p2[0], curve.p3[0]);
    let (mut d1y, mut d2y, d3y) = init(curve.p0[1], curve.p1[1], curve.p2[1], curve.p3[1]);
    let mut out = Vec::with_capacity(segs + 1);
    let (mut x, mut y) = (curve.p0[0], curve.p0[1]);
    out.push([x, y]);
    for _ in 0..segs {
        x += d1x;
        y += d1y;
        d1x += d2x;
        d1y += d2y;
        d2x += d3x;
        d2y += d3y;
        out.push([x, y]);
    }
    *out.last_mut().unwrap() = curve.eval(tau_max);
    out
}

/// 自适应细分：不够平就 de Casteljau 在 τ=0.5 处切成两段递归
///
/// 平直判据用**沿时间轴的 y 偏差**：控制点与弦在同一 x 处的 y 差 < eps。
/// 它直接约束"时间→值"分段线性插值的函数误差（比 2D 垂距更贴合动画查询；
/// 陡峭段垂距天然小，但 y 偏差能正确触发细分）
fn adaptive(curve: &Bezier, eps: f32) -> Vec<Pt2> {
    const MAX_DEPTH: usize = 24;
    fn mid(a: Pt2, b: Pt2) -> Pt2 {
        [(a[0] + b[0]) * 0.5, (a[1] + b[1]) * 0.5]
    }
    fn lerp_pt(a: Pt2, b: Pt2, t: f32) -> Pt2 {
        [a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t]
    }
    fn subdivide(p0: Pt2, p1: Pt2, p2: Pt2, p3: Pt2, eps: f32, depth: usize, out: &mut Vec<Pt2>) {
        let span_x = p3[0] - p0[0];
        let flat = if span_x.abs() <= 1e-12 {
            // x 方向退化（垂直段），按垂距近似判据
            let d = |p: Pt2, a: Pt2, b: Pt2| {
                let (abx, aby) = (b[0] - a[0], b[1] - a[1]);
                let len2 = abx * abx + aby * aby;
                if len2 <= 1e-20 {
                    return ((p[0] - a[0]).powi(2) + (p[1] - a[1]).powi(2)).sqrt();
                }
                let t = ((p[0] - a[0]) * abx + (p[1] - a[1]) * aby) / len2;
                let (qx, qy) = (a[0] + t * abx, a[1] + t * aby);
                ((p[0] - qx).powi(2) + (p[1] - qy).powi(2)).sqrt()
            };
            d(p1, p0, p3) <= eps && d(p2, p0, p3) <= eps
        } else {
            let chord_y = |x: f32| p0[1] + (x - p0[0]) / span_x * (p3[1] - p0[1]);
            (p1[1] - chord_y(p1[0])).abs() <= eps && (p2[1] - chord_y(p2[0])).abs() <= eps
        };
        if flat || depth >= MAX_DEPTH {
            return;
        }
        let l1 = mid(p0, p1);
        let h = mid(p1, p2);
        let r1 = mid(p2, p3);
        let l2 = mid(l1, h);
        let r2 = mid(h, r1);
        let m = mid(l2, r2);
        subdivide(p0, l1, l2, m, eps, depth + 1, out);
        out.push(m);
        subdivide(m, r2, r1, p3, eps, depth + 1, out);
    }
    let mut out = Vec::new();
    // 反求 t1 对应的上升分支 τ，de Casteljau 切出左段 [0, τ_max]（回绕曲线不采 τ=1 的 p3）
    let tau_max = curve.solve_tau_upper(curve.p3[0]);
    let l1 = lerp_pt(curve.p0, curve.p1, tau_max);
    let h = lerp_pt(curve.p1, curve.p2, tau_max);
    let r1 = lerp_pt(curve.p2, curve.p3, tau_max);
    let l2 = lerp_pt(l1, h, tau_max);
    let r2 = lerp_pt(h, r1, tau_max);
    let m = lerp_pt(l2, r2, tau_max);
    let (q0, q1, q2, q3) = (curve.p0, l1, l2, m);
    out.push(q0);
    subdivide(q0, q1, q2, q3, eps, 0, &mut out);
    out.push(q3);
    out
}

/// 时间轴均匀采样：x 从 p0.x 到 p3.x 等间隔，每个 x 反求 τ（官方上升分支语义）再取 y，
/// 最后点 = 反求 t1 对应的 τ 后的曲线点（回绕曲线取上升分支，不强制 p3 终点分支）
fn time_uniform(curve: &Bezier, segs: usize) -> Vec<Pt2> {
    let (x0, x1) = (curve.p0[0], curve.p3[0]);
    let span = x1 - x0;
    let mut out = Vec::with_capacity(segs + 1);
    for i in 0..=segs {
        let x = if span == 0.0 {
            x0
        } else {
            x0 + span * (i as f32 / segs as f32)
        };
        let tau = curve.solve_tau_upper(x);
        let pt = curve.eval(tau);
        out.push([x, pt[1]]);
    }
    out[0] = curve.p0;
    out
}

/// 在折线点表上按 x 二分定位 + 线性插值采样
pub fn sample_curve(points: &[Pt2], x: f32) -> f32 {
    let Some(first) = points.first() else {
        return 0.0;
    };
    if x <= first[0] {
        return first[1];
    }
    let last = points.last().unwrap();
    if x >= last[0] {
        return last[1];
    }
    let idx = points.partition_point(|p| p[0] < x).clamp(1, points.len() - 1);
    let a = points[idx - 1];
    let b = points[idx];
    let span = b[0] - a[0];
    if span <= 0.0 {
        return b[1];
    }
    a[1] + (x - a[0]) / span * (b[1] - a[1])
}

/// 输出折线点表为 CSV 文本（表头 x,y，小数点后 6 位）
pub fn points_to_csv(points: &[Pt2]) -> String {
    let mut s = String::from("x,y\n");
    for p in points {
        s.push_str(&format!("{:.6},{:.6}\n", p[0], p[1]));
    }
    s
}

/// 把多个折线点表合并为一个对齐的宽表 CSV（表头 x{name},y{name}，行数 = 最长的点数，不足留空）
pub fn points_to_csv_columns(sets: &[(&str, Vec<Pt2>)]) -> String {
    let n = sets.iter().map(|(_, p)| p.len()).max().unwrap_or(0);
    let mut headers: Vec<String> = Vec::new();
    for (name, _) in sets {
        headers.push(format!("x{name}"));
        headers.push(format!("y{name}"));
    }
    let mut s = String::new();
    s.push_str(&headers.join(","));
    s.push('\n');
    for r in 0..n {
        let mut cells: Vec<String> = Vec::new();
        for (_, pts) in sets {
            match pts.get(r) {
                Some(pt) => {
                    cells.push(format!("{:.6}", pt[0]));
                    cells.push(format!("{:.6}", pt[1]));
                }
                None => {
                    cells.push(String::new());
                    cells.push(String::new());
                }
            }
        }
        s.push_str(&cells.join(","));
        s.push('\n');
    }
    s
}

#[cfg(test)]
mod tests {
    use super::*;

    fn test_curve() -> Bezier {
        Bezier {
            p0: [0.0, 10.0],
            p1: [0.3, 30.0],
            p2: [0.7, 70.0],
            p3: [1.0, 100.0],
        }
    }

    fn check_shape(points: &[Pt2]) {
        assert_eq!(points[0], [0.0, 10.0], "起点必须是 P0");
        assert_eq!(*points.last().unwrap(), [1.0, 100.0], "终点必须是 P3");
        for w in points.windows(2) {
            assert!(w[1][0] >= w[0][0], "x 必须非降: {points:?}");
        }
    }

    #[test]
    fn differenced_generates_segs_plus_one() {
        let pts = prebuild_curve(&test_curve(), CurveSampling::Differenced { segs: 10 });
        assert_eq!(pts.len(), 11);
        check_shape(&pts);
    }

    #[test]
    fn adaptive_respects_eps() {
        let curve = test_curve();
        let loose = prebuild_curve(&curve, CurveSampling::Adaptive { eps: 5.0 });
        let tight = prebuild_curve(&curve, CurveSampling::Adaptive { eps: 0.01 });
        assert!(tight.len() > loose.len(), "容差越小点数越多");
        check_shape(&tight);
    }

    #[test]
    fn time_uniform_x_is_exact() {
        let pts = prebuild_curve(&test_curve(), CurveSampling::TimeUniform { segs: 4 });
        assert_eq!(pts.len(), 5);
        check_shape(&pts);
        assert_eq!(pts[0][0], 0.0);
        assert_eq!(pts[2][0], 0.5);
        assert_eq!(pts[4][0], 1.0);
    }

    #[test]
    fn sampling_interpolates_between_points() {
        let pts = prebuild_curve(&test_curve(), CurveSampling::TimeUniform { segs: 4 });
        let v = sample_curve(&pts, 0.375);
        let expect = 0.5 * (pts[1][1] + pts[2][1]);
        assert!((v - expect).abs() < 1e-4);
    }

    #[test]
    fn three_methods_agree_within_tolerance() {
        let curve = test_curve();
        let exact = |x: f32| {
            let tau = curve.solve_tau(x);
            curve.eval(tau)[1]
        };
        let diff = prebuild_curve(&curve, CurveSampling::Differenced { segs: 256 });
        let adapt = prebuild_curve(&curve, CurveSampling::Adaptive { eps: 1e-3 });
        let uni = prebuild_curve(&curve, CurveSampling::TimeUniform { segs: 256 });
        for x in [0.1, 0.25, 0.5, 0.75, 0.9] {
            let e = exact(x);
            for (name, pts) in [("diff", &diff), ("adapt", &adapt), ("uni", &uni)] {
                let v = sample_curve(pts, x);
                assert!(
                    (v - e).abs() < 5e-3,
                    "{name} x={x} got={v} exact={e}"
                );
            }
        }
    }

    #[test]
    fn adaptive_geometry_error_within_eps() {
        let curve = test_curve();
        let eps = 0.05;
        let pts = prebuild_curve(&curve, CurveSampling::Adaptive { eps });
        let dist = |p: Pt2, a: Pt2, b: Pt2| {
            let (abx, aby) = (b[0] - a[0], b[1] - a[1]);
            let len2 = abx * abx + aby * aby;
            if len2 <= 1e-20 {
                return ((p[0] - a[0]).powi(2) + (p[1] - a[1]).powi(2)).sqrt();
            }
            let t = ((p[0] - a[0]) * abx + (p[1] - a[1]) * aby) / len2;
            let (qx, qy) = (a[0] + t * abx, a[1] + t * aby);
            ((p[0] - qx).powi(2) + (p[1] - qy).powi(2)).sqrt()
        };
        let (x0, x1) = (pts[0][0], pts.last().unwrap()[0]);
        let mut max_dev = 0.0f32;
        for i in 0..200 {
            let p = curve.eval(i as f32 / 199.0);
            if p[0] < x0 || p[0] > x1 {
                continue;
            }
            let idx = pts.partition_point(|q| q[0] < p[0]).clamp(1, pts.len() - 1);
            max_dev = max_dev.max(dist(p, pts[idx - 1], pts[idx]));
        }
        assert!(max_dev <= eps, "折线到曲线最大垂距 {max_dev} > eps {eps}");
    }

    #[test]
    fn csv_dump_all_methods() {
        let curve = test_curve();
        let mut out = String::new();
        for (name, m) in [
            ("differenced_segs10", CurveSampling::Differenced { segs: 10 }),
            ("adaptive_eps0.1", CurveSampling::Adaptive { eps: 0.1 }),
            ("time_uniform_segs10", CurveSampling::TimeUniform { segs: 10 }),
        ] {
            let pts = prebuild_curve(&curve, m);
            out.push_str(&format!("== {name} ({} points) ==\n", pts.len()));
            out.push_str(&points_to_csv(&pts));
        }
        eprintln!("{out}");
    }

    /// 三种方式 + 50 点时间均匀采样，合并到同一个 CSV 写盘（对比用）
    #[test]
    fn dump_three_methods_combined_csv() {
        let curve = test_curve();
        let sets = [
            (
                "diff_segs10",
                prebuild_curve(&curve, CurveSampling::Differenced { segs: 10 }),
            ),
            (
                "adapt_eps0.1",
                prebuild_curve(&curve, CurveSampling::Adaptive { eps: 0.1 }),
            ),
            (
                "uni_segs10",
                prebuild_curve(&curve, CurveSampling::TimeUniform { segs: 10 }),
            ),
            (
                "uni_segs50",
                prebuild_curve(&curve, CurveSampling::TimeUniform { segs: 50 }),
            ),
        ];
        let csv = points_to_csv_columns(&sets);
        // uni_segs50 共 51 点，行数 = 表头 1 + 51
        assert_eq!(csv.lines().count(), 52);
        let dir = std::path::Path::new(env!("CARGO_MANIFEST_DIR"))
            .join("..")
            .join("target");
        std::fs::create_dir_all(&dir).unwrap();
        let path = dir.join("curve_sampling_compare.csv");
        if let Err(e) = std::fs::write(&path, &csv) {
            panic!("write {} failed: {e}", path.display());
        }
        eprintln!("written: {}\n{csv}", path.display());
    }

    /// 画线段（DDA 步进）
    fn draw_line(rgba: &mut [u8], w: usize, h: usize, x0: f32, y0: f32, x1: f32, y1: f32, c: [u8; 4]) {
        let steps = (((x1 - x0).abs()).max((y1 - y0).abs()) as usize).max(1);
        for i in 0..=steps {
            let t = i as f32 / steps as f32;
            let x = (x0 + (x1 - x0) * t).round() as isize;
            let y = (y0 + (y1 - y0) * t).round() as isize;
            if x >= 0 && x < w as isize && y >= 0 && y < h as isize {
                let idx = (y as usize * w + x as usize) * 4;
                rgba[idx] = c[0];
                rgba[idx + 1] = c[1];
                rgba[idx + 2] = c[2];
                rgba[idx + 3] = c[3];
            }
        }
    }

    /// 画实心圆点（半径 r）
    fn draw_dot(rgba: &mut [u8], w: usize, h: usize, cx: f32, cy: f32, r: f32, c: [u8; 4]) {
        let r = r as isize;
        for dy in -r..=r {
            for dx in -r..=r {
                if dx * dx + dy * dy <= r * r {
                    let x = cx.round() as isize + dx;
                    let y = cy.round() as isize + dy;
                    if x >= 0 && x < w as isize && y >= 0 && y < h as isize {
                        let idx = (y as usize * w + x as usize) * 4;
                        rgba[idx] = c[0];
                        rgba[idx + 1] = c[1];
                        rgba[idx + 2] = c[2];
                        rgba[idx + 3] = c[3];
                    }
                }
            }
        }
    }

    /// 把四条曲线（差分/自适应/时间均匀10段/时间均匀50段）绘制到一张 BMP（RGBA，白底网格）
    /// 前三条把采样点画成明显圆点；50 点那条过密只画折线
    #[test]
    fn draw_four_curves_to_bmp() {
        let (w, h) = (3840usize, 2160usize);
        let mut rgba = vec![255u8; w * h * 4];
        let (xmin, xmax) = (0.0f32, 1.0f32);
        let (ymin, ymax) = (10.0f32, 100.0f32);
        let margin = 24.0f32;
        let to_px = |x: f32, y: f32| {
            let px = margin + (x - xmin) / (xmax - xmin) * (w as f32 - 2.0 * margin);
            let py = margin + (ymax - y) / (ymax - ymin) * (h as f32 - 2.0 * margin);
            (px, py)
        };
        // 网格（浅灰）
        let grid: [u8; 4] = [230, 230, 230, 255];
        for i in 0..=4 {
            let x = xmin + (xmax - xmin) * i as f32 / 4.0;
            let (px, _) = to_px(x, ymin);
            draw_line(&mut rgba, w, h, px, margin, px, h as f32 - margin, grid);
            let y = ymin + (ymax - ymin) * i as f32 / 4.0;
            let (_, py) = to_px(xmin, y);
            draw_line(&mut rgba, w, h, margin, py, w as f32 - margin, py, grid);
        }
        // 四条曲线
        let curve = test_curve();
        let sets = [
            (
                "diff_segs10",
                prebuild_curve(&curve, CurveSampling::Differenced { segs: 10 }),
                [255, 0, 0, 255],
            ),
            (
                "adapt_eps0.1",
                prebuild_curve(&curve, CurveSampling::Adaptive { eps: 0.1 }),
                [0, 200, 0, 255],
            ),
            (
                "uni_segs10",
                prebuild_curve(&curve, CurveSampling::TimeUniform { segs: 10 }),
                [0, 0, 255, 255],
            ),
            (
                "uni_segs50",
                prebuild_curve(&curve, CurveSampling::TimeUniform { segs: 50 }),
                [255, 128, 0, 255],
            ),
        ];
        for (name, pts, color) in &sets {
            for seg in pts.windows(2) {
                let (x0, y0) = to_px(seg[0][0], seg[0][1]);
                let (x1, y1) = to_px(seg[1][0], seg[1][1]);
                draw_line(&mut rgba, w, h, x0, y0, x1, y1, *color);
            }
            // 前三条把采样点画成明显圆点；50 点那条过密只画折线
            if *name != "uni_segs50" {
                for p in pts {
                    let (px, py) = to_px(p[0], p[1]);
                    draw_dot(&mut rgba, w, h, px, py, 14.0, *color);
                }
            }
        }
        let bmp = crate::utils::write_bmp_rgba(w as u32, h as u32, &rgba);
        let dir = std::path::Path::new(env!("CARGO_MANIFEST_DIR"))
            .join("..")
            .join("target");
        std::fs::create_dir_all(&dir).unwrap();
        let path = dir.join("four_curves.bmp");
        std::fs::write(&path, &bmp).unwrap();
        eprintln!("written: {}", path.display());
    }

    /// 手工构造的 10 条贝塞尔曲线，覆盖常见形状：ease-in/out、overshoot、回绕（x 非单调）、
    /// 陡升缓降、平坦段、拱形、大值域等（非任何游戏资产数据）
    fn aya_curves() -> [(&'static str, Bezier); 10] {
        [
            ("ease-in", Bezier { p0: [0.0, 0.0], p1: [0.25, 0.0], p2: [0.6, 1.0], p3: [1.0, 1.0] }),
            ("wraparound", Bezier { p0: [0.667, 5.0], p1: [1.59, 10.0], p2: [1.75, 5.0], p3: [1.33, 0.0] }),
            ("ease-out", Bezier { p0: [0.0, 1.0], p1: [0.3, 1.0], p2: [0.75, 0.0], p3: [1.0, 0.0] }),
            ("overshoot", Bezier { p0: [0.0, 2.0], p1: [0.3, 4.0], p2: [0.7, -2.0], p3: [1.0, 1.0] }),
            ("ease-in-out", Bezier { p0: [0.0, 0.0], p1: [0.42, 0.0], p2: [0.58, 1.0], p3: [1.0, 1.0] }),
            ("steep-rise", Bezier { p0: [0.0, -4.0], p1: [0.25, -8.0], p2: [0.6, 0.0], p3: [1.0, 4.0] }),
            ("near-linear", Bezier { p0: [0.0, 0.0], p1: [0.4, 1.0], p2: [0.6, 2.0], p3: [1.0, 3.0] }),
            ("plateau", Bezier { p0: [0.0, 2.0], p1: [0.3, 2.0], p2: [0.5, 0.5], p3: [1.0, 0.0] }),
            ("bounce", Bezier { p0: [0.0, 1.0], p1: [0.2, 3.0], p2: [0.8, 3.0], p3: [1.0, 1.0] }),
            ("long-fall", Bezier { p0: [0.0, 10.0], p1: [0.3, 3.0], p2: [0.7, 1.0], p3: [1.0, 0.0] }),
        ]
    }

    /// 点表 pts 相对参考点表 ref_pts，在 [x0,x1] 区间内 19 个均匀查询点的平均绝对误差
    fn avg_err(pts: &[Pt2], ref_pts: &[Pt2], x0: f32, x1: f32) -> f32 {
        let n = 20;
        let mut sum = 0.0f32;
        for i in 1..n {
            let x = x0 + (x1 - x0) * i as f32 / n as f32;
            sum += (sample_curve(pts, x) - sample_curve(ref_pts, x)).abs();
        }
        sum / (n - 1) as f32
    }

    /// 同上，但取 19 个查询点的最大绝对误差
    fn max_err(pts: &[Pt2], ref_pts: &[Pt2], x0: f32, x1: f32) -> f32 {
        let n = 20;
        let mut m = 0.0f32;
        for i in 1..n {
            let x = x0 + (x1 - x0) * i as f32 / n as f32;
            m = m.max((sample_curve(pts, x) - sample_curve(ref_pts, x)).abs());
        }
        m
    }

    /// 手工构造的 10 条贝塞尔曲线，
    /// 每条渲染到一个子图（网格 + 控制多边形 + 高精度曲线 + 控制点），拼成 4K 大图
    /// 自适应采样按"该曲线 diff10 相对 uni50 的平均误差"配置 eps
    #[test]
    fn draw_aya_bezier_curves() {
        let curves = aya_curves();
        let (cols, rows) = (5usize, 2usize);
        let (w, h) = (3840usize, 2160usize);
        let (cell_w, cell_h) = (w / cols, h / rows);
        let mut rgba = vec![255u8; w * h * 4];

        for (i, (_, curve)) in curves.iter().enumerate() {
            let (cx0, cy0) = (i % cols * cell_w, i / cols * cell_h);
            // bbox（含控制点）
            let mut xmin = f32::MAX;
            let mut xmax = f32::MIN;
            let mut ymin = f32::MAX;
            let mut ymax = f32::MIN;
            let cp = [curve.p0, curve.p1, curve.p2, curve.p3];
            for p in cp {
                xmin = xmin.min(p[0]);
                xmax = xmax.max(p[0]);
                ymin = ymin.min(p[1]);
                ymax = ymax.max(p[1]);
            }
            let (sx, sy) = ((xmax - xmin).max(1e-6), (ymax - ymin).max(1e-6));
            let m = (cell_w as f32 * 0.06).min(cell_h as f32 * 0.1);
            let to_px = |x: f32, y: f32| {
                let px = cx0 as f32 + m + (x - xmin) / sx * (cell_w as f32 - 2.0 * m);
                let py = cy0 as f32 + m + (ymax - y) / sy * (cell_h as f32 - 2.0 * m);
                (px, py)
            };
            // 子图边框
            let border: [u8; 4] = [80, 80, 80, 255];
            draw_line(&mut rgba, w, h, cx0 as f32, cy0 as f32, (cx0 + cell_w) as f32, cy0 as f32, border);
            draw_line(&mut rgba, w, h, cx0 as f32, cy0 as f32, cx0 as f32, (cy0 + cell_h) as f32, border);
            draw_line(&mut rgba, w, h, (cx0 + cell_w) as f32, (cy0 + cell_h) as f32, (cx0 + cell_w) as f32, cy0 as f32, border);
            draw_line(&mut rgba, w, h, (cx0 + cell_w) as f32, (cy0 + cell_h) as f32, cx0 as f32, (cy0 + cell_h) as f32, border);
            // 网格（浅灰）
            let grid: [u8; 4] = [235, 235, 235, 255];
            for k in 1..4 {
                let x = xmin + sx * k as f32 / 4.0;
                let (px, _) = to_px(x, ymin);
                draw_line(&mut rgba, w, h, px, cy0 as f32 + m, px, cy0 as f32 + cell_h as f32 - m, grid);
                let y = ymin + sy * k as f32 / 4.0;
                let (_, py) = to_px(xmin, y);
                draw_line(&mut rgba, w, h, cx0 as f32 + m, py, cx0 as f32 + cell_w as f32 - m, py, grid);
            }
            // 控制多边形（浅灰）
            let ctrl_line: [u8; 4] = [200, 200, 200, 255];
            for seg in cp.windows(2) {
                let (a, b) = (to_px(seg[0][0], seg[0][1]), to_px(seg[1][0], seg[1][1]));
                draw_line(&mut rgba, w, h, a.0, a.1, b.0, b.1, ctrl_line);
            }
            // 基准曲线（深灰，200 个 τ 点，作为采样线参照）
            let base_col: [u8; 4] = [90, 90, 90, 255];
            let mut prev = to_px(curve.eval(0.0)[0], curve.eval(0.0)[1]);
            for k in 1..=200 {
                let pt = curve.eval(k as f32 / 200.0);
                let cur = to_px(pt[0], pt[1]);
                draw_line(&mut rgba, w, h, prev.0, prev.1, cur.0, cur.1, base_col);
                prev = cur;
            }
            // 精确参考求值（端点语义：回绕曲线取上升分支，而非终点 value2）
            let ref_y = |t: f32| {
                let tau = curve.solve_tau_upper(t);
                curve.eval(tau)[1]
            };
            // 三种采样折线：diff 红 / adapt 绿 / uni 蓝（按参数采样，x 非单调也能绘制）
            // 自适应 eps 按"该曲线 diff9 相对 uni50 的最大误差"配置
            let (x0, x1) = (
                curve.p0[0].min(curve.p3[0]),
                curve.p0[0].max(curve.p3[0]),
            );
            let ref_pts = prebuild_curve(curve, CurveSampling::TimeUniform { segs: 50 });
            let e_max = max_err(
                &prebuild_curve(curve, CurveSampling::Differenced { segs: 9 }),
                &ref_pts,
                x0,
                x1,
            );
            let eps = e_max.max(sy * 1e-3);
            // 三种采样点表直接来自 prebuild：端点已用反求 τ 取上升分支，x 单调覆盖 [t0,t1]
            let diff_pts = prebuild_curve(curve, CurveSampling::Differenced { segs: 10 });
            let adapt_pts = prebuild_curve(curve, CurveSampling::Adaptive { eps });
            let uni_pts = prebuild_curve(curve, CurveSampling::TimeUniform { segs: 10 });
            let samplings = [
                ("diff", diff_pts, [220, 20, 20, 255]),
                ("adapt", adapt_pts, [20, 170, 20, 255]),
                ("uni", uni_pts, [20, 20, 220, 255]),
            ];
            for (_, pts, color) in &samplings {
                for seg in pts.windows(2) {
                    let (a, b) = (
                        to_px(seg[0][0], seg[0][1]),
                        to_px(seg[1][0], seg[1][1]),
                    );
                    draw_line(&mut rgba, w, h, a.0, a.1, b.0, b.1, *color);
                }
                // 采样点（折线顶点）也画成同色圆点
                for p in pts {
                    let (px, py) = to_px(p[0], p[1]);
                    draw_dot(&mut rgba, w, h, px, py, 6.0, *color);
                }
            }
            // 精确参考逐点采样线（最上层，黄色）：对 [t0,t1] 内 100 个 time 各反求一次上升分支 τ 连线
            let off_col: [u8; 4] = [255, 255, 0, 255];
            let (t0, t1) = (
                curve.p0[0].min(curve.p3[0]),
                curve.p0[0].max(curve.p3[0]),
            );
            let n = 100;
            let mut prev_off = to_px(t0, ref_y(t0));
            for k in 1..=n {
                let t = t0 + (t1 - t0) * k as f32 / n as f32;
                let cur = to_px(t, ref_y(t));
                // 黑色描边，保证亮黄在白色背景上清晰可见
                for (dx, dy) in [(1.0, 0.0), (-1.0, 0.0), (0.0, 1.0), (0.0, -1.0)] {
                    draw_line(
                        &mut rgba,
                        w,
                        h,
                        prev_off.0 + dx,
                        prev_off.1 + dy,
                        cur.0 + dx,
                        cur.1 + dy,
                        [0, 0, 0, 255],
                    );
                }
                draw_line(&mut rgba, w, h, prev_off.0, prev_off.1, cur.0, cur.1, off_col);
                prev_off = cur;
            }
            // 控制点：P1/P2 红点，P0/P3 蓝点
            let (p1, p2, p0, p3) = (
                to_px(cp[1][0], cp[1][1]),
                to_px(cp[2][0], cp[2][1]),
                to_px(cp[0][0], cp[0][1]),
                to_px(cp[3][0], cp[3][1]),
            );
            draw_dot(&mut rgba, w, h, p1.0, p1.1, 9.0, [230, 30, 30, 255]);
            draw_dot(&mut rgba, w, h, p2.0, p2.1, 9.0, [230, 30, 30, 255]);
            draw_dot(&mut rgba, w, h, p0.0, p0.1, 6.0, [30, 30, 230, 255]);
            draw_dot(&mut rgba, w, h, p3.0, p3.1, 6.0, [30, 30, 230, 255]);
        }

        let bmp = crate::utils::write_bmp_rgba(w as u32, h as u32, &rgba);
        let dir = std::path::Path::new(env!("CARGO_MANIFEST_DIR"))
            .join("..")
            .join("target");
        std::fs::create_dir_all(&dir).unwrap();
        let path = dir.join("aya_bezier_curves.bmp");
        std::fs::write(&path, &bmp).unwrap();
        eprintln!("written: {}", path.display());
        // 统计黄色（官方线）像素，验证确实画上
        let yellow = rgba
            .chunks_exact(4)
            .filter(|p| p[0] > 200 && p[1] > 180 && p[2] < 100)
            .count();
        eprintln!("official-line yellow pixels: {yellow}");
    }

    /// 诊断：回绕曲线（x 非单调）三种采样点表在"采样成数值"时的行为
    #[test]
    fn diagnose_non_monotonic_sampling() {
        let curve = Bezier {
            p0: [0.667, 5.0],
            p1: [1.59, 10.0],
            p2: [1.75, 5.0],
            p3: [1.33, 0.0],
        };
        let diff = prebuild_curve(&curve, CurveSampling::Differenced { segs: 10 });
        let uni = prebuild_curve(&curve, CurveSampling::TimeUniform { segs: 10 });
        let adapt = prebuild_curve(&curve, CurveSampling::Adaptive { eps: 0.05 });
        eprintln!("diff x: {:?}", diff.iter().map(|p| p[0]).collect::<Vec<_>>());
        eprintln!("uni  x: {:?}", uni.iter().map(|p| p[0]).collect::<Vec<_>>());
        eprintln!("adapt x: {:?}", adapt.iter().map(|p| p[0]).collect::<Vec<_>>());
        let exact = {
            let tau = curve.solve_tau(1.0);
            curve.eval(tau)[1]
        };
        eprintln!("exact(1.0)       = {exact}");
        eprintln!("sample diff(1.0) = {}", sample_curve(&diff, 1.0));
        eprintln!("sample uni(1.0)  = {}", sample_curve(&uni, 1.0));
        eprintln!("sample adapt(1.0)= {}", sample_curve(&adapt, 1.0));
        // 采样点表在线性扫描语义下取"首个 x >= time"的分支，与官方一致
        let non_mono = diff.windows(2).any(|w| w[1][0] < w[0][0]);
        eprintln!("diff 点表 x 非单调 = {non_mono}（sample_curve 二分结果未定义）");
    }

    /// 诊断：回绕曲线（x 非单调）精确参考线逐点反求，找异常值/跳变
    #[test]
    fn diagnose_reference_line_termination() {
        let curve = Bezier {
            p0: [0.667, 5.0],
            p1: [1.59, 10.0],
            p2: [1.75, 5.0],
            p3: [1.33, 0.0],
        };
        let ref_y = |t: f32| {
            let tau = curve.solve_tau_upper(t);
            curve.eval(tau)[1]
        };
        let (t0, t1) = (0.667f32, 1.33f32);
        let n = 100;
        let mut prev_y = ref_y(t0);
        let mut jumps = 0;
        let mut bad = 0;
        for k in 1..=n {
            let t = t0 + (t1 - t0) * k as f32 / n as f32;
            let y = ref_y(t);
            if y.is_nan() || y.is_infinite() {
                bad += 1;
                eprintln!("t={t}: y={y}");
            } else if (y - prev_y).abs() > 20.0 {
                jumps += 1;
                eprintln!("t={t}: y={y} (jump from {prev_y})");
            }
            prev_y = y;
        }
        eprintln!("bad={bad} jumps={jumps} (100 samples)");
        // 打印关键区间
        for t in [1.0, 1.1, 1.2, 1.25, 1.3, 1.32, 1.33] {
            eprintln!("t={t} -> y={}", ref_y(t));
        }
    }

    /// 回绕曲线（x 非单调）三种 prebuild 点表：端点用反求 τ 取上升分支、
    /// x 单调覆盖 [t0,t1]、内部查询与精确值一致
    #[test]
    fn prebuild_wraparound_endpoint_on_upper_branch() {
        let curve = Bezier {
            p0: [0.667, 5.0],
            p1: [1.59, 10.0],
            p2: [1.75, 5.0],
            p3: [1.33, 0.0],
        };
        // 精确参考求值（回绕曲线在 t=time2 处取上升分支，而非终点 value2）
        let ref_y = |t: f32| {
            let tau = curve.solve_tau_upper(t);
            curve.eval(tau)[1]
        };
        let (t0, t1) = (curve.p0[0], curve.p3[0]);
        let methods = [
            ("diff10", CurveSampling::Differenced { segs: 10 }),
            ("adapt0.01", CurveSampling::Adaptive { eps: 0.01 }),
            ("uni10", CurveSampling::TimeUniform { segs: 10 }),
        ];
        for (name, method) in methods {
            let pts = prebuild_curve(&curve, method);
            for w in pts.windows(2) {
                assert!(w[1][0] >= w[0][0], "{name} 点表 x 必须非降");
            }
            assert!(pts[0][0] <= t0 + 1e-4, "{name} 覆盖左端点");
            assert!(pts.last().unwrap()[0] >= t1 - 1e-4, "{name} 覆盖右端点");
            // 端点 = 反求 τ 后的上升分支值（远离 value2=0 的跳变）
            let end_y = pts.last().unwrap()[1];
            assert!(
                (end_y - ref_y(t1)).abs() < 0.15,
                "{name} 终点 y={end_y} 应接近精确上升分支值 {}",
                ref_y(t1)
            );
            // 内部查询与精确上升分支值一致
            for t in [0.8, 1.0, 1.2] {
                let v = sample_curve(&pts, t);
                let exact = curve.value_at_time(t);
                assert!((v - exact).abs() < 2e-2, "{name} t={t} v={v} exact={exact}");
            }
            eprintln!(
                "{name}: {} pts, end=({:.4},{:.4}), x {:?}",
                pts.len(),
                pts.last().unwrap()[0],
                pts.last().unwrap()[1],
                pts.iter().map(|p| (p[0] * 100.0).round() / 100.0).collect::<Vec<_>>()
            );
        }
    }

    /// 计算 10 段采样（diff/uni/adapt1%）相对 50 段时间均匀采样的平均误差，
    /// 并按每条曲线 diff10 的误差配置自适应采样，对比点数
    #[test]
    fn report_10pts_error_and_adaptive() {
        let curves = aya_curves();
        let mut sum_diff = 0.0f32;
        let mut sum_uni = 0.0f32;
        let mut sum_adapt = 0.0f32;
        for (name, curve) in &curves {
            let (x0, x1) = (
                curve.p0[0].min(curve.p3[0]),
                curve.p0[0].max(curve.p3[0]),
            );
            let (mut ymin, mut ymax) = (f32::MAX, f32::MIN);
            for p in [curve.p0, curve.p1, curve.p2, curve.p3] {
                ymin = ymin.min(p[1]);
                ymax = ymax.max(p[1]);
            }
            let ref_pts = prebuild_curve(curve, CurveSampling::TimeUniform { segs: 50 });
            let diff9 = prebuild_curve(curve, CurveSampling::Differenced { segs: 9 });
            let uni = prebuild_curve(curve, CurveSampling::TimeUniform { segs: 10 });
            let adapt = prebuild_curve(curve, CurveSampling::Adaptive { eps: (ymax - ymin) * 0.01 });
            let e_avg = avg_err(&diff9, &ref_pts, x0, x1);
            let e_max = max_err(&diff9, &ref_pts, x0, x1);
            let e_uni = avg_err(&uni, &ref_pts, x0, x1);
            let e_adapt = avg_err(&adapt, &ref_pts, x0, x1);
            sum_diff += e_avg;
            sum_uni += e_uni;
            sum_adapt += e_adapt;
            // 自适应 eps = 该曲线 diff9 的最大误差
            let eps = e_max.max((ymax - ymin) * 1e-3);
            let adapt_by_err = prebuild_curve(curve, CurveSampling::Adaptive { eps });
            eprintln!(
                "{name}\n  diff9 avg={e_avg:.4} max={e_max:.4}  uni10 avg={e_uni:.4}  adapt1% avg={e_adapt:.4}  adapt(eps={eps:.4}) -> {} pts",
                adapt_by_err.len()
            );
        }
        let n = curves.len() as f32;
        eprintln!(
            "== avg over {n} curves: diff10={:.4}  uni10={:.4}  adapt1%={:.4} ==",
            sum_diff / n,
            sum_uni / n,
            sum_adapt / n
        );
    }
}