pxdct 0.3.0

/*
 * // Copyright (c) Radzivon Bartoshyk 11/2025. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
#![cfg_attr(
    all(feature = "fcma", target_arch = "aarch64"),
    feature(stdarch_neon_fcma)
)]
#![allow(clippy::too_many_arguments, clippy::needless_range_loop)]
#[cfg(all(target_arch = "x86_64", feature = "avx"))]
mod avx;
mod bidirectional;
mod butterflies;
mod dct2;
mod dct3;
mod dct4;
mod dst2;
mod dst3;
mod dst3_butterfly;
mod factory_dct2;
mod factory_dct3;
mod factory_dct4;
mod factory_dst2;
mod factory_scaled_dct2;
mod identity;
mod mla;
#[cfg(all(target_arch = "aarch64", feature = "neon"))]
mod neon;
mod prime_factors;
mod pxdct_error;
mod scaling;
mod spectrum_mul;
mod transpose;
mod twiddles;
mod two_dims;
mod util;

use crate::dct2::Dct2Fft;
use crate::dct3::SplitRadixDst3;
use crate::dst2::Dst2Fft;
use crate::dst3::Dst3Fft;
use crate::dst3_butterfly::{
    Dst3Butterfly2, Dst3Butterfly3, Dst3Butterfly4, Dst3Butterfly5, Dst3Butterfly8, Dst3Butterfly16,
};
use crate::factory_dct2::Dct2Factory;
use crate::factory_dct3::Dct3Factory;
use crate::factory_dct4::Dct4Factory;
use crate::factory_dst2::Dst2Factory;
use crate::factory_scaled_dct2::ScaledDct2Factory;
use crate::identity::DctIdentity;
use crate::prime_factors::PrimeFactors;
use crate::scaling::ScalingInterceptor;
use crate::transpose::TransposeFactory;
use crate::two_dims::TwoDimensionalDct;
use crate::util::DctSample;
use dct2::power2_butterflies::{Dst2Butterfly2, Dst2Butterfly4};
use num_traits::AsPrimitive;
pub use pxdct_error::PxdctError;
use std::sync::{Arc, OnceLock};
pub use two_dims::MultidimensionalDctExecutor;
use zaft::FftDirection;

/// The main entry point for creating DCT (Discrete Cosine Transform) executors.
///
/// `Pxdct` provides convenient factory methods to construct optimized
/// executors for DCT-II and DCT-III transforms using single (`f32`) or
/// double (`f64`) precision. Each executor implements the [`PxdctExecutor`]
/// trait and can be used to perform an in-place DCT transform on a data slice.
pub struct Pxdct {}

macro_rules! make_dct2_butterflies {
    ($length: expr, $f_type: ident) => {
        if $length == 2 {
            return Ok($f_type::dct2_butterfly2());
        } else if $length == 3 {
            return Ok($f_type::dct2_butterfly3());
        } else if $length == 4 {
            return Ok($f_type::dct2_butterfly4());
        } else if $length == 5 {
            return Ok($f_type::dct2_butterfly5());
        } else if $length == 6 {
            return Ok($f_type::dct2_butterfly6());
        } else if $length == 7 {
            return Ok($f_type::dct2_butterfly7());
        } else if $length == 8 {
            return Ok($f_type::dct2_butterfly8());
        } else if $length == 9 {
            return Ok($f_type::dct2_butterfly9());
        } else if $length == 10 {
            return Ok($f_type::dct2_butterfly10());
        } else if $length == 11 {
            return Ok($f_type::dct2_butterfly11());
        } else if $length == 12 {
            return Ok($f_type::dct2_butterfly12());
        } else if $length == 13 {
            return Ok($f_type::dct2_butterfly13());
        } else if $length == 14 {
            return Ok($f_type::dct2_butterfly14());
        } else if $length == 15 {
            return Ok($f_type::dct2_butterfly15());
        } else if $length == 16 {
            return Ok($f_type::dct2_butterfly16());
        } else if $length == 17 {
            return Ok($f_type::dct2_butterfly17());
        } else if $length == 18 {
            return Ok($f_type::dct2_butterfly18());
        } else if $length == 19 {
            return Ok($f_type::dct2_butterfly19());
        } else if $length == 20 {
            return Ok($f_type::dct2_butterfly20());
        } else if $length == 21 {
            return Ok($f_type::dct2_butterfly21());
        } else if $length == 22 {
            return Ok($f_type::dct2_butterfly22());
        } else if $length == 23 {
            return Ok($f_type::dct2_butterfly23());
        } else if $length == 24 {
            return Ok($f_type::dct2_butterfly24());
        } else if $length == 25 {
            return Ok($f_type::dct2_butterfly25());
        } else if $length == 26 {
            return Ok($f_type::dct2_butterfly26());
        } else if $length == 27 {
            return Ok($f_type::dct2_butterfly27());
        } else if $length == 29 {
            return Ok($f_type::dct2_butterfly29());
        } else if $length == 30 {
            return Ok($f_type::dct2_butterfly30());
        } else if $length == 31 {
            return Ok($f_type::dct2_butterfly31());
        } else if $length == 32 {
            return Ok($f_type::dct2_butterfly32());
        } else if $length == 35 {
            return Ok($f_type::dct2_butterfly35());
        } else if $length == 36 {
            return Ok($f_type::dct2_butterfly36());
        } else if $length == 37 {
            return Ok($f_type::dct2_butterfly37());
        } else if $length == 42 {
            return Ok($f_type::dct2_butterfly42());
        } else if $length == 48 {
            return Ok($f_type::dct2_butterfly48());
        } else if $length == 49 {
            return Ok($f_type::dct2_butterfly49());
        } else if $length == 64 {
            return Ok($f_type::dct2_butterfly64());
        } else if $length == 81 {
            return Ok($f_type::dct2_butterfly81());
        } else if $length == 128 {
            return Ok($f_type::dct2_butterfly128());
        } else if $length == 216 {
            return Ok($f_type::dct2_butterfly216());
        } else if $length == 243 {
            return Ok($f_type::dct2_butterfly243());
        } else if $length == 256 {
            return Ok($f_type::dct2_butterfly256());
        } else if $length == 512 {
            return Ok($f_type::dct2_butterfly512());
        }
    };
}

impl Pxdct {
    fn dct2_strategy<T: DctSample + Dct2Factory + Send + Sync>(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<T> + Send + Sync>, PxdctError>
    where
        f64: AsPrimitive<T>,
    {
        if length == 0 {
            return Err(PxdctError::ZeroSizedDct);
        }
        if length == 1 {
            return Ok(Arc::new(DctIdentity::default()));
        }

        make_dct2_butterflies!(length, T);

        if length.is_power_of_two() && length > 2 {
            return T::split_radix(
                length,
                Pxdct::dct2_strategy(length / 2)?,
                Pxdct::dct2_strategy(length / 4)?,
            );
        }

        if length.is_multiple_of(6) {
            return T::mixed_radix6(length, Pxdct::dct2_strategy(length / 6)?);
        }

        if length.is_multiple_of(9) {
            return T::mixed_radix9(length, Pxdct::dct2_strategy(length / 9)?);
        }

        if length.is_multiple_of(3) {
            return T::mixed_radix3(length, Pxdct::dct2_strategy(length / 3)?);
        }

        if length.is_multiple_of(5) {
            return T::mixed_radix5(length, Pxdct::dct2_strategy(length / 5)?);
        }

        if length.is_multiple_of(7) {
            return T::mixed_radix7(length, Pxdct::dct2_strategy(length / 7)?);
        }

        if length.is_multiple_of(11) {
            return T::mixed_radix11(length, Pxdct::dct2_strategy(length / 11)?);
        }

        if length.is_multiple_of(13) {
            return T::mixed_radix13(length, Pxdct::dct2_strategy(length / 13)?);
        }

        let prime_factors = PrimeFactors::from_number(length as u64);

        // if number is too big and 2^k*Q is big this is not super effective to use mixed-radix,
        // go with FFT instead
        let factor2 = prime_factors.factor_of_2();

        if (length.is_multiple_of(2) && factor2 < 5) || (length.is_multiple_of(2) && length < 1000)
        {
            return T::mixed_radix2(length, Pxdct::dct2_strategy(length / 2)?);
        }

        if length == 29 * 11 {
            return T::dct2_relatively_prime(Pxdct::dct2_strategy(29)?, Pxdct::dct2_strategy(11)?);
        }

        Dct2Fft::new(length).map(|x| Arc::new(x) as Arc<dyn PxdctExecutor<T> + Send + Sync>)
    }

    fn strategy_scaled_dct2<T: DctSample + ScaledDct2Factory + Dct2Factory>(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<T> + Send + Sync>, PxdctError>
    where
        f64: AsPrimitive<T>,
    {
        if length == 0 {
            return Err(PxdctError::ZeroSizedDct);
        }
        if length == 1 {
            return Ok(Arc::new(DctIdentity::default()));
        }
        if length.is_power_of_two() {
            return match length {
                2 => Ok(T::scaled_dct2_butterfly2()),
                4 => Ok(T::scaled_dct2_butterfly4()),
                8 => Ok(T::scaled_dct2_butterfly8()),
                16 => Ok(T::scaled_dct2_butterfly16()),
                32 => Ok(T::scaled_dct2_butterfly32()),
                64 => Ok(T::scaled_dct2_butterfly64()),
                128 => Ok(T::scaled_dct2_butterfly128()),
                256 => Ok(T::scaled_dct2_butterfly256()),
                512 => Ok(T::scaled_dct2_butterfly512()),
                _ => T::scaled_split_radix(
                    length,
                    Pxdct::dct2_strategy(length / 2)?,
                    Pxdct::dct2_strategy(length / 4)?,
                ),
            };
        }
        Ok(Arc::new(ScalingInterceptor {
            interceptor: Pxdct::dct2_strategy(length)?,
        }))
    }

    /// Creates a single-precision (f32) DCT-II executor.
    pub fn make_dct2_f32(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f32> + Send + Sync>, PxdctError> {
        Pxdct::dct2_strategy(length)
    }

    /// Creates a scaled single-precision (f32) DCT-II executor.
    ///
    /// Results are scaled by SQRT(2/N).
    /// Scaling routine is completely incorporated only in power of 2 executors.
    /// For anything else if absolute performance is preferred consider to
    /// incorporate scaling factor into your processing routine.
    pub fn make_scaled_dct2_f32(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f32> + Send + Sync>, PxdctError> {
        Pxdct::strategy_scaled_dct2(length)
    }

    /// Creates a double-precision (f64) DCT-II executor.
    pub fn make_dct2_f64(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f64> + Send + Sync>, PxdctError> {
        Pxdct::dct2_strategy(length)
    }

    /// Creates a scaled double-precision (f64) DCT-II executor.
    ///
    /// Results are scaled by SQRT(2/N).
    /// Scaling routine is completely incorporated only in power of 2 executors.
    /// For anything else if absolute performance is preferred consider to
    /// incorporate scaling factor into your processing routine.
    pub fn make_scaled_dct2_f64(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f64> + Send + Sync>, PxdctError> {
        Pxdct::strategy_scaled_dct2(length)
    }

    fn dct3_strategy<T: Copy + Dct3Factory>(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<T> + Send + Sync>, PxdctError> {
        if length == 0 {
            return Err(PxdctError::ZeroSizedDct);
        }

        match length {
            1 => return Ok(T::dct3_identity()),
            2 => return Ok(T::dct3_butterfly2()),
            3 => return Ok(T::dct3_butterfly3()),
            4 => return Ok(T::dct3_butterfly4()),
            5 => return Ok(T::dct3_butterfly5()),
            6 => return Ok(T::dct3_butterfly6()),
            7 => return Ok(T::dct3_butterfly7()),
            8 => return Ok(T::dct3_butterfly8()),
            9 => return Ok(T::dct3_butterfly9()),
            10 => return Ok(T::dct3_butterfly10()),
            11 => return Ok(T::dct3_butterfly11()),
            12 => return Ok(T::dct3_butterfly12()),
            13 => return Ok(T::dct3_butterfly13()),
            14 => return Ok(T::dct3_butterfly14()),
            15 => return Ok(T::dct3_butterfly15()),
            16 => return Ok(T::dct3_butterfly16()),
            18 => return Ok(T::dct3_butterfly18()),
            20 => return Ok(T::dct3_butterfly20()),
            21 => return Ok(T::dct3_butterfly21()),
            24 => return Ok(T::dct3_butterfly24()),
            26 => return Ok(T::dct3_butterfly26()),
            28 => return Ok(T::dct3_butterfly28()),
            30 => return Ok(T::dct3_butterfly30()),
            32 => return Ok(T::dct3_butterfly32()),
            35 => return Ok(T::dct3_butterfly35()),
            36 => return Ok(T::dct3_butterfly36()),
            64 => return Ok(T::dct3_butterfly64()),
            _ => {}
        }

        if length.is_power_of_two() && length > 2 {
            return T::dct3_split_radix(
                length,
                Pxdct::dct3_strategy(length / 2)?,
                Pxdct::dct3_strategy(length / 4)?,
            );
        }

        T::dct3_fft(length)
    }

    /// Creates a single-precision (f32) DCT-III executor.
    pub fn make_dct3_f32(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f32> + Send + Sync>, PxdctError> {
        Pxdct::dct3_strategy(length)
    }

    /// Creates a single-precision (f32) DCT-III executor.
    /// Results are scaled by SQRT(2/N).
    /// Scaling routine is completely incorporated only in power of 2 executors.
    /// For anything else if absolute performance is preferred consider to
    /// incorporate scaling factor into your processing routine.
    pub fn make_scaled_dct3_f32(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f32> + Send + Sync>, PxdctError> {
        Ok(Arc::new(ScalingInterceptor {
            interceptor: Pxdct::dct3_strategy(length)?,
        }))
    }

    /// Creates a double-precision (f64) DCT-III executor.
    pub fn make_dct3_f64(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f64> + Send + Sync>, PxdctError> {
        Pxdct::dct3_strategy(length)
    }

    /// Creates a double-precision (f64) DCT-III executor.
    pub fn make_scaling_dct3_f64(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f64> + Send + Sync>, PxdctError> {
        Ok(Arc::new(ScalingInterceptor {
            interceptor: Pxdct::dct3_strategy(length)?,
        }))
    }

    /// Creates a single-precision (f32) DST-II executor.
    pub fn make_dst2_f32(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f32> + Send + Sync>, PxdctError> {
        if length == 0 {
            return Err(PxdctError::ZeroSizedDct);
        }
        #[cfg(all(target_arch = "x86_64", feature = "avx"))]
        if std::arch::is_x86_feature_detected!("avx2") && std::arch::is_x86_feature_detected!("fma")
        {
            #[allow(clippy::collapsible_if)]
            if length == 4 {
                use crate::avx::AvxDst2Butterfly4;
                static Q: OnceLock<Arc<dyn PxdctExecutor<f32> + Send + Sync>> = OnceLock::new();
                return Ok(Q
                    .get_or_init(|| {
                        Arc::new(AvxDst2Butterfly4::default())
                            as Arc<dyn PxdctExecutor<f32> + Send + Sync>
                    })
                    .clone());
            }
        }
        if length == 2 {
            static Q: OnceLock<Arc<dyn PxdctExecutor<f32> + Send + Sync>> = OnceLock::new();
            return Ok(Q
                .get_or_init(|| {
                    Arc::new(Dst2Butterfly2::default()) as Arc<dyn PxdctExecutor<f32> + Send + Sync>
                })
                .clone());
        } else if length == 4 {
            static Q: OnceLock<Arc<dyn PxdctExecutor<f32> + Send + Sync>> = OnceLock::new();
            return Ok(Q
                .get_or_init(|| {
                    Arc::new(Dst2Butterfly4::default()) as Arc<dyn PxdctExecutor<f32> + Send + Sync>
                })
                .clone());
        }
        Dst2Fft::new(length).map(|x| Arc::new(x) as Arc<dyn PxdctExecutor<f32> + Send + Sync>)
    }

    /// Creates a double-precision (f64) DST-II executor.
    pub fn make_dst2_f64(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f64> + Send + Sync>, PxdctError> {
        if length == 0 {
            return Err(PxdctError::ZeroSizedDct);
        }
        #[cfg(all(target_arch = "x86_64", feature = "avx"))]
        if std::arch::is_x86_feature_detected!("avx2") && std::arch::is_x86_feature_detected!("fma")
        {
            #[allow(clippy::collapsible_if)]
            if length == 4 {
                use crate::avx::AvxDst2Butterfly4;
                static Q: OnceLock<Arc<dyn PxdctExecutor<f64> + Send + Sync>> = OnceLock::new();
                return Ok(Q
                    .get_or_init(|| {
                        Arc::new(AvxDst2Butterfly4::default())
                            as Arc<dyn PxdctExecutor<f64> + Send + Sync>
                    })
                    .clone());
            }
        }
        if length == 2 {
            static Q: OnceLock<Arc<dyn PxdctExecutor<f64> + Send + Sync>> = OnceLock::new();
            return Ok(Q
                .get_or_init(|| {
                    Arc::new(Dst2Butterfly2::default()) as Arc<dyn PxdctExecutor<f64> + Send + Sync>
                })
                .clone());
        } else if length == 4 {
            static Q: OnceLock<Arc<dyn PxdctExecutor<f64> + Send + Sync>> = OnceLock::new();
            return Ok(Q
                .get_or_init(|| {
                    Arc::new(Dst2Butterfly4::default()) as Arc<dyn PxdctExecutor<f64> + Send + Sync>
                })
                .clone());
        }

        if length.is_power_of_two() && length > 2 {
            return f64::dst2_split_radix(
                length,
                Pxdct::make_dct2_f64(length / 2)?,
                Pxdct::make_dct2_f64(length / 4)?,
            );
        }

        Dst2Fft::new(length).map(|x| Arc::new(x) as Arc<dyn PxdctExecutor<f64> + Send + Sync>)
    }

    /// Creates a single-precision (f32) DST-III executor.
    pub fn make_dst3_f32(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f32> + Send + Sync>, PxdctError> {
        if length == 0 {
            return Err(PxdctError::ZeroSizedDct);
        }

        if length == 2 {
            static Q: OnceLock<Arc<dyn PxdctExecutor<f32> + Send + Sync>> = OnceLock::new();
            return Ok(Q
                .get_or_init(|| {
                    Arc::new(Dst3Butterfly2::default()) as Arc<dyn PxdctExecutor<f32> + Send + Sync>
                })
                .clone());
        } else if length == 3 {
            static Q: OnceLock<Arc<dyn PxdctExecutor<f32> + Send + Sync>> = OnceLock::new();
            return Ok(Q
                .get_or_init(|| {
                    Arc::new(Dst3Butterfly3::default()) as Arc<dyn PxdctExecutor<f32> + Send + Sync>
                })
                .clone());
        } else if length == 4 {
            static Q: OnceLock<Arc<dyn PxdctExecutor<f32> + Send + Sync>> = OnceLock::new();
            return Ok(Q
                .get_or_init(|| {
                    Arc::new(Dst3Butterfly4::default()) as Arc<dyn PxdctExecutor<f32> + Send + Sync>
                })
                .clone());
        } else if length == 5 {
            static Q: OnceLock<Arc<dyn PxdctExecutor<f32> + Send + Sync>> = OnceLock::new();
            return Ok(Q
                .get_or_init(|| {
                    Arc::new(Dst3Butterfly5::default()) as Arc<dyn PxdctExecutor<f32> + Send + Sync>
                })
                .clone());
        } else if length == 8 {
            static Q: OnceLock<Arc<dyn PxdctExecutor<f32> + Send + Sync>> = OnceLock::new();
            return Ok(Q
                .get_or_init(|| {
                    Arc::new(Dst3Butterfly8::default()) as Arc<dyn PxdctExecutor<f32> + Send + Sync>
                })
                .clone());
        } else if length == 16 {
            static Q: OnceLock<Arc<dyn PxdctExecutor<f32> + Send + Sync>> = OnceLock::new();
            return Ok(Q
                .get_or_init(|| {
                    Arc::new(Dst3Butterfly16::default())
                        as Arc<dyn PxdctExecutor<f32> + Send + Sync>
                })
                .clone());
        }

        if length.is_power_of_two() && length > 2 {
            return Ok(Arc::new(SplitRadixDst3::new(
                length,
                Pxdct::make_dct3_f32(length / 2)?,
                Pxdct::make_dct3_f32(length / 4)?,
            )?));
        }

        Dst3Fft::new(length).map(|x| Arc::new(x) as Arc<dyn PxdctExecutor<f32> + Send + Sync>)
    }

    /// Creates a double-precision (f64) DST-III executor.
    pub fn make_dst3_f64(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f64> + Send + Sync>, PxdctError> {
        if length == 0 {
            return Err(PxdctError::ZeroSizedDct);
        }

        if length.is_power_of_two() && length > 2 {
            return Ok(Arc::new(SplitRadixDst3::new(
                length,
                Pxdct::make_dct3_f64(length / 2)?,
                Pxdct::make_dct3_f64(length / 4)?,
            )?));
        }

        Dst3Fft::new(length).map(|x| Arc::new(x) as Arc<dyn PxdctExecutor<f64> + Send + Sync>)
    }

    fn strategy_dct4<T: DctSample + Dct4Factory + Dct2Factory>(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<T> + Send + Sync>, PxdctError>
    where
        f64: AsPrimitive<T>,
    {
        if length == 0 {
            return Err(PxdctError::ZeroSizedDct);
        }

        if length == 1 {
            return Ok(T::dct4_identity());
        } else if length == 2 {
            return Ok(T::dct4_butterfly2());
        } else if length == 3 {
            return Ok(T::dct4_butterfly3());
        } else if length == 4 {
            return Ok(T::dct4_butterfly4());
        } else if length == 5 {
            return Ok(T::dct4_butterfly5());
        } else if length == 6 {
            return Ok(T::dct4_butterfly6());
        } else if length == 7 {
            return Ok(T::dct4_butterfly7());
        } else if length == 8 {
            return Ok(T::dct4_butterfly8());
        } else if length == 9 {
            return Ok(T::dct4_butterfly9());
        } else if length == 10 {
            return Ok(T::dct4_butterfly10());
        } else if length == 11 {
            return Ok(T::dct4_butterfly11());
        } else if length == 12 {
            return Ok(T::dct4_butterfly12());
        } else if length == 13 {
            return Ok(T::dct4_butterfly13());
        } else if length == 14 {
            return Ok(T::dct4_butterfly14());
        } else if length == 16 {
            return Ok(T::dct4_butterfly16());
        } else if length == 17 {
            return Ok(T::dct4_butterfly17());
        } else if length == 18 {
            return Ok(T::dct4_butterfly18());
        } else if length == 19 {
            return Ok(T::dct4_butterfly19());
        } else if length == 20 {
            return Ok(T::dct4_butterfly20());
        } else if length == 22 {
            return Ok(T::dct4_butterfly22());
        } else if length == 23 {
            return Ok(T::dct4_butterfly23());
        } else if length == 24 {
            return Ok(T::dct4_butterfly24());
        } else if length == 26 {
            return Ok(T::dct4_butterfly26());
        } else if length == 27 {
            return Ok(T::dct4_butterfly27());
        } else if length == 28 {
            return Ok(T::dct4_butterfly28());
        } else if length == 29 {
            return Ok(T::dct4_butterfly29());
        } else if length == 30 {
            return Ok(T::dct4_butterfly30());
        } else if length == 32 {
            return Ok(T::dct4_butterfly32());
        }

        if length.is_multiple_of(9) {
            let half_length = length / 9;
            return T::dct4_mixed_radix9(length, Pxdct::strategy_dct4(half_length)?);
        }

        if length.is_multiple_of(7) {
            let half_length = length / 7;
            return T::dct4_mixed_radix7(length, Pxdct::strategy_dct4(half_length)?);
        }

        if length.is_multiple_of(5) {
            let half_length = length / 5;
            return T::dct4_mixed_radix5(length, Pxdct::strategy_dct4(half_length)?);
        }

        if length.is_multiple_of(11) {
            let half_length = length / 11;
            return T::dct4_mixed_radix11(length, Pxdct::strategy_dct4(half_length)?);
        }

        if length.is_multiple_of(13) {
            let half_length = length / 13;
            return T::dct4_mixed_radix13(length, Pxdct::strategy_dct4(half_length)?);
        }

        if length.is_multiple_of(17) {
            let half_length = length / 17;
            return T::dct4_mixed_radix17(length, Pxdct::strategy_dct4(half_length)?);
        }

        if length.is_multiple_of(19) {
            let half_length = length / 19;
            return T::dct4_mixed_radix19(length, Pxdct::strategy_dct4(half_length)?);
        }

        if length.is_multiple_of(3) {
            let half_length = length / 3;
            return T::dct4_mixed_radix3(length, Pxdct::strategy_dct4(half_length)?);
        }

        if length.is_power_of_two() {
            return T::dct4_radix2(length, Pxdct::dct2_strategy(length / 2)?);
        }

        if length.is_multiple_of(2) {
            let half_length = length / 2;
            return T::dct4_mixed_radix2(length, Pxdct::dct2_strategy(half_length)?);
        }

        T::dct4_fft(
            T::make_fft(length, FftDirection::Forward)
                .map_err(|x| PxdctError::FftError(x.to_string()))?,
        )
    }

    /// Creates a single-precision (f32) DCT-IV executor.
    pub fn make_dct4_f32(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f32> + Send + Sync>, PxdctError> {
        Pxdct::strategy_dct4(length)
    }

    /// Creates a double-precision (f32) DCT-IV executor.
    pub fn make_dct4_f64(
        length: usize,
    ) -> Result<Arc<dyn PxdctExecutor<f64> + Send + Sync>, PxdctError> {
        Pxdct::strategy_dct4(length)
    }

    /// Creates 2D DCT executor.
    ///
    /// For matrix WxH to get an inverse use H as width and W as height.
    pub fn make_2d_dct_f32(
        width_dct: Arc<dyn PxdctExecutor<f32> + Send + Sync>,
        height_dct: Arc<dyn PxdctExecutor<f32> + Send + Sync>,
    ) -> Arc<dyn MultidimensionalDctExecutor<f32> + Send + Sync> {
        let width = width_dct.length();
        let height = height_dct.length();
        let width_scratch_size = width_dct.scratch_size();
        let height_scratch_size = height_dct.scratch_size();
        Arc::new(TwoDimensionalDct {
            width,
            height,
            height_executor: height_dct,
            width_executor: width_dct,
            transpose_width_to_height: f32::make_transpose(width, height),
            width_scratch_size,
            height_scratch_size,
        })
    }

    /// Creates 2D DCT executor.
    ///
    /// For matrix WxH to get an inverse use H as width and W as height.
    pub fn make_2d_dct_f64(
        width_dct: Arc<dyn PxdctExecutor<f64> + Send + Sync>,
        height_dct: Arc<dyn PxdctExecutor<f64> + Send + Sync>,
    ) -> Arc<dyn MultidimensionalDctExecutor<f64> + Send + Sync> {
        let width = width_dct.length();
        let height = height_dct.length();
        let width_scratch_size = width_dct.scratch_size();
        let height_scratch_size = height_dct.scratch_size();
        Arc::new(TwoDimensionalDct {
            width,
            height,
            height_executor: height_dct,
            width_executor: width_dct,
            transpose_width_to_height: f64::make_transpose(width, height),
            width_scratch_size,
            height_scratch_size,
        })
    }
}

/// Trait implemented by all PXDCT executors.
///
/// This trait defines the common interface for performing an in-place
/// DCT (Discrete Cosine Transform) on a data slice.
pub trait PxdctExecutor<T> {
    /// Executes the DCT transform in-place on the given data buffer.
    fn execute(&self, data: &mut [T]) -> Result<(), PxdctError>;
    fn execute_with_scratch(&self, data: &mut [T], scratch: &mut [T]) -> Result<(), PxdctError>;
    fn execute_into(&self, input: &[T], output: &mut [T]) -> Result<(), PxdctError>;
    fn execute_into_with_scratch(
        &self,
        input: &[T],
        output: &mut [T],
        scratch: &mut [T],
    ) -> Result<(), PxdctError>;
    /// Returns the length of the transform supported by this executor.
    fn length(&self) -> usize;
    fn scratch_size(&self) -> usize;
}

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

    pub(crate) fn naive_dct2(input: &[f64]) -> Vec<f64> {
        let mut result = Vec::new();

        for output_index in 0..input.len() {
            let mut entry = 0.0;
            for input_index in 0..input.len() {
                let cos_inner =
                    (output_index as f64) * (input_index as f64 + 0.5) * std::f64::consts::PI
                        / (input.len() as f64);
                let twiddle = cos_inner.cos();
                entry += input[input_index] * twiddle;
            }
            result.push(entry);
        }

        result
    }

    pub(crate) fn naive_scaled_dct2(input: &[f64]) -> Vec<f64> {
        let mut result = Vec::new();

        let scale = (2f64 / input.len() as f64).sqrt();

        for output_index in 0..input.len() {
            let mut entry = 0.0;
            for input_index in 0..input.len() {
                let cos_inner =
                    (output_index as f64) * (input_index as f64 + 0.5) * std::f64::consts::PI
                        / (input.len() as f64);
                let twiddle = cos_inner.cos();
                entry += input[input_index] * twiddle;
            }
            result.push(entry * scale);
        }

        result
    }

    #[allow(unused)]
    pub(crate) fn naive_scaled_dct2_f32(input: &[f32]) -> Vec<f32> {
        let mut result = Vec::new();

        let scale = (2f32 / input.len() as f32).sqrt();

        for output_index in 0..input.len() {
            let mut entry = 0.0;
            for input_index in 0..input.len() {
                let cos_inner =
                    (output_index as f32) * (input_index as f32 + 0.5) * std::f32::consts::PI
                        / (input.len() as f32);
                let twiddle = cos_inner.cos();
                entry += input[input_index] * twiddle;
            }
            result.push(entry * scale);
        }

        result
    }
    #[allow(unused)]
    pub(crate) fn naive_dct2_f32(input: &[f32]) -> Vec<f32> {
        let mut result = Vec::new();

        for output_index in 0..input.len() {
            let mut entry = 0.0;
            for input_index in 0..input.len() {
                let cos_inner =
                    (output_index as f32) * (input_index as f32 + 0.5) * std::f32::consts::PI
                        / (input.len() as f32);
                let twiddle = cos_inner.cos();
                entry += input[input_index] * twiddle;
            }
            result.push(entry);
        }

        result
    }

    pub(crate) fn naive_dst2(input: &[f64]) -> Vec<f64> {
        let mut result = Vec::new();
        for output_index in 0..input.len() {
            let mut entry = 0.0;
            for input_index in 0..input.len() {
                let sin_inner =
                    (output_index as f64 + 1.0) * (input_index as f64 + 0.5) * std::f64::consts::PI
                        / (input.len() as f64);
                let twiddle = sin_inner.sin();
                entry += input[input_index] * twiddle;
            }
            result.push(entry);
        }
        result
    }

    pub fn naive_dct3(input: &[f64]) -> Vec<f64> {
        let mut result = Vec::new();

        for output_index in 0..input.len() {
            let mut entry = 0.0;
            for input_index in 0..input.len() {
                let multiplier = if input_index == 0 { 0.5 } else { 1.0 };
                let cos_inner =
                    (output_index as f64 + 0.5) * (input_index as f64) * std::f64::consts::PI
                        / (input.len() as f64);
                let twiddle = cos_inner.cos();
                entry += input[input_index] * twiddle * multiplier;
            }
            result.push(entry);
        }

        result
    }

    pub fn naive_dct3_f32(input: &[f32]) -> Vec<f32> {
        let mut result = Vec::new();

        for output_index in 0..input.len() {
            let mut entry = 0.0;
            for input_index in 0..input.len() {
                let multiplier = if input_index == 0 { 0.5 } else { 1.0 };
                let cos_inner =
                    (output_index as f32 + 0.5) * (input_index as f32) * std::f32::consts::PI
                        / (input.len() as f32);
                let twiddle = cos_inner.cos();
                entry += input[input_index] * twiddle * multiplier;
            }
            result.push(entry);
        }

        result
    }

    pub fn naive_dst3(input: &[f64]) -> Vec<f64> {
        let mut result = Vec::new();
        for output_index in 0..input.len() {
            let mut entry = 0.0;
            for input_index in 0..input.len() {
                let multiplier = if input_index == input.len() - 1 {
                    0.5
                } else {
                    1.0
                };
                let sin_inner =
                    (output_index as f64 + 0.5) * (input_index as f64 + 1.0) * std::f64::consts::PI
                        / (input.len() as f64);
                let twiddle = sin_inner.sin();
                entry += input[input_index] * twiddle * multiplier;
            }
            result.push(entry);
        }
        result
    }

    pub fn naive_dct4(input: &[f64]) -> Vec<f64> {
        let mut result = Vec::new();

        for output_index in 0..input.len() {
            let mut entry = 0.0;
            for input_index in 0..input.len() {
                let cos_inner =
                    (output_index as f64 + 0.5) * (input_index as f64 + 0.5) * std::f64::consts::PI
                        / (input.len() as f64);
                let twiddle = cos_inner.cos();
                entry += input[input_index] * twiddle;
            }
            result.push(entry);
        }

        result
    }

    pub fn naive_dct4_f32(input: &[f32]) -> Vec<f32> {
        let mut result = Vec::new();

        for output_index in 0..input.len() {
            let mut entry = 0.0;
            for input_index in 0..input.len() {
                let cos_inner =
                    (output_index as f32 + 0.5) * (input_index as f32 + 0.5) * std::f32::consts::PI
                        / (input.len() as f32);
                let twiddle = cos_inner.cos();
                entry += input[input_index] * twiddle;
            }
            result.push(entry);
        }

        result
    }

    #[test]
    fn dct2_roundtrip() {
        for i in 1..250 {
            let mut array = vec![0f32; i];
            for i in 1..i + 1 {
                array[i - 1] = i as f32;
            }
            let mut working_array = array.clone();
            let dct_forward = Pxdct::make_dct2_f32(array.len()).unwrap();
            let dct_inverse = Pxdct::make_dct3_f32(array.len()).unwrap();

            dct_forward.execute(&mut working_array).unwrap();
            dct_inverse.execute(&mut working_array).unwrap();

            for k in working_array.iter_mut() {
                *k = *k / (i as f32) * 2.;
            }

            working_array.iter().zip(array.iter()).enumerate().for_each(|(k, (&x, &c))| {
                assert!((x - c).abs() < 0.01, "Difference to control values exceeded 0.01 when it shouldn't, value {x}, control {c} at {k} for size {i}");
            });
        }
    }

    #[test]
    fn dct2_roundtrip_oof() {
        for i in 1..250 {
            let mut array = vec![0f32; i];
            for i in 1..i + 1 {
                array[i - 1] = i as f32;
            }
            let mut working_array = array.clone();

            let mut transient = vec![0f32; i];

            let dct_forward = Pxdct::make_dct2_f32(array.len()).unwrap();
            let dct_inverse = Pxdct::make_dct3_f32(array.len()).unwrap();

            dct_forward
                .execute_into(&working_array, &mut transient)
                .expect(&format!("Failed to execute DCT on size {i}"));
            dct_inverse
                .execute_into(&transient, &mut working_array)
                .expect(&format!("Failed to execute DCT on size {i}"));

            for k in working_array.iter_mut() {
                *k = *k / (i as f32) * 2.;
            }

            working_array.iter().zip(array.iter()).enumerate().for_each(|(k, (&x, &c))| {
                assert!((x - c).abs() < 0.01, "Difference to control values exceeded 0.01 when it shouldn't, value {x}, control {c} at {k} for size {i}");
            });
        }
    }

    #[test]
    fn dct2_roundtrip_f64() {
        for i in 1..250 {
            let mut array = vec![0f64; i];
            for i in 1..i + 1 {
                array[i - 1] = i as f64;
            }
            let mut working_array = array.clone();
            let dct_forward = Pxdct::make_dct2_f64(array.len()).unwrap();
            let dct_inverse = Pxdct::make_dct3_f64(array.len()).unwrap();

            dct_forward.execute(&mut working_array).unwrap();
            dct_inverse.execute(&mut working_array).unwrap();

            for k in working_array.iter_mut() {
                *k = *k / (i as f64) * 2.;
            }

            working_array.iter().zip(array.iter()).enumerate().for_each(|(k, (&x, &c))| {
                assert!((x - c).abs() < 0.01, "Difference to control values exceeded 0.01 when it shouldn't, value {x}, control {c} at {k} for size {i}");
            });
        }
    }

    #[test]
    fn dct2_roundtrip_oof_f64() {
        for i in 1..250 {
            let mut array = vec![0f64; i];
            for i in 1..i + 1 {
                array[i - 1] = i as f64;
            }
            let mut working_array = array.clone();
            let mut transient = vec![0f64; i];

            let dct_forward = Pxdct::make_dct2_f64(array.len()).unwrap();
            let dct_inverse = Pxdct::make_dct3_f64(array.len()).unwrap();

            dct_forward
                .execute_into(&working_array, &mut transient)
                .expect(&format!("Failed to execute DCT on size {i}"));
            dct_inverse
                .execute_into(&transient, &mut working_array)
                .expect(&format!("Failed to execute DCT on size {i}"));

            for k in working_array.iter_mut() {
                *k = *k / (i as f64) * 2.;
            }

            working_array.iter().zip(array.iter()).enumerate().for_each(|(k, (&x, &c))| {
                assert!((x - c).abs() < 0.01, "Difference to control values exceeded 0.01 when it shouldn't, value {x}, control {c} at {k} for size {i}");
            });
        }
    }

    #[test]
    fn dst2_roundtrip() {
        for i in 1..250 {
            let mut array = vec![0f32; i];
            for i in 1..i + 1 {
                array[i - 1] = i as f32;
            }
            let mut working_array = array.clone();
            let dct_forward = Pxdct::make_dst2_f32(array.len()).unwrap();
            let dct_inverse = Pxdct::make_dst3_f32(array.len()).unwrap();

            dct_forward.execute(&mut working_array).unwrap();
            dct_inverse.execute(&mut working_array).unwrap();

            for k in working_array.iter_mut() {
                *k = *k / (i as f32) * 2.;
            }

            working_array.iter().zip(array.iter()).enumerate().for_each(|(k, (&x, &c))| {
                assert!((x - c).abs() < 0.01, "Difference to control values exceeded 0.01 when it shouldn't, value {x}, control {c} at {k} for size {i}");
            });
        }
    }

    #[test]
    fn dst2_roundtrip_oof() {
        for i in 1..250 {
            let mut array = vec![0f32; i];
            for i in 1..i + 1 {
                array[i - 1] = i as f32;
            }
            let mut working_array = array.clone();
            let mut transient = vec![0f32; i];

            let dct_forward = Pxdct::make_dct2_f32(array.len()).unwrap();
            let dct_inverse = Pxdct::make_dct3_f32(array.len()).unwrap();

            dct_forward
                .execute_into(&working_array, &mut transient)
                .expect(&format!("Failed to execute DST-II on size {i}"));
            dct_inverse
                .execute_into(&transient, &mut working_array)
                .expect(&format!("Failed to execute DST-III on size {i}"));

            for k in working_array.iter_mut() {
                *k = *k / (i as f32) * 2.;
            }

            working_array.iter().zip(array.iter()).enumerate().for_each(|(k, (&x, &c))| {
                assert!((x - c).abs() < 0.01, "Difference to control values exceeded 0.01 when it shouldn't, value {x}, control {c} at {k} for size {i}");
            });
        }
    }

    #[test]
    fn dct4_roundtrip() {
        for i in 1..250 {
            let mut array = vec![0f32; i];
            for i in 1..i + 1 {
                array[i - 1] = i as f32;
            }
            let mut working_array = array.clone();
            let dct_forward = Pxdct::make_dct4_f32(array.len()).unwrap();
            let dct_inverse = Pxdct::make_dct4_f32(array.len()).unwrap();

            dct_forward.execute(&mut working_array).unwrap();
            dct_inverse.execute(&mut working_array).unwrap();

            for k in working_array.iter_mut() {
                *k = *k / (i as f32) * 2.;
            }

            working_array.iter().zip(array.iter()).enumerate().for_each(|(k, (&x, &c))| {
                assert!((x - c).abs() < 0.01, "Difference to control values exceeded 0.01 when it shouldn't, value {x}, control {c} at {k} for size {i}");
            });
        }
    }

    #[test]
    fn dct4_roundtrip_oof() {
        for i in 1..250 {
            let mut array = vec![0f32; i];
            for i in 1..i + 1 {
                array[i - 1] = i as f32;
            }
            let mut working_array = array.clone();

            let mut transient = vec![0f32; i];

            let dct_forward = Pxdct::make_dct4_f32(array.len()).unwrap();
            let dct_inverse = Pxdct::make_dct4_f32(array.len()).unwrap();

            dct_forward
                .execute_into(&working_array, &mut transient)
                .expect(&format!("Failed to execute DCT-IV on size {i}"));
            dct_inverse
                .execute_into(&transient, &mut working_array)
                .expect(&format!("Failed to execute DCT-IV on size {i}"));

            for k in working_array.iter_mut() {
                *k = *k / (i as f32) * 2.;
            }

            working_array.iter().zip(array.iter()).enumerate().for_each(|(k, (&x, &c))| {
                assert!((x - c).abs() < 0.01, "Difference to control values exceeded 0.01 when it shouldn't, value {x}, control {c} at {k} for size {i}");
            });
        }
    }

    #[test]
    fn dct4_roundtrip_f64() {
        for i in 1..300 {
            let mut array = vec![0.; i];
            for i in 1..i + 1 {
                array[i - 1] = i as f64;
            }
            let mut working_array = array.clone();
            let dct_forward = Pxdct::make_dct4_f64(array.len()).unwrap();
            let dct_inverse = Pxdct::make_dct4_f64(array.len()).unwrap();

            dct_forward.execute(&mut working_array).unwrap();
            dct_inverse.execute(&mut working_array).unwrap();

            for k in working_array.iter_mut() {
                *k = *k / (i as f64) * 2.;
            }

            working_array.iter().zip(array.iter()).enumerate().for_each(|(k, (&x, &c))| {
                assert!((x - c).abs() < 0.00001, "Difference to control values exceeded 0.01 when it shouldn't, value {x}, control {c} at {k} for size {i}");
            });
        }
    }

    #[test]
    fn dct4_roundtrip_oof_f64() {
        for i in 1..300 {
            let mut array = vec![0.; i];
            for i in 1..i + 1 {
                array[i - 1] = i as f64;
            }
            let mut working_array = array.clone();

            let mut transient = vec![0f64; i];

            let dct_forward = Pxdct::make_dct4_f64(array.len()).unwrap();
            let dct_inverse = Pxdct::make_dct4_f64(array.len()).unwrap();

            dct_forward
                .execute_into(&working_array, &mut transient)
                .expect(&format!("Failed to execute DCT-IV on size {i}"));
            dct_inverse
                .execute_into(&transient, &mut working_array)
                .expect(&format!("Failed to execute DCT-IV on size {i}"));

            for k in working_array.iter_mut() {
                *k = *k / (i as f64) * 2.;
            }

            working_array.iter().zip(array.iter()).enumerate().for_each(|(k, (&x, &c))| {
                assert!((x - c).abs() < 0.00001, "Difference to control values exceeded 0.01 when it shouldn't, value {x}, control {c} at {k} for size {i}");
            });
        }
    }
}