cddlib_rs/

lib.rs

//! Safe, idiomatic Rust bindings on top of the raw `cddlib-sys` FFI.
//!
//! This crate owns the unsafe glue for working with cddlib matrices, polyhedra,
//! and LP solutions. Callers should be able to use these types without touching
//! raw pointers or the `cddlib-sys` bindings directly.

use std::convert::TryFrom;
use std::ffi::c_void;
use std::marker::PhantomData;
use std::mem;
use std::os::raw::c_long;
use std::ptr::{self, NonNull};
use std::rc::Rc;
use std::sync::Once;

#[cfg(any(feature = "gmp", feature = "gmprational"))]
use std::mem::MaybeUninit;

use thiserror::Error;

pub use cddlib_sys as sys;

#[cfg(all(
    not(feature = "gmprational"),
    not(feature = "gmp"),
    not(feature = "f64")
))]
compile_error!("cddlib-rs requires at least one backend feature: f64, gmp, or gmprational");

#[cfg(feature = "gmp")]
pub struct CddFloat {
    inner: sys::gmpfloat::mytype,
}

#[cfg(feature = "gmprational")]
pub struct CddRational {
    inner: sys::gmprational::mytype,
}

#[cfg(feature = "f64")]
pub type DefaultNumber = f64;

#[cfg(all(not(feature = "f64"), feature = "gmprational"))]
pub type DefaultNumber = CddRational;

#[cfg(all(not(feature = "f64"), not(feature = "gmprational"), feature = "gmp"))]
pub type DefaultNumber = CddFloat;

mod sealed {
    pub trait Sealed {}
}

pub trait CddNumber: sealed::Sealed + 'static {
    const MYTYPE_SIZE: usize;
    const DEFAULT_NUMBER_TYPE: NumberType;

    fn ensure_initialized();

    #[doc(hidden)]
    fn dd_no_error() -> u32;

    #[doc(hidden)]
    unsafe fn dd_create_matrix(rows: c_long, cols: c_long) -> *mut c_void;

    #[doc(hidden)]
    unsafe fn dd_copy_matrix(matrix: *mut c_void) -> *mut c_void;

    #[doc(hidden)]
    unsafe fn dd_free_matrix(matrix: *mut c_void);

    #[doc(hidden)]
    unsafe fn dd_matrix_append_to(dst: *mut *mut c_void, src: *mut c_void) -> i32;

    #[doc(hidden)]
    unsafe fn dd_append_matrix(a: *mut c_void, b: *mut c_void) -> *mut c_void;

    #[doc(hidden)]
    unsafe fn dd_matrix_row_remove(matrix: *mut *mut c_void, row: c_long) -> i32;

    #[doc(hidden)]
    unsafe fn dd_initialize_arow(cols: c_long, out: *mut *mut c_void);

    #[doc(hidden)]
    unsafe fn dd_free_arow(cols: c_long, arow: *mut c_void);

    #[doc(hidden)]
    unsafe fn dd_redundant(
        matrix: *mut c_void,
        row: c_long,
        cert: *mut c_void,
        err: *mut u32,
    ) -> i32;

    #[doc(hidden)]
    unsafe fn dd_redundant_rows(matrix: *mut c_void, err: *mut u32) -> *mut libc::c_ulong;

    #[doc(hidden)]
    unsafe fn dd_matrix_canonicalize(
        matrix: *mut *mut c_void,
        impl_lin: *mut *mut libc::c_ulong,
        redset: *mut *mut libc::c_ulong,
        newpos: *mut *mut c_long,
        err: *mut u32,
    ) -> i32;

    #[doc(hidden)]
    unsafe fn set_groundsize(set: *mut libc::c_ulong) -> c_long;

    #[doc(hidden)]
    unsafe fn set_member(elem: c_long, set: *mut libc::c_ulong) -> i32;

    #[doc(hidden)]
    unsafe fn set_free(set: *mut libc::c_ulong);

    #[doc(hidden)]
    unsafe fn dd_matrix2poly(matrix: *mut c_void, err: *mut u32) -> *mut c_void;

    #[doc(hidden)]
    unsafe fn dd_free_polyhedra(poly: *mut c_void);

    #[doc(hidden)]
    unsafe fn dd_copy_inequalities(poly: *mut c_void) -> *mut c_void;

    #[doc(hidden)]
    unsafe fn dd_copy_generators(poly: *mut c_void) -> *mut c_void;

    #[doc(hidden)]
    unsafe fn dd_copy_adjacency(poly: *mut c_void) -> *mut c_void;

    #[doc(hidden)]
    unsafe fn dd_copy_input_adjacency(poly: *mut c_void) -> *mut c_void;

    #[doc(hidden)]
    unsafe fn dd_copy_incidence(poly: *mut c_void) -> *mut c_void;

    #[doc(hidden)]
    unsafe fn dd_copy_input_incidence(poly: *mut c_void) -> *mut c_void;

    #[doc(hidden)]
    unsafe fn dd_append_matrix2poly(poly: *mut *mut c_void, rows: *mut c_void) -> i32;

    #[doc(hidden)]
    unsafe fn dd_free_set_family(family: *mut c_void);

    #[doc(hidden)]
    unsafe fn dd_matrix2lp(matrix: *mut c_void, err: *mut u32) -> *mut c_void;

    #[doc(hidden)]
    unsafe fn dd_free_lp_data(lp: *mut c_void);

    #[doc(hidden)]
    unsafe fn dd_lp_solve_dual_simplex(lp: *mut c_void, err: *mut u32) -> i32;

    #[doc(hidden)]
    unsafe fn lp_status_raw(lp: *mut c_void) -> u32;

    #[doc(hidden)]
    unsafe fn dd_copy_lp_solution(lp: *mut c_void) -> *mut c_void;

    #[doc(hidden)]
    unsafe fn dd_free_lp_solution(sol: *mut c_void);

    #[doc(hidden)]
    unsafe fn lp_solution_optvalue_ptr(sol: *mut c_void) -> *const c_void;

    #[doc(hidden)]
    unsafe fn write_mytype_real(target: *mut c_void, value: f64);

    #[doc(hidden)]
    unsafe fn write_mytype_int(target: *mut c_void, value: c_long);

    #[doc(hidden)]
    unsafe fn read_mytype_real(source: *const c_void) -> f64;

    #[doc(hidden)]
    unsafe fn write_mytype(target: *mut c_void, value: &Self);

    #[doc(hidden)]
    unsafe fn read_mytype(source: *const c_void) -> Self;
}

#[cfg(feature = "f64")]
impl sealed::Sealed for f64 {}

#[cfg(feature = "f64")]
impl CddNumber for f64 {
    const MYTYPE_SIZE: usize = mem::size_of::<sys::f64::mytype>();
    const DEFAULT_NUMBER_TYPE: NumberType = NumberType::Real;

    fn ensure_initialized() {
        static INIT: Once = Once::new();
        INIT.call_once(|| unsafe {
            sys::f64::dd_set_global_constants();
        });
    }

    fn dd_no_error() -> u32 {
        sys::f64::dd_ErrorType_dd_NoError
    }

    unsafe fn dd_create_matrix(rows: c_long, cols: c_long) -> *mut c_void {
        unsafe { sys::f64::dd_CreateMatrix(rows, cols).cast() }
    }

    unsafe fn dd_copy_matrix(matrix: *mut c_void) -> *mut c_void {
        unsafe { sys::f64::dd_CopyMatrix(matrix.cast()).cast() }
    }

    unsafe fn dd_free_matrix(matrix: *mut c_void) {
        unsafe {
            sys::f64::dd_FreeMatrix(matrix.cast());
        }
    }

    unsafe fn dd_matrix_append_to(dst: *mut *mut c_void, src: *mut c_void) -> i32 {
        unsafe { sys::f64::dd_MatrixAppendTo(dst.cast(), src.cast()) }
    }

    unsafe fn dd_append_matrix(a: *mut c_void, b: *mut c_void) -> *mut c_void {
        unsafe { sys::f64::dd_AppendMatrix(a.cast(), b.cast()).cast() }
    }

    unsafe fn dd_matrix_row_remove(matrix: *mut *mut c_void, row: c_long) -> i32 {
        unsafe { sys::f64::dd_MatrixRowRemove(matrix.cast(), row) }
    }

    unsafe fn dd_initialize_arow(cols: c_long, out: *mut *mut c_void) {
        unsafe {
            sys::f64::dd_InitializeArow(cols, out.cast());
        }
    }

    unsafe fn dd_free_arow(cols: c_long, arow: *mut c_void) {
        unsafe {
            sys::f64::dd_FreeArow(cols, arow.cast());
        }
    }

    unsafe fn dd_redundant(
        matrix: *mut c_void,
        row: c_long,
        cert: *mut c_void,
        err: *mut u32,
    ) -> i32 {
        unsafe { sys::f64::dd_Redundant(matrix.cast(), row, cert.cast(), err.cast()) }
    }

