runmat_runtime/builtins/math/poly/

polyder.rs

//! MATLAB-compatible `polyder` builtin with GPU-aware semantics for RunMat.

use log::trace;
use num_complex::Complex64;
use runmat_builtins::{ComplexTensor, Tensor, Value};
use runmat_macros::runtime_builtin;

use crate::builtins::common::random_args::complex_tensor_into_value;
use crate::builtins::common::spec::{
    BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
    ProviderHook, ReductionNaN, ResidencyPolicy, ScalarType, ShapeRequirements,
};
use crate::builtins::common::{tensor, tensor::tensor_into_value};
use crate::builtins::math::poly::type_resolvers::polyder_type;
use crate::dispatcher;
use crate::{build_runtime_error, BuiltinResult, RuntimeError};

const EPS: f64 = 1.0e-12;
const BUILTIN_NAME: &str = "polyder";

#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::math::poly::polyder")]
pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
    name: "polyder",
    op_kind: GpuOpKind::Custom("polynomial-derivative"),
    supported_precisions: &[ScalarType::F32, ScalarType::F64],
    broadcast: BroadcastSemantics::None,
    provider_hooks: &[
        ProviderHook::Custom("polyder-single"),
        ProviderHook::Custom("polyder-product"),
        ProviderHook::Custom("polyder-quotient"),
    ],
    constant_strategy: ConstantStrategy::InlineLiteral,
    residency: ResidencyPolicy::NewHandle,
    nan_mode: ReductionNaN::Include,
    two_pass_threshold: None,
    workgroup_size: None,
    accepts_nan_mode: false,
    notes: "Runs on-device when providers expose polyder hooks; falls back to the host for complex coefficients or unsupported shapes.",
};

fn polyder_error(message: impl Into<String>) -> RuntimeError {
    build_runtime_error(message)
        .with_builtin(BUILTIN_NAME)
        .build()
}

#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::math::poly::polyder")]
pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
    name: "polyder",
    shape: ShapeRequirements::Any,
    constant_strategy: ConstantStrategy::InlineLiteral,
    elementwise: None,
    reduction: None,
    emits_nan: false,
    notes: "Symbolic operation on coefficient vectors; fusion bypasses this builtin.",
};

#[runtime_builtin(
    name = "polyder",
    category = "math/poly",
    summary = "Differentiate polynomials, products, and ratios with MATLAB-compatible coefficient vectors.",
    keywords = "polyder,polynomial,derivative,product,quotient",
    type_resolver(polyder_type),
    builtin_path = "crate::builtins::math::poly::polyder"
)]
async fn polyder_builtin(first: Value, rest: Vec<Value>) -> crate::BuiltinResult<Value> {
    if let Some(out_count) = crate::output_count::current_output_count() {
        if out_count <= 1 {
            let result = match rest.len() {
                0 => derivative_single(first).await,
                1 => derivative_product(first, rest.into_iter().next().unwrap()).await,
                _ => Err(polyder_error("polyder: too many input arguments")),
            }?;
            if out_count == 0 {
                return Ok(Value::OutputList(Vec::new()));
            }
            return Ok(Value::OutputList(vec![result]));
        }
        if rest.len() != 1 {
            return Err(polyder_error(
                "Not enough input arguments for quotient form.",
            ));
        }
        let eval = evaluate_quotient(first, rest.into_iter().next().unwrap()).await?;
        let outputs = vec![eval.numerator(), eval.denominator()];
        return Ok(crate::output_count::output_list_with_padding(
            out_count, outputs,
        ));
    }
    match rest.len() {
        0 => derivative_single(first).await,
        1 => derivative_product(first, rest.into_iter().next().unwrap()).await,
        _ => Err(polyder_error("polyder: too many input arguments")),
    }
}

