1#include <string.h>
  2
  3#include "../coeffs/coeffs_arithmetic.h"
  4#include "vec_rnx_arithmetic_private.h"
  5
  6void rnx_add_ref(uint64_t nn, double* res, const double* a, const double* b) {
  7  for (uint64_t i = 0; i < nn; ++i) {
  8    res[i] = a[i] + b[i];
  9  }
 10}
 11
 12void rnx_sub_ref(uint64_t nn, double* res, const double* a, const double* b) {
 13  for (uint64_t i = 0; i < nn; ++i) {
 14    res[i] = a[i] - b[i];
 15  }
 16}
 17
 18void rnx_negate_ref(uint64_t nn, double* res, const double* a) {
 19  for (uint64_t i = 0; i < nn; ++i) {
 20    res[i] = -a[i];
 21  }
 22}
 23
 24/** @brief sets res = a + b */
 25EXPORT void vec_rnx_add_ref(                          //
 26    const MOD_RNX* module,                            // N
 27    double* res, uint64_t res_size, uint64_t res_sl,  // res
 28    const double* a, uint64_t a_size, uint64_t a_sl,  // a
 29    const double* b, uint64_t b_size, uint64_t b_sl   // b
 30) {
 31  const uint64_t nn = module->n;
 32  if (a_size < b_size) {
 33    const uint64_t msize = res_size < a_size ? res_size : a_size;
 34    const uint64_t nsize = res_size < b_size ? res_size : b_size;
 35    for (uint64_t i = 0; i < msize; ++i) {
 36      rnx_add_ref(nn, res + i * res_sl, a + i * a_sl, b + i * b_sl);
 37    }
 38    for (uint64_t i = msize; i < nsize; ++i) {
 39      memcpy(res + i * res_sl, b + i * b_sl, nn * sizeof(double));
 40    }
 41    for (uint64_t i = nsize; i < res_size; ++i) {
 42      memset(res + i * res_sl, 0, nn * sizeof(double));
 43    }
 44  } else {
 45    const uint64_t msize = res_size < b_size ? res_size : b_size;
 46    const uint64_t nsize = res_size < a_size ? res_size : a_size;
 47    for (uint64_t i = 0; i < msize; ++i) {
 48      rnx_add_ref(nn, res + i * res_sl, a + i * a_sl, b + i * b_sl);
 49    }
 50    for (uint64_t i = msize; i < nsize; ++i) {
 51      memcpy(res + i * res_sl, a + i * a_sl, nn * sizeof(double));
 52    }
 53    for (uint64_t i = nsize; i < res_size; ++i) {
 54      memset(res + i * res_sl, 0, nn * sizeof(double));
 55    }
 56  }
 57}
 58
 59/** @brief sets res = 0 */
 60EXPORT void vec_rnx_zero_ref(                        //
 61    const MOD_RNX* module,                           // N
 62    double* res, uint64_t res_size, uint64_t res_sl  // res
 63) {
 64  const uint64_t nn = module->n;
 65  for (uint64_t i = 0; i < res_size; ++i) {
 66    memset(res + i * res_sl, 0, nn * sizeof(double));
 67  }
 68}
 69
 70/** @brief sets res = a */
 71EXPORT void vec_rnx_copy_ref(                         //
 72    const MOD_RNX* module,                            // N
 73    double* res, uint64_t res_size, uint64_t res_sl,  // res
 74    const double* a, uint64_t a_size, uint64_t a_sl   // a
 75) {
 76  const uint64_t nn = module->n;
 77
 78  const uint64_t rot_end_idx = res_size < a_size ? res_size : a_size;
 79  // rotate up to the smallest dimension
 80  for (uint64_t i = 0; i < rot_end_idx; ++i) {
 81    double* res_ptr = res + i * res_sl;
 82    const double* a_ptr = a + i * a_sl;
 83    memcpy(res_ptr, a_ptr, nn * sizeof(double));
 84  }
 85  // then extend with zeros
 86  for (uint64_t i = rot_end_idx; i < res_size; ++i) {
 87    memset(res + i * res_sl, 0, nn * sizeof(double));
 88  }
 89}
 90
 91/** @brief sets res = -a */
 92EXPORT void vec_rnx_negate_ref(                       //
 93    const MOD_RNX* module,                            // N
 94    double* res, uint64_t res_size, uint64_t res_sl,  // res
 95    const double* a, uint64_t a_size, uint64_t a_sl   // a
 96) {
 97  const uint64_t nn = module->n;
 98
 99  const uint64_t rot_end_idx = res_size < a_size ? res_size : a_size;
100  // rotate up to the smallest dimension
101  for (uint64_t i = 0; i < rot_end_idx; ++i) {
102    double* res_ptr = res + i * res_sl;
103    const double* a_ptr = a + i * a_sl;
104    rnx_negate_ref(nn, res_ptr, a_ptr);
105  }
106  // then extend with zeros
107  for (uint64_t i = rot_end_idx; i < res_size; ++i) {
108    memset(res + i * res_sl, 0, nn * sizeof(double));
109  }
110}
111
112/** @brief sets res = a - b */
113EXPORT void vec_rnx_sub_ref(                          //
114    const MOD_RNX* module,                            // N
115    double* res, uint64_t res_size, uint64_t res_sl,  // res
116    const double* a, uint64_t a_size, uint64_t a_sl,  // a
117    const double* b, uint64_t b_size, uint64_t b_sl   // b
118) {
119  const uint64_t nn = module->n;
120  if (a_size < b_size) {
121    const uint64_t msize = res_size < a_size ? res_size : a_size;
122    const uint64_t nsize = res_size < b_size ? res_size : b_size;
123    for (uint64_t i = 0; i < msize; ++i) {
124      rnx_sub_ref(nn, res + i * res_sl, a + i * a_sl, b + i * b_sl);
125    }
126    for (uint64_t i = msize; i < nsize; ++i) {
127      rnx_negate_ref(nn, res + i * res_sl, b + i * b_sl);
128    }
129    for (uint64_t i = nsize; i < res_size; ++i) {
130      memset(res + i * res_sl, 0, nn * sizeof(double));
131    }
132  } else {
133    const uint64_t msize = res_size < b_size ? res_size : b_size;
134    const uint64_t nsize = res_size < a_size ? res_size : a_size;
135    for (uint64_t i = 0; i < msize; ++i) {
136      rnx_sub_ref(nn, res + i * res_sl, a + i * a_sl, b + i * b_sl);
137    }
138    for (uint64_t i = msize; i < nsize; ++i) {
139      memcpy(res + i * res_sl, a + i * a_sl, nn * sizeof(double));
140    }
141    for (uint64_t i = nsize; i < res_size; ++i) {
142      memset(res + i * res_sl, 0, nn * sizeof(double));
143    }
144  }
145}
146
147/** @brief sets res = a . X^p */
148EXPORT void vec_rnx_rotate_ref(                       //
149    const MOD_RNX* module,                            // N
150    const int64_t p,                                  // rotation value
151    double* res, uint64_t res_size, uint64_t res_sl,  // res
152    const double* a, uint64_t a_size, uint64_t a_sl   // a
153) {
154  const uint64_t nn = module->n;
155
156  const uint64_t rot_end_idx = res_size < a_size ? res_size : a_size;
157  // rotate up to the smallest dimension
158  for (uint64_t i = 0; i < rot_end_idx; ++i) {
159    double* res_ptr = res + i * res_sl;
160    const double* a_ptr = a + i * a_sl;
161    if (res_ptr == a_ptr) {
162      rnx_rotate_inplace_f64(nn, p, res_ptr);
163    } else {
164      rnx_rotate_f64(nn, p, res_ptr, a_ptr);
165    }
166  }
167  // then extend with zeros
168  for (uint64_t i = rot_end_idx; i < res_size; ++i) {
169    memset(res + i * res_sl, 0, nn * sizeof(double));
170  }
171}
172
173/** @brief sets res = a(X^p) */
174EXPORT void vec_rnx_automorphism_ref(                 //
175    const MOD_RNX* module,                            // N
176    int64_t p,                                        // X -> X^p
177    double* res, uint64_t res_size, uint64_t res_sl,  // res
178    const double* a, uint64_t a_size, uint64_t a_sl   // a
179) {
180  const uint64_t nn = module->n;
181
182  const uint64_t rot_end_idx = res_size < a_size ? res_size : a_size;
183  // rotate up to the smallest dimension
184  for (uint64_t i = 0; i < rot_end_idx; ++i) {
185    double* res_ptr = res + i * res_sl;
186    const double* a_ptr = a + i * a_sl;
187    if (res_ptr == a_ptr) {
188      rnx_automorphism_inplace_f64(nn, p, res_ptr);
189    } else {
190      rnx_automorphism_f64(nn, p, res_ptr, a_ptr);
191    }
192  }
193  // then extend with zeros
194  for (uint64_t i = rot_end_idx; i < res_size; ++i) {
195    memset(res + i * res_sl, 0, nn * sizeof(double));
196  }
197}
198
199/** @brief sets res = a . (X^p - 1) */
200EXPORT void vec_rnx_mul_xp_minus_one_ref(             //
201    const MOD_RNX* module,                            // N
202    const int64_t p,                                  // rotation value
203    double* res, uint64_t res_size, uint64_t res_sl,  // res
204    const double* a, uint64_t a_size, uint64_t a_sl   // a
205) {
206  const uint64_t nn = module->n;
207
208  const uint64_t rot_end_idx = res_size < a_size ? res_size : a_size;
209  // rotate up to the smallest dimension
210  for (uint64_t i = 0; i < rot_end_idx; ++i) {
211    double* res_ptr = res + i * res_sl;
212    const double* a_ptr = a + i * a_sl;
213    if (res_ptr == a_ptr) {
214      rnx_mul_xp_minus_one_inplace_f64(nn, p, res_ptr);
215    } else {
216      rnx_mul_xp_minus_one_f64(nn, p, res_ptr, a_ptr);
217    }
218  }
219  // then extend with zeros
220  for (uint64_t i = rot_end_idx; i < res_size; ++i) {
221    memset(res + i * res_sl, 0, nn * sizeof(double));
222  }
223}