    unsafe fn dd_redundant_rows(matrix: *mut c_void, err: *mut u32) -> *mut libc::c_ulong {
        unsafe { sys::f64::dd_RedundantRows(matrix.cast(), err.cast()).cast() }
    }

    unsafe fn dd_matrix_canonicalize(
        matrix: *mut *mut c_void,
        impl_lin: *mut *mut libc::c_ulong,
        redset: *mut *mut libc::c_ulong,
        newpos: *mut *mut c_long,
        err: *mut u32,
    ) -> i32 {
        unsafe {
            sys::f64::dd_MatrixCanonicalize(
                matrix.cast(),
                impl_lin.cast(),
                redset.cast(),
                newpos.cast(),
                err.cast(),
            )
        }
    }

    unsafe fn set_groundsize(set: *mut libc::c_ulong) -> c_long {
        unsafe { sys::f64::set_groundsize(set.cast()) }
    }

    unsafe fn set_member(elem: c_long, set: *mut libc::c_ulong) -> i32 {
        unsafe { sys::f64::set_member(elem, set.cast()) }
    }

    unsafe fn set_free(set: *mut libc::c_ulong) {
        unsafe {
            sys::f64::set_free(set.cast());
        }
    }

    unsafe fn dd_matrix2poly(matrix: *mut c_void, err: *mut u32) -> *mut c_void {
        unsafe { sys::f64::dd_DDMatrix2Poly(matrix.cast(), err.cast()).cast() }
    }

    unsafe fn dd_free_polyhedra(poly: *mut c_void) {
        unsafe {
            sys::f64::dd_FreePolyhedra(poly.cast());
        }
    }

    unsafe fn dd_copy_inequalities(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::f64::dd_CopyInequalities(poly.cast()).cast() }
    }

    unsafe fn dd_copy_generators(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::f64::dd_CopyGenerators(poly.cast()).cast() }
    }

    unsafe fn dd_copy_adjacency(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::f64::dd_CopyAdjacency(poly.cast()).cast() }
    }

    unsafe fn dd_copy_input_adjacency(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::f64::dd_CopyInputAdjacency(poly.cast()).cast() }
    }

    unsafe fn dd_copy_incidence(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::f64::dd_CopyIncidence(poly.cast()).cast() }
    }

    unsafe fn dd_copy_input_incidence(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::f64::dd_CopyInputIncidence(poly.cast()).cast() }
    }

    unsafe fn dd_append_matrix2poly(poly: *mut *mut c_void, rows: *mut c_void) -> i32 {
        unsafe { sys::f64::dd_AppendMatrix2Poly(poly.cast(), rows.cast()) }
    }

    unsafe fn dd_free_set_family(family: *mut c_void) {
        unsafe {
            sys::f64::dd_FreeSetFamily(family.cast());
        }
    }

    unsafe fn dd_matrix2lp(matrix: *mut c_void, err: *mut u32) -> *mut c_void {
        unsafe { sys::f64::dd_Matrix2LP(matrix.cast(), err.cast()).cast() }
    }

    unsafe fn dd_free_lp_data(lp: *mut c_void) {
        unsafe {
            sys::f64::dd_FreeLPData(lp.cast());
        }
    }

    unsafe fn dd_lp_solve_dual_simplex(lp: *mut c_void, err: *mut u32) -> i32 {
        unsafe {
            sys::f64::dd_LPSolve(
                lp.cast(),
                sys::f64::dd_LPSolverType_dd_DualSimplex,
                err.cast(),
            )
        }
    }

    unsafe fn lp_status_raw(lp: *mut c_void) -> u32 {
        unsafe { (*lp.cast::<sys::f64::dd_lpdata>()).LPS }
    }

    unsafe fn dd_copy_lp_solution(lp: *mut c_void) -> *mut c_void {
        unsafe { sys::f64::dd_CopyLPSolution(lp.cast()).cast() }
    }

    unsafe fn dd_free_lp_solution(sol: *mut c_void) {
        unsafe {
            sys::f64::dd_FreeLPSolution(sol.cast());
        }
    }

    unsafe fn lp_solution_optvalue_ptr(sol: *mut c_void) -> *const c_void {
        unsafe {
            let sol = sol.cast::<sys::f64::dd_lpsolution>();
            (&(*sol).optvalue as *const sys::f64::mytype).cast()
        }
    }

    unsafe fn write_mytype_real(target: *mut c_void, value: f64) {
        unsafe {
            (*target.cast::<sys::f64::mytype>())[0] = value;
        }
    }

    unsafe fn write_mytype_int(target: *mut c_void, value: c_long) {
        unsafe {
            (*target.cast::<sys::f64::mytype>())[0] = value as f64;
        }
    }

    unsafe fn read_mytype_real(source: *const c_void) -> f64 {
        unsafe { (*source.cast::<sys::f64::mytype>())[0] }
    }

    unsafe fn write_mytype(target: *mut c_void, value: &Self) {
        unsafe {
            Self::write_mytype_real(target, *value);
        }
    }

    unsafe fn read_mytype(source: *const c_void) -> Self {
        unsafe { Self::read_mytype_real(source) }
    }
}

#[cfg(feature = "gmp")]
impl sealed::Sealed for CddFloat {}

#[cfg(feature = "gmp")]
impl CddNumber for CddFloat {
    const MYTYPE_SIZE: usize = mem::size_of::<sys::gmpfloat::mytype>();
    const DEFAULT_NUMBER_TYPE: NumberType = NumberType::Real;

    fn ensure_initialized() {
        static INIT: Once = Once::new();
        INIT.call_once(|| unsafe {
            sys::gmpfloat::dd_set_global_constants();
        });
    }

    fn dd_no_error() -> u32 {
        sys::gmpfloat::dd_ErrorType_dd_NoError
    }

    unsafe fn dd_create_matrix(rows: c_long, cols: c_long) -> *mut c_void {
        unsafe { sys::gmpfloat::dd_CreateMatrix(rows, cols).cast() }
    }

    unsafe fn dd_copy_matrix(matrix: *mut c_void) -> *mut c_void {
        unsafe { sys::gmpfloat::dd_CopyMatrix(matrix.cast()).cast() }
    }

    unsafe fn dd_free_matrix(matrix: *mut c_void) {
        unsafe {
            sys::gmpfloat::dd_FreeMatrix(matrix.cast());
        }
    }

    unsafe fn dd_matrix_append_to(dst: *mut *mut c_void, src: *mut c_void) -> i32 {
        unsafe { sys::gmpfloat::dd_MatrixAppendTo(dst.cast(), src.cast()) }
    }

    unsafe fn dd_append_matrix(a: *mut c_void, b: *mut c_void) -> *mut c_void {
        unsafe { sys::gmpfloat::dd_AppendMatrix(a.cast(), b.cast()).cast() }
    }

    unsafe fn dd_matrix_row_remove(matrix: *mut *mut c_void, row: c_long) -> i32 {
        unsafe { sys::gmpfloat::dd_MatrixRowRemove(matrix.cast(), row) }
    }

    unsafe fn dd_initialize_arow(cols: c_long, out: *mut *mut c_void) {
        unsafe {
            sys::gmpfloat::dd_InitializeArow(cols, out.cast());
        }
    }

    unsafe fn dd_free_arow(cols: c_long, arow: *mut c_void) {
        unsafe {
            sys::gmpfloat::dd_FreeArow(cols, arow.cast());
        }
    }

    unsafe fn dd_redundant(
        matrix: *mut c_void,
        row: c_long,
        cert: *mut c_void,
        err: *mut u32,
    ) -> i32 {
        unsafe { sys::gmpfloat::dd_Redundant(matrix.cast(), row, cert.cast(), err.cast()) }
    }

    unsafe fn dd_redundant_rows(matrix: *mut c_void, err: *mut u32) -> *mut libc::c_ulong {
        unsafe { sys::gmpfloat::dd_RedundantRows(matrix.cast(), err.cast()).cast() }
    }

    unsafe fn dd_matrix_canonicalize(
        matrix: *mut *mut c_void,
        impl_lin: *mut *mut libc::c_ulong,
        redset: *mut *mut libc::c_ulong,
        newpos: *mut *mut c_long,
        err: *mut u32,
    ) -> i32 {
        unsafe {
            sys::gmpfloat::dd_MatrixCanonicalize(
                matrix.cast(),
                impl_lin.cast(),
                redset.cast(),
                newpos.cast(),
                err.cast(),
            )
        }
    }

    unsafe fn set_groundsize(set: *mut libc::c_ulong) -> c_long {
        unsafe { sys::gmpfloat::set_groundsize(set.cast()) }
    }

    unsafe fn set_member(elem: c_long, set: *mut libc::c_ulong) -> i32 {
        unsafe { sys::gmpfloat::set_member(elem, set.cast()) }
    }

    unsafe fn set_free(set: *mut libc::c_ulong) {
        unsafe {
            sys::gmpfloat::set_free(set.cast());
        }
    }