async fn try_gpu_derivative_single(value: &Value) -> BuiltinResult<Option<Value>> {
    let Value::GpuTensor(handle) = value else {
        return Ok(None);
    };
    let Some(provider) = runmat_accelerate_api::provider() else {
        return Ok(None);
    };
    match provider.polyder_single(handle).await {
        Ok(out) => Ok(Some(Value::GpuTensor(out))),
        Err(err) => {
            trace!("polyder: provider polyder_single fallback: {err}");
            Ok(None)
        }
    }
}

async fn try_gpu_derivative_product(first: &Value, second: &Value) -> BuiltinResult<Option<Value>> {
    match (first, second) {
        (Value::GpuTensor(p), Value::GpuTensor(q)) => {
            let Some(provider) = runmat_accelerate_api::provider() else {
                return Ok(None);
            };
            match provider.polyder_product(p, q).await {
                Ok(out) => Ok(Some(Value::GpuTensor(out))),
                Err(err) => {
                    trace!("polyder: provider polyder_product fallback: {err}");
                    Ok(None)
                }
            }
        }
        _ => Ok(None),
    }
}

async fn try_gpu_quotient(u: &Value, v: &Value) -> BuiltinResult<Option<PolyderEval>> {
    match (u, v) {
        (Value::GpuTensor(uh), Value::GpuTensor(vh)) => {
            let Some(provider) = runmat_accelerate_api::provider() else {
                return Ok(None);
            };
            match provider.polyder_quotient(uh, vh).await {
                Ok(result) => Ok(Some(PolyderEval {
                    numerator: Value::GpuTensor(result.numerator),
                    denominator: Value::GpuTensor(result.denominator),
                })),
                Err(err) => {
                    trace!("polyder: provider polyder_quotient fallback: {err}");
                    Ok(None)
                }
            }
        }
        _ => Ok(None),
    }
}

/// Evaluate the quotient rule derivative `[num, den] = polyder(u, v)`.
pub async fn evaluate_quotient(u: Value, v: Value) -> BuiltinResult<PolyderEval> {
    if let Some(eval) = try_gpu_quotient(&u, &v).await? {
        return Ok(eval);
    }
    let u_poly = parse_polynomial("polyder", "U", u).await?;
    let v_poly = parse_polynomial("polyder", "V", v).await?;
    let numerator = quotient_numerator(&u_poly, &v_poly)?;
    let denominator = quotient_denominator(&v_poly)?;
    Ok(PolyderEval {
        numerator,
        denominator,
    })
}

/// Differentiated outputs for the quotient rule.
#[derive(Clone)]
pub struct PolyderEval {
    numerator: Value,
    denominator: Value,
}

impl PolyderEval {
    /// Numerator coefficients of the derivative `(u' * v - u * v')`.
    pub fn numerator(&self) -> Value {
        self.numerator.clone()
    }

    /// Denominator coefficients `v^2`.
    pub fn denominator(&self) -> Value {
        self.denominator.clone()
    }
}

pub async fn derivative_single(value: Value) -> BuiltinResult<Value> {
    if let Some(out) = try_gpu_derivative_single(&value).await? {
        return Ok(out);
    }
    let poly = parse_polynomial("polyder", "P", value).await?;
    differentiate_polynomial(&poly)
}

pub async fn derivative_product(first: Value, second: Value) -> BuiltinResult<Value> {
    if let Some(out) = try_gpu_derivative_product(&first, &second).await? {
        return Ok(out);
    }
    let p = parse_polynomial("polyder", "P", first).await?;
    let q = parse_polynomial("polyder", "A", second).await?;
    product_derivative(&p, &q)
}

fn quotient_numerator(u: &Polynomial, v: &Polynomial) -> BuiltinResult<Value> {
    let du = raw_derivative(&u.coeffs);
    let dv = raw_derivative(&v.coeffs);
    let term1 = poly_convolve(&du, &v.coeffs);
    let term2 = poly_convolve(&u.coeffs, &dv);
    let mut numerator = poly_sub(&term1, &term2);
    numerator = trim_leading_zeros(&numerator);
    coeffs_to_value(&numerator, u.orientation)
}

