1#ifndef SPQLIOS_NEGACYCLIC_POLYNOMIAL_IMPL_H
  2#define SPQLIOS_NEGACYCLIC_POLYNOMIAL_IMPL_H
  3
  4#include "negacyclic_polynomial.h"
  5
  6template <typename T>
  7polynomial<T>::polynomial(uint64_t N, const T* c) : coeffs(N) {
  8  for (uint64_t i = 0; i < N; ++i) coeffs[i] = c[i];
  9}
 10/** @brief zero polynomial of dimension N */
 11template <typename T>
 12polynomial<T>::polynomial(uint64_t N) : coeffs(N, 0) {}
 13/** @brief empty polynomial (dim 0) */
 14template <typename T>
 15polynomial<T>::polynomial() {}
 16
 17/** @brief ring dimension */
 18template <typename T>
 19uint64_t polynomial<T>::nn() const {
 20  uint64_t n = coeffs.size();
 21  REQUIRE_DRAMATICALLY(is_pow2(n), "polynomial dim is not a pow of 2");
 22  return n;
 23}
 24
 25/** @brief special setter (accept any indexes, and does the negacyclic translation) */
 26template <typename T>
 27void polynomial<T>::set_coeff(int64_t i, T v) {
 28  const uint64_t n = nn();
 29  const uint64_t _2nm = 2 * n - 1;
 30  uint64_t pos = uint64_t(i) & _2nm;
 31  if (pos < n) {
 32    coeffs[pos] = v;
 33  } else {
 34    coeffs[pos - n] = -v;
 35  }
 36}
 37/** @brief special getter (accept any indexes, and does the negacyclic translation) */
 38template <typename T>
 39T polynomial<T>::get_coeff(int64_t i) const {
 40  const uint64_t n = nn();
 41  const uint64_t _2nm = 2 * n - 1;
 42  uint64_t pos = uint64_t(i) & _2nm;
 43  if (pos < n) {
 44    return coeffs[pos];
 45  } else {
 46    return -coeffs[pos - n];
 47  }
 48}
 49/** @brief returns the coefficient layout */
 50template <typename T>
 51T* polynomial<T>::data() {
 52  return coeffs.data();
 53}
 54
 55template <typename T>
 56void polynomial<T>::save_as(T* dest) const {
 57  const uint64_t n = nn();
 58  for (uint64_t i = 0; i < n; ++i) {
 59    dest[i] = coeffs[i];
 60  }
 61}
 62
 63/** @brief returns the coefficient layout (const version) */
 64template <typename T>
 65const T* polynomial<T>::data() const {
 66  return coeffs.data();
 67}
 68
 69/** @brief returns the coefficient layout (const version) */
 70template <typename T>
 71polynomial<T> polynomial<T>::zero(uint64_t n) {
 72  return polynomial<T>(n);
 73}
 74
 75/** @brief equality operator (used during tests) */
 76template <typename T>
 77bool operator==(const polynomial<T>& a, const polynomial<T>& b) {
 78  uint64_t n = a.nn();
 79  REQUIRE_DRAMATICALLY(b.nn() == n, "wrong dimensions");
 80  for (uint64_t i = 0; i < n; ++i) {
 81    if (a.get_coeff(i) != b.get_coeff(i)) return false;
 82  }
 83  return true;
 84}
 85
 86/** @brief addition operator (used during tests) */
 87template <typename T>
 88polynomial<T> operator+(const polynomial<T>& a, const polynomial<T>& b) {
 89  uint64_t n = a.nn();
 90  REQUIRE_DRAMATICALLY(b.nn() == n, "wrong dimensions");
 91  polynomial<T> res(n);
 92  for (uint64_t i = 0; i < n; ++i) {
 93    res.set_coeff(i, a.get_coeff(i) + b.get_coeff(i));
 94  }
 95  return res;
 96}
 97
 98/** @brief subtraction operator (used during tests) */
 99template <typename T>