    unsafe fn dd_matrix2poly(matrix: *mut c_void, err: *mut u32) -> *mut c_void {
        unsafe { sys::gmpfloat::dd_DDMatrix2Poly(matrix.cast(), err.cast()).cast() }
    }

    unsafe fn dd_free_polyhedra(poly: *mut c_void) {
        unsafe {
            sys::gmpfloat::dd_FreePolyhedra(poly.cast());
        }
    }

    unsafe fn dd_copy_inequalities(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::gmpfloat::dd_CopyInequalities(poly.cast()).cast() }
    }

    unsafe fn dd_copy_generators(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::gmpfloat::dd_CopyGenerators(poly.cast()).cast() }
    }

    unsafe fn dd_copy_adjacency(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::gmpfloat::dd_CopyAdjacency(poly.cast()).cast() }
    }

    unsafe fn dd_copy_input_adjacency(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::gmpfloat::dd_CopyInputAdjacency(poly.cast()).cast() }
    }

    unsafe fn dd_copy_incidence(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::gmpfloat::dd_CopyIncidence(poly.cast()).cast() }
    }

    unsafe fn dd_copy_input_incidence(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::gmpfloat::dd_CopyInputIncidence(poly.cast()).cast() }
    }

    unsafe fn dd_append_matrix2poly(poly: *mut *mut c_void, rows: *mut c_void) -> i32 {
        unsafe { sys::gmpfloat::dd_AppendMatrix2Poly(poly.cast(), rows.cast()) }
    }

    unsafe fn dd_free_set_family(family: *mut c_void) {
        unsafe {
            sys::gmpfloat::dd_FreeSetFamily(family.cast());
        }
    }

    unsafe fn dd_matrix2lp(matrix: *mut c_void, err: *mut u32) -> *mut c_void {
        unsafe { sys::gmpfloat::dd_Matrix2LP(matrix.cast(), err.cast()).cast() }
    }

    unsafe fn dd_free_lp_data(lp: *mut c_void) {
        unsafe {
            sys::gmpfloat::dd_FreeLPData(lp.cast());
        }
    }

    unsafe fn dd_lp_solve_dual_simplex(lp: *mut c_void, err: *mut u32) -> i32 {
        unsafe {
            sys::gmpfloat::dd_LPSolve(
                lp.cast(),
                sys::gmpfloat::dd_LPSolverType_dd_DualSimplex,
                err.cast(),
            )
        }
    }

    unsafe fn lp_status_raw(lp: *mut c_void) -> u32 {
        unsafe { (*lp.cast::<sys::gmpfloat::dd_lpdata>()).LPS }
    }

    unsafe fn dd_copy_lp_solution(lp: *mut c_void) -> *mut c_void {
        unsafe { sys::gmpfloat::dd_CopyLPSolution(lp.cast()).cast() }
    }

    unsafe fn dd_free_lp_solution(sol: *mut c_void) {
        unsafe {
            sys::gmpfloat::dd_FreeLPSolution(sol.cast());
        }
    }

    unsafe fn lp_solution_optvalue_ptr(sol: *mut c_void) -> *const c_void {
        unsafe {
            let sol = sol.cast::<sys::gmpfloat::dd_lpsolution>();
            (&(*sol).optvalue as *const sys::gmpfloat::mytype).cast()
        }
    }

    unsafe fn write_mytype_real(target: *mut c_void, value: f64) {
        unsafe {
            sys::gmpfloat::__gmpf_set_d(target.cast(), value);
        }
    }

    unsafe fn write_mytype_int(target: *mut c_void, value: c_long) {
        unsafe {
            sys::gmpfloat::__gmpf_set_si(target.cast(), value);
        }
    }

    unsafe fn read_mytype_real(source: *const c_void) -> f64 {
        unsafe { sys::gmpfloat::__gmpf_get_d(source.cast()) }
    }

    unsafe fn write_mytype(target: *mut c_void, value: &Self) {
        unsafe {
            sys::gmpfloat::ddd_set(target.cast(), value.inner.as_ptr().cast_mut());
        }
    }

    unsafe fn read_mytype(source: *const c_void) -> Self {
        unsafe {
            Self::ensure_initialized();
            let mut inner = MaybeUninit::<sys::gmpfloat::mytype>::uninit();
            sys::gmpfloat::ddd_init(inner.as_mut_ptr().cast());
            sys::gmpfloat::ddd_set(
                inner.as_mut_ptr().cast(),
                source.cast::<sys::gmpfloat::__mpf_struct>().cast_mut(),
            );
            CddFloat {
                inner: inner.assume_init(),
            }
        }
    }
}

#[cfg(feature = "gmprational")]
impl sealed::Sealed for CddRational {}

#[cfg(feature = "gmprational")]
impl CddNumber for CddRational {
    const MYTYPE_SIZE: usize = mem::size_of::<sys::gmprational::mytype>();
    const DEFAULT_NUMBER_TYPE: NumberType = NumberType::Rational;

    fn ensure_initialized() {
        static INIT: Once = Once::new();
        INIT.call_once(|| unsafe {
            sys::gmprational::dd_set_global_constants();
        });
    }

    fn dd_no_error() -> u32 {
        sys::gmprational::dd_ErrorType_dd_NoError
    }

    unsafe fn dd_create_matrix(rows: c_long, cols: c_long) -> *mut c_void {
        unsafe { sys::gmprational::dd_CreateMatrix(rows, cols).cast() }
    }

    unsafe fn dd_copy_matrix(matrix: *mut c_void) -> *mut c_void {
        unsafe { sys::gmprational::dd_CopyMatrix(matrix.cast()).cast() }
    }

    unsafe fn dd_free_matrix(matrix: *mut c_void) {
        unsafe {
            sys::gmprational::dd_FreeMatrix(matrix.cast());
        }
    }

    unsafe fn dd_matrix_append_to(dst: *mut *mut c_void, src: *mut c_void) -> i32 {
        unsafe { sys::gmprational::dd_MatrixAppendTo(dst.cast(), src.cast()) }
    }

    unsafe fn dd_append_matrix(a: *mut c_void, b: *mut c_void) -> *mut c_void {
        unsafe { sys::gmprational::dd_AppendMatrix(a.cast(), b.cast()).cast() }
    }

    unsafe fn dd_matrix_row_remove(matrix: *mut *mut c_void, row: c_long) -> i32 {
        unsafe { sys::gmprational::dd_MatrixRowRemove(matrix.cast(), row) }
    }

    unsafe fn dd_initialize_arow(cols: c_long, out: *mut *mut c_void) {
        unsafe {
            sys::gmprational::dd_InitializeArow(cols, out.cast());
        }
    }

    unsafe fn dd_free_arow(cols: c_long, arow: *mut c_void) {
        unsafe {
            sys::gmprational::dd_FreeArow(cols, arow.cast());
        }
    }

    unsafe fn dd_redundant(
        matrix: *mut c_void,
        row: c_long,
        cert: *mut c_void,
        err: *mut u32,
    ) -> i32 {
        unsafe { sys::gmprational::dd_Redundant(matrix.cast(), row, cert.cast(), err.cast()) }
    }

    unsafe fn dd_redundant_rows(matrix: *mut c_void, err: *mut u32) -> *mut libc::c_ulong {
        unsafe { sys::gmprational::dd_RedundantRows(matrix.cast(), err.cast()).cast() }
    }

    unsafe fn dd_matrix_canonicalize(
        matrix: *mut *mut c_void,
        impl_lin: *mut *mut libc::c_ulong,
        redset: *mut *mut libc::c_ulong,
        newpos: *mut *mut c_long,
        err: *mut u32,
    ) -> i32 {
        unsafe {
            sys::gmprational::dd_MatrixCanonicalize(
                matrix.cast(),
                impl_lin.cast(),
                redset.cast(),
                newpos.cast(),
                err.cast(),
            )
        }
    }

    unsafe fn set_groundsize(set: *mut libc::c_ulong) -> c_long {
        unsafe { sys::gmprational::set_groundsize(set.cast()) }
    }

    unsafe fn set_member(elem: c_long, set: *mut libc::c_ulong) -> i32 {
        unsafe { sys::gmprational::set_member(elem, set.cast()) }
    }

    unsafe fn set_free(set: *mut libc::c_ulong) {
        unsafe {
            sys::gmprational::set_free(set.cast());
        }
    }

    unsafe fn dd_matrix2poly(matrix: *mut c_void, err: *mut u32) -> *mut c_void {
        unsafe { sys::gmprational::dd_DDMatrix2Poly(matrix.cast(), err.cast()).cast() }
    }

    unsafe fn dd_free_polyhedra(poly: *mut c_void) {
        unsafe {
            sys::gmprational::dd_FreePolyhedra(poly.cast());
        }
    }

