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
125#[inline(always)]
126pub fn reim4_vec_mat2cols_product_ref(
127    nrows: usize,
128    dst: &mut [f64], // 16 doubles: [re1(4), im1(4), re2(4), im2(4)]
129    u: &[f64],       // nrows * 8 doubles: [ur(4) | ui(4)] per row
130    v: &[f64],       // nrows * 16 doubles: [ar(4) | ai(4) | br(4) | bi(4)] per row
131) {
132    #[cfg(debug_assertions)]
133    {
134        assert_eq!(dst.len(), 16, "dst must have 16 doubles");
135        assert!(u.len() >= nrows * 8, "u must be at least nrows * 8 doubles");
136        assert!(
137            v.len() >= nrows * 16,
138            "v must be at least nrows * 16 doubles"
139        );
140    }
141
142    // zero accumulators
143    let mut acc_0: [f64; 8] = [0f64; 8];
144    let mut acc_1: [f64; 8] = [0f64; 8];
145    for i in 0..nrows {
146        let _1j: usize = i << 3;
147        let _2j: usize = i << 4;
148        let u_j: &[f64; 8] = as_arr(&u[_1j..]);
149        reim4_add_mul(&mut acc_0, u_j, as_arr(&v[_2j..]));
150        reim4_add_mul(&mut acc_1, u_j, as_arr(&v[_2j + 8..]));
151    }
152    dst[0..8].copy_from_slice(&acc_0);
153    dst[8..16].copy_from_slice(&acc_1);
154}
155
156#[inline(always)]
157pub fn reim4_vec_mat2cols_2ndcol_product_ref(
158    nrows: usize,
159    dst: &mut [f64], // 8 doubles: [re1(4), im1(4), re2(4), im2(4)]
160    u: &[f64],       // nrows * 8 doubles: [ur(4) | ui(4)] per row
161    v: &[f64],       // nrows * 16 doubles: [x | x | br(4) | bi(4)] per row
162) {
163    #[cfg(debug_assertions)]
164    {
165        assert!(
166            dst.len() >= 8,
167            "dst must be at least 8 doubles but is {}",
168            dst.len()
169        );
170        assert!(
171            u.len() >= nrows * 8,
172            "u must be at least nrows={} * 8 doubles but is {}",
173            nrows,
174            u.len()
175        );
176        assert!(
177            v.len() >= nrows * 16,
178            "v must be at least nrows={} * 16 doubles but is {}",
179            nrows,
180            v.len()
181        );
182    }
183
184    // zero accumulators
185    let mut acc: [f64; 8] = [0f64; 8];
186    for i in 0..nrows {
187        let _1j: usize = i << 3;
188        let _2j: usize = i << 4;
189        reim4_add_mul(&mut acc, as_arr(&u[_1j..]), as_arr(&v[_2j + 8..]));
190    }
191    dst[0..8].copy_from_slice(&acc);
192}
193
194#[inline(always)]
195pub fn reim4_add_mul(dst: &mut [f64; 8], a: &[f64; 8], b: &[f64; 8]) {
196    for k in 0..4 {
197        let ar: f64 = a[k];
198        let br: f64 = b[k];
199        let ai: f64 = a[k + 4];
200        let bi: f64 = b[k + 4];
201        dst[k] += ar * br - ai * bi;
202        dst[k + 4] += ar * bi + ai * br;
203    }
204}