fn quotient_denominator(v: &Polynomial) -> BuiltinResult<Value> {
    let mut denominator = poly_convolve(&v.coeffs, &v.coeffs);
    denominator = trim_leading_zeros(&denominator);
    coeffs_to_value(&denominator, v.orientation)
}

fn differentiate_polynomial(poly: &Polynomial) -> BuiltinResult<Value> {
    let mut coeffs = raw_derivative(&poly.coeffs);
    coeffs = trim_leading_zeros(&coeffs);
    coeffs_to_value(&coeffs, poly.orientation)
}

fn product_derivative(p: &Polynomial, q: &Polynomial) -> BuiltinResult<Value> {
    let dp = raw_derivative(&p.coeffs);
    let dq = raw_derivative(&q.coeffs);
    let term1 = poly_convolve(&dp, &q.coeffs);
    let term2 = poly_convolve(&p.coeffs, &dq);
    let mut result = poly_add(&term1, &term2);
    result = trim_leading_zeros(&result);
    coeffs_to_value(&result, p.orientation)
}

fn raw_derivative(coeffs: &[Complex64]) -> Vec<Complex64> {
    if coeffs.len() <= 1 {
        return vec![Complex64::new(0.0, 0.0)];
    }
    let mut output = Vec::with_capacity(coeffs.len() - 1);
    let mut power = coeffs.len() - 1;
    for coeff in coeffs.iter().take(coeffs.len() - 1) {
        output.push(*coeff * (power as f64));
        power -= 1;
    }
    output
}

fn poly_convolve(a: &[Complex64], b: &[Complex64]) -> Vec<Complex64> {
    if a.is_empty() || b.is_empty() {
        return Vec::new();
    }
    let mut result = vec![Complex64::new(0.0, 0.0); a.len() + b.len() - 1];
    for (i, &ai) in a.iter().enumerate() {
        for (j, &bj) in b.iter().enumerate() {
            result[i + j] += ai * bj;
        }
    }
    result
}

fn poly_add(a: &[Complex64], b: &[Complex64]) -> Vec<Complex64> {
    let len = a.len().max(b.len());
    let mut result = vec![Complex64::new(0.0, 0.0); len];
    for (idx, &value) in a.iter().enumerate() {
        result[len - a.len() + idx] += value;
    }
    for (idx, &value) in b.iter().enumerate() {
        result[len - b.len() + idx] += value;
    }
    result
}

fn poly_sub(a: &[Complex64], b: &[Complex64]) -> Vec<Complex64> {
    let len = a.len().max(b.len());
    let mut result = vec![Complex64::new(0.0, 0.0); len];
    for (idx, &value) in a.iter().enumerate() {
        result[len - a.len() + idx] += value;
    }
    for (idx, &value) in b.iter().enumerate() {
        result[len - b.len() + idx] -= value;
    }
    result
}

fn trim_leading_zeros(coeffs: &[Complex64]) -> Vec<Complex64> {
    let mut first = None;
    for (idx, coeff) in coeffs.iter().enumerate() {
        if coeff.norm() > EPS {
            first = Some(idx);
            break;
        }
    }
    match first {
        Some(idx) => coeffs[idx..].to_vec(),
        None => vec![Complex64::new(0.0, 0.0)],
    }
}

fn coeffs_to_value(coeffs: &[Complex64], orientation: Orientation) -> BuiltinResult<Value> {
    if coeffs.iter().all(|c| c.im.abs() <= EPS) {
        let data: Vec<f64> = coeffs.iter().map(|c| c.re).collect();
        let shape = orientation.shape_for_len(data.len());
        let tensor =
            Tensor::new(data, shape).map_err(|e| polyder_error(format!("polyder: {e}")))?;
        Ok(tensor_into_value(tensor))
    } else {
        let data: Vec<(f64, f64)> = coeffs.iter().map(|c| (c.re, c.im)).collect();
        let shape = orientation.shape_for_len(data.len());
        let tensor =
            ComplexTensor::new(data, shape).map_err(|e| polyder_error(format!("polyder: {e}")))?;
        Ok(complex_tensor_into_value(tensor))
    }
}