100polynomial<T> operator-(const polynomial<T>& a, const polynomial<T>& b) {
101  uint64_t n = a.nn();
102  REQUIRE_DRAMATICALLY(b.nn() == n, "wrong dimensions");
103  polynomial<T> res(n);
104  for (uint64_t i = 0; i < n; ++i) {
105    res.set_coeff(i, a.get_coeff(i) - b.get_coeff(i));
106  }
107  return res;
108}
109
110/** @brief subtraction operator (used during tests) */
111template <typename T>
112polynomial<T> operator-(const polynomial<T>& a) {
113  uint64_t n = a.nn();
114  polynomial<T> res(n);
115  for (uint64_t i = 0; i < n; ++i) {
116    res.set_coeff(i, -a.get_coeff(i));
117  }
118  return res;
119}
120
121/** @brief random polynomial */
122template <typename T>
123polynomial<T> random_polynomial(uint64_t n);
124
125/** @brief random int64 polynomial */
126template <>
127polynomial<int64_t> random_polynomial(uint64_t n) {
128  polynomial<int64_t> res(n);
129  for (uint64_t i = 0; i < n; ++i) {
130    res.set_coeff(i, uniform_i64());
131  }
132  return res;
133}
134
135/** @brief random float64 gaussian polynomial */
136template <>
137polynomial<double> random_polynomial(uint64_t n) {
138  polynomial<double> res(n);
139  for (uint64_t i = 0; i < n; ++i) {
140    res.set_coeff(i, random_f64_gaussian());
141  }
142  return res;
143}
144
145template <typename T>
146polynomial<T> random_polynomial_bounds(uint64_t n, const T lb, const T ub);
147
148/** @brief random int64 polynomial */
149template <>
150polynomial<int64_t> random_polynomial_bounds(uint64_t n, const int64_t lb, const int64_t ub) {
151  polynomial<int64_t> res(n);
152  for (uint64_t i = 0; i < n; ++i) {
153    res.set_coeff(i, uniform_i64_bounds(lb, ub));
154  }
155  return res;
156}
157
158/** @brief random float64 gaussian polynomial */
159template <>
160polynomial<double> random_polynomial_bounds(uint64_t n, const double lb, const double ub) {
161  polynomial<double> res(n);
162  for (uint64_t i = 0; i < n; ++i) {
163    res.set_coeff(i, uniform_f64_bounds(lb, ub));
164  }
165  return res;
166}
167
168/** @brief random int64 polynomial */
169template <>
170polynomial<__int128_t> random_polynomial_bounds(uint64_t n, const __int128_t lb, const __int128_t ub) {
171  polynomial<__int128_t> res(n);
172  for (uint64_t i = 0; i < n; ++i) {
173    res.set_coeff(i, uniform_i128_bounds(lb, ub));
174  }
175  return res;
176}
177
178template <typename T>
179polynomial<T> random_polynomial_bits(uint64_t n, const uint64_t bits) {
180  T b = UINT64_C(1) << bits;
181  return random_polynomial_bounds(n, -b, b);
182}
183
184template <>
185polynomial<int64_t> polynomial<int64_t>::random_log2bound(uint64_t n, uint64_t log2bound) {
186  polynomial<int64_t> res(n);
187  for (uint64_t i = 0; i < n; ++i) {
188    res.set_coeff(i, uniform_i64_bits(log2bound));
189  }
190  return res;
191}
192
193template <>
194polynomial<int64_t> polynomial<int64_t>::random(uint64_t n) {
195  polynomial<int64_t> res(n);
196  for (uint64_t i = 0; i < n; ++i) {
197    res.set_coeff(i, uniform_u64());
198  }
199  return res;
200}
201
202template <>
203polynomial<double> polynomial<double>::random_log2bound(uint64_t n, uint64_t log2bound) {
204  polynomial<double> res(n);
205  double bound = pow(2., log2bound);
206  for (uint64_t i = 0; i < n; ++i) {
207    res.set_coeff(i, uniform_f64_bounds(-bound, bound));
208  }
209  return res;
210}
211
212template <>
213polynomial<double> polynomial<double>::random(uint64_t n) {
214  polynomial<double> res(n);
215  double bound = 2.;
216  for (uint64_t i = 0; i < n; ++i) {
217    res.set_coeff(i, uniform_f64_bounds(-bound, bound));
218  }
219  return res;
220}
221
222template <typename T>
223polynomial<T> naive_product(const polynomial<T>& a, const polynomial<T>& b) {
224  const int64_t nn = a.nn();
225  REQUIRE_DRAMATICALLY(b.nn() == uint64_t(nn), "dimension mismatch!");
226  polynomial<T> res(nn);
227  for (int64_t i = 0; i < nn; ++i) {
228    T ri = 0;
229    for (int64_t j = 0; j < nn; ++j) {
230      ri += a.get_coeff(j) * b.get_coeff(i - j);
231    }
232    res.set_coeff(i, ri);
233  }
234  return res;
235}
236
237#define EXPLICIT_INSTANTIATE_POLYNOMIAL(TYPE)                                                    \
238  template class polynomial<TYPE>;                                                               \
239  template bool operator==(const polynomial<TYPE>& a, const polynomial<TYPE>& b);                \
240  template polynomial<TYPE> operator+(const polynomial<TYPE>& a, const polynomial<TYPE>& b);     \
241  template polynomial<TYPE> operator-(const polynomial<TYPE>& a, const polynomial<TYPE>& b);     \
242  template polynomial<TYPE> operator-(const polynomial<TYPE>& a);                                \
243  template polynomial<TYPE> random_polynomial_bits(uint64_t n, const uint64_t bits);             \
244  template polynomial<TYPE> naive_product(const polynomial<TYPE>& a, const polynomial<TYPE>& b); \
245  // template polynomial<TYPE> random_polynomial(uint64_t n);
246
247#endif  // SPQLIOS_NEGACYCLIC_POLYNOMIAL_IMPL_H