1#ifndef SPQLIOS_REIM4_ARITHMETIC_H
  2#define SPQLIOS_REIM4_ARITHMETIC_H
  3
  4#include "../commons.h"
  5
  6// the reim4 structure represent 4 complex numbers,
  7// represented as re0,re1,re2,re3,im0,im1,im2,im3
  8
  9// TODO implement these basic arithmetic functions. (only ref is needed)
 10
 11/** @brief dest = 0 */
 12EXPORT void reim4_zero(double* dest);
 13/** @brief dest = a + b */
 14EXPORT void reim4_add(double* dest, const double* u, const double* v);
 15/** @brief dest = a * b */
 16EXPORT void reim4_mul(double* dest, const double* u, const double* v);
 17/** @brief dest += a * b */
 18EXPORT void reim4_add_mul(double* dest, const double* a, const double* b);
 19
 20// TODO add the dot products: vec x mat1col, vec x mat2cols functions here (ref and avx)
 21
 22// TODO implement the convolution functions here (ref and avx)
 23/**
 24 * @brief k-th coefficient of the convolution product a * b.
 25 *
 26 * The k-th coefficient is defined as sum_{i+j=k} a[i].b[j].
 27 * (in this sum, i and j must remain within the bounds [0,sizea[ and [0,sizeb[)
 28 *
 29 * In practice, accounting for these bounds, the convolution function boils down to
 30 * ```
 31 * res := 0
 32 * if (k<sizeb+sizea-1) then
 33 *   for j from max(0,k+1-sizea) incl to min(sizeb,k+1) excl do
 34 *     res += a[k-j] * b[j]
 35 * return res
 36 * ```
 37 */
 38EXPORT void reim4_convolution_1coeff_ref(uint64_t k, double* dest, const double* a, uint64_t sizea, const double* b,
 39                                         uint64_t sizeb);
 40EXPORT void reim4_convolution_1coeff_avx(uint64_t k, double* dest, const double* a, uint64_t sizea, const double* b,
 41                                         uint64_t sizeb);
 42/** @brief returns two consecutive convolution coefficients: k and k+1 */
 43EXPORT void reim4_convolution_2coeff_ref(uint64_t k, double* dest, const double* a, uint64_t sizea, const double* b,
 44                                         uint64_t sizeb);
 45EXPORT void reim4_convolution_2coeff_avx(uint64_t k, double* dest, const double* a, uint64_t sizea, const double* b,
 46                                         uint64_t sizeb);
 47
 48/**
 49 * @brief From the convolution a * b, return the coefficients between offset and offset + size
 50 * For the full convolution, use offset=0 and size=sizea+sizeb-1.
 51 */
 52EXPORT void reim4_convolution_ref(double* dest, uint64_t dest_size, uint64_t dest_offset, const double* a,
 53                                  uint64_t sizea, const double* b, uint64_t sizeb);
 54EXPORT void reim4_convolution_avx(double* dest, uint64_t dest_size, uint64_t dest_offset, const double* a,
 55                                  uint64_t sizea, const double* b, uint64_t sizeb);
 56
 57// The reim fft layout that encodes m = nn/2 complex evaluations
 58// can be decomposed into m/4 reim4 blocks
 59// the i-th block of v is noted v(i);
 60// it contains the 4 evaluations (e[4i],e[4i+1],e[4i+2],e[4i+3])
 61
 62/**
 63 * @brief extract 1 reim4 block from one reim vectors of m complexes
 64 * @param m the size of each reim
 65 * @param blk the block id to extract (<m/4)
 66 * @param dst the output: nrows reim4's  dst[i] = src[i](blk)
 67 * @param src the input: nrows reim's
 68 */
 69EXPORT void reim4_extract_1blk_from_reim_ref(uint64_t m, uint64_t blk,
 70                                             double* const dst,       // 8 doubles
 71                                             const double* const src  // a reim vector
 72);
 73EXPORT void reim4_extract_1blk_from_reim_avx(uint64_t m, uint64_t blk,
 74                                             double* const dst,       // 8 doubles
 75                                             const double* const src  // a reim vectors
 76);
 77
 78/**
 79 * @brief extract 1 reim vectors from one reim4 block of m complexes
 80 * @param m the size of each reim
 81 * @param blk the block id to extract (<m/4)
 82 * @param dst the output: nrows reim4's  dst[i](blk) = src[i]
 83 * @param src the input: nrows reim's
 84 */
 85EXPORT void reim4_extract_reim_from_1blk_ref(uint64_t m, uint64_t blk,
 86                                             double* const dst,       // 8 doubles
 87                                             const double* const src  // a reim vector
 88);
 89
 90/**
 91 * @brief extract 1 reim4 block from nrows reim vectors of m complexes
 92 * @param nrows the number of reim (fft) vectors
 93 * @param m the size of each reim
 94 * @param blk the block id to extract (<m/4)
 95 * @param dst the output: nrows reim4's  dst[i] = src[i](blk)
 96 * @param src the input: nrows reim's
 97 */
 98EXPORT void reim4_extract_1blk_from_contiguous_reim_ref(uint64_t m, uint64_t nrows, uint64_t blk, double* const dst,
 99                                                        const double* const src);
