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::{
    BuiltinCompletionPolicy, BuiltinDescriptor, BuiltinErrorDescriptor, BuiltinExtensionDescriptor,
    BuiltinExtensionMode, BuiltinIntegerBackendRule, BuiltinIntegerCapabilityDescriptor,
    BuiltinIntegerComputationDomain, BuiltinIntegerInputAvailability,
    BuiltinIntegerInputCapability, BuiltinIntegerOutputClassRule, BuiltinIntegerOverflowRule,
    BuiltinIntegerOverloadKind, BuiltinIntegerScalarDoubleRule, BuiltinOutputMode,
    BuiltinParamArity, BuiltinParamDescriptor, BuiltinParamType, BuiltinSignatureDescriptor,
};
use runmat_macros::runtime_builtin;
use runmat_value::{ComplexTensor, Tensor, Value};

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";

const POLYDER_OUTPUT_D: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
    name: "d",
    ty: BuiltinParamType::Any,
    arity: BuiltinParamArity::Required,
    default: None,
    description: "Derivative coefficient vector.",
}];

const POLYDER_OUTPUT_NUM_DEN: [BuiltinParamDescriptor; 2] = [
    BuiltinParamDescriptor {
        name: "num",
        ty: BuiltinParamType::Any,
        arity: BuiltinParamArity::Required,
        default: None,
        description: "Quotient-rule numerator coefficients.",
    },
    BuiltinParamDescriptor {
        name: "den",
        ty: BuiltinParamType::Any,
        arity: BuiltinParamArity::Required,
        default: None,
        description: "Quotient-rule denominator coefficients.",
    },
];

const POLYDER_INPUTS_SINGLE: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
    name: "p",
    ty: BuiltinParamType::Any,
    arity: BuiltinParamArity::Required,
    default: None,
    description: "Polynomial coefficient vector.",
}];

const POLYDER_INPUTS_BINARY: [BuiltinParamDescriptor; 2] = [
    BuiltinParamDescriptor {
        name: "a",
        ty: BuiltinParamType::Any,
        arity: BuiltinParamArity::Required,
        default: None,
        description: "First polynomial coefficient vector.",
    },
    BuiltinParamDescriptor {
        name: "b",
        ty: BuiltinParamType::Any,
        arity: BuiltinParamArity::Required,
        default: None,
        description: "Second polynomial coefficient vector.",
    },
];

const POLYDER_SIGNATURES: [BuiltinSignatureDescriptor; 3] = [
    BuiltinSignatureDescriptor {
        label: "d = polyder(p)",
        inputs: &POLYDER_INPUTS_SINGLE,
        outputs: &POLYDER_OUTPUT_D,
    },
    BuiltinSignatureDescriptor {
        label: "d = polyder(a, b)",
        inputs: &POLYDER_INPUTS_BINARY,
        outputs: &POLYDER_OUTPUT_D,
    },
    BuiltinSignatureDescriptor {
        label: "[num, den] = polyder(u, v)",
        inputs: &POLYDER_INPUTS_BINARY,
        outputs: &POLYDER_OUTPUT_NUM_DEN,
    },
];

const POLYDER_ERROR_INVALID_ARGUMENT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
    code: "RM.POLYDER.INVALID_ARGUMENT",
    identifier: Some("RunMat:polyder:InvalidArgument"),
    when: "Input arity/output mode combination is invalid.",
    message: "polyder: invalid argument",
};

const POLYDER_ERROR_INVALID_INPUT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
    code: "RM.POLYDER.INVALID_INPUT",
    identifier: Some("RunMat:polyder:InvalidInput"),
    when: "Inputs cannot be interpreted as numeric coefficient vectors.",
    message: "polyder: invalid input",
};

const POLYDER_ERROR_INTERNAL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
    code: "RM.POLYDER.INTERNAL",
    identifier: Some("RunMat:polyder:Internal"),
    when: "Runtime fails while building derivative outputs or provider fallback paths.",
    message: "polyder: internal runtime failure",
};