    unsafe fn dd_copy_inequalities(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::gmprational::dd_CopyInequalities(poly.cast()).cast() }
    }

    unsafe fn dd_copy_generators(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::gmprational::dd_CopyGenerators(poly.cast()).cast() }
    }

    unsafe fn dd_copy_adjacency(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::gmprational::dd_CopyAdjacency(poly.cast()).cast() }
    }

    unsafe fn dd_copy_input_adjacency(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::gmprational::dd_CopyInputAdjacency(poly.cast()).cast() }
    }

    unsafe fn dd_copy_incidence(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::gmprational::dd_CopyIncidence(poly.cast()).cast() }
    }

    unsafe fn dd_copy_input_incidence(poly: *mut c_void) -> *mut c_void {
        unsafe { sys::gmprational::dd_CopyInputIncidence(poly.cast()).cast() }
    }

    unsafe fn dd_append_matrix2poly(poly: *mut *mut c_void, rows: *mut c_void) -> i32 {
        unsafe { sys::gmprational::dd_AppendMatrix2Poly(poly.cast(), rows.cast()) }
    }

    unsafe fn dd_free_set_family(family: *mut c_void) {
        unsafe {
            sys::gmprational::dd_FreeSetFamily(family.cast());
        }
    }

    unsafe fn dd_matrix2lp(matrix: *mut c_void, err: *mut u32) -> *mut c_void {
        unsafe { sys::gmprational::dd_Matrix2LP(matrix.cast(), err.cast()).cast() }
    }

    unsafe fn dd_free_lp_data(lp: *mut c_void) {
        unsafe {
            sys::gmprational::dd_FreeLPData(lp.cast());
        }
    }

    unsafe fn dd_lp_solve_dual_simplex(lp: *mut c_void, err: *mut u32) -> i32 {
        unsafe {
            sys::gmprational::dd_LPSolve(
                lp.cast(),
                sys::gmprational::dd_LPSolverType_dd_DualSimplex,
                err.cast(),
            )
        }
    }

    unsafe fn lp_status_raw(lp: *mut c_void) -> u32 {
        unsafe { (*lp.cast::<sys::gmprational::dd_lpdata>()).LPS }
    }

    unsafe fn dd_copy_lp_solution(lp: *mut c_void) -> *mut c_void {
        unsafe { sys::gmprational::dd_CopyLPSolution(lp.cast()).cast() }
    }

    unsafe fn dd_free_lp_solution(sol: *mut c_void) {
        unsafe {
            sys::gmprational::dd_FreeLPSolution(sol.cast());
        }
    }

    unsafe fn lp_solution_optvalue_ptr(sol: *mut c_void) -> *const c_void {
        unsafe {
            let sol = sol.cast::<sys::gmprational::dd_lpsolution>();
            (&(*sol).optvalue as *const sys::gmprational::mytype).cast()
        }
    }

    unsafe fn write_mytype_real(target: *mut c_void, value: f64) {
        unsafe {
            sys::gmprational::__gmpq_set_d(target.cast(), value);
        }
    }

    unsafe fn write_mytype_int(target: *mut c_void, value: c_long) {
        unsafe {
            sys::gmprational::__gmpq_set_si(target.cast(), value, 1);
        }
    }

    unsafe fn read_mytype_real(source: *const c_void) -> f64 {
        unsafe { sys::gmprational::__gmpq_get_d(source.cast()) }
    }

    unsafe fn write_mytype(target: *mut c_void, value: &Self) {
        unsafe {
            sys::gmprational::ddd_set(target.cast(), value.inner.as_ptr().cast_mut());
        }
    }

    unsafe fn read_mytype(source: *const c_void) -> Self {
        unsafe {
            Self::ensure_initialized();
            let mut inner = MaybeUninit::<sys::gmprational::mytype>::uninit();
            sys::gmprational::ddd_init(inner.as_mut_ptr().cast());
            sys::gmprational::ddd_set(
                inner.as_mut_ptr().cast(),
                source.cast::<sys::gmprational::__mpq_struct>().cast_mut(),
            );
            CddRational {
                inner: inner.assume_init(),
            }
        }
    }
}

#[cfg(feature = "gmp")]
impl CddFloat {
    pub fn to_f64(&self) -> f64 {
        <Self as CddNumber>::ensure_initialized();
        unsafe { sys::gmpfloat::__gmpf_get_d(self.inner.as_ptr()) }
    }
}

#[cfg(feature = "gmp")]
impl Clone for CddFloat {
    fn clone(&self) -> Self {
        <Self as CddNumber>::ensure_initialized();
        unsafe {
            let mut inner = MaybeUninit::<sys::gmpfloat::mytype>::uninit();
            sys::gmpfloat::ddd_init(inner.as_mut_ptr().cast());
            sys::gmpfloat::ddd_set(inner.as_mut_ptr().cast(), self.inner.as_ptr().cast_mut());
            CddFloat {
                inner: inner.assume_init(),
            }
        }
    }
}

#[cfg(feature = "gmp")]
impl Drop for CddFloat {
    fn drop(&mut self) {
        <Self as CddNumber>::ensure_initialized();
        unsafe {
            sys::gmpfloat::ddd_clear(self.inner.as_mut_ptr());
        }
    }
}

#[cfg(feature = "gmp")]
impl From<f64> for CddFloat {
    fn from(value: f64) -> Self {
        <Self as CddNumber>::ensure_initialized();
        unsafe {
            let mut inner = MaybeUninit::<sys::gmpfloat::mytype>::uninit();
            sys::gmpfloat::ddd_init(inner.as_mut_ptr().cast());
            sys::gmpfloat::__gmpf_set_d(inner.as_mut_ptr().cast(), value);
            CddFloat {
                inner: inner.assume_init(),
            }
        }
    }
}

#[cfg(feature = "gmprational")]
impl CddRational {
    pub fn to_f64(&self) -> f64 {
        <Self as CddNumber>::ensure_initialized();
        unsafe { sys::gmprational::__gmpq_get_d(self.inner.as_ptr()) }
    }
}

#[cfg(feature = "gmprational")]
impl Clone for CddRational {
    fn clone(&self) -> Self {
        <Self as CddNumber>::ensure_initialized();
        unsafe {
            let mut inner = MaybeUninit::<sys::gmprational::mytype>::uninit();
            sys::gmprational::ddd_init(inner.as_mut_ptr().cast());
            sys::gmprational::ddd_set(inner.as_mut_ptr().cast(), self.inner.as_ptr().cast_mut());
            CddRational {
                inner: inner.assume_init(),
            }
        }
    }
}

#[cfg(feature = "gmprational")]
impl Drop for CddRational {
    fn drop(&mut self) {
        <Self as CddNumber>::ensure_initialized();
        unsafe {
            sys::gmprational::ddd_clear(self.inner.as_mut_ptr());
        }
    }
}

#[cfg(feature = "gmprational")]
impl From<f64> for CddRational {
    fn from(value: f64) -> Self {
        <Self as CddNumber>::ensure_initialized();
        unsafe {
            let mut inner = MaybeUninit::<sys::gmprational::mytype>::uninit();
            sys::gmprational::ddd_init(inner.as_mut_ptr().cast());
            sys::gmprational::__gmpq_set_d(inner.as_mut_ptr().cast(), value);
            CddRational {
                inner: inner.assume_init(),
            }
        }
    }
}

#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub enum CddErrorCode {
    DimensionTooLarge,
    ImproperInputFormat,
    NegativeMatrixSize,
    EmptyVRepresentation,
    EmptyHRepresentation,
    EmptyRepresentation,
    InputFileNotFound,
    OutputFileNotOpen,
    NoLpObjective,
    NoRealNumberSupport,
    NotAvailForH,
    NotAvailForV,
    CannotHandleLinearity,
    RowIndexOutOfRange,
    ColIndexOutOfRange,
    LpCycling,
    NumericallyInconsistent,
    NoError,
}

impl CddErrorCode {
    pub fn from_raw(raw: u32) -> Self {
        match raw {
            0 => Self::DimensionTooLarge,
            1 => Self::ImproperInputFormat,
            2 => Self::NegativeMatrixSize,
            3 => Self::EmptyVRepresentation,
            4 => Self::EmptyHRepresentation,
            5 => Self::EmptyRepresentation,
            6 => Self::InputFileNotFound,
            7 => Self::OutputFileNotOpen,
            8 => Self::NoLpObjective,
            9 => Self::NoRealNumberSupport,
            10 => Self::NotAvailForH,
            11 => Self::NotAvailForV,
            12 => Self::CannotHandleLinearity,
            13 => Self::RowIndexOutOfRange,
            14 => Self::ColIndexOutOfRange,
            15 => Self::LpCycling,
            16 => Self::NumericallyInconsistent,
            17 => Self::NoError,
            other => panic!("unknown dd_ErrorType value {other}"),
        }
    }

