1#include <assert.h>
  2#include <string.h>
  3
  4#include "reim4_arithmetic.h"
  5
  6// Stores the first 4 values (RE) of src + blk*4 and he first 4 values (IM) of src + blk*4 + m
  7// contiguously into dst
  8void reim4_extract_1blk_from_reim_ref(uint64_t m, uint64_t blk,
  9                                      double* const dst,       // 8 doubles
 10                                      const double* const src  // one reim4 vector
 11) {
 12  assert(blk < (m >> 2));
 13  const double* src_ptr = src + (blk << 2);
 14  // copy the real parts
 15  dst[0] = src_ptr[0];
 16  dst[1] = src_ptr[1];
 17  dst[2] = src_ptr[2];
 18  dst[3] = src_ptr[3];
 19  src_ptr += m;
 20  // copy the imaginary parts
 21  dst[4] = src_ptr[0];
 22  dst[5] = src_ptr[1];
 23  dst[6] = src_ptr[2];
 24  dst[7] = src_ptr[3];
 25}
 26
 27void reim4_extract_reim_from_1blk_ref(uint64_t m, uint64_t blk,
 28                                      double* const dst,       // 8 doubles
 29                                      const double* const src  // one reim4 vector
 30) {
 31  assert(blk < (m >> 2));
 32  double* dst_ptr = dst + (blk << 2);
 33  // copy the real parts
 34  dst_ptr[0] = src[0];
 35  dst_ptr[1] = src[1];
 36  dst_ptr[2] = src[2];
 37  dst_ptr[3] = src[3];
 38  dst_ptr += m;
 39  // copy the imaginary parts
 40  dst_ptr[0] = src[4];
 41  dst_ptr[1] = src[5];
 42  dst_ptr[2] = src[6];
 43  dst_ptr[3] = src[7];
 44}
 45
 46EXPORT void reim4_extract_1blk_from_contiguous_reim_ref(uint64_t m, uint64_t nrows, uint64_t blk, double* const dst,
 47                                                        const double* const src) {
 48  assert(blk < (m >> 2));
 49
 50  const double* src_ptr = src + (blk << 2);
 51  double* dst_ptr = dst;
 52  for (uint64_t i = 0; i < nrows * 2; ++i) {
 53    dst_ptr[0] = src_ptr[0];
 54    dst_ptr[1] = src_ptr[1];
 55    dst_ptr[2] = src_ptr[2];
 56    dst_ptr[3] = src_ptr[3];
 57    dst_ptr += 4;
 58    src_ptr += m;
 59  }
 60}
 61
 62EXPORT void reim4_extract_1blk_from_contiguous_reim_sl_ref(uint64_t m, uint64_t sl, uint64_t nrows, uint64_t blk,
 63                                                           double* const dst, const double* const src) {
 64  assert(blk < (m >> 2));
 65
 66  const double* src_ptr = src + (blk << 2);
 67  double* dst_ptr = dst;
 68  const uint64_t sl_minus_m = sl - m;
 69  for (uint64_t i = 0; i < nrows; ++i) {
 70    dst_ptr[0] = src_ptr[0];
 71    dst_ptr[1] = src_ptr[1];
 72    dst_ptr[2] = src_ptr[2];
 73    dst_ptr[3] = src_ptr[3];
 74    dst_ptr += 4;
 75    src_ptr += m;
 76    dst_ptr[0] = src_ptr[0];
 77    dst_ptr[1] = src_ptr[1];
 78    dst_ptr[2] = src_ptr[2];
 79    dst_ptr[3] = src_ptr[3];
 80    dst_ptr += 4;
 81    src_ptr += sl_minus_m;
 82  }
 83}
 84
 85// dest(i)=src
 86// use scalar or sse or avx? Code
 87// should be the inverse of reim4_extract_1col_from_reim
 88void reim4_save_1blk_to_reim_ref(uint64_t m, uint64_t blk,
 89                                 double* dst,       // 1 reim vector of length m
 90                                 const double* src  // 8 doubles
 91) {
 92  assert(blk < (m >> 2));
 93  double* dst_ptr = dst + (blk << 2);
 94  // save the real part
 95  dst_ptr[0] = src[0];
 96  dst_ptr[1] = src[1];
 97  dst_ptr[2] = src[2];
 98  dst_ptr[3] = src[3];
 99  dst_ptr += m;
100  // save the imag part
101  dst_ptr[0] = src[4];
102  dst_ptr[1] = src[5];
103  dst_ptr[2] = src[6];
104  dst_ptr[3] = src[7];
105}
106
107void reim4_add_1blk_to_reim_ref(uint64_t m, uint64_t blk,
108                                double* dst,       // 1 reim vector of length m
109                                const double* src  // 8 doubles
110) {
111  assert(blk < (m >> 2));
112  double* dst_ptr = dst + (blk << 2);
113  // add the real part
114  dst_ptr[0] += src[0];
115  dst_ptr[1] += src[1];
116  dst_ptr[2] += src[2];
117  dst_ptr[3] += src[3];
118  dst_ptr += m;
119  // add the imag part
120  dst_ptr[0] += src[4];
121  dst_ptr[1] += src[5];
122  dst_ptr[2] += src[6];
123  dst_ptr[3] += src[7];
124}
125
126// dest = 0
127void reim4_zero(double* const dst  // 8 doubles
128) {
129  for (uint64_t i = 0; i < 8; ++i) dst[i] = 0;
130}
131
132/** @brief dest = a + b */
133void reim4_add(double* const dst,      // 8 doubles
134               const double* const u,  // nrows * 8 doubles
135               const double* const v   // nrows * 8 doubles
136) {
137  for (uint64_t k = 0; k < 4; ++k) {
138    const double a = u[k];
139    const double c = v[k];
140    const double b = u[k + 4];
141    const double d = v[k + 4];
142
143    dst[k] = a + c;
144    dst[k + 4] = b + d;
145  }
146}
147
148/** @brief dest = a * b */
149void reim4_mul(double* const dst,      // 8 doubles
150               const double* const u,  // 8 doubles
151               const double* const v   // 8 doubles
152) {
153  for (uint64_t k = 0; k < 4; ++k) {
154    const double a = u[k];
155    const double c = v[k];
156    const double b = u[k + 4];
157    const double d = v[k + 4];
158
159    dst[k] = a * c - b * d;
160    dst[k + 4] = a * d + b * c;
161  }
162}
163
164/** @brief dest += a * b */
165void reim4_add_mul(double* const dst,      // 8 doubles
166                   const double* const u,  // 8 doubles
167                   const double* const v   // 8 doubles
168) {
169  for (uint64_t k = 0; k < 4; ++k) {
170    const double a = u[k];
171    const double c = v[k];
172    const double b = u[k + 4];
173    const double d = v[k + 4];
174
175    dst[k] += a * c - b * d;
176    dst[k + 4] += a * d + b * c;
177  }
178}
179
180/** dest = uT * v  where u is a vector of size nrows, and v is a nrows x 1 matrix */
181void reim4_vec_mat1col_product_ref(const uint64_t nrows,
182                                   double* const dst,      // 8 doubles
183                                   const double* const u,  // nrows * 8 doubles
184                                   const double* const v   // nrows * 8 doubles
185) {
186  reim4_zero(dst);
187
188  for (uint64_t i = 0, j = 0; i < nrows; ++i, j += 8) {
189    reim4_add_mul(dst, u + j, v + j);
190  }
191}
192
193/** dest = uT * v  where u is a vector of size nrows, and v is a nrows x 2 matrix */
194void reim4_vec_mat2cols_product_ref(const uint64_t nrows,
195                                    double* const dst,      // 16 doubles
196                                    const double* const u,  // nrows * 8 doubles
197                                    const double* const v   // nrows * 16 doubles
198) {
199  reim4_zero(dst);
200  reim4_zero(dst + 8);
201
202  double* dst1 = dst;
203  double* dst2 = dst + 8;
204
205  for (uint64_t i = 0, j = 0; i < nrows; ++i, j += 8) {
206    uint64_t double_j = j << 1;
207    reim4_add_mul(dst1, u + j, v + double_j);
208    reim4_add_mul(dst2, u + j, v + double_j + 8);
209  }
210}
211
212/**
213 * @brief k-th coefficient of the convolution product a * b.
214 *
215 * The k-th coefficient is defined as sum_{i+j=k} a[i].b[j].
216 * (in this sum, i and j must remain within the bounds [0,sizea[ and [0,sizeb[)
217 *
218 * In practice, accounting for these bounds, the convolution function boils down to
219 * ```
220 * res := 0
221 * if (k<sizeb+sizea-1) then
222 *   for j from max(0,k+1-sizea) incl to min(sizeb,k+1) excl do
223 *     res += a[k-j] * b[j]
224 * return res
225 * ```
226 */
227EXPORT void reim4_convolution_1coeff_ref(uint64_t k, double* dest, const double* a, uint64_t sizea, const double* b,
228                                         uint64_t sizeb) {
229  reim4_zero(dest);
230  if (k >= sizea + sizeb) return;
231  uint64_t jmin = k >= sizea ? k + 1 - sizea : 0;
232  uint64_t jmax = k < sizeb ? k + 1 : sizeb;
233  for (uint64_t j = jmin; j < jmax; ++j) {
234    reim4_add_mul(dest, a + 8 * (k - j), b + 8 * j);
235  }
236}
237
238/** @brief returns two consecutive convolution coefficients: k and k+1 */
239EXPORT void reim4_convolution_2coeff_ref(uint64_t k, double* dest, const double* a, uint64_t sizea, const double* b,
240                                         uint64_t sizeb) {
241  reim4_convolution_1coeff_ref(k, dest, a, sizea, b, sizeb);
242  reim4_convolution_1coeff_ref(k + 1, dest + 8, a, sizea, b, sizeb);
243}
244
245/**
246 * @brief From the convolution a * b, return the coefficients between offset and offset + size
247 * For the full convolution, use offset=0 and size=sizea+sizeb-1.
248 */
249EXPORT void reim4_convolution_ref(double* dest, uint64_t dest_size, uint64_t dest_offset, const double* a,
250                                  uint64_t sizea, const double* b, uint64_t sizeb) {
251  for (uint64_t k = 0; k < dest_size; ++k) {
252    reim4_convolution_1coeff_ref(k + dest_offset, dest + 8 * k, a, sizea, b, sizeb);
253  }
254}