async fn parse_polynomial(context: &str, label: &str, value: Value) -> BuiltinResult<Polynomial> {
    let gathered = dispatcher::gather_if_needed_async(&value).await?;
    let (coeffs, orientation) = match gathered {
        Value::Tensor(tensor) => {
            ensure_vector_shape(context, label, &tensor.shape)?;
            let orientation = orientation_from_shape(&tensor.shape);
            if tensor.data.is_empty() {
                (vec![Complex64::new(0.0, 0.0)], orientation)
            } else {
                (
                    tensor
                        .data
                        .into_iter()
                        .map(|re| Complex64::new(re, 0.0))
                        .collect(),
                    orientation,
                )
            }
        }
        Value::ComplexTensor(tensor) => {
            ensure_vector_shape(context, label, &tensor.shape)?;
            let orientation = orientation_from_shape(&tensor.shape);
            if tensor.data.is_empty() {
                (vec![Complex64::new(0.0, 0.0)], orientation)
            } else {
                (
                    tensor
                        .data
                        .into_iter()
                        .map(|(re, im)| Complex64::new(re, im))
                        .collect(),
                    orientation,
                )
            }
        }
        Value::LogicalArray(logical) => {
            let tensor = tensor::logical_to_tensor(&logical).map_err(polyder_error)?;
            ensure_vector_shape(context, label, &tensor.shape)?;
            let orientation = orientation_from_shape(&tensor.shape);
            if tensor.data.is_empty() {
                (vec![Complex64::new(0.0, 0.0)], orientation)
            } else {
                (
                    tensor
                        .data
                        .into_iter()
                        .map(|re| Complex64::new(re, 0.0))
                        .collect(),
                    orientation,
                )
            }
        }
        Value::Num(n) => (vec![Complex64::new(n, 0.0)], Orientation::Scalar),
        Value::Int(i) => (vec![Complex64::new(i.to_f64(), 0.0)], Orientation::Scalar),
        Value::Bool(b) => (
            vec![Complex64::new(if b { 1.0 } else { 0.0 }, 0.0)],
            Orientation::Scalar,
        ),
        Value::Complex(re, im) => (vec![Complex64::new(re, im)], Orientation::Scalar),
        other => {
            return Err(polyder_error(format!(
                "{context}: expected {label} to be a numeric vector, got {other:?}"
            )));
        }
    };

    Ok(Polynomial {
        coeffs,
        orientation,
    })
}

fn ensure_vector_shape(context: &str, label: &str, shape: &[usize]) -> BuiltinResult<()> {
    let non_unit = shape.iter().copied().filter(|&dim| dim > 1).count();
    if non_unit <= 1 {
        Ok(())
    } else {
        Err(polyder_error(format!(
            "{context}: {label} must be a vector of coefficients"
        )))
    }
}

#[derive(Clone, Copy)]
enum Orientation {
    Scalar,
    Row,
    Column,
}

impl Orientation {
    fn shape_for_len(self, len: usize) -> Vec<usize> {
        if len <= 1 {
            return vec![1, 1];
        }
        match self {
            Orientation::Scalar => vec![1, len],
            Orientation::Row => vec![1, len],
            Orientation::Column => vec![len, 1],
        }
    }
}

fn orientation_from_shape(shape: &[usize]) -> Orientation {
    for (idx, &dim) in shape.iter().enumerate() {
        if dim != 1 {
            return match idx {
                0 => Orientation::Column,
                1 => Orientation::Row,
                _ => Orientation::Column,
            };
        }
    }
    Orientation::Scalar
}

#[derive(Clone)]
struct Polynomial {
    coeffs: Vec<Complex64>,
    orientation: Orientation,
}