    pub fn as_raw(self) -> u32 {
        match self {
            Self::DimensionTooLarge => 0,
            Self::ImproperInputFormat => 1,
            Self::NegativeMatrixSize => 2,
            Self::EmptyVRepresentation => 3,
            Self::EmptyHRepresentation => 4,
            Self::EmptyRepresentation => 5,
            Self::InputFileNotFound => 6,
            Self::OutputFileNotOpen => 7,
            Self::NoLpObjective => 8,
            Self::NoRealNumberSupport => 9,
            Self::NotAvailForH => 10,
            Self::NotAvailForV => 11,
            Self::CannotHandleLinearity => 12,
            Self::RowIndexOutOfRange => 13,
            Self::ColIndexOutOfRange => 14,
            Self::LpCycling => 15,
            Self::NumericallyInconsistent => 16,
            Self::NoError => 17,
        }
    }

    fn cddlib_name(self) -> &'static str {
        match self {
            Self::DimensionTooLarge => "dd_DimensionTooLarge",
            Self::ImproperInputFormat => "dd_ImproperInputFormat",
            Self::NegativeMatrixSize => "dd_NegativeMatrixSize",
            Self::EmptyVRepresentation => "dd_EmptyVrepresentation",
            Self::EmptyHRepresentation => "dd_EmptyHrepresentation",
            Self::EmptyRepresentation => "dd_EmptyRepresentation",
            Self::InputFileNotFound => "dd_IFileNotFound",
            Self::OutputFileNotOpen => "dd_OFileNotOpen",
            Self::NoLpObjective => "dd_NoLPObjective",
            Self::NoRealNumberSupport => "dd_NoRealNumberSupport",
            Self::NotAvailForH => "dd_NotAvailForH",
            Self::NotAvailForV => "dd_NotAvailForV",
            Self::CannotHandleLinearity => "dd_CannotHandleLinearity",
            Self::RowIndexOutOfRange => "dd_RowIndexOutOfRange",
            Self::ColIndexOutOfRange => "dd_ColIndexOutOfRange",
            Self::LpCycling => "dd_LPCycling",
            Self::NumericallyInconsistent => "dd_NumericallyInconsistent",
            Self::NoError => "dd_NoError",
        }
    }
}

impl std::fmt::Display for CddErrorCode {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(f, "{} (code={})", self.cddlib_name(), self.as_raw())
    }
}

#[derive(Debug, Error)]
pub enum CddError {
    #[error("cddlib returned {0}")]
    Cdd(CddErrorCode),
    #[error("cddlib returned a null pointer")]
    NullPointer,
    #[error("LP error")]
    LpError,
    #[error("LP did not have an optimal solution (status={0:?})")]
    LpStatus(LpStatus),
    #[error("cddlib operation returned failure")]
    OpFailed,
}

#[derive(Debug, Error)]
pub enum CddWrapperError {
    #[error(transparent)]
    Cdd(#[from] CddError),
    #[error("invalid matrix dimensions (rows={rows}, cols={cols})")]
    InvalidMatrix { rows: usize, cols: usize },
}

pub type CddResult<T> = std::result::Result<T, CddWrapperError>;

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum Representation {
    Inequality,
    Generator,
}

impl Representation {
    fn to_raw(self) -> u32 {
        match self {
            Representation::Inequality => 1,
            Representation::Generator => 2,
        }
    }

    fn from_raw(raw: u32) -> Self {
        match raw {
            1 => Representation::Inequality,
            2 => Representation::Generator,
            other => panic!("unknown dd_RepresentationType value {other}"),
        }
    }
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum NumberType {
    Real,
    Rational,
}

impl NumberType {
    fn to_raw(self) -> u32 {
        match self {
            NumberType::Real => 1,
            NumberType::Rational => 2,
        }
    }

    fn from_raw(raw: u32) -> Self {
        match raw {
            1 => NumberType::Real,
            2 => NumberType::Rational,
            other => panic!("unknown dd_NumberType value {other}"),
        }
    }
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum LpObjective {
    None,
    Maximize,
    Minimize,
}

impl LpObjective {
    fn to_raw(self) -> u32 {
        match self {
            LpObjective::None => 0,
            LpObjective::Maximize => 1,
            LpObjective::Minimize => 2,
        }
    }
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum LpStatus {
    Undecided,
    Optimal,
    Inconsistent,
    DualInconsistent,
    StructuralInconsistent,
    StructuralDualInconsistent,
    Unbounded,
    DualUnbounded,
}

impl LpStatus {
    fn from_raw(raw: u32) -> Self {
        match raw {
            0 => LpStatus::Undecided,
            1 => LpStatus::Optimal,
            2 => LpStatus::Inconsistent,
            3 => LpStatus::DualInconsistent,
            4 => LpStatus::StructuralInconsistent,
            5 => LpStatus::StructuralDualInconsistent,
            6 => LpStatus::Unbounded,
            7 => LpStatus::DualUnbounded,
            other => panic!("unknown dd_LPStatusType value {other}"),
        }
    }
}

#[repr(C)]
struct RawMatrixData {
    rowsize: c_long,
    linset: *mut libc::c_ulong,
    colsize: c_long,
    representation: u32,
    numbtype: u32,
    matrix: *mut *mut c_void,
    objective: u32,
    rowvec: *mut c_void,
}

#[repr(C)]
struct RawSetFamily {
    famsize: c_long,
    setsize: c_long,
    set: *mut *mut libc::c_ulong,
}

#[derive(Debug)]
pub struct Matrix<N: CddNumber = DefaultNumber> {
    ptr: NonNull<RawMatrixData>,
    _marker: PhantomData<N>,
    _no_send_sync: PhantomData<Rc<()>>,
}

#[derive(Debug)]
pub struct CanonicalForm<N: CddNumber = DefaultNumber> {
    pub matrix: Matrix<N>,
    pub implicit_linearity: Vec<usize>,
    pub redundant_rows: Vec<usize>,
    pub positions: Vec<isize>,
}

impl<N: CddNumber> Matrix<N> {
    pub fn new(
        rows: usize,
        cols: usize,
        repr: Representation,
        num_type: NumberType,
    ) -> CddResult<Self> {
        N::ensure_initialized();

        if rows == 0 || cols == 0 {
            return Err(CddWrapperError::InvalidMatrix { rows, cols });
        }

        let rowrange =
            c_long::try_from(rows).map_err(|_| CddWrapperError::InvalidMatrix { rows, cols })?;
        let colrange =
            c_long::try_from(cols).map_err(|_| CddWrapperError::InvalidMatrix { rows, cols })?;

        let ptr = unsafe { N::dd_create_matrix(rowrange, colrange) };
        let mut ptr = NonNull::new(ptr.cast::<RawMatrixData>()).ok_or(CddError::NullPointer)?;

        unsafe {
            ptr.as_mut().representation = repr.to_raw();
            ptr.as_mut().numbtype = num_type.to_raw();
        }

        Ok(Matrix {
            ptr,
            _marker: PhantomData,
            _no_send_sync: PhantomData,
        })
    }

    pub fn rows(&self) -> usize {
        unsafe { self.ptr.as_ref().rowsize as usize }
    }

    pub fn cols(&self) -> usize {
        unsafe { self.ptr.as_ref().colsize as usize }
    }

    pub fn representation(&self) -> Representation {
        unsafe { Representation::from_raw(self.ptr.as_ref().representation) }
    }

    pub fn number_type(&self) -> NumberType {
        unsafe { NumberType::from_raw(self.ptr.as_ref().numbtype) }
    }

    pub fn as_raw(&self) -> *mut c_void {
        self.ptr.as_ptr().cast()
    }

    /// # Safety
    ///
    /// - `ptr` must be a non-null `dd_MatrixPtr` compatible with the backend `N`.
    /// - The pointer must be valid for the lifetime of this wrapper and must be
    ///   owned by this wrapper (it will be freed on drop).
    pub unsafe fn from_raw(ptr: *mut c_void) -> Self {
        Matrix {
            ptr: NonNull::new(ptr.cast()).expect("from_raw called with null dd_MatrixPtr"),
            _marker: PhantomData,
            _no_send_sync: PhantomData,
        }
    }

    pub fn clone_cdd(&self) -> CddResult<Self> {
        N::ensure_initialized();

        let ptr = unsafe { N::dd_copy_matrix(self.as_raw()) };
        let ptr = NonNull::new(ptr).ok_or(CddError::NullPointer)?;
        Ok(unsafe { Self::from_raw(ptr.as_ptr()) })
    }

    pub fn set(&mut self, row: usize, col: usize, value: &N) {
        assert!(row < self.rows());
        assert!(col < self.cols());
        let cell = unsafe { matrix_cell::<N>(self.ptr, row, col) };
        unsafe {
            N::write_mytype(cell, value);
        }
    }

    pub fn get(&self, row: usize, col: usize) -> N {
        assert!(row < self.rows());
        assert!(col < self.cols());
        let cell = unsafe { matrix_cell::<N>(self.ptr, row, col) };
        unsafe { N::read_mytype(cell.cast_const()) }
    }

