1use crate::reference::fft64::reim::as_arr;
  2
  3#[inline(always)]
  4pub fn reim4_extract_1blk_from_reim_ref(m: usize, rows: usize, blk: usize, dst: &mut [f64], src: &[f64]) {
  5    let mut offset: usize = blk << 2;
  6
  7    debug_assert!(blk < (m >> 2));
  8    debug_assert!(dst.len() >= 2 * rows * 4);
  9
 10    for chunk in dst.chunks_exact_mut(4).take(2 * rows) {
 11        chunk.copy_from_slice(&src[offset..offset + 4]);
 12        offset += m
 13    }
 14}
 15
 16#[inline(always)]
 17pub fn reim4_save_1blk_to_reim_ref<const OVERWRITE: bool>(m: usize, blk: usize, dst: &mut [f64], src: &[f64]) {
 18    let mut offset: usize = blk << 2;
 19
 20    debug_assert!(blk < (m >> 2));
 21    debug_assert!(dst.len() >= offset + m + 4);
 22    debug_assert!(src.len() >= 8);
 23
 24    let dst_off = &mut dst[offset..offset + 4];
 25
 26    if OVERWRITE {
 27        dst_off.copy_from_slice(&src[0..4]);
 28    } else {
 29        dst_off[0] += src[0];
 30        dst_off[1] += src[1];
 31        dst_off[2] += src[2];
 32        dst_off[3] += src[3];
 33    }
 34
 35    offset += m;
 36
 37    let dst_off = &mut dst[offset..offset + 4];
 38    if OVERWRITE {
 39        dst_off.copy_from_slice(&src[4..8]);
 40    } else {
 41        dst_off[0] += src[4];
 42        dst_off[1] += src[5];
 43        dst_off[2] += src[6];
 44        dst_off[3] += src[7];
 45    }
 46}
 47
 48#[inline(always)]
 49pub fn reim4_save_2blk_to_reim_ref<const OVERWRITE: bool>(m: usize, blk: usize, dst: &mut [f64], src: &[f64]) {
 50    let mut offset: usize = blk << 2;
 51
 52    debug_assert!(blk < (m >> 2));
 53    debug_assert!(dst.len() >= offset + 3 * m + 4);
 54    debug_assert!(src.len() >= 16);
 55
 56    let dst_off = &mut dst[offset..offset + 4];
 57    if OVERWRITE {
 58        dst_off.copy_from_slice(&src[0..4]);
 59    } else {
 60        dst_off[0] += src[0];
 61        dst_off[1] += src[1];
 62        dst_off[2] += src[2];
 63        dst_off[3] += src[3];
 64    }
 65
 66    offset += m;
 67    let dst_off = &mut dst[offset..offset + 4];
 68    if OVERWRITE {
 69        dst_off.copy_from_slice(&src[4..8]);
 70    } else {
 71        dst_off[0] += src[4];
 72        dst_off[1] += src[5];
 73        dst_off[2] += src[6];
 74        dst_off[3] += src[7];
 75    }
 76
 77    offset += m;
 78
 79    let dst_off = &mut dst[offset..offset + 4];
 80    if OVERWRITE {
 81        dst_off.copy_from_slice(&src[8..12]);
 82    } else {
 83        dst_off[0] += src[8];
 84        dst_off[1] += src[9];
 85        dst_off[2] += src[10];
 86        dst_off[3] += src[11];
 87    }
 88
 89    offset += m;
 90    let dst_off = &mut dst[offset..offset + 4];
 91    if OVERWRITE {
 92        dst_off.copy_from_slice(&src[12..16]);
 93    } else {
 94        dst_off[0] += src[12];
 95        dst_off[1] += src[13];
 96        dst_off[2] += src[14];
 97        dst_off[3] += src[15];
 98    }
 99}
