4023 <T>(a.get_n(), k, GRS_G(a, d, k)), a(a), d(d) {
4024 const size_t n =
this->n;
4026 for (size_t i = 0; i < n; ++i) {
4027 if (d[i] == T(0))
throw std::invalid_argument(
"GRS codes must have nonzero column multipliers!");
4030 this->set_dmin(n - k + 1);
4032 catch (
const std::invalid_argument& e) {
4033 throw std::invalid_argument(std::string(
"Cannot construct GRS code: ") + e.what());
4053 if (
this->dmin.has_value()) res.set_dmin(*(
this->dmin));
4054 if constexpr (FiniteFieldType<T>) {
4055 if (
this->weight_enumerator.has_value()) res.set_weight_enumerator(*(
this->weight_enumerator));
4056 if (
this->p_ary_image_weight_enumerator.has_value())
4057 res.p_ary_image_weight_enumerator =
this->p_ary_image_weight_enumerator;
4059 if (
this->minimal_trellis.has_value()) res.minimal_trellis =
this->minimal_trellis;
4060 if (
this->codewords.has_value()) res.codewords =
this->codewords;
4065 if (
this->k ==
this->n)
4066 throw std::logic_error(
4067 "The dual of a full-dimension GRS code is the zero code, construct a ZeroCode instead!");
4068 return GRSCode<T>(a, dual_multipliers(),
this->n -
this->k);
4072 if constexpr (!FiniteFieldType<T>) {
4073 throw std::logic_error(
"Cannot calculate weight enumerator of code over infinite field!");
4075 this->weight_enumerator.call_once([
this] {
4076 if (
this->weight_enumerator.has_value())
return;
4077 constexpr size_t q = T::get_size();
4078 const size_t n =
this->n;
4079 const size_t k =
this->k;
4081 Polynomial<InfInt> enumerator;
4082 enumerator.set_coefficient(0, 1);
4083 for (size_t i = n - k + 1; i <= n; ++i) {
4085 for (size_t j = 0; j <= i - (n - k + 1); ++j) {
4086 InfInt s = j % 2 ? -1 : 1;
4087 sum += s * bin<
InfInt>(i - 1, j) * sqm<
InfInt>(q, i - j - (n - k + 1));
4089 enumerator.set_coefficient(i, bin<
InfInt>(n, i) * (q - 1) * sum);
4091 this->weight_enumerator.emplace(
std::move(enumerator));
4093 return this->weight_enumerator.value();
4098 for (size_t i = 0; i <
this->n; ++i) {
4099 if (a[i] == T(0))
return true;
4105 if constexpr (!FiniteFieldType<T>) {
4108 const size_t n =
this->n;
4109 if (n != T::get_size() - 1)
return false;
4110 for (size_t i = 0; i < n; ++i) {
4111 if (a[i] == T(0))
return false;
4118 for (size_t i = 0; i <
this->n; ++i) {
4119 if (d[i] != T(1))
return false;
4127 if (
this->
template info_prelude<
LinearCode<T>>(os)) {
4128 os <<
BOLD(
"GRS code")
" with properties: { a = " << a <<
", d = " << d;
4130 os <<
" singly-extended";
4136 os <<
" narrow-sense";
4139 os <<
" normalized";
4146#ifdef CECCO_ERASURE_SUPPORT
4147 if (LinearCode<T>::erasures_present(r))
return dec_WBA_EE(r);
4149 return dec_WBA_impl(r, {});
4153#ifdef CECCO_ERASURE_SUPPORT
4154 if (LinearCode<T>::erasures_present(r))
return dec_WBA_impl(r,
this->erasure_positions(r),
true);
4156 return dec_WBA_impl(r, {},
true);
4160#ifdef CECCO_ERASURE_SUPPORT
4161 if (LinearCode<T>::erasures_present(r))
return dec_BMA_EE(r);
4163 return dec_BMA_impl(r, {});
4169#ifdef CECCO_ERASURE_SUPPORT
4176 return dec_GSA_list(r, 1, list_size);
4180 return dec_GSA_list_optimized(r, 1, list_size);
4194#ifdef CECCO_ERASURE_SUPPORT
4195 return dec_GSA_list_impl(r,
this->erasure_positions(r), multiplicity, list_size,
false);
4197 return dec_GSA_list_impl(r, {}, multiplicity, list_size,
false);
4203#ifdef CECCO_ERASURE_SUPPORT
4204 return dec_GSA_list_impl(r,
this->erasure_positions(r), multiplicity, list_size,
true);
4206 return dec_GSA_list_impl(r, {}, multiplicity, list_size,
true);
4212 if constexpr (!FiniteFieldType<T>) {
4213 throw std::logic_error(
"Koetter-Vardy decoding only available for codes over finite fields!");
4214 }
else if constexpr (T::get_size() != 2) {
4215 throw std::logic_error(
4216 "Soft-input KVA LLR vector decoding only available for binary codes; use the Matrix<double> overload!");
4218 return dec_KVA_list(Matrix<
double>(llrs), strength, list_size);
4230 if constexpr (!FiniteFieldType<T>) {
4231 throw std::logic_error(
"Koetter-Vardy decoding only available for codes over finite fields!");
4233 this->validate_length(llrs);
4234 if (strength == 0)
throw std::invalid_argument(
"KVA strength must be positive!");
4235 if (strength > 1000000)
throw std::invalid_argument(
"KVA strength must be at most 10^6!");
4236 if (list_size == 0)
throw std::invalid_argument(
"KVA list size must be positive!");
4238 constexpr size_t q = T::get_size();
4239 const size_t n =
this->n;
4243 for (size_t beta = 1; beta < q; ++beta)
4244 for (size_t i = 0; i < n; ++i)
4245 if (
std::isnan(llrs(beta - 1, i)))
throw std::invalid_argument(
"KVA LLRs must not be NaN!");
4247 auto Theta = ZeroMatrix<
int>(q, n);
4248 bool nonzero =
false;
4252 std::vector<
double> weights(q - 1);
4253 for (size_t i = 0; i < n; ++i) {
4254 double denominator = 1.0;
4255 for (size_t beta = 1; beta < q; ++beta) {
4256 weights[beta - 1] = std::exp(-std::clamp(llrs(beta - 1, i), -300.0, 300.0));
4257 denominator += weights[beta - 1];
4259 for (size_t beta = 0; beta < q; ++beta) {
4260 const double Pi = (beta == 0 ? 1.0 : weights[beta - 1]) / denominator;
4261 const int multiplicity =
static_cast<
int>(strength * Pi);
4262 if (multiplicity == 0)
continue;
4263 Theta.set_component((T(beta) / d[i]).get_label(), i, multiplicity);
4267 if (!nonzero)
return {};
4269 return GSA_rank_candidates(dec_KVA_impl(a, Theta, list_size), [&llrs](
const Vector<T>& c) {
4272 double finite_sum = 0.0;
4273 bool pos_inf =
false;
4274 bool neg_inf =
false;
4275 for (size_t i = 0; i < c.get_n(); ++i) {
4276 const size_t label = c[i].get_label();
4277 if (label == 0)
continue;
4278 const double l = llrs(label - 1, i);
4280 (l > 0.0 ? pos_inf : neg_inf) =
true;
4284 if (pos_inf)
return std::numeric_limits<
double>::infinity();
4285 if (neg_inf)
return -std::numeric_limits<
double>::infinity();
4296 const auto& Gp = C.get_G();
4297 const size_t n = C.get_n();
4298 const size_t k = C.get_k();
4299 if (k == 0)
throw std::invalid_argument(
"Linear code of dimension 0 cannot be converted into GRS code!");
4302 for (size_t i = 0; i < n; ++i) {
4303 if (Gp(0, i).is_zero())
throw std::invalid_argument(
"Linear code cannot be converted into GRS code!");
4304 d.set_component(i, Gp(0, i));
4310 if constexpr (FiniteFieldType<T>) {
4311 if (n > T::get_size())
throw std::invalid_argument(
"Linear code cannot be converted into GRS code!");
4313 for (size_t i = 0; i < n; ++i) a.set_component(i, T(i));
4314 return {
std::move(a),
std::move(d)};
4317 for (size_t i = 0; i < n; ++i) a.set_component(i, Gp(1, i) / Gp(0, i));
4318 for (size_t j = 2; j < k; ++j)
4319 for (size_t i = 0; i < n; ++i)
4320 if (Gp(j, i) != Gp(j - 1, i) * a[i])
4321 throw std::invalid_argument(
"Linear code cannot be converted into GRS code!");
4322 if (!a.is_pairwise_distinct())
throw std::invalid_argument(
"Linear code cannot be converted into GRS code!");
4324 return {
std::move(a),
std::move(d)};
4337 Vector<T> dec_WBA_impl(
const Vector<T>& r,
const std::vector<size_t>& X,
bool reencode =
false)
const {
4338 this->validate_length(r);
4340 const size_t tau = X.size();
4341 if (tau >
this->get_dmin() - 1)
4342 throw decoding_failure(
"GRS code WBA error/erasure decoder failed (too many erasures)!");
4344 const size_t k =
this->k;
4346 const auto ap = delete_components(a, X);
4347 const auto dp = delete_components(d, X);
4348 const size_t np =
this->n - tau;
4349 const size_t tmaxp = (np - k) / 2;
4350 const auto r_punct = delete_components(r, X);
4352 G_canonical.call_once([
this] {
4353 if (G_canonical.has_value())
return;
4354 G_canonical.emplace(VandermondeMatrix<T>(a,
this->k) * DiagonalMatrix<T>(d));
4358 auto Q0 = ZeroPolynomial<T>();
4359 auto Q1 = ZeroPolynomial<T>();
4362 c_reenc = reencoding_projection(ap, dp, r_punct, rp);
4364 if (dH(r, c_reenc) > tmaxp)
4365 throw decoding_failure(
4366 "GRS code WBA error/erasure decoder failed (identified error beyond BD radius)!");
4372 reencoding_prefactors(pf, ap, 1, 1);
4374 auto U0 = ZeroPolynomial<T>();
4375 if (np ==
this->n && tmaxp > 0) {
4379 R_canonical.call_once([
this, &pf] {
4380 if (R_canonical.has_value())
return;
4381 const size_t k =
this->k;
4382 const size_t n =
this->n;
4383 auto S = ZeroMatrix<T>(n - k, n - k);
4384 for (size_t i = k; i < n; ++i) {
4385 const T scale = pf.F[0](a[i]);
4387 for (size_t j = 0; j < n - k; ++j) {
4388 S.set_component(i - k, j, scale * pow);
4393 R_canonical.emplace(
std::move(S));
4396 auto M1 = ZeroMatrix<T>(np - k, np - tmaxp - k + 1);
4397 for (size_t i = k; i < np; ++i) {
4398 const T y = rp[i] / dp[i];
4400 for (size_t j = 0; j <= np - tmaxp - k; ++j) {
4401 M1.set_component(i - k, j, y * pow);
4405 const auto& R = R_canonical.value();
4406 const auto B = (R.get_submatrix(np - tmaxp - k, 0, tmaxp, np - k) * M1).basis_of_kernel();
4407 for (size_t j = 0; j <= np - tmaxp - k; ++j) Q1.set_coefficient(j, B(0, j));
4411 for (size_t i = k; i < np; ++i) w.set_component(i - k, M1(i - k, 0) * Q1(ap[i]));
4412 for (size_t i = 0; i < np - tmaxp - k; ++i) {
4414 for (size_t j = 0; j < np - k; ++j) s += R(i, j) * w[j];
4415 U0.set_coefficient(i, T(0) - s);
4418 Vector<T> ap_act(np - k), y_act(np - k), F0_act(np - k);
4419 for (size_t i = k; i < np; ++i) {
4420 ap_act.set_component(i - k, ap[i]);
4421 y_act.set_component(i - k, rp[i] / dp[i]);
4422 F0_act.set_component(i - k, pf.F[0](ap[i]));
4425 const auto M0 = DiagonalMatrix(F0_act) * transpose(VandermondeMatrix<T>(ap_act, np - tmaxp - k));
4426 const auto M1 = DiagonalMatrix(y_act) * transpose(VandermondeMatrix<T>(ap_act, np - tmaxp - k + 1));
4427 const auto M = horizontal_join(M0, M1);
4429 const auto B = M.basis_of_kernel();
4432 for (size_t j = 0; j < np - tmaxp - k; ++j, ++i) U0.set_coefficient(j, B(0, i));
4433 for (size_t j = 0; j <= np - tmaxp - k; ++j, ++i) Q1.set_coefficient(j, B(0, i));
4437 c_reenc =
Vector<T>(
this->n);
4439 for (size_t i = 0; i < np; ++i) rp_norm.set_component(i, r_punct[i] / dp[i]);
4441 const auto M0 = transpose(VandermondeMatrix<T>(ap, np - tmaxp));
4442 const auto M1 = DiagonalMatrix(rp_norm) * transpose(VandermondeMatrix<T>(ap, np - tmaxp - k + 1));
4443 const auto M = horizontal_join(M0, M1);
4445 const auto B = M.basis_of_kernel();
4448 for (size_t j = 0; j <= np - tmaxp - 1; ++j, ++i) Q0.set_coefficient(j, B(0, i));
4449 for (size_t j = 0; j <= np - tmaxp - k; ++j, ++i) Q1.set_coefficient(j, B(0, i));
4452 const auto temp = poly_long_div(-Q0, Q1);
4453 const auto quotient = temp.first;
4454 const auto remainder = temp.second;
4455 if (!remainder.is_zero() || (!quotient.is_zero() && quotient.degree() >= k))
4456 throw decoding_failure(
"GRS code WBA error/erasure decoder failed (true error beyond BD radius)!");
4458 const auto u_est = pad_back(
Vector<T>(quotient), k);
4459 const auto c_est = u_est * G_canonical.value() + c_reenc;
4460 if (dH(r, c_est) > tmaxp)
4461 throw decoding_failure(
"GRS code WBA error/erasure decoder failed (identified error beyond BD radius)!");
4466 Vector<T> dec_BMA_impl(
const Vector<T>& r,
const std::vector<size_t>& X)
const {
4467 this->validate_length(r);
4469 const size_t n =
this->n;
4470 const size_t redundancy = n -
this->k;
4471 const size_t tau = X.size();
4472 if (tau > redundancy)
throw decoding_failure(
"GRS code BMA error/erasure decoder failed (too many erasures)!");
4473 if (redundancy == 0)
return r;
4476 for (size_t i = 0; i < tau; ++i) rp.set_component(X[i], T(0));
4478 HT_canonical.call_once([
this] {
4479 if (HT_canonical.has_value())
return;
4480 const auto D = DiagonalMatrix<T>(dual_multipliers());
4481 const auto V = VandermondeMatrix<T>(a,
this->n -
this->k);
4482 HT_canonical.emplace(transpose(V * D));
4484 const auto& HT = HT_canonical.value();
4486 const auto s = rp * HT;
4487 if (s.is_zero() && tau == 0)
return r;
4490 auto Psi = OnePolynomial<T>();
4491 for (size_t i = 0; i < tau; ++i) Psi *=
Polynomial<T>({-a[X[i]], T(1)});
4494 auto sxp = ZeroPolynomial<T>();
4495 for (size_t i = 0; i < redundancy - tau; ++i) {
4497 for (size_t j = 0; j <= tau; ++j) coeff += Psi[j] * s[i + j];
4498 sxp.set_coefficient(i, coeff);
4502 auto Lambda_e = OnePolynomial<T>();
4504 auto B_ref = OnePolynomial<T>();
4506 T discrepancy_ref(1);
4507 size_t last_update = 0;
4509 for (size_t i = 0; i < redundancy - tau; ++i) {
4511 const size_t j0 = (t > i ? t - i : 0);
4512 for (size_t j = j0; j <= t; ++j) discrepancy += Lambda_e[j] * sxp[i + j - t];
4514 if (!discrepancy.is_zero()) {
4515 const auto ratio = discrepancy / discrepancy_ref;
4516 const size_t shift = i - last_update + 1;
4518 if (shift + t_ref <= t) {
4519 Lambda_e -= shift_up(ratio * B_ref, t - shift - t_ref);
4522 const auto temp = t;
4523 Lambda_e.shift_up(shift + t_ref - t);
4524 Lambda_e -= ratio * B_ref;
4526 B_ref =
std::move(B);
4528 discrepancy_ref = discrepancy;
4529 last_update = i + 1;
4534 if (2 * t + tau > redundancy)
4535 throw decoding_failure(
"GRS code BMA error/erasure decoder failed (identified error beyond BD radius)!");
4539 std::vector<size_t> E(X);
4541 std::vector<
bool> erased(n,
false);
4542 for (size_t i = 0; i < tau; ++i) erased[X[i]] =
true;
4543 for (size_t i = 0; i < n && E.size() < tau + t; ++i)
4544 if (!erased[i] && Lambda_e(a[i]).is_zero()) E.push_back(i);
4546 if (E.size() != tau + t)
4547 throw decoding_failure(
"GRS code BMA error/erasure decoder failed (true error beyond BD radius)!");
4549 auto Lambda = Psi * Lambda_e;
4552 auto Omega = ZeroPolynomial<T>();
4553 for (size_t j = 0; j < tau + t; ++j) {
4555 for (size_t m = j + 1; m <= tau + t; ++m) coeff += Lambda[m] * sx[m - j - 1];
4556 Omega.set_coefficient(j, coeff);
4559 Lambda.differentiate(1);
4562 for (size_t i = 0; i < E.size(); ++i) {
4563 const auto locator = a[E[i]];
4564 c_est.set_component(E[i], c_est[E[i]] - Omega(locator) / (HT(E[i], 0) * Lambda(locator)));
4567 if (2 * dH(r, c_est) + tau > redundancy)
4568 throw decoding_failure(
"GRS code BMA error/erasure decoder failed (result outside the decoding radius)!");
4575 static size_t GSA_block_size(size_t k, size_t h, size_t mu) {
return h > mu * (k - 1) ? h - mu * (k - 1) : 0; }
4579 static size_t GSA_degree_bound(size_t k, size_t ell, size_t cost) {
4581 size_t n_unknowns = 0;
4585 for (size_t mu = 0; mu <= ell; ++mu) n_unknowns += GSA_block_size(k, h, mu);
4586 }
while (n_unknowns <= cost);
4590 std::vector<Vector<T>> dec_GSA_list_impl(
const Vector<T>& r,
const std::vector<size_t>& X, size_t multiplicity,
4591 size_t list_size,
bool optimized =
true)
const {
4592 if constexpr (!FiniteFieldType<T>) {
4593 throw std::logic_error(
"Guruswami-Sudan decoding only available for codes over finite fields!");
4595 this->validate_length(r);
4596 if (multiplicity == 0)
throw std::invalid_argument(
"GSA multiplicity must be positive!");
4597 if (list_size == 0)
throw std::invalid_argument(
"GSA list size must be positive!");
4599 const size_t tau = X.size();
4600 if (tau >
this->get_dmin() - 1)
return {};
4602 const auto ap = delete_components(a, X);
4603 const auto dp = delete_components(d, X);
4604 const auto r_punct = delete_components(r, X);
4606 return GSA_rank_candidates(dec_GSA_impl(ap, dp, r_punct, multiplicity, list_size, optimized, optimized),
4607 [&r_punct, &X](
const Vector<T>& c) {
4608 return static_cast<
double>(dH(r_punct, delete_components(c, X)));
4617 size_t reencoded = 0;
4618 std::vector<Polynomial<T>> F;
4619 std::vector<size_t> deg_F;
4623 const size_t k =
this->k;
4624 const size_t np = ap.get_n();
4626 for (size_t i = 0; i < k; ++i) rJ.set_component(i, r[i]);
4628 if (np ==
this->n) {
4629 if (np == k)
return r;
4630 W_canonical.call_once([
this] {
4631 if (W_canonical.has_value())
return;
4632 W_canonical.emplace(
this->MI *
this->G.get_submatrix(0,
this->k,
this->k,
this->n -
this->k));
4634 const auto p = rJ * W_canonical.value();
4635 for (size_t i = k; i < np; ++i) rp.set_component(i, r[i] - p[i - k]);
4638 const auto Gp = VandermondeMatrix<T>(ap, k) * DiagonalMatrix<T>(dp);
4639 const auto u_reenc = rJ * Gp.get_submatrix(0, 0, k, k).invert();
4641 const auto p = u_reenc * Gp.get_submatrix(0, k, k, np - k);
4642 for (size_t i = k; i < np; ++i) rp.set_component(i, r[i] - p[i - k]);
4644 return u_reenc * canonical_G();
4647 static Polynomial<T> linear_factor(
const T& root) {
4648 auto res = Monomial<T>(1, T(1));
4649 res.set_coefficient(0, T(0) - root);
4653 Vector<T> dual_multipliers()
const {
4654 const size_t n =
this->n;
4656 for (size_t i = 0; i < n; ++i) {
4658 for (size_t j = 0; j < n; ++j)
4659 if (j != i) prod *= a[i] - a[j];
4660 dperp.set_component(i, T(1) / prod);
4665 const Matrix<T>& canonical_G()
const {
4666 G_canonical.call_once([
this] {
4667 if (G_canonical.has_value())
return;
4668 G_canonical.emplace(VandermondeMatrix<T>(a,
this->k) * DiagonalMatrix<T>(d));
4670 return G_canonical.value();
4673 void reencoding_prefactors(Prefactors& pf,
const Vector<T>& ap, size_t multiplicity, size_t list_size)
const {
4674 const size_t k =
this->k;
4675 const size_t r = multiplicity;
4676 const size_t nf = std::min(r, list_size + 1);
4678 pf.F.assign(nf, OnePolynomial<T>());
4679 pf.deg_F.assign(nf, 0);
4683 if (ap.get_n() ==
this->n) {
4684 P_canonical.call_once([
this] {
4685 if (P_canonical.has_value())
return;
4686 auto poly = OnePolynomial<T>();
4687 for (size_t i = 0; i <
this->k; ++i) poly *= linear_factor(a[i]);
4688 P_canonical.emplace(
std::move(poly));
4690 P = P_canonical.value();
4692 P = OnePolynomial<T>();
4693 for (size_t i = 0; i < k; ++i) P *= linear_factor(ap[i]);
4695 std::vector<Polynomial<T>> P_pow(r + 1, OnePolynomial<T>());
4696 for (size_t e = 1; e <= r; ++e) P_pow[e] = P_pow[e - 1] * P;
4698 for (size_t mu = 0; mu < nf; ++mu) {
4699 pf.F[mu] *= P_pow[r - mu];
4700 pf.deg_F[mu] += k * (r - mu);
4705 void Sierpinski_prefactors(Prefactors& pf,
const Vector<T>& ap, size_t multiplicity, size_t list_size,
4707 const size_t np = ap.get_n();
4708 const size_t r = multiplicity;
4709 const size_t ell = list_size;
4711 const size_t nf = std::min(r, ell + 1);
4715 std::vector<size_t> g(r, r);
4716 std::vector<
bool> resolvable(r,
false);
4717 for (size_t t = r; t-- > 0;) {
4719 for (size_t v = t + 1; v <= ell && ok; ++v)
4720 if (triangle(v, t) != 0) ok = v < r && resolvable[v];
4723 for (size_t t = 0; t < r; ++t) {
4724 if (!resolvable[t])
continue;
4726 for (size_t v = t + 1; v <= ell; ++v)
4727 if (triangle(v, t) != 0) g[t] = v;
4731 for (size_t mu = 0; mu < nf; ++mu)
4732 if (g[mu] < r) e_max =
std::max(e_max, r - g[mu]);
4733 if (e_max == 0)
return;
4736 pf.F.assign(nf, OnePolynomial<T>());
4737 pf.deg_F.assign(nf, 0);
4741 if (np ==
this->n && pf.reencoded ==
this->k) {
4742 PN_canonical.call_once([
this] {
4743 if (PN_canonical.has_value())
return;
4744 auto poly = OnePolynomial<T>();
4745 for (size_t i =
this->k; i <
this->n; ++i) poly *= linear_factor(a[i]);
4746 PN_canonical.emplace(
std::move(poly));
4748 PN = PN_canonical.value();
4750 PN = OnePolynomial<T>();
4751 for (size_t i = pf.reencoded; i < np; ++i) PN *= linear_factor(ap[i]);
4753 std::vector<Polynomial<T>> PN_pow(e_max + 1, OnePolynomial<T>());
4754 for (size_t e = 1; e <= e_max; ++e) PN_pow[e] = PN_pow[e - 1] * PN;
4756 for (size_t mu = 0; mu < nf; ++mu) {
4757 const size_t e = g[mu] < r ? r - g[mu] : 0;
4758 pf.F[mu] *= PN_pow[e];
4759 pf.deg_F[mu] += (np - pf.reencoded) * e;
4767 std::vector<Vector<T>> dec_GSA_impl(
const Vector<T>& ap,
const Vector<T>& dp,
const Vector<T>& r,
4768 size_t multiplicity, size_t list_size,
bool reencode,
bool Sierpinski)
const {
4769 const size_t k =
this->k;
4770 const size_t np = ap.get_n();
4771 const size_t ell = list_size;
4772 constexpr size_t p = T::get_characteristic();
4773 if (k == 0)
return {
Vector<T>(
this->n)};
4778 c_reenc = reencoding_projection(ap, dp, r, rp);
4779 if (np == k)
return {c_reenc};
4782 c_reenc =
Vector<T>(
this->n);
4785 const size_t cost = np * bin<size_t>(multiplicity + 1, 2);
4786 const size_t h = GSA_degree_bound(k, ell, cost);
4792 if (reencode) reencoding_prefactors(pf, ap, multiplicity, ell);
4793 if (Sierpinski) Sierpinski_prefactors(pf, ap, multiplicity, ell, triangle);
4794 const auto has_prefactor = [&pf](size_t mu) {
return mu < pf.F.size() && pf.deg_F[mu] != 0; };
4796 std::vector<size_t> width(ell + 1);
4797 std::vector<size_t> offset(ell + 1);
4798 size_t n_unknowns = 0;
4799 for (size_t mu = 0; mu <= ell; ++mu) {
4800 const size_t bs = GSA_block_size(k, h, mu);
4801 const size_t dF = mu < pf.F.size() ? pf.deg_F[mu] : 0;
4802 width[mu] = bs > dF ? bs - dF : 0;
4803 offset[mu] = n_unknowns;
4804 n_unknowns += width[mu];
4809 std::vector<Matrix<T>> FD(pf.F.size());
4810 for (size_t mu = 0; mu < pf.F.size(); ++mu) {
4811 if (width[mu] == 0 || !has_prefactor(mu))
continue;
4812 FD[mu] = ZeroMatrix<T>(multiplicity, np - pf.reencoded);
4813 for (size_t j = 0; j < multiplicity; ++j) {
4814 const auto Fd = Hasse_derivative(pf.F[mu], j);
4815 for (size_t i = pf.reencoded; i < np; ++i) FD[mu].set_component(j, i - pf.reencoded, Fd(ap[i]));
4821 auto M = ZeroMatrix<T>((np - pf.reencoded) * bin<size_t>(multiplicity + 1, 2), n_unknowns);
4822 std::vector<T> xpows(multiplicity, T(1));
4824 for (size_t i = pf.reencoded; i < np; ++i) {
4825 const T y = rp[i] / dp[i];
4826 for (size_t s = 0; s < multiplicity; ++s) {
4827 for (size_t t = 0; s + t < multiplicity; ++t) {
4829 for (size_t mu = t; mu <= ell; ++mu) {
4830 const size_t ct = triangle(mu, t);
4831 if (ct != 0 && has_prefactor(mu)) {
4832 for (size_t j = 0; j <= s; ++j) xpows[j] = T(1);
4833 for (size_t nu = 0; nu < width[mu]; ++nu) {
4835 for (size_t j = 0; j <=
std::min(s, nu); ++j) {
4836 const size_t binomials = (ct * triangle(nu, j)) % p;
4838 entry += binomials * (xpows[j] * FD[mu](s - j, i - pf.reencoded));
4840 if (!entry.is_zero()) M.set_component(row, offset[mu] + nu, entry * ypow);
4841 for (size_t j = s; j > 0; --j) xpows[j] = xpows[j - 1];
4844 }
else if (ct != 0) {
4846 for (size_t nu = s; nu < width[mu]; ++nu) {
4847 const size_t binomials = (ct * triangle(nu, s)) % p;
4848 if (binomials != 0) M.set_component(row, offset[mu] + nu, binomials * (xpow * ypow));
4859 const auto B = M.basis_of_kernel();
4860 if (B.get_m() == 0)
return {};
4861 std::vector<Polynomial<T>> Q(ell + 1, ZeroPolynomial<T>());
4862 for (size_t mu = 0; mu <= ell; ++mu) {
4863 if (width[mu] == 0)
continue;
4864 auto G = ZeroPolynomial<T>();
4865 for (size_t nu = 0; nu < width[mu]; ++nu) G.set_coefficient(nu, B(0, offset[mu] + nu));
4866 Q[mu] = has_prefactor(mu) ? G * pf.F[mu] : G;
4870 std::vector<Polynomial<T>> roots;
4871 auto w = ZeroPolynomial<T>();
4872 find_y_roots(Q, k, 0, w, roots, triangle);
4874 std::vector<Vector<T>> res;
4875 for (size_t j = 0; j < roots.size(); ++j)
4876 res.push_back(pad_back(
Vector<T>(roots[j]), k) * canonical_G() + c_reenc);
4883 static void find_y_roots(
const std::vector<Polynomial<T>>& Q, size_t k, size_t lambda,
Polynomial<T>& g,
4885 constexpr size_t q = T::get_size();
4888 size_t m = std::numeric_limits<size_t>::max();
4889 for (size_t mu = 0; mu < Q.size(); ++mu)
4890 if (!Q[mu].is_zero()) m =
std::min(m, Q[mu].trailing_degree());
4891 if (m == std::numeric_limits<size_t>::max())
return;
4894 std::vector<Polynomial<T>> shifted;
4897 for (size_t mu = 0; mu < shifted.size(); ++mu) shifted[mu].shift_down(m);
4899 const auto& R = m > 0 ? shifted : Q;
4903 for (size_t mu = 0; mu < R.size(); ++mu) f.set_coefficient(mu, R[mu][0]);
4905 std::vector<T> gamma_pow(R.size(), T(1));
4906 for (size_t label = 0; label < q; ++label) {
4907 const T gamma(label);
4908 if (!f(gamma).is_zero())
continue;
4909 for (size_t e = 1; e < gamma_pow.size(); ++e) gamma_pow[e] = gamma_pow[e - 1] * gamma;
4910 g.set_coefficient(lambda, gamma);
4911 if (lambda == k - 1) {
4912 auto acc = ZeroPolynomial<T>();
4913 for (size_t mu = 0; mu < R.size(); ++mu) {
4914 if (gamma_pow[mu].is_zero())
continue;
4915 acc += gamma_pow[mu] * R[mu];
4917 if (acc.is_zero()) U.push_back(g);
4920 std::vector<Polynomial<T>> S(R.size());
4921 for (size_t t = 0; t < R.size(); ++t) {
4922 auto acc = ZeroPolynomial<T>();
4923 for (size_t mu = t; mu < R.size(); ++mu) {
4924 const T scale = triangle(mu, t) * gamma_pow[mu - t];
4925 if (scale.is_zero())
continue;
4926 acc += scale * R[mu];
4929 S[t] =
std::move(acc);
4931 find_y_roots(S, k, lambda + 1, g, U, triangle);
4939 static std::vector<Vector<T>> GSA_rank_candidates(std::vector<Vector<T>> candidates,
const M& metric) {
4940 if (candidates.size() < 2)
return candidates;
4942 std::vector<std::pair<
double, Vector<T>>> ranked;
4943 ranked.reserve(candidates.size());
4944 for (size_t j = 0; j < candidates.size(); ++j) {
4945 ranked.emplace_back(metric(candidates[j]),
std::move(candidates[j]));
4950 for (size_t j = 1; j < ranked.size(); ++j) {
4951 if (ranked[j].first < ranked[choice].first) {
4954 }
else if (ranked[j].first == ranked[choice].first) {
4955 if (details::reservoir_accept(++ties)) choice = j;
4958 std::swap(ranked[0], ranked[choice]);
4959 std::stable_sort(ranked.begin() + 1, ranked.end(),
4960 [](
const auto& lhs,
const auto& rhs) {
return lhs.first < rhs.first; });
4963 for (size_t j = 0; j < ranked.size(); ++j) candidates.push_back(
std::move(ranked[j].second));
4967 std::vector<Vector<T>> dec_KVA_impl(
const Vector<T>& ap,
const Matrix<
int>& Theta, size_t list_size)
const {
4968 const size_t k =
this->k;
4969 const size_t np = ap.get_n();
4970 const size_t q = Theta.get_m();
4971 const size_t ell = list_size;
4972 constexpr size_t p = T::get_characteristic();
4973 if (k == 0)
return {
Vector<T>(
this->n)};
4977 size_t max_multiplicity = 0;
4978 for (size_t beta = 0; beta < q; ++beta)
4979 for (size_t i = 0; i < np; ++i) {
4980 const size_t multiplicity =
static_cast<size_t>(Theta(beta, i));
4981 cost += bin<size_t>(multiplicity + 1, 2);
4982 max_multiplicity =
std::max(max_multiplicity, multiplicity);
4985 const size_t h = GSA_degree_bound(k, ell, cost);
4986 std::vector<size_t> width(ell + 1);
4987 std::vector<size_t> offset(ell + 1);
4988 size_t n_unknowns = 0;
4989 for (size_t mu = 0; mu <= ell; ++mu) {
4990 width[mu] = GSA_block_size(k, h, mu);
4991 offset[mu] = n_unknowns;
4992 n_unknowns += width[mu];
4997 std::max(max_multiplicity > 0 ? max_multiplicity - 1 : 0, ell), p);
4998 auto M = ZeroMatrix<T>(cost, n_unknowns);
5000 for (size_t i = 0; i < np; ++i) {
5001 for (size_t beta = 0; beta < q; ++beta) {
5002 const size_t multiplicity =
static_cast<size_t>(Theta(beta, i));
5003 if (multiplicity == 0)
continue;
5004 const T y(
static_cast<
int>(beta));
5005 for (size_t s = 0; s < multiplicity; ++s) {
5006 for (size_t t = 0; s + t < multiplicity; ++t) {
5008 for (size_t mu = t; mu <= ell; ++mu) {
5009 const size_t ct = triangle(mu, t);
5012 for (size_t nu = s; nu < width[mu]; ++nu) {
5013 const size_t binomials = (ct * triangle(nu, s)) % p;
5015 M.set_component(row, offset[mu] + nu, binomials * (xpow * ypow));
5027 const auto B = M.basis_of_kernel();
5028 if (B.get_m() == 0)
return {};
5029 std::vector<Polynomial<T>> Q(ell + 1, ZeroPolynomial<T>());
5030 for (size_t mu = 0; mu <= ell; ++mu)
5031 for (size_t nu = 0; nu < width[mu]; ++nu) Q[mu].set_coefficient(nu, B(0, offset[mu] + nu));
5033 std::vector<Polynomial<T>> roots;
5034 auto g = ZeroPolynomial<T>();
5035 find_y_roots(Q, k, 0, g, roots, triangle);
5037 std::vector<Vector<T>> res;
5038 for (size_t j = 0; j < roots.size(); ++j) res.push_back(pad_back(
Vector<T>(roots[j]), k) * canonical_G());
5043 if (a.get_n() != d.get_n())
5044 throw std::invalid_argument(
"code locators a and column multipliers d must have the same length");
5046 throw std::invalid_argument(
5047 "code dimension must not exceed the number of code locators/column multipliers");
5048 if (!a.is_pairwise_distinct())
5049 throw std::invalid_argument(
"GRS codes must have pairwise distinct code locators!");
5050 return VandermondeMatrix<T>(a, k) * DiagonalMatrix<T>(d);
5053 GRSCode(
std::pair<
Vector<T>,
Vector<T>> params, size_t k) :
GRSCode<T>(params.first, params.second, k) {}