    pub fn set_real(&mut self, row: usize, col: usize, value: f64) {
        assert!(row < self.rows());
        assert!(col < self.cols());
        let cell = unsafe { matrix_cell::<N>(self.ptr, row, col) };
        unsafe {
            N::write_mytype_real(cell, value);
        }
    }

    pub fn get_real(&self, row: usize, col: usize) -> f64 {
        assert!(row < self.rows());
        assert!(col < self.cols());
        let cell = unsafe { matrix_cell::<N>(self.ptr, row, col) };
        unsafe { N::read_mytype_real(cell.cast_const()) }
    }

    pub fn set_generator_type(&mut self, row: usize, is_vertex: bool) {
        assert!(row < self.rows());
        let cell = unsafe { matrix_cell::<N>(self.ptr, row, 0) };
        let value = if is_vertex { 1 } else { 0 };
        unsafe {
            N::write_mytype_int(cell, value);
        }
    }

    pub fn set_objective_real(&mut self, coeffs: &[f64]) {
        assert_eq!(coeffs.len(), self.cols());
        let rowvec = unsafe { self.ptr.as_ref().rowvec };
        assert!(!rowvec.is_null(), "matrix objective rowvec is null");
        for (j, &c) in coeffs.iter().enumerate() {
            let cell = unsafe { rowvec_cell::<N>(rowvec, j) };
            unsafe {
                N::write_mytype_real(cell, c);
            }
        }
    }

    pub fn from_vertices<const D: usize>(vertices: &[[N; D]]) -> CddResult<Self> {
        if vertices.is_empty() {
            return Err(CddWrapperError::InvalidMatrix {
                rows: 0,
                cols: D + 1,
            });
        }

        let mut m = Self::new(
            vertices.len(),
            D + 1,
            Representation::Generator,
            N::DEFAULT_NUMBER_TYPE,
        )?;
        for (i, v) in vertices.iter().enumerate() {
            m.set_generator_type(i, true);
            for (j, coord) in v.iter().enumerate() {
                m.set(i, j + 1, coord);
            }
        }
        Ok(m)
    }

    /// Build a generator matrix from a dynamic list of vertices.
    ///
    /// The ambient dimension is determined from the first vertex. All vertices
    /// must have the same length.
    pub fn from_vertex_rows(vertices: &[Vec<N>]) -> CddResult<Self> {
        let Some(first) = vertices.first() else {
            return Err(CddWrapperError::InvalidMatrix { rows: 0, cols: 0 });
        };

        let dim = first.len();
        if dim == 0 {
            return Err(CddWrapperError::InvalidMatrix {
                rows: vertices.len(),
                cols: 0,
            });
        }

        if vertices.iter().skip(1).any(|v| v.len() != dim) {
            return Err(CddWrapperError::InvalidMatrix {
                rows: vertices.len(),
                cols: dim + 1,
            });
        }

        let mut m = Self::new(
            vertices.len(),
            dim + 1,
            Representation::Generator,
            N::DEFAULT_NUMBER_TYPE,
        )?;
        for (row, coords) in vertices.iter().enumerate() {
            m.set_generator_type(row, true);
            for (col, coord) in coords.iter().enumerate() {
                m.set(row, col + 1, coord);
            }
        }
        Ok(m)
    }

    pub fn append_rows_in_place(&mut self, rows: &Matrix<N>) -> CddResult<()> {
        N::ensure_initialized();

        if self.cols() != rows.cols() {
            return Err(CddWrapperError::InvalidMatrix {
                rows: self.rows() + rows.rows(),
                cols: rows.cols(),
            });
        }

        if self.representation() != rows.representation()
            || self.number_type() != rows.number_type()
        {
            return Err(CddWrapperError::InvalidMatrix {
                rows: self.rows() + rows.rows(),
                cols: rows.cols(),
            });
        }

        let mut raw = self.as_raw();
        let ok = unsafe { N::dd_matrix_append_to(&mut raw, rows.as_raw()) };
        if ok == 0 {
            return Err(CddError::OpFailed.into());
        }

        let raw = NonNull::new(raw.cast()).ok_or(CddError::NullPointer)?;
        if raw != self.ptr {
            self.ptr = raw;
        }
        Ok(())
    }

    pub fn append_rows(&self, rows: &Matrix<N>) -> CddResult<Self> {
        N::ensure_initialized();

        if self.cols() != rows.cols() {
            return Err(CddWrapperError::InvalidMatrix {
                rows: self.rows() + rows.rows(),
                cols: rows.cols(),
            });
        }

        if self.representation() != rows.representation()
            || self.number_type() != rows.number_type()
        {
            return Err(CddWrapperError::InvalidMatrix {
                rows: self.rows() + rows.rows(),
                cols: rows.cols(),
            });
        }

        let ptr = unsafe { N::dd_append_matrix(self.as_raw(), rows.as_raw()) };
        let ptr = NonNull::new(ptr).ok_or(CddError::NullPointer)?;
        Ok(unsafe { Self::from_raw(ptr.as_ptr()) })
    }

    pub fn append_row(&self, coords: &[N], is_vertex: bool) -> CddResult<Self> {
        if coords.len() + 1 != self.cols() {
            return Err(CddWrapperError::InvalidMatrix {
                rows: self.rows() + 1,
                cols: coords.len() + 1,
            });
        }

        let mut row = Self::new(1, self.cols(), self.representation(), self.number_type())?;
        row.set_generator_type(0, is_vertex);
        for (col, val) in coords.iter().enumerate() {
            row.set(0, col + 1, val);
        }

        self.append_rows(&row)
    }

    pub fn remove_row(&mut self, row: usize) -> CddResult<()> {
        N::ensure_initialized();

        if row >= self.rows() {
            return Err(CddWrapperError::InvalidMatrix {
                rows: self.rows(),
                cols: self.cols(),
            });
        }

        let row_idx = c_long::try_from(row + 1).map_err(|_| CddWrapperError::InvalidMatrix {
            rows: self.rows(),
            cols: self.cols(),
        })?;

        let mut raw = self.as_raw();
        let ok = unsafe { N::dd_matrix_row_remove(&mut raw, row_idx) };
        if ok == 0 {
            return Err(CddError::OpFailed.into());
        }

        let raw = NonNull::new(raw.cast()).ok_or(CddError::NullPointer)?;
        if raw != self.ptr {
            self.ptr = raw;
        }
        Ok(())
    }

    pub fn is_row_redundant(&self, row: usize) -> CddResult<bool> {
        N::ensure_initialized();

        if row >= self.rows() {
            return Err(CddWrapperError::InvalidMatrix {
                rows: self.rows(),
                cols: self.cols(),
            });
        }

        let row_idx = c_long::try_from(row + 1).map_err(|_| CddWrapperError::InvalidMatrix {
            rows: self.rows(),
            cols: self.cols(),
        })?;

        let mut err = N::dd_no_error();
        let mut cert: *mut c_void = ptr::null_mut();

        let cert_cols = match self.representation() {
            Representation::Generator => self.cols() + 1,
            Representation::Inequality => self.cols(),
        };

        let colrange = c_long::try_from(cert_cols).map_err(|_| CddWrapperError::InvalidMatrix {
            rows: self.rows(),
            cols: self.cols(),
        })?;

        unsafe {
            N::dd_initialize_arow(colrange, &mut cert);
        }

        let redundant = unsafe { N::dd_redundant(self.as_raw(), row_idx, cert, &mut err) };

        unsafe {
            N::dd_free_arow(colrange, cert);
        }

        let err = CddErrorCode::from_raw(err);
        if err != CddErrorCode::NoError {
            return Err(CddError::Cdd(err).into());
        }
        Ok(redundant != 0)
    }

    pub fn redundant_rows(&self) -> CddResult<Vec<usize>> {
        N::ensure_initialized();

        let mut err = N::dd_no_error();
        let set = unsafe { N::dd_redundant_rows(self.as_raw(), &mut err) };
        let err = CddErrorCode::from_raw(err);
        if err != CddErrorCode::NoError {
            return Err(CddError::Cdd(err).into());
        }
        let set = NonNull::new(set).ok_or(CddError::NullPointer)?;

        let ground = unsafe { N::set_groundsize(set.as_ptr()) };
        let ground = usize::try_from(ground).map_err(|_| CddWrapperError::InvalidMatrix {
            rows: self.rows(),
            cols: self.cols(),
        })?;

        let mut rows = Vec::new();
        for i in 1..=ground {
            let member = unsafe { N::set_member(i as c_long, set.as_ptr()) };
            if member != 0 {
                rows.push(i - 1);
            }
        }

        unsafe {
            N::set_free(set.as_ptr());
        }

        Ok(rows)
    }