const POLYDER_ERRORS: [BuiltinErrorDescriptor; 3] = [
    POLYDER_ERROR_INVALID_ARGUMENT,
    POLYDER_ERROR_INVALID_INPUT,
    POLYDER_ERROR_INTERNAL,
];

pub const POLYDER_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
    signatures: &POLYDER_SIGNATURES,
    output_mode: BuiltinOutputMode::ByRequestedOutputCount,
    completion_policy: BuiltinCompletionPolicy::Public,
    errors: &POLYDER_ERRORS,
};

const POLYDER_INTEGER_COEFFICIENTS_EXTENSION: BuiltinExtensionDescriptor =
    BuiltinExtensionDescriptor {
        id: "polyder-integer-coefficients",
        mode: BuiltinExtensionMode::RunMatOnly,
        description: "polyder accepts typed-integer coefficient vectors as a RunMat extension",
        error_identifier: Some("RunMat:compatibility:PolyderIntegerCoefficientsExtension"),
    };
pub const POLYDER_EXTENSIONS: [BuiltinExtensionDescriptor; 1] =
    [POLYDER_INTEGER_COEFFICIENTS_EXTENSION];
const POLYDER_INTEGER_INPUTS: [BuiltinIntegerInputCapability; 1] =
    [BuiltinIntegerInputCapability {
        name: "p, a, or b",
        classes: &crate::builtins::common::integer_capability::ALL_INTEGER_CLASSES,
        availability: BuiltinIntegerInputAvailability::RunMatOnly,
        scalar_double: BuiltinIntegerScalarDoubleRule::NotApplicable,
        notes: "The compatibility target documents single and double polynomial coefficients. RunMat admits typed integers only after exact conversion to the floating polynomial domain.",
    }];
pub const POLYDER_INTEGER_CAPABILITIES: [BuiltinIntegerCapabilityDescriptor; 1] =
    [BuiltinIntegerCapabilityDescriptor {
        form: "d = polyder(integer_p) or d = polyder(integer_a,b) or [q,d] = polyder(integer_a,b)",
        inputs: &POLYDER_INTEGER_INPUTS,
        computation_domain: BuiltinIntegerComputationDomain::FloatingPoint,
        output_class: BuiltinIntegerOutputClassRule::FunctionSpecific,
        overflow: BuiltinIntegerOverflowRule::Error,
        backend: BuiltinIntegerBackendRule::GatherFallback,
        overload: BuiltinIntegerOverloadKind::Multiple,
        notes: "Every integer coefficient operand is checked before provider dispatch; derivative, convolution, and quotient arithmetic then use the selected floating precision domain.",
    }];

#[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 {
    polyder_error_with(message, &POLYDER_ERROR_INVALID_INPUT)
}

fn polyder_argument_error(message: impl Into<String>) -> RuntimeError {
    polyder_error_with(message, &POLYDER_ERROR_INVALID_ARGUMENT)
}

