use nalgebra::{*, allocator::Allocator};
use std::f64::consts::{PI, E};
/*use std::ops::Sub;*/
/*use typenum::{B1, UInt, UTerm};*/

/* dynamic matrices */
pub fn rand_eigval_series<N>(time_res: usize) -> Vec<OVector<Complex<f64>, Dyn>>
    where
        N: ToTypenum + DimName + DimSub<U1>,
        DefaultAllocator:
            Allocator<N> +
            Allocator<N, N> +
            Allocator<<N as DimSub<U1>>::Output> +
            Allocator<N, <N as DimSub<U1>>::Output>
{
    // initialize the random matrix
    let dim = N::try_to_usize().unwrap();
    let mut rand_mat = DMatrix::<f64>::from_fn(dim, dim, |j, k| {
        let n = j*dim + k;
        E*((n*n) as f64) % 2.0 - 1.0
    }) * (3.0 / (dim as f64)).sqrt();
    
    // initialize the rotation step
    let mut rot_step = DMatrix::<f64>::identity(dim, dim);
    let max_freq = 4;
    for n in (0..dim).step_by(2) {
        let ang = PI * ((n % max_freq) as f64) / (time_res as f64);
        let ang_cos = ang.cos();
        let ang_sin = ang.sin();
        rot_step[(n, n)] = ang_cos;
        rot_step[(n+1, n)] = ang_sin;
        rot_step[(n, n+1)] = -ang_sin;
        rot_step[(n+1, n+1)] = ang_cos;
    }
    
    // find the eigenvalues
    let mut eigval_series = Vec::<OVector<Complex<f64>, Dyn>>::with_capacity(time_res);
    eigval_series.push(rand_mat.complex_eigenvalues());
    for _ in 1..time_res {
        rand_mat = &rot_step * rand_mat;
        eigval_series.push(rand_mat.complex_eigenvalues());
    }
    eigval_series
}

/* dynamic single float matrices */
/*pub fn rand_eigval_series<N>(time_res: usize) -> Vec<OVector<Complex<f32>, Dyn>>
    where
        N: ToTypenum + DimName + DimSub<U1>,
        DefaultAllocator:
            Allocator<N> +
            Allocator<N, N> +
            Allocator<<N as DimSub<U1>>::Output> +
            Allocator<N, <N as DimSub<U1>>::Output>
{
    // initialize the random matrix
    let dim = N::try_to_usize().unwrap();
    let mut rand_mat = DMatrix::<f32>::from_fn(dim, dim, |j, k| {
        let n = j*dim + k;
        (E as f32)*((n*n) as f32) % 2.0_f32 - 1.0_f32
    }) * (3.0_f32 / (dim as f32)).sqrt();
    
    // initialize the rotation step
    let mut rot_step = DMatrix::<f32>::identity(dim, dim);
    let max_freq = 4;
    for n in (0..dim).step_by(2) {
        let ang = (PI as f32) * ((n % max_freq) as f32) / (time_res as f32);
        let ang_cos = ang.cos();
        let ang_sin = ang.sin();
        rot_step[(n, n)] = ang_cos;
        rot_step[(n+1, n)] = ang_sin;
        rot_step[(n, n+1)] = -ang_sin;
        rot_step[(n+1, n+1)] = ang_cos;
    }
    
    // find the eigenvalues
    let mut eigval_series = Vec::<OVector<Complex<f32>, Dyn>>::with_capacity(time_res);
    eigval_series.push(rand_mat.complex_eigenvalues());
    for _ in 1..time_res {
        rand_mat = &rot_step * rand_mat;
        eigval_series.push(rand_mat.complex_eigenvalues());
    }
    eigval_series
}*/

/* static matrices. should only be used when the dimension is really small */
/*pub fn rand_eigval_series<N>(time_res: usize) -> Vec<OVector<Complex<f64>, N>>
    where
        N: ToTypenum + DimName + DimSub<U1>,
        DefaultAllocator:
            Allocator<N> +
            Allocator<N, N> +
            Allocator<<N as DimSub<U1>>::Output> +
            Allocator<N, <N as DimSub<U1>>::Output>
{
    // initialize the random matrix
    let dim = N::try_to_usize().unwrap();
    let mut rand_mat = OMatrix::<f64, N, N>::from_fn(|j, k| {
        let n = j*dim + k;
        E*((n*n) as f64) % 2.0 - 1.0
    }) * (3.0 / (dim as f64)).sqrt();
    /*let mut rand_mat = OMatrix::<f64, N, N>::identity();*/
    
    // initialize the rotation step
    let mut rot_step = OMatrix::<f64, N, N>::identity();
    let max_freq = 4;
    for n in (0..dim).step_by(2) {
        let ang = PI * ((n % max_freq) as f64) / (time_res as f64);
        let ang_cos = ang.cos();
        let ang_sin = ang.sin();
        rot_step[(n, n)] = ang_cos;
        rot_step[(n+1, n)] = ang_sin;
        rot_step[(n, n+1)] = -ang_sin;
        rot_step[(n+1, n+1)] = ang_cos;
    }
    
    // find the eigenvalues
    let mut eigval_series = Vec::<OVector<Complex<f64>, N>>::with_capacity(time_res);
    eigval_series.push(rand_mat.complex_eigenvalues());
    for _ in 1..time_res {
        rand_mat = &rot_step * rand_mat;
        eigval_series.push(rand_mat.complex_eigenvalues());
    }
    eigval_series
}*/

/* another attempt at static matrices. i couldn't get the types to work out */
/*pub fn random_eigval_series<const N: usize>(time_res: usize) -> Vec<OVector<Complex<f64>, Const<N>>>
    where
        Const<N>: ToTypenum,
        <Const<N> as ToTypenum>::Typenum: Sub<UInt<UTerm, B1>>,
        <<Const<N> as ToTypenum>::Typenum as Sub<UInt<UTerm, B1>>>::Output: ToConst
{
    // initialize the random matrix
    /*let mut rand_mat = SMatrix::<f64, N, N>::zeros();
    for n in 0..N*N {
        rand_mat[n] = E*((n*n) as f64) % 2.0 - 1.0;
    }*/
    let rand_mat = OMatrix::<f64, Const<N>, Const<N>>::from_fn(|j, k| {
        let n = j*N + k;
        E*((n*n) as f64) % 2.0 - 1.0
    });
    
    // initialize the rotation step
    let mut rot_step = OMatrix::<f64, Const<N>, Const<N>>::identity();
    let max_freq = 4;
    for n in (0..N).step_by(2) {
        let ang = PI * ((n % max_freq) as f64) / (time_res as f64);
        let ang_cos = ang.cos();
        let ang_sin = ang.sin();
        rot_step[(n, n)] = ang_cos;
        rot_step[(n+1, n)] = ang_sin;
        rot_step[(n, n+1)] = -ang_sin;
        rot_step[(n+1, n+1)] = ang_cos;
    }
    
    // find the eigenvalues
    let mut eigvals = Vec::<OVector<Complex<f64>, Const<N>>>::with_capacity(time_res);
    unsafe { eigvals.set_len(time_res); }
    for t in 0..time_res {
        eigvals[t] = rand_mat.complex_eigenvalues();
    }
    eigvals
}*/