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