    pub fn canonicalize(&self) -> CddResult<CanonicalForm<N>> {
        N::ensure_initialized();

        let clone = self.clone_cdd()?;
        let mut raw = clone.as_raw();
        let mut impl_lin: *mut libc::c_ulong = ptr::null_mut();
        let mut redset: *mut libc::c_ulong = ptr::null_mut();
        let mut newpos: *mut c_long = ptr::null_mut();
        let mut err = N::dd_no_error();

        let ok = unsafe {
            N::dd_matrix_canonicalize(&mut raw, &mut impl_lin, &mut redset, &mut newpos, &mut err)
        };
        let err = CddErrorCode::from_raw(err);
        if ok == 0 || err != CddErrorCode::NoError {
            unsafe {
                if !impl_lin.is_null() {
                    N::set_free(impl_lin);
                }
                if !redset.is_null() {
                    N::set_free(redset);
                }
                if !newpos.is_null() {
                    libc::free(newpos.cast());
                }
            }
            if ok == 0 && err == CddErrorCode::NoError {
                return Err(CddError::OpFailed.into());
            }
            return Err(CddError::Cdd(err).into());
        }

        let canon = unsafe { Self::from_raw(raw) };
        let rows = canon.rows();
        let implicit = unsafe { read_rowset(impl_lin, rows) };
        let redundant = unsafe { read_rowset(redset, rows) };
        let positions = unsafe { read_rowindex(newpos, rows) };

        unsafe {
            N::set_free(impl_lin);
            N::set_free(redset);
            libc::free(newpos.cast());
        }
        mem::forget(clone);

        Ok(CanonicalForm {
            matrix: canon,
            implicit_linearity: implicit,
            redundant_rows: redundant,
            positions,
        })
    }
}

impl<N: CddNumber> Drop for Matrix<N> {
    fn drop(&mut self) {
        unsafe {
            N::dd_free_matrix(self.as_raw());
        }
    }
}

#[derive(Debug)]
pub struct SetFamily<N: CddNumber = DefaultNumber> {
    ptr: NonNull<RawSetFamily>,
    _marker: PhantomData<N>,
    _no_send_sync: PhantomData<Rc<()>>,
}

impl<N: CddNumber> SetFamily<N> {
    /// # Safety
    ///
    /// - `ptr` must be a non-null `dd_SetFamilyPtr` compatible with the backend `N`.
    /// - The pointer must be valid for the lifetime of this wrapper and must be
    ///   owned by this wrapper (it will be freed on drop).
    pub unsafe fn from_raw(ptr: *mut c_void) -> Self {
        SetFamily {
            ptr: NonNull::new(ptr.cast()).expect("from_raw called with null dd_SetFamilyPtr"),
            _marker: PhantomData,
            _no_send_sync: PhantomData,
        }
    }

    pub fn as_raw(&self) -> *mut c_void {
        self.ptr.as_ptr().cast()
    }

    pub fn len(&self) -> usize {
        unsafe { self.ptr.as_ref().famsize as usize }
    }

    pub fn is_empty(&self) -> bool {
        self.len() == 0
    }

    pub fn universe_size(&self) -> usize {
        unsafe { self.ptr.as_ref().setsize as usize }
    }

    pub fn to_adjacency_lists(&self) -> Vec<Vec<usize>> {
        let n = self.len();
        let universe = self.universe_size();
        let mut result = Vec::with_capacity(n);

        for i in 0..n {
            let set = unsafe { *self.ptr.as_ref().set.add(i) };
            let neighbors = unsafe { set_members_bounded(set, universe) };
            result.push(neighbors);
        }

        result
    }
}

impl<N: CddNumber> Drop for SetFamily<N> {
    fn drop(&mut self) {
        unsafe {
            N::dd_free_set_family(self.as_raw());
        }
    }
}

#[derive(Debug)]
pub struct Polyhedron<N: CddNumber = DefaultNumber> {
    ptr: NonNull<c_void>,
    _marker: PhantomData<N>,
    _no_send_sync: PhantomData<Rc<()>>,
}

impl<N: CddNumber> Polyhedron<N> {
    pub fn from_matrix(m: &Matrix<N>) -> CddResult<Self> {
        N::ensure_initialized();
        let mut err = N::dd_no_error();
        let ptr = unsafe { N::dd_matrix2poly(m.as_raw(), &mut err) };
        let err = CddErrorCode::from_raw(err);
        if err != CddErrorCode::NoError {
            return Err(CddError::Cdd(err).into());
        }
        let ptr = NonNull::new(ptr).ok_or(CddError::NullPointer)?;
        Ok(Polyhedron {
            ptr,
            _marker: PhantomData,
            _no_send_sync: PhantomData,
        })
    }

    pub fn from_generators_matrix(m: &Matrix<N>) -> CddResult<Self> {
        Self::from_matrix(m)
    }

    pub fn from_inequalities_matrix(m: &Matrix<N>) -> CddResult<Self> {
        Self::from_matrix(m)
    }

    pub fn from_vertices<const D: usize>(vertices: &[[N; D]]) -> CddResult<Self> {
        let m = Matrix::<N>::from_vertices(vertices)?;
        Self::from_generators_matrix(&m)
    }

    pub fn from_vertex_rows(vertices: &[Vec<N>]) -> CddResult<Self> {
        let m = Matrix::<N>::from_vertex_rows(vertices)?;
        Self::from_generators_matrix(&m)
    }

    pub fn facets(&self) -> CddResult<Matrix<N>> {
        let ptr = unsafe { N::dd_copy_inequalities(self.ptr.as_ptr()) };
        let ptr = NonNull::new(ptr).ok_or(CddError::NullPointer)?;
        Ok(unsafe { Matrix::<N>::from_raw(ptr.as_ptr()) })
    }

    pub fn generators(&self) -> CddResult<Matrix<N>> {
        let ptr = unsafe { N::dd_copy_generators(self.ptr.as_ptr()) };
        let ptr = NonNull::new(ptr).ok_or(CddError::NullPointer)?;
        Ok(unsafe { Matrix::<N>::from_raw(ptr.as_ptr()) })
    }

    pub fn adjacency(&self) -> CddResult<SetFamily<N>> {
        let ptr = unsafe { N::dd_copy_adjacency(self.ptr.as_ptr()) };
        let ptr = NonNull::new(ptr).ok_or(CddError::NullPointer)?;
        Ok(unsafe { SetFamily::<N>::from_raw(ptr.as_ptr()) })
    }

    pub fn input_adjacency(&self) -> CddResult<SetFamily<N>> {
        let ptr = unsafe { N::dd_copy_input_adjacency(self.ptr.as_ptr()) };
        let ptr = NonNull::new(ptr).ok_or(CddError::NullPointer)?;
        Ok(unsafe { SetFamily::<N>::from_raw(ptr.as_ptr()) })
    }

    pub fn incidence(&self) -> CddResult<SetFamily<N>> {
        let ptr = unsafe { N::dd_copy_incidence(self.ptr.as_ptr()) };
        let ptr = NonNull::new(ptr).ok_or(CddError::NullPointer)?;
        Ok(unsafe { SetFamily::<N>::from_raw(ptr.as_ptr()) })
    }

    pub fn input_incidence(&self) -> CddResult<SetFamily<N>> {
        let ptr = unsafe { N::dd_copy_input_incidence(self.ptr.as_ptr()) };
        let ptr = NonNull::new(ptr).ok_or(CddError::NullPointer)?;
        Ok(unsafe { SetFamily::<N>::from_raw(ptr.as_ptr()) })
    }

    pub fn append_input_rows(&mut self, rows: &Matrix<N>) -> CddResult<()> {
        let mut raw = self.ptr.as_ptr();
        let ok = unsafe { N::dd_append_matrix2poly(&mut raw, rows.as_raw()) };
        if ok == 0 {
            return Err(CddError::OpFailed.into());
        }
        if raw != self.ptr.as_ptr() {
            self.ptr = NonNull::new(raw).ok_or(CddError::NullPointer)?;
        }
        Ok(())
    }
}

impl<N: CddNumber> Drop for Polyhedron<N> {
    fn drop(&mut self) {
        unsafe {
            N::dd_free_polyhedra(self.ptr.as_ptr());
        }
    }
}

#[derive(Debug)]
pub struct Lp<N: CddNumber = DefaultNumber> {
    ptr: NonNull<c_void>,
    _marker: PhantomData<N>,
    _no_send_sync: PhantomData<Rc<()>>,
}

#[derive(Debug)]
pub struct LpSolution<N: CddNumber = DefaultNumber> {
    ptr: NonNull<c_void>,
    _marker: PhantomData<N>,
    _no_send_sync: PhantomData<Rc<()>>,
}

impl<N: CddNumber> Lp<N> {
    pub fn from_matrix(matrix: &mut Matrix<N>, objective: LpObjective) -> CddResult<Self> {
        N::ensure_initialized();
        unsafe {
            matrix.ptr.as_mut().objective = objective.to_raw();
        }
        let mut err = N::dd_no_error();
        let ptr = unsafe { N::dd_matrix2lp(matrix.as_raw(), &mut err) };
        let err = CddErrorCode::from_raw(err);
        if err != CddErrorCode::NoError {
            return Err(CddError::Cdd(err).into());
        }
        let ptr = NonNull::new(ptr).ok_or(CddError::NullPointer)?;
        Ok(Lp {
            ptr,
            _marker: PhantomData,
            _no_send_sync: PhantomData,
        })
    }