100EXPORT void reim4_extract_1blk_from_contiguous_reim_avx(uint64_t m, uint64_t nrows, uint64_t blk, double* const dst,
101                                                        const double* const src);
102
103/**
104 * @brief extract 1 reim4 block from nrows reim vectors of m complexes with slice sl
105 * @param m the size of each reim
106 * @param sl the slice size
107 * @param nrows the number of reim (fft) vectors to extract
108 * @param blk the block id to extract (<m/4)
109 * @param dst the output: nrows reim4's  dst[i] = src[i](blk)
110 * @param src the input: nrows reim's
111 */
112EXPORT void reim4_extract_1blk_from_contiguous_reim_sl_ref(uint64_t m, uint64_t sl, uint64_t nrows, uint64_t blk,
113                                                           double* const dst, const double* const src);
114EXPORT void reim4_extract_1blk_from_contiguous_reim_sl_avx(uint64_t m, uint64_t sl, uint64_t nrows, uint64_t blk,
115                                                           double* const dst, const double* const src);
116
117/**
118 * @brief saves 1 single reim4 block in a reim vectors of m complexes
119 * @param m the size of each reim
120 * @param blk the block id to save (<m/4)
121 * @param dest the output reim: dst(blk) = src
122 * @param src the input reim4
123 */
124EXPORT void reim4_save_1blk_to_reim_ref(uint64_t m, uint64_t blk,
125                                        double* dest,      // 1 reim vector of length m
126                                        const double* src  // 8 doubles
127);
128EXPORT void reim4_save_1blk_to_reim_avx(uint64_t m, uint64_t blk,
129                                        double* dest,      // 1 reim vector of length m
130                                        const double* src  // 8 doubles
131);
132
133/**
134 * @brief Adds 1 single reim4 block in a reim vectors of m complexes
135 * @param m the size of each reim
136 * @param blk the block id to save (<m/4)
137 * @param dest the output reim: dst(blk) += src
138 * @param src the input reim4
139 */
140EXPORT void reim4_add_1blk_to_reim_ref(uint64_t m, uint64_t blk,
141                                       double* dest,      // 1 reim vector of length m
142                                       const double* src  // 8 doubles
143);
144EXPORT void reim4_add_1blk_to_reim_avx(uint64_t m, uint64_t blk,
145                                       double* dest,      // 1 reim vector of length m
146                                       const double* src  // 8 doubles
147);
148
149// dest = sum u[i].v[i]
150// make ref and avx2 implem, needs to be fast
151// can accumulate the real and imaginary parts on 2 separate registers, and combine at the end
152EXPORT void reim4_vec_mat1col_product_ref(const uint64_t nrows,
153                                          double* const dst,      // 8 doubles
154                                          const double* const u,  // nrows * 8 doubles
155                                          const double* const v   // nrows * 8 doubles
156);
157
158EXPORT void reim4_vec_mat1col_product_avx2(const uint64_t nrows,
159                                           double* const dst,      // 8 doubles
160                                           const double* const u,  // nrows * 8 doubles
161                                           const double* const v   // nrows * 8 doubles
162);
163
164// dest = sum u[i].v[i]
165// make ref and avx2 implem, needs to be fast
166// can accumulate the real and imaginary parts on 2 separate registers, and recombine at the end
167EXPORT void reim4_vec_mat2cols_product_ref(uint64_t nrows, double* const dst,  // 16 doubles
168                                           const double* const u,              // nrows * 8 doubles
169                                           const double* const v               // nrows * 16 doubles
170);
171
172EXPORT void reim4_vec_mat2cols_product_avx2(uint64_t nrows, double* const dst,  // 16 doubles
173                                            const double* const u,              // nrows * 8 doubles
174                                            const double* const v               // nrows * 16 doubles
175);
176
177#endif  // SPQLIOS_REIM4_ARITHMETIC_H