100
101#[inline(always)]
102pub fn reim4_vec_mat1col_product_ref(
103    nrows: usize,
104    dst: &mut [f64], // 8 doubles: [re1(4), im1(4)]
105    u: &[f64],       // nrows * 8 doubles: [ur(4) | ui(4)] per row
106    v: &[f64],       // nrows * 8 doubles: [ar(4) | ai(4)] per row
107) {
108    #[cfg(debug_assertions)]
109    {
110        assert!(dst.len() >= 8, "dst must have at least 8 doubles");
111        assert!(u.len() >= nrows * 8, "u must be at least nrows * 8 doubles");
112        assert!(v.len() >= nrows * 8, "v must be at least nrows * 8 doubles");
113    }
114
115    let mut acc: [f64; 8] = [0f64; 8];
116    let mut j = 0;
117    for _ in 0..nrows {
118        reim4_add_mul(&mut acc, as_arr(&u[j..]), as_arr(&v[j..]));
119        j += 8;
120    }
121    dst[0..8].copy_from_slice(&acc);
122}
123
124#[inline(always)]
125pub fn reim4_vec_mat2cols_product_ref(
126    nrows: usize,
127    dst: &mut [f64], // 16 doubles: [re1(4), im1(4), re2(4), im2(4)]
128    u: &[f64],       // nrows * 8 doubles: [ur(4) | ui(4)] per row
129    v: &[f64],       // nrows * 16 doubles: [ar(4) | ai(4) | br(4) | bi(4)] per row
130) {
131    #[cfg(debug_assertions)]
132    {
133        assert_eq!(dst.len(), 16, "dst must have 16 doubles");
134        assert!(u.len() >= nrows * 8, "u must be at least nrows * 8 doubles");
135        assert!(
136            v.len() >= nrows * 16,
137            "v must be at least nrows * 16 doubles"
138        );
139    }
140
141    // zero accumulators
142    let mut acc_0: [f64; 8] = [0f64; 8];
143    let mut acc_1: [f64; 8] = [0f64; 8];
144    for i in 0..nrows {
145        let _1j: usize = i << 3;
146        let _2j: usize = i << 4;
147        let u_j: &[f64; 8] = as_arr(&u[_1j..]);
148        reim4_add_mul(&mut acc_0, u_j, as_arr(&v[_2j..]));
149        reim4_add_mul(&mut acc_1, u_j, as_arr(&v[_2j + 8..]));
150    }
151    dst[0..8].copy_from_slice(&acc_0);
152    dst[8..16].copy_from_slice(&acc_1);
153}
154
155#[inline(always)]
156pub fn reim4_vec_mat2cols_2ndcol_product_ref(
157    nrows: usize,
158    dst: &mut [f64], // 8 doubles: [re1(4), im1(4), re2(4), im2(4)]
159    u: &[f64],       // nrows * 8 doubles: [ur(4) | ui(4)] per row
160    v: &[f64],       // nrows * 16 doubles: [x | x | br(4) | bi(4)] per row
161) {
162    #[cfg(debug_assertions)]
163    {
164        assert!(
165            dst.len() >= 8,
166            "dst must be at least 8 doubles but is {}",
167            dst.len()
168        );
169        assert!(
170            u.len() >= nrows * 8,
171            "u must be at least nrows={} * 8 doubles but is {}",
172            nrows,
173            u.len()
174        );
175        assert!(
176            v.len() >= nrows * 16,
177            "v must be at least nrows={} * 16 doubles but is {}",
178            nrows,
179            v.len()
180        );
181    }
182
183    // zero accumulators
184    let mut acc: [f64; 8] = [0f64; 8];
185    for i in 0..nrows {
186        let _1j: usize = i << 3;
187        let _2j: usize = i << 4;
188        reim4_add_mul(&mut acc, as_arr(&u[_1j..]), as_arr(&v[_2j + 8..]));
189    }
190    dst[0..8].copy_from_slice(&acc);
191}
192
193#[inline(always)]
194pub fn reim4_add_mul(dst: &mut [f64; 8], a: &[f64; 8], b: &[f64; 8]) {
195    for k in 0..4 {
196        let ar: f64 = a[k];
197        let br: f64 = b[k];
198        let ai: f64 = a[k + 4];
199        let bi: f64 = b[k + 4];
200        dst[k] += ar * br - ai * bi;
201        dst[k + 4] += ar * bi + ai * br;
202    }
203}