#[cfg(test)]
pub(crate) mod tests {
    use super::*;
    use crate::builtins::common::test_support;
    use futures::executor::block_on;
    use runmat_builtins::{IntValue, Tensor};

    fn assert_error_contains(err: crate::RuntimeError, needle: &str) {
        assert!(
            err.message().contains(needle),
            "expected error containing '{needle}', got '{}'",
            err.message()
        );
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn derivative_of_cubic_polynomial_is_correct() {
        let tensor = Tensor::new(vec![3.0, -2.0, 5.0, 7.0], vec![1, 4]).unwrap();
        let result = derivative_single(Value::Tensor(tensor)).expect("polyder");
        match result {
            Value::Tensor(t) => {
                assert_eq!(t.shape, vec![1, 3]);
                assert!(t
                    .data
                    .iter()
                    .zip([9.0, -4.0, 5.0])
                    .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
            }
            other => panic!("expected tensor result, got {other:?}"),
        }
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn derivative_of_product_matches_manual_rule() {
        let p = Tensor::new(vec![1.0, 0.0, -2.0], vec![1, 3]).unwrap();
        let a = Tensor::new(vec![1.0, 1.0], vec![1, 2]).unwrap();
        let result =
            derivative_product(Value::Tensor(p), Value::Tensor(a)).expect("polyder product");
        match result {
            Value::Tensor(t) => {
                assert_eq!(t.shape, vec![1, 3]);
                assert!(t
                    .data
                    .iter()
                    .zip([3.0, 2.0, -2.0])
                    .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
            }
            other => panic!("expected tensor result, got {other:?}"),
        }
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn quotient_rule_produces_expected_num_and_den() {
        let u = Tensor::new(vec![1.0, 0.0, -4.0], vec![1, 3]).unwrap();
        let v = Tensor::new(vec![1.0, -1.0], vec![1, 2]).unwrap();
        let eval = evaluate_quotient(Value::Tensor(u), Value::Tensor(v)).expect("polyder quotient");
        match eval.numerator() {
            Value::Tensor(t) => {
                assert_eq!(t.shape, vec![1, 3]);
                assert!(t
                    .data
                    .iter()
                    .zip([1.0, -2.0, 4.0])
                    .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
            }
            other => panic!("expected tensor numerator, got {other:?}"),
        }
        match eval.denominator() {
            Value::Tensor(t) => {
                assert_eq!(t.shape, vec![1, 3]);
                assert!(t
                    .data
                    .iter()
                    .zip([1.0, -2.0, 1.0])
                    .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
            }
            other => panic!("expected tensor denominator, got {other:?}"),
        }
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn column_vector_orientation_is_preserved() {
        let tensor = Tensor::new(vec![1.0, 0.0, -3.0], vec![3, 1]).unwrap();
        let result = derivative_single(Value::Tensor(tensor)).expect("polyder column");
        match result {
            Value::Tensor(t) => {
                assert_eq!(t.shape, vec![2, 1]);
                assert!(t
                    .data
                    .iter()
                    .zip([2.0, 0.0])
                    .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
            }
            other => panic!("expected column tensor, got {other:?}"),
        }
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn complex_coefficients_are_supported() {
        let tensor =
            ComplexTensor::new(vec![(1.0, 2.0), (-3.0, 0.0), (0.0, 4.0)], vec![1, 3]).unwrap();
        let result = derivative_single(Value::ComplexTensor(tensor)).expect("polyder complex");
        match result {
            Value::ComplexTensor(t) => {
                assert_eq!(t.shape, vec![1, 2]);
                let expected = [Complex64::new(2.0, 4.0), Complex64::new(-3.0, 0.0)];
                assert!(t
                    .data
                    .iter()
                    .zip(expected.iter())
                    .all(|((re, im), expected)| {
                        (re - expected.re).abs() < 1e-12 && (im - expected.im).abs() < 1e-12
                    }));
            }
            other => panic!("expected complex tensor, got {other:?}"),
        }
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn empty_polynomial_returns_zero() {
        let tensor = Tensor::new(Vec::new(), vec![1, 0]).unwrap();
        let result = derivative_single(Value::Tensor(tensor)).expect("polyder empty");
        assert_eq!(result, Value::Num(0.0));
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn rejects_matrix_input() {
        let tensor = Tensor::new(vec![1.0, 2.0, 3.0, 4.0], vec![2, 2]).unwrap();
        let err = derivative_single(Value::Tensor(tensor)).unwrap_err();
        assert_error_contains(err, "vector of coefficients");
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn rejects_string_input() {
        let err = derivative_single(Value::String("abc".into())).unwrap_err();
        assert_error_contains(err, "numeric vector");
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn mixed_gpu_cpu_product_falls_back_to_host() {
        test_support::with_test_provider(|provider| {
            let p = Tensor::new(vec![1.0, 0.0, -2.0], vec![1, 3]).unwrap();
            let q = Tensor::new(vec![1.0, 1.0], vec![1, 2]).unwrap();
            let cpu_expected =
                derivative_product(Value::Tensor(p.clone()), Value::Tensor(q.clone()))
                    .expect("cpu product");
            let Value::Tensor(cpu_tensor) = cpu_expected else {
                panic!("expected tensor result");
            };

            let view_p = runmat_accelerate_api::HostTensorView {
                data: &p.data,
                shape: &p.shape,
            };
            let handle_p = provider.upload(&view_p).expect("upload p");
            let result = derivative_product(Value::GpuTensor(handle_p), Value::Tensor(q))
                .expect("mixed product");
            let Value::Tensor(host_tensor) = result else {
                panic!("expected host tensor result");
            };
            assert_eq!(host_tensor.shape, cpu_tensor.shape);
            assert!(host_tensor
                .data
                .iter()
                .zip(cpu_tensor.data.iter())
                .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
        });
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn builtin_rejects_too_many_inputs() {
        let err = futures::executor::block_on(super::polyder_builtin(
            Value::Num(1.0),
            vec![Value::Num(2.0), Value::Num(3.0)],
        ))
        .unwrap_err();
        assert_error_contains(err, "too many input arguments");
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn gpu_inputs_remain_on_device() {
        test_support::with_test_provider(|provider| {
            let tensor = Tensor::new(vec![2.0, 0.0, -5.0, 4.0], vec![1, 4]).unwrap();
            let view = runmat_accelerate_api::HostTensorView {
                data: &tensor.data,
                shape: &tensor.shape,
            };
            let handle = provider.upload(&view).expect("upload");
            let result = derivative_single(Value::GpuTensor(handle)).expect("polyder gpu");
            let Value::GpuTensor(out_handle) = result else {
                panic!("expected GPU tensor result");
            };
            let gathered = test_support::gather(Value::GpuTensor(out_handle)).expect("gather");
            assert_eq!(gathered.shape, vec![1, 3]);
            assert!(gathered
                .data
                .iter()
                .zip([6.0, 0.0, -5.0])
                .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
        });
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn gpu_product_matches_cpu() {
        test_support::with_test_provider(|provider| {
            let p = Tensor::new(vec![1.0, 0.0, -2.0], vec![1, 3]).unwrap();
            let q = Tensor::new(vec![1.0, 1.0], vec![1, 2]).unwrap();
            let expected = derivative_product(Value::Tensor(p.clone()), Value::Tensor(q.clone()))
                .expect("cpu product");
            let Value::Tensor(expected_tensor) = expected else {
                panic!("expected tensor output");
            };

            let view_p = runmat_accelerate_api::HostTensorView {
                data: &p.data,
                shape: &p.shape,
            };
            let view_q = runmat_accelerate_api::HostTensorView {
                data: &q.data,
                shape: &q.shape,
            };
            let handle_p = provider.upload(&view_p).expect("upload p");
            let handle_q = provider.upload(&view_q).expect("upload q");
            let gpu_result =
                derivative_product(Value::GpuTensor(handle_p), Value::GpuTensor(handle_q))
                    .expect("gpu product");
            let Value::GpuTensor(gpu_handle) = gpu_result else {
                panic!("expected GPU tensor");
            };
            let gathered = test_support::gather(Value::GpuTensor(gpu_handle)).expect("gather");
            assert_eq!(gathered.shape, expected_tensor.shape);
            assert!(gathered
                .data
                .iter()
                .zip(expected_tensor.data.iter())
                .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
        });
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn gpu_quotient_matches_cpu() {
        test_support::with_test_provider(|provider| {
            let u = Tensor::new(vec![1.0, 0.0, -4.0], vec![1, 3]).unwrap();
            let v = Tensor::new(vec![1.0, -1.0], vec![1, 2]).unwrap();
            let expected = evaluate_quotient(Value::Tensor(u.clone()), Value::Tensor(v.clone()))
                .expect("cpu quotient");
            let Value::Tensor(expected_num) = expected.numerator() else {
                panic!("expected tensor numerator");
            };
            let Value::Tensor(expected_den) = expected.denominator() else {
                panic!("expected tensor denominator");
            };

            let view_u = runmat_accelerate_api::HostTensorView {
                data: &u.data,
                shape: &u.shape,
            };
            let view_v = runmat_accelerate_api::HostTensorView {
                data: &v.data,
                shape: &v.shape,
            };
            let handle_u = provider.upload(&view_u).expect("upload u");
            let handle_v = provider.upload(&view_v).expect("upload v");
            let gpu_eval =
                evaluate_quotient(Value::GpuTensor(handle_u), Value::GpuTensor(handle_v))
                    .expect("gpu quotient");
            let gpu_num = test_support::gather(gpu_eval.numerator()).expect("gather num");
            let gpu_den = test_support::gather(gpu_eval.denominator()).expect("gather den");
            assert_eq!(gpu_num.shape, expected_num.shape);
            assert_eq!(gpu_den.shape, expected_den.shape);
            assert!(gpu_num
                .data
                .iter()
                .zip(expected_num.data.iter())
                .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
            assert!(gpu_den
                .data
                .iter()
                .zip(expected_den.data.iter())
                .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
        });
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    #[cfg(feature = "wgpu")]
    fn wgpu_polyder_single_matches_cpu() {
        let _ = runmat_accelerate::backend::wgpu::provider::register_wgpu_provider(
            runmat_accelerate::backend::wgpu::provider::WgpuProviderOptions::default(),
        );
        let provider = runmat_accelerate_api::provider().expect("wgpu provider");
        let tensor = Tensor::new(vec![3.0, -2.0, 5.0, 7.0], vec![1, 4]).unwrap();
        let expected = derivative_single(Value::Tensor(tensor.clone())).expect("cpu polyder");
        let Value::Tensor(expected_tensor) = expected else {
            panic!("expected tensor");
        };
        let view = runmat_accelerate_api::HostTensorView {
            data: &tensor.data,
            shape: &tensor.shape,
        };
        let handle = provider.upload(&view).expect("upload");
        let gpu_result = derivative_single(Value::GpuTensor(handle)).expect("gpu polyder");
        let gathered = test_support::gather(gpu_result).expect("gather");
        assert_eq!(gathered.shape, expected_tensor.shape);
        assert!(gathered
            .data
            .iter()
            .zip(expected_tensor.data.iter())
            .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    #[cfg(feature = "wgpu")]
    fn wgpu_polyder_product_matches_cpu() {
        let _ = runmat_accelerate::backend::wgpu::provider::register_wgpu_provider(
            runmat_accelerate::backend::wgpu::provider::WgpuProviderOptions::default(),
        );
        let provider = runmat_accelerate_api::provider().expect("wgpu provider");
        let p = Tensor::new(vec![1.0, 0.0, -2.0], vec![1, 3]).unwrap();
        let q = Tensor::new(vec![1.0, 1.0], vec![1, 2]).unwrap();
        let expected = derivative_product(Value::Tensor(p.clone()), Value::Tensor(q.clone()))
            .expect("cpu product");
        let Value::Tensor(expected_tensor) = expected else {
            panic!("expected tensor");
        };
        let view_p = runmat_accelerate_api::HostTensorView {
            data: &p.data,
            shape: &p.shape,
        };
        let view_q = runmat_accelerate_api::HostTensorView {
            data: &q.data,
            shape: &q.shape,
        };
        let handle_p = provider.upload(&view_p).expect("upload p");
        let handle_q = provider.upload(&view_q).expect("upload q");
        let gpu_result = derivative_product(Value::GpuTensor(handle_p), Value::GpuTensor(handle_q))
            .expect("gpu product");
        let gathered = test_support::gather(gpu_result).expect("gather");
        assert_eq!(gathered.shape, expected_tensor.shape);
        assert!(gathered
            .data
            .iter()
            .zip(expected_tensor.data.iter())
            .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    #[cfg(feature = "wgpu")]
    fn wgpu_polyder_quotient_matches_cpu() {
        let _ = runmat_accelerate::backend::wgpu::provider::register_wgpu_provider(
            runmat_accelerate::backend::wgpu::provider::WgpuProviderOptions::default(),
        );
        let provider = runmat_accelerate_api::provider().expect("wgpu provider");
        let u = Tensor::new(vec![1.0, 0.0, -4.0], vec![1, 3]).unwrap();
        let v = Tensor::new(vec![1.0, -1.0], vec![1, 2]).unwrap();
        let expected = evaluate_quotient(Value::Tensor(u.clone()), Value::Tensor(v.clone()))
            .expect("cpu quotient");
        let expected_num = match expected.numerator() {
            Value::Tensor(t) => t,
            other => panic!("expected tensor numerator, got {other:?}"),
        };
        let expected_den = match expected.denominator() {
            Value::Tensor(t) => t,
            other => panic!("expected tensor denominator, got {other:?}"),
        };
        let view_u = runmat_accelerate_api::HostTensorView {
            data: &u.data,
            shape: &u.shape,
        };
        let view_v = runmat_accelerate_api::HostTensorView {
            data: &v.data,
            shape: &v.shape,
        };
        let handle_u = provider.upload(&view_u).expect("upload u");
        let handle_v = provider.upload(&view_v).expect("upload v");
        let gpu_eval = evaluate_quotient(Value::GpuTensor(handle_u), Value::GpuTensor(handle_v))
            .expect("gpu quotient");
        let gpu_num = test_support::gather(gpu_eval.numerator()).expect("gather num");
        let gpu_den = test_support::gather(gpu_eval.denominator()).expect("gather den");
        assert_eq!(gpu_num.shape, expected_num.shape);
        assert_eq!(gpu_den.shape, expected_den.shape);
        assert!(gpu_num
            .data
            .iter()
            .zip(expected_num.data.iter())
            .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
        assert!(gpu_den
            .data
            .iter()
            .zip(expected_den.data.iter())
            .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn derivative_promotes_integers() {
        let value = Value::Int(IntValue::I32(5));
        let result = derivative_single(value).expect("polyder int");
        assert_eq!(result, Value::Num(0.0));
    }

    fn derivative_single(value: Value) -> BuiltinResult<Value> {
        block_on(super::derivative_single(value))
    }

    fn derivative_product(first: Value, second: Value) -> BuiltinResult<Value> {
        block_on(super::derivative_product(first, second))
    }

    fn evaluate_quotient(u: Value, v: Value) -> BuiltinResult<PolyderEval> {
        block_on(super::evaluate_quotient(u, v))
    }
}