1use crate::reference::fft64::reim::{as_arr, as_arr_mut, reim_zero_ref};
  2
  3#[inline(always)]
  4pub fn reim4_extract_1blk_from_reim_contiguous_ref(m: usize, rows: usize, blk: usize, dst: &mut [f64], src: &[f64]) {
  5    debug_assert!(blk < (m >> 2));
  6    debug_assert!(dst.len() >= 2 * rows * 4);
  7
  8    let offset: usize = blk << 2;
  9
 10    // src = 4-values chunks spaced by m, dst = sequential 4-values chunks
 11    let src_rows = src.chunks_exact(m).take(2 * rows);
 12    let dst_chunks = dst.chunks_exact_mut(4).take(2 * rows);
 13
 14    for (dst_chunk, src_row) in dst_chunks.zip(src_rows) {
 15        dst_chunk.copy_from_slice(&src_row[offset..offset + 4]);
 16    }
 17}
 18
 19#[inline(always)]
 20pub fn reim4_save_1blk_to_reim_contiguous_ref(m: usize, rows: usize, blk: usize, dst: &mut [f64], src: &[f64]) {
 21    debug_assert!(blk < (m >> 2));
 22    debug_assert!(src.len() >= 2 * rows * 4);
 23
 24    let offset: usize = blk << 2;
 25
 26    // dst = 4-values chunks spaced by m, src = sequential 4-values chunks
 27    let dst_rows = dst.chunks_exact_mut(m).take(2 * rows);
 28    let src_chunks = src.chunks_exact(4).take(2 * rows);
 29
 30    for (dst_row, src_chunk) in dst_rows.zip(src_chunks) {
 31        dst_row[offset..offset + 4].copy_from_slice(src_chunk);
 32    }
 33}
 34
 35#[inline(always)]
 36pub fn reim4_save_1blk_to_reim_ref<const OVERWRITE: bool>(m: usize, blk: usize, dst: &mut [f64], src: &[f64]) {
 37    let mut offset: usize = blk << 2;
 38
 39    debug_assert!(blk < (m >> 2));
 40    debug_assert!(dst.len() >= offset + m + 4);
 41    debug_assert!(src.len() >= 8);
 42
 43    let dst_off = &mut dst[offset..offset + 4];
 44
 45    if OVERWRITE {
 46        dst_off.copy_from_slice(&src[0..4]);
 47    } else {
 48        dst_off[0] += src[0];
 49        dst_off[1] += src[1];
 50        dst_off[2] += src[2];
 51        dst_off[3] += src[3];
 52    }
 53
 54    offset += m;
 55
 56    let dst_off = &mut dst[offset..offset + 4];
 57    if OVERWRITE {
 58        dst_off.copy_from_slice(&src[4..8]);
 59    } else {
 60        dst_off[0] += src[4];
 61        dst_off[1] += src[5];
 62        dst_off[2] += src[6];
 63        dst_off[3] += src[7];
 64    }
 65}
 66
 67#[inline(always)]
 68pub fn reim4_save_2blk_to_reim_ref<const OVERWRITE: bool>(m: usize, blk: usize, dst: &mut [f64], src: &[f64]) {
 69    let mut offset: usize = blk << 2;
 70
 71    debug_assert!(blk < (m >> 2));
 72    debug_assert!(dst.len() >= offset + 3 * m + 4);
 73    debug_assert!(src.len() >= 16);
 74
 75    let dst_off: &mut [f64] = &mut dst[offset..offset + 4];
 76    if OVERWRITE {
 77        dst_off.copy_from_slice(&src[0..4]);
 78    } else {
 79        dst_off[0] += src[0];
 80        dst_off[1] += src[1];
 81        dst_off[2] += src[2];
 82        dst_off[3] += src[3];
 83    }
 84
 85    offset += m;
 86    let dst_off: &mut [f64] = &mut dst[offset..offset + 4];
 87    if OVERWRITE {
 88        dst_off.copy_from_slice(&src[4..8]);
 89    } else {
 90        dst_off[0] += src[4];
 91        dst_off[1] += src[5];
 92        dst_off[2] += src[6];
 93        dst_off[3] += src[7];
 94    }
 95
 96    offset += m;
 97
 98    let dst_off: &mut [f64] = &mut dst[offset..offset + 4];
 99    if OVERWRITE {
100        dst_off.copy_from_slice(&src[8..12]);
101    } else {
102        dst_off[0] += src[8];
103        dst_off[1] += src[9];
104        dst_off[2] += src[10];
105        dst_off[3] += src[11];
106    }
107
108    offset += m;
109    let dst_off: &mut [f64] = &mut dst[offset..offset + 4];
110    if OVERWRITE {
111        dst_off.copy_from_slice(&src[12..16]);
112    } else {
113        dst_off[0] += src[12];
114        dst_off[1] += src[13];
115        dst_off[2] += src[14];
116        dst_off[3] += src[15];
117    }
118}
119
120#[inline(always)]
121pub fn reim4_vec_mat1col_product_ref(
122    nrows: usize,
123    dst: &mut [f64], // 8 doubles: [re1(4), im1(4)]
124    u: &[f64],       // nrows * 8 doubles: [ur(4) | ui(4)] per row
125    v: &[f64],       // nrows * 8 doubles: [ar(4) | ai(4)] per row
126) {
127    #[cfg(debug_assertions)]
128    {
129        assert!(dst.len() >= 8, "dst must have at least 8 doubles");
130        assert!(u.len() >= nrows * 8, "u must be at least nrows * 8 doubles");
131        assert!(v.len() >= nrows * 8, "v must be at least nrows * 8 doubles");
132    }
133
134    let mut acc: [f64; 8] = [0f64; 8];
135    let mut j = 0;
136    for _ in 0..nrows {
137        reim4_add_mul(&mut acc, as_arr(&u[j..]), as_arr(&v[j..]));
138        j += 8;
139    }
140    dst[0..8].copy_from_slice(&acc);
141}
142
143#[inline(always)]
144pub fn reim4_vec_mat2cols_product_ref(
145    nrows: usize,
146    dst: &mut [f64], // 16 doubles: [re1(4), im1(4), re2(4), im2(4)]
147    u: &[f64],       // nrows * 8 doubles: [ur(4) | ui(4)] per row
148    v: &[f64],       // nrows * 16 doubles: [ar(4) | ai(4) | br(4) | bi(4)] per row
149) {
150    #[cfg(debug_assertions)]
151    {
152        assert_eq!(dst.len(), 16, "dst must have 16 doubles");
153        assert!(u.len() >= nrows * 8, "u must be at least nrows * 8 doubles");
154        assert!(v.len() >= nrows * 16, "v must be at least nrows * 16 doubles");
155    }
156
157    // zero accumulators
158    let mut acc_0: [f64; 8] = [0f64; 8];
159    let mut acc_1: [f64; 8] = [0f64; 8];
160    for i in 0..nrows {
161        let _1j: usize = i << 3;
162        let _2j: usize = i << 4;
163        let u_j: &[f64; 8] = as_arr(&u[_1j..]);
164        reim4_add_mul(&mut acc_0, u_j, as_arr(&v[_2j..]));
165        reim4_add_mul(&mut acc_1, u_j, as_arr(&v[_2j + 8..]));
166    }
167    dst[0..8].copy_from_slice(&acc_0);
168    dst[8..16].copy_from_slice(&acc_1);
169}
170
171#[inline(always)]
172pub fn reim4_vec_mat2cols_2ndcol_product_ref(
173    nrows: usize,
174    dst: &mut [f64], // 8 doubles: [re1(4), im1(4), re2(4), im2(4)]
175    u: &[f64],       // nrows * 8 doubles: [ur(4) | ui(4)] per row
176    v: &[f64],       // nrows * 16 doubles: [x | x | br(4) | bi(4)] per row
177) {
178    #[cfg(debug_assertions)]
179    {
180        assert!(dst.len() >= 8, "dst must be at least 8 doubles but is {}", dst.len());
181        assert!(
182            u.len() >= nrows * 8,
183            "u must be at least nrows={} * 8 doubles but is {}",
184            nrows,
185            u.len()
186        );
187        assert!(
188            v.len() >= nrows * 16,
189            "v must be at least nrows={} * 16 doubles but is {}",
190            nrows,
191            v.len()
192        );
193    }
194
195    // zero accumulators
196    let mut acc: [f64; 8] = [0f64; 8];
197    for i in 0..nrows {
198        let _1j: usize = i << 3;
199        let _2j: usize = i << 4;
200        reim4_add_mul(&mut acc, as_arr(&u[_1j..]), as_arr(&v[_2j + 8..]));
201    }
202    dst[0..8].copy_from_slice(&acc);
203}
204
205#[inline(always)]
206pub fn reim4_add_mul(dst: &mut [f64; 8], a: &[f64; 8], b: &[f64; 8]) {
207    for k in 0..4 {
208        let ar: f64 = a[k];
209        let br: f64 = b[k];
210        let ai: f64 = a[k + 4];
211        let bi: f64 = b[k + 4];
212        dst[k] += ar * br - ai * bi;
213        dst[k + 4] += ar * bi + ai * br;
214    }
215}
216
217#[inline(always)]
218pub fn reim4_convolution_1coeff_ref(k: usize, dst: &mut [f64; 8], a: &[f64], a_size: usize, b: &[f64], b_size: usize) {
219    reim_zero_ref(dst);
220
221    if k >= a_size + b_size {
222        return;
223    }
224    let j_min: usize = k.saturating_sub(a_size - 1);
225    let j_max: usize = (k + 1).min(b_size);
226
227    for j in j_min..j_max {
228        reim4_add_mul(dst, as_arr(&a[8 * (k - j)..]), as_arr(&b[8 * j..]));
229    }
230}
231
232#[inline(always)]
233pub fn reim4_convolution_2coeffs_ref(k: usize, dst: &mut [f64; 16], a: &[f64], a_size: usize, b: &[f64], b_size: usize) {
234    reim4_convolution_1coeff_ref(k, as_arr_mut(dst), a, a_size, b, b_size);
235    reim4_convolution_1coeff_ref(k + 1, as_arr_mut(&mut dst[8..]), a, a_size, b, b_size);
236}
237
238#[inline(always)]
239pub fn reim4_add_mul_b_real_const(dst: &mut [f64; 8], a: &[f64; 8], b: f64) {
240    for k in 0..4 {
241        let ar: f64 = a[k];
242        let ai: f64 = a[k + 4];
243        dst[k] += ar * b;
244        dst[k + 4] += ai * b;
245    }
246}
247
248#[inline(always)]
249pub fn reim4_convolution_by_real_const_1coeff_ref(k: usize, dst: &mut [f64; 8], a: &[f64], a_size: usize, b: &[f64]) {
250    reim_zero_ref(dst);
251
252    let b_size: usize = b.len();
253
254    if k >= a_size + b_size {
255        return;
256    }
257    let j_min: usize = k.saturating_sub(a_size - 1);
258    let j_max: usize = (k + 1).min(b_size);
259
260    for j in j_min..j_max {
261        reim4_add_mul_b_real_const(dst, as_arr(&a[8 * (k - j)..]), b[j]);
262    }
263}
264
265#[inline(always)]
266pub fn reim4_convolution_by_real_const_2coeffs_ref(k: usize, dst: &mut [f64; 16], a: &[f64], a_size: usize, b: &[f64]) {
267    reim4_convolution_by_real_const_1coeff_ref(k, as_arr_mut(dst), a, a_size, b);
268    reim4_convolution_by_real_const_1coeff_ref(k + 1, as_arr_mut(&mut dst[8..]), a, a_size, b);
269}