    pub fn solve(&self) -> CddResult<LpSolution<N>> {
        N::ensure_initialized();

        let mut err = N::dd_no_error();
        let ok = unsafe { N::dd_lp_solve_dual_simplex(self.ptr.as_ptr(), &mut err) };
        let err = CddErrorCode::from_raw(err);
        if err != CddErrorCode::NoError {
            return Err(CddError::Cdd(err).into());
        }
        if ok == 0 {
            return Err(CddError::LpError.into());
        }
        let status = unsafe { N::lp_status_raw(self.ptr.as_ptr()) };
        let status = LpStatus::from_raw(status);
        if status != LpStatus::Optimal {
            return Err(CddError::LpStatus(status).into());
        }
        let sol_ptr = unsafe { N::dd_copy_lp_solution(self.ptr.as_ptr()) };
        let sol_ptr = NonNull::new(sol_ptr).ok_or(CddError::NullPointer)?;
        Ok(LpSolution {
            ptr: sol_ptr,
            _marker: PhantomData,
            _no_send_sync: PhantomData,
        })
    }
}

impl<N: CddNumber> Drop for Lp<N> {
    fn drop(&mut self) {
        unsafe {
            N::dd_free_lp_data(self.ptr.as_ptr());
        }
    }
}

impl<N: CddNumber> LpSolution<N> {
    pub fn opt_value_real(&self) -> f64 {
        unsafe {
            let opt_ptr = N::lp_solution_optvalue_ptr(self.ptr.as_ptr());
            N::read_mytype_real(opt_ptr)
        }
    }
}

impl<N: CddNumber> Drop for LpSolution<N> {
    fn drop(&mut self) {
        unsafe {
            N::dd_free_lp_solution(self.ptr.as_ptr());
        }
    }
}

pub fn width_in_direction_real<N: CddNumber>(
    poly: &Polyhedron<N>,
    direction: &[f64],
) -> CddResult<f64> {
    let mut h_min = poly.facets()?;
    let mut h_max = h_min.clone_cdd()?;

    if direction.len() + 1 != h_max.cols() {
        return Err(CddWrapperError::InvalidMatrix {
            rows: h_max.rows(),
            cols: h_max.cols(),
        });
    }

    let cols = h_max.cols();
    let mut coeffs = vec![0.0f64; cols];

    for (i, &u_i) in direction.iter().enumerate() {
        coeffs[i + 1] = u_i;
    }

    h_max.set_objective_real(&coeffs);
    let lp_max = Lp::<N>::from_matrix(&mut h_max, LpObjective::Maximize)?;
    let sol_max = lp_max.solve()?;
    let max_val = sol_max.opt_value_real();

    h_min.set_objective_real(&coeffs);
    let lp_min = Lp::<N>::from_matrix(&mut h_min, LpObjective::Minimize)?;
    let sol_min = lp_min.solve()?;
    let min_val = sol_min.opt_value_real();

    Ok(max_val - min_val)
}

#[cfg(feature = "f64")]
pub type MatrixF64 = Matrix<f64>;
#[cfg(feature = "f64")]
pub type PolyhedronF64 = Polyhedron<f64>;

#[cfg(feature = "gmp")]
pub type MatrixGmpFloat = Matrix<CddFloat>;
#[cfg(feature = "gmp")]
pub type PolyhedronGmpFloat = Polyhedron<CddFloat>;

#[cfg(feature = "gmprational")]
pub type MatrixGmpRational = Matrix<CddRational>;
#[cfg(feature = "gmprational")]
pub type PolyhedronGmpRational = Polyhedron<CddRational>;

impl<N: CddNumber> Matrix<N> {
    pub fn convert<M: CddNumber>(&self) -> CddResult<Matrix<M>> {
        convert_matrix_via_real::<N, M>(self)
    }
}

#[cfg(all(feature = "f64", feature = "gmp"))]
impl From<Matrix<f64>> for Matrix<CddFloat> {
    fn from(value: Matrix<f64>) -> Self {
        value
            .convert()
            .expect("Matrix<f64> -> Matrix<CddFloat> conversion failed")
    }
}

#[cfg(all(feature = "f64", feature = "gmprational"))]
impl From<Matrix<f64>> for Matrix<CddRational> {
    fn from(value: Matrix<f64>) -> Self {
        value
            .convert()
            .expect("Matrix<f64> -> Matrix<CddRational> conversion failed")
    }
}

#[cfg(all(feature = "f64", feature = "gmp"))]
impl From<Matrix<CddFloat>> for Matrix<f64> {
    fn from(value: Matrix<CddFloat>) -> Self {
        value
            .convert()
            .expect("Matrix<CddFloat> -> Matrix<f64> conversion failed")
    }
}

#[cfg(all(feature = "f64", feature = "gmprational"))]
impl From<Matrix<CddRational>> for Matrix<f64> {
    fn from(value: Matrix<CddRational>) -> Self {
        value
            .convert()
            .expect("Matrix<CddRational> -> Matrix<f64> conversion failed")
    }
}

#[cfg(all(feature = "gmp", feature = "gmprational"))]
impl From<Matrix<CddFloat>> for Matrix<CddRational> {
    fn from(value: Matrix<CddFloat>) -> Self {
        value
            .convert()
            .expect("Matrix<CddFloat> -> Matrix<CddRational> conversion failed")
    }
}

#[cfg(all(feature = "gmp", feature = "gmprational"))]
impl From<Matrix<CddRational>> for Matrix<CddFloat> {
    fn from(value: Matrix<CddRational>) -> Self {
        value
            .convert()
            .expect("Matrix<CddRational> -> Matrix<CddFloat> conversion failed")
    }
}

fn convert_matrix_via_real<Src: CddNumber, Dst: CddNumber>(
    src: &Matrix<Src>,
) -> CddResult<Matrix<Dst>> {
    let rows = src.rows();
    let cols = src.cols();
    let mut out = Matrix::<Dst>::new(rows, cols, src.representation(), Dst::DEFAULT_NUMBER_TYPE)?;

    for i in 0..rows {
        for j in 0..cols {
            out.set_real(i, j, src.get_real(i, j));
        }
    }

    unsafe {
        out.ptr.as_mut().objective = src.ptr.as_ref().objective;
    }

    let rowvec = unsafe { src.ptr.as_ref().rowvec };
    if !rowvec.is_null() {
        let mut coeffs = Vec::with_capacity(cols);
        for j in 0..cols {
            let cell = unsafe { rowvec_cell::<Src>(rowvec, j) };
            coeffs.push(unsafe { Src::read_mytype_real(cell.cast_const()) });
        }
        out.set_objective_real(&coeffs);
    }

    Ok(out)
}

unsafe fn matrix_cell<N: CddNumber>(
    matrix: NonNull<RawMatrixData>,
    row: usize,
    col: usize,
) -> *mut c_void {
    unsafe {
        let matrix_rows = matrix.as_ref().matrix;
        assert!(!matrix_rows.is_null(), "matrix rows pointer is null");
        let row_ptr = *matrix_rows.add(row);
        assert!(!row_ptr.is_null(), "matrix row pointer is null");
        row_ptr
            .cast::<u8>()
            .add(col * N::MYTYPE_SIZE)
            .cast::<c_void>()
    }
}

unsafe fn rowvec_cell<N: CddNumber>(rowvec: *mut c_void, col: usize) -> *mut c_void {
    unsafe {
        rowvec
            .cast::<u8>()
            .add(col * N::MYTYPE_SIZE)
            .cast::<c_void>()
    }
}

unsafe fn read_rowset(set: *mut libc::c_ulong, rows: usize) -> Vec<usize> {
    unsafe { set_members_bounded(set, rows) }
}

unsafe fn set_members_bounded(set: *mut libc::c_ulong, max_elems: usize) -> Vec<usize> {
    assert!(!set.is_null(), "set pointer is null");

    let ground = unsafe { *set } as usize;
    let bound = max_elems.min(ground);
    if bound == 0 {
        return Vec::new();
    }

    let bits_per_block = mem::size_of::<libc::c_ulong>() * 8;
    let blocks = (bound - 1) / bits_per_block + 1;

    let mut out = Vec::new();
    for block_idx in 0..blocks {
        let mut word = unsafe { *set.add(1 + block_idx) };
        let base = block_idx * bits_per_block;
        while word != 0 {
            let tz = word.trailing_zeros() as usize;
            let idx = base + tz;
            if idx < bound {
                out.push(idx);
            }
            word &= word - 1;
        }
    }

    out
}

unsafe fn read_rowindex(index: *mut c_long, rows: usize) -> Vec<isize> {
    unsafe {
        let mut out = Vec::with_capacity(rows);
        for i in 1..=rows {
            out.push(*index.add(i) as isize);
        }
        out
    }
}