fn polyder_error_with(
    message: impl Into<String>,
    error: &'static BuiltinErrorDescriptor,
) -> RuntimeError {
    let mut builder = build_runtime_error(message).with_builtin(BUILTIN_NAME);
    if let Some(identifier) = error.identifier {
        builder = builder.with_identifier(identifier);
    }
    builder.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.",
    keywords = "polyder,polynomial,derivative,product,quotient",
    type_resolver(polyder_type),
    descriptor(crate::builtins::math::poly::polyder::POLYDER_DESCRIPTOR),
    extensions(crate::builtins::math::poly::polyder::POLYDER_EXTENSIONS),
    integer_capabilities(crate::builtins::math::poly::polyder::POLYDER_INTEGER_CAPABILITIES),
    builtin_path = "crate::builtins::math::poly::polyder"
)]
async fn polyder_builtin(first: Value, rest: Vec<Value>) -> crate::BuiltinResult<Value> {
    if rest.len() > 1 {
        return Err(polyder_argument_error("polyder: too many input arguments"));
    }
    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,
                _ => unreachable!("input count validated above"),
            }?;
            if out_count == 0 {
                return Ok(Value::OutputList(Vec::new()));
            }
            return Ok(Value::OutputList(vec![result]));
        }
        if rest.len() != 1 {
            return Err(polyder_argument_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,
        _ => unreachable!("input count validated above"),
    }
}

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> {
    gate_coefficient(&u).await?;
    gate_coefficient(&v).await?;
    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> {
    gate_coefficient(&value).await?;
    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> {
    gate_coefficient(&first).await?;
    gate_coefficient(&second).await?;
    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)
}

async fn gate_coefficient(value: &Value) -> BuiltinResult<()> {
    crate::builtins::common::validation::reject_typed_complex_integer(value, BUILTIN_NAME)?;
    crate::builtins::common::validation::ensure_runmat_integer_f64_boundary(
        value,
        &POLYDER_INTEGER_COEFFICIENTS_EXTENSION,
        BUILTIN_NAME,
        "coefficient",
    )
    .await
}

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::tensor_element_len(&tensor) == 0 {
                (vec![Complex64::new(0.0, 0.0)], orientation)
            } else {
                (
                    tensor::tensor_values_f64(&tensor)
                        .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::complex_tensor_element_len(&tensor) == 0 {
                (vec![Complex64::new(0.0, 0.0)], orientation)
            } else {
                (
                    tensor::complex_tensor_into_values_complex64(tensor),
                    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.is_empty() {
                (vec![Complex64::new(0.0, 0.0)], orientation)
            } else {
                (
                    tensor::tensor_values_f64(&tensor)
                        .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_value::{IntValue, IntegerComplexStorage, IntegerStorage, Tensor};

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

    #[test]
    fn polyder_descriptor_signatures_cover_core_forms() {
        let labels: Vec<&str> = POLYDER_DESCRIPTOR
            .signatures
            .iter()
            .map(|signature| signature.label)
            .collect();
        assert!(labels.contains(&"d = polyder(p)"));
        assert!(labels.contains(&"d = polyder(a, b)"));
        assert!(labels.contains(&"[num, den] = polyder(u, v)"));
    }

    #[test]
    fn polyder_descriptor_errors_have_stable_codes() {
        let codes: Vec<&str> = POLYDER_DESCRIPTOR
            .errors
            .iter()
            .map(|error| error.code)
            .collect();
        assert!(codes.contains(&"RM.POLYDER.INVALID_ARGUMENT"));
        assert!(codes.contains(&"RM.POLYDER.INVALID_INPUT"));
        assert!(codes.contains(&"RM.POLYDER.INTERNAL"));
    }

    #[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
                    .materialize_f64()
                    .iter()
                    .zip([9.0, -4.0, 5.0])
                    .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
            }
            other => panic!("expected tensor result, got {other:?}"),
        }
    }

    #[test]
    fn derivative_typed_integer_coefficients_cross_double_boundary_exactly() {
        let _extensions = crate::compatibility::push_runmat_extensions_enabled(true);
        let tensor =
            Tensor::new_integer(IntegerStorage::I16(vec![3, -2, 5, 7]), 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
                    .materialize_f64()
                    .iter()
                    .zip([9.0, -4.0, 5.0])
                    .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
                assert!(t.integer_storage().is_none());
            }
            other => panic!("expected tensor result, got {other:?}"),
        }
    }

    #[test]
    fn derivative_typed_complex_integer_coefficients_reject_before_conversion() {
        let _extensions = crate::compatibility::push_runmat_extensions_enabled(true);
        let tensor = ComplexTensor::new_integer(
            IntegerComplexStorage::new(
                IntegerStorage::I16(vec![3, -2, 5, 7]),
                IntegerStorage::I16(vec![1, 0, -1, 2]),
            )
            .expect("complex integer storage"),
            vec![1, 4],
        )
        .expect("complex integer tensor");

        let error = derivative_single(Value::ComplexTensor(tensor))
            .expect_err("typed complex integer arithmetic must reject");
        assert!(
            error
                .message()
                .contains("complex numbers with integer types are not supported"),
            "{error:?}"
        );
    }

    #[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
                    .materialize_f64()
                    .iter()
                    .zip([3.0, 2.0, -2.0])
                    .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
            }
            other => panic!("expected tensor result, got {other:?}"),
        }
    }

    #[test]
    fn product_typed_integer_coefficients_cross_double_boundary_exactly() {
        let _extensions = crate::compatibility::push_runmat_extensions_enabled(true);
        let p = Tensor::new_integer(IntegerStorage::I16(vec![1, 0, -2]), vec![1, 3]).unwrap();
        let a = Tensor::new_integer(IntegerStorage::U16(vec![1, 1]), 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
                    .materialize_f64()
                    .iter()
                    .zip([3.0, 2.0, -2.0])
                    .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
                assert!(t.integer_storage().is_none());
            }
            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
                    .materialize_f64()
                    .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
                    .materialize_f64()
                    .iter()
                    .zip([1.0, -2.0, 1.0])
                    .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
            }
            other => panic!("expected tensor denominator, got {other:?}"),
        }
    }

    #[test]
    fn quotient_typed_integer_coefficients_cross_double_boundary_exactly() {
        let _extensions = crate::compatibility::push_runmat_extensions_enabled(true);
        let u = Tensor::new_integer(IntegerStorage::I16(vec![1, 0, -4]), vec![1, 3]).unwrap();
        let v = Tensor::new_integer(IntegerStorage::I16(vec![1, -1]), 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
                    .materialize_f64()
                    .iter()
                    .zip([1.0, -2.0, 4.0])
                    .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
                assert!(t.integer_storage().is_none());
            }
            other => panic!("expected tensor numerator, got {other:?}"),
        }
        match eval.denominator() {
            Value::Tensor(t) => {
                assert_eq!(t.shape, vec![1, 3]);
                assert!(t
                    .materialize_f64()
                    .iter()
                    .zip([1.0, -2.0, 1.0])
                    .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
                assert!(t.integer_storage().is_none());
            }
            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
                    .materialize_f64()
                    .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.materialize_f64().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.materialize_f64(),
                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
                .materialize_f64()
                .iter()
                .zip(cpu_tensor.materialize_f64().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_eq!(err.identifier(), POLYDER_ERROR_INVALID_ARGUMENT.identifier);
        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.materialize_f64(),
                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
                .materialize_f64()
                .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.materialize_f64(),
                shape: &p.shape,
            };
            let view_q = runmat_accelerate_api::HostTensorView {
                data: &q.materialize_f64(),
                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
                .materialize_f64()
                .iter()
                .zip(expected_tensor.materialize_f64().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.materialize_f64(),
                shape: &u.shape,
            };
            let view_v = runmat_accelerate_api::HostTensorView {
                data: &v.materialize_f64(),
                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
                .materialize_f64()
                .iter()
                .zip(expected_num.materialize_f64().iter())
                .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
            assert!(gpu_den
                .materialize_f64()
                .iter()
                .zip(expected_den.materialize_f64().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.materialize_f64(),
            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
            .materialize_f64()
            .iter()
            .zip(expected_tensor.materialize_f64().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.materialize_f64(),
            shape: &p.shape,
        };
        let view_q = runmat_accelerate_api::HostTensorView {
            data: &q.materialize_f64(),
            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
            .materialize_f64()
            .iter()
            .zip(expected_tensor.materialize_f64().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.materialize_f64(),
            shape: &u.shape,
        };
        let view_v = runmat_accelerate_api::HostTensorView {
            data: &v.materialize_f64(),
            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
            .materialize_f64()
            .iter()
            .zip(expected_num.materialize_f64().iter())
            .all(|(lhs, rhs)| (lhs - rhs).abs() < 1e-12));
        assert!(gpu_den
            .materialize_f64()
            .iter()
            .zip(expected_den.materialize_f64().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 _extensions = crate::compatibility::push_runmat_extensions_enabled(true);
        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))
    }
}