1#include "mod_q120.h"
  2
  3#include <cstdint>
  4#include <random>
  5
  6int64_t centermod(int64_t v, int64_t q) {
  7  int64_t t = v % q;
  8  if (t >= (q + 1) / 2) return t - q;
  9  if (t < -q / 2) return t + q;
 10  return t;
 11}
 12
 13int64_t centermod(uint64_t v, int64_t q) {
 14  int64_t t = int64_t(v % uint64_t(q));
 15  if (t >= q / 2) return t - q;
 16  return t;
 17}
 18
 19mod_q120::mod_q120() {
 20  for (uint64_t i = 0; i < 4; ++i) {
 21    a[i] = 0;
 22  }
 23}
 24
 25mod_q120::mod_q120(int64_t a0, int64_t a1, int64_t a2, int64_t a3) {
 26  a[0] = centermod(a0, Qi[0]);
 27  a[1] = centermod(a1, Qi[1]);
 28  a[2] = centermod(a2, Qi[2]);
 29  a[3] = centermod(a3, Qi[3]);
 30}
 31
 32mod_q120 operator+(const mod_q120& x, const mod_q120& y) {
 33  mod_q120 r;
 34  for (uint64_t i = 0; i < 4; ++i) {
 35    r.a[i] = centermod(x.a[i] + y.a[i], mod_q120::Qi[i]);
 36  }
 37  return r;
 38}
 39
 40mod_q120 operator-(const mod_q120& x, const mod_q120& y) {
 41  mod_q120 r;
 42  for (uint64_t i = 0; i < 4; ++i) {
 43    r.a[i] = centermod(x.a[i] - y.a[i], mod_q120::Qi[i]);
 44  }
 45  return r;
 46}
 47
 48mod_q120 operator*(const mod_q120& x, const mod_q120& y) {
 49  mod_q120 r;
 50  for (uint64_t i = 0; i < 4; ++i) {
 51    r.a[i] = centermod(x.a[i] * y.a[i], mod_q120::Qi[i]);
 52  }
 53  return r;
 54}
 55
 56mod_q120& operator+=(mod_q120& x, const mod_q120& y) {
 57  for (uint64_t i = 0; i < 4; ++i) {
 58    x.a[i] = centermod(x.a[i] + y.a[i], mod_q120::Qi[i]);
 59  }
 60  return x;
 61}
 62
 63mod_q120& operator-=(mod_q120& x, const mod_q120& y) {
 64  for (uint64_t i = 0; i < 4; ++i) {
 65    x.a[i] = centermod(x.a[i] - y.a[i], mod_q120::Qi[i]);
 66  }
 67  return x;
 68}
 69
 70mod_q120& operator*=(mod_q120& x, const mod_q120& y) {
 71  for (uint64_t i = 0; i < 4; ++i) {
 72    x.a[i] = centermod(x.a[i] * y.a[i], mod_q120::Qi[i]);
 73  }
 74  return x;
 75}
 76
 77int64_t modq_pow(int64_t x, int32_t k, int64_t q) {
 78  k = (k % (q - 1) + q - 1) % (q - 1);
 79
 80  int64_t res = 1;
 81  int64_t x_pow = centermod(x, q);
 82  while (k != 0) {
 83    if (k & 1) res = centermod(res * x_pow, q);
 84    x_pow = centermod(x_pow * x_pow, q);
 85    k >>= 1;
 86  }
 87  return res;
 88}
 89
 90mod_q120 pow(const mod_q120& x, int32_t k) {
 91  const int64_t r0 = modq_pow(x.a[0], k, x.Qi[0]);
 92  const int64_t r1 = modq_pow(x.a[1], k, x.Qi[1]);
 93  const int64_t r2 = modq_pow(x.a[2], k, x.Qi[2]);
 94  const int64_t r3 = modq_pow(x.a[3], k, x.Qi[3]);
 95  return mod_q120{r0, r1, r2, r3};
 96}
 97
 98static int64_t half_modq(int64_t x, int64_t q) {
 99  // q must be odd in this function
100  if (x % 2 == 0) return x / 2;
101  return centermod((x + q) / 2, q);
102}
103
104mod_q120 half(const mod_q120& x) {
105  const int64_t r0 = half_modq(x.a[0], x.Qi[0]);
106  const int64_t r1 = half_modq(x.a[1], x.Qi[1]);
107  const int64_t r2 = half_modq(x.a[2], x.Qi[2]);
108  const int64_t r3 = half_modq(x.a[3], x.Qi[3]);
109  return mod_q120{r0, r1, r2, r3};
110}
111
112bool operator==(const mod_q120& x, const mod_q120& y) {
113  for (uint64_t i = 0; i < 4; ++i) {
114    if (x.a[i] != y.a[i]) return false;
115  }
116  return true;
117}
118
119std::ostream& operator<<(std::ostream& out, const mod_q120& x) {
120  return out << "q120{" << x.a[0] << "," << x.a[1] << "," << x.a[2] << "," << x.a[3] << "}";
121}
122
123mod_q120 mod_q120::from_q120a(const void* addr) {
124  static const uint64_t _2p32 = UINT64_C(1) << 32;
125  const uint64_t* in = (const uint64_t*)addr;
126  mod_q120 r;
127  for (uint64_t i = 0; i < 4; ++i) {
128    REQUIRE_DRAMATICALLY(in[i] < _2p32, "invalid layout a q120");
129    r.a[i] = centermod(in[i], mod_q120::Qi[i]);
130  }
131  return r;
132}
133
134mod_q120 mod_q120::from_q120b(const void* addr) {
135  const uint64_t* in = (const uint64_t*)addr;
136  mod_q120 r;
137  for (uint64_t i = 0; i < 4; ++i) {
138    r.a[i] = centermod(in[i], mod_q120::Qi[i]);
139  }
140  return r;
141}
142
143mod_q120 mod_q120::from_q120c(const void* addr) {
144  // static const uint64_t _mask_2p32 = (uint64_t(1) << 32) - 1;
145  const uint32_t* in = (const uint32_t*)addr;
146  mod_q120 r;
147  for (uint64_t i = 0, k = 0; i < 8; i += 2, ++k) {
148    const uint64_t q = mod_q120::Qi[k];
149    uint64_t u = in[i];
150    uint64_t w = in[i + 1];
151    REQUIRE_DRAMATICALLY(((u << 32) % q) == (w % q),
152                         "invalid layout q120c: " << u << ".2^32 != " << (w >> 32) << " mod " << q);
153    r.a[k] = centermod(u, q);
154  }
155  return r;
156}
157__int128_t mod_q120::to_int128() const {
158  static const __int128_t qm[] = {(__int128_t(Qi[1]) * Qi[2]) * Qi[3], (__int128_t(Qi[0]) * Qi[2]) * Qi[3],
159                                  (__int128_t(Qi[0]) * Qi[1]) * Qi[3], (__int128_t(Qi[0]) * Qi[1]) * Qi[2]};
160  static const int64_t CRTi[] = {Q1_CRT_CST, Q2_CRT_CST, Q3_CRT_CST, Q4_CRT_CST};
161  static const __int128_t q = qm[0] * Qi[0];
162  static const __int128_t qs2 = q / 2;
163  __int128_t res = 0;
164  for (uint64_t i = 0; i < 4; ++i) {
165    res += (a[i] * CRTi[i] % Qi[i]) * qm[i];
166  }
167  res = (((res % q) + q + qs2) % q) - qs2;  // centermod
168  return res;
169}
170void mod_q120::save_as_q120a(void* dest) const {
171  int64_t* d = (int64_t*)dest;
172  for (uint64_t i = 0; i < 4; ++i) {
173    d[i] = a[i] + Qi[i];
174  }
175}
176void mod_q120::save_as_q120b(void* dest) const {
177  int64_t* d = (int64_t*)dest;
178  for (uint64_t i = 0; i < 4; ++i) {
179    d[i] = a[i] + (Qi[i] * (1 + uniform_u64_bits(32)));
180  }
181}
182void mod_q120::save_as_q120c(void* dest) const {
183  int32_t* d = (int32_t*)dest;
184  for (uint64_t i = 0; i < 4; ++i) {
185    d[2 * i] = a[i] + 3 * Qi[i];
186    d[2 * i + 1] = (uint64_t(d[2 * i]) << 32) % uint64_t(Qi[i]);
187  }
188}
189
190mod_q120 uniform_q120() {
191  test_rng& gen = randgen();
192  std::uniform_int_distribution<uint64_t> dista(0, mod_q120::Qi[0]);
193  std::uniform_int_distribution<uint64_t> distb(0, mod_q120::Qi[1]);
194  std::uniform_int_distribution<uint64_t> distc(0, mod_q120::Qi[2]);
195  std::uniform_int_distribution<uint64_t> distd(0, mod_q120::Qi[3]);
196  return mod_q120(dista(gen), distb(gen), distc(gen), distd(gen));
197}
198
199void uniform_q120a(void* dest) {
200  uint64_t* res = (uint64_t*)dest;
201  for (uint64_t i = 0; i < 4; ++i) {
202    res[i] = uniform_u64_bits(32);
203  }
204}
205
206void uniform_q120b(void* dest) {
207  uint64_t* res = (uint64_t*)dest;
208  for (uint64_t i = 0; i < 4; ++i) {
209    res[i] = uniform_u64();
210  }
211}
212
213void uniform_q120c(void* dest) {
214  uint32_t* res = (uint32_t*)dest;
215  static const uint64_t _2p32 = uint64_t(1) << 32;
216  for (uint64_t i = 0, k = 0; i < 8; i += 2, ++k) {
217    const uint64_t q = mod_q120::Qi[k];
218    const uint64_t z = uniform_u64_bits(32);
219    const uint64_t z_pow_red = (z << 32) % q;
220    const uint64_t room = (_2p32 - z_pow_red) / q;
221    const uint64_t z_pow = z_pow_red + (uniform_u64() % room) * q;
222    REQUIRE_DRAMATICALLY(z < _2p32, "bug!");
223    REQUIRE_DRAMATICALLY(z_pow < _2p32, "bug!");
224    REQUIRE_DRAMATICALLY(z_pow % q == (z << 32) % q, "bug!");
225
226    res[i] = (uint32_t)z;
227    res[i + 1] = (uint32_t)z_pow;
228  }
229}