Skip to content

const-correct geometric lagrange #150

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
spaceships merged 1 commit into from
Aug 27, 2026
Merged

const-correct geometric lagrange #150

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 34 additions & 44 deletions libs/math/include/nil/crypto3/math/domains/arithmetic_sequence_domain.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,21 +56,28 @@ namespace nil {
std::vector<field_value_type> arithmetic_sequence;
field_value_type arithmetic_generator;

private:
// Queries derive domain points directly so they do not populate the transform precomputation.
field_value_type arithmetic_sequence_element(const std::size_t idx) const {
return arithmetic_generator * field_value_type(idx);
}

public:
void do_precomputation() {
compute_subproduct_tree<FieldType>(this->subproduct_tree, log2(this->m));

arithmetic_generator = field_value_type(fields::arithmetic_params<FieldType>::arithmetic_generator);

arithmetic_sequence = std::vector<field_value_type>(this->m);

for (std::size_t i = 0; i < arithmetic_sequence.size(); ++i) {
this->arithmetic_sequence[i] = this->arithmetic_generator * field_value_type(i);
this->arithmetic_sequence[i] = arithmetic_sequence_element(i);
}

precomputation_sentinel = true;
}

arithmetic_sequence_domain(const std::size_t m) : evaluation_domain<FieldType, ValueType>(m) {
arithmetic_sequence_domain(const std::size_t m) :
evaluation_domain<FieldType, ValueType>(m), precomputation_sentinel(false),
arithmetic_generator(field_value_type(fields::arithmetic_params<FieldType>::arithmetic_generator)) {
if (m <= 1) {
throw std::invalid_argument("arithmetic(): expected m > 1");
}
Expand All @@ -80,8 +87,6 @@ namespace nil {
"arithmetic(): expected arithmetic_params<FieldType>::arithmetic_generator.is_zero() "
"!= true");
}

precomputation_sentinel = false;
}

void fft(std::vector<value_type> &a) override {
Expand Down Expand Up @@ -163,21 +168,18 @@ namespace nil {
throw std::logic_error {"Not implemented yet"};
}

std::vector<field_value_type> evaluate_all_lagrange_polynomials(const field_value_type &t) override {
std::vector<field_value_type>
evaluate_all_lagrange_polynomials(const field_value_type &t) const override {
/* Compute Lagrange polynomial of size m, with m+1 points (x_0, y_0), ... ,(x_m, y_m) */
/* Evaluate for x = t */
/* Return coeffs for each l_j(x) = (l / l_i[j]) * w[j] */

if (!precomputation_sentinel)
do_precomputation();

/**
* If t equals one of the arithmetic progression values,
* then output 1 at the right place, and 0 elsewhere.
*/
for (std::size_t i = 0; i < this->m; ++i) {
if (arithmetic_sequence[i] == t) // i.e., t equals this->arithmetic_sequence[i]
{
if (arithmetic_sequence_element(i) == t) {
std::vector<field_value_type> res(this->m, field_value_type::zero());
res[i] = field_value_type::one();
return res;
Expand All @@ -189,15 +191,15 @@ namespace nil {
* then compute each Lagrange coefficient.
*/
std::vector<field_value_type> l(this->m);
l[0] = t - this->arithmetic_sequence[0];
l[0] = t - arithmetic_sequence_element(0);

field_value_type l_vanish = l[0];
field_value_type g_vanish = field_value_type::one();

for (std::size_t i = 1; i < this->m; i++) {
l[i] = t - this->arithmetic_sequence[i];
l[i] = t - arithmetic_sequence_element(i);
l_vanish *= l[i];
g_vanish *= -this->arithmetic_sequence[i];
g_vanish *= -arithmetic_sequence_element(i);
}

std::vector<field_value_type> w(this->m);
Expand All @@ -206,8 +208,8 @@ namespace nil {
l[0] = l_vanish * l[0].inversed() * w[0];
for (std::size_t i = 1; i < this->m; i++) {
field_value_type num =
this->arithmetic_sequence[i - 1] - this->arithmetic_sequence[this->m - 1];
w[i] = w[i - 1] * num * this->arithmetic_sequence[i].inversed();
arithmetic_sequence_element(i - 1) - arithmetic_sequence_element(this->m - 1);
w[i] = w[i - 1] * num * arithmetic_sequence_element(i).inversed();
l[i] = l_vanish * l[i].inversed() * w[i];
}

Expand All @@ -216,7 +218,7 @@ namespace nil {

std::vector<value_type> evaluate_all_lagrange_polynomials(
const typename std::vector<value_type>::const_iterator &t_powers_begin,
const typename std::vector<value_type>::const_iterator &t_powers_end) override {
const typename std::vector<value_type>::const_iterator &t_powers_end) const override {
if (std::size_t(std::distance(t_powers_begin, t_powers_end)) < this->m) {
throw std::invalid_argument(
"arithmetic_sequence_radix2: expected std::distance(t_powers_begin, t_powers_end) >= "
Expand All @@ -227,16 +229,12 @@ namespace nil {
/* Evaluate for x = t */
/* Return coeffs for each l_j(x) = (l / l_i[j]) * w[j] */

if (!precomputation_sentinel)
do_precomputation();

/**
* If t equals one of the arithmetic progression values,
* then output 1 at the right place, and 0 elsewhere.
*/
for (std::size_t i = 0; i < this->m; ++i) {
if (arithmetic_sequence[i] * t_powers_begin[0] == t_powers_begin[1]) // i.e., t equals a[i]
{
if (arithmetic_sequence_element(i) * t_powers_begin[0] == t_powers_begin[1]) {
std::vector<value_type> res(this->m, value_type::zero());
res[i] = t_powers_begin[0];
return res;
Expand All @@ -248,16 +246,16 @@ namespace nil {
* then compute each Lagrange coefficient.
*/
std::vector<polynomial<field_value_type>> l(this->m);
l[0] = polynomial<field_value_type>({-arithmetic_sequence[0], field_value_type::one()});
l[0] = polynomial<field_value_type>({-arithmetic_sequence_element(0), field_value_type::one()});
;

polynomial<field_value_type> l_vanish = l[0];
field_value_type g_vanish = field_value_type::one();

for (std::size_t i = 1; i < this->m; i++) {
l[i] = polynomial<field_value_type>({-arithmetic_sequence[i], field_value_type::one()});
l[i] = polynomial<field_value_type>({-arithmetic_sequence_element(i), field_value_type::one()});
l_vanish = l_vanish * l[i];
g_vanish *= -this->arithmetic_sequence[i];
g_vanish *= -arithmetic_sequence_element(i);
}

std::vector<field_value_type> w(this->m);
Expand All @@ -276,8 +274,8 @@ namespace nil {

for (std::size_t i = 1; i < this->m; i++) {
field_value_type num =
this->arithmetic_sequence[i - 1] - this->arithmetic_sequence[this->m - 1];
w[i] = w[i - 1] * num * this->arithmetic_sequence[i].inversed();
arithmetic_sequence_element(i - 1) - arithmetic_sequence_element(this->m - 1);
w[i] = w[i - 1] * num * arithmetic_sequence_element(i).inversed();

for (std::size_t j = 0; j < l[i].size(); ++j) {
result[i] = result[i] + t_powers_begin[j] * l[i][j];
Expand All @@ -289,36 +287,28 @@ namespace nil {
}

// This one is not the unity root actually, but it's ok for our purposes.
const field_value_type &get_unity_root() override {
const field_value_type &get_unity_root() const override {
return arithmetic_generator;
}

field_value_type get_domain_element(const std::size_t idx) override {
if (!this->precomputation_sentinel)
do_precomputation();

return this->arithmetic_sequence[idx];
field_value_type get_domain_element(const std::size_t idx) const override {
return arithmetic_sequence_element(idx);
}

field_value_type compute_vanishing_polynomial(const field_value_type &t) override {
if (!this->precomputation_sentinel)
do_precomputation();

field_value_type compute_vanishing_polynomial(const field_value_type &t) const override {
/* Notes: Z = prod_{i = 0 to m} (t - a[i]) */
field_value_type Z = field_value_type::one();
for (std::size_t i = 0; i < this->m; i++) {
Z *= (t - this->arithmetic_sequence[i]);
Z *= (t - arithmetic_sequence_element(i));
}
return Z;
}

polynomial<field_value_type> get_vanishing_polynomial() override {
if (!precomputation_sentinel)
do_precomputation();

polynomial<field_value_type> get_vanishing_polynomial() const override {
polynomial<field_value_type> z({field_value_type::one()});
for (std::size_t i = 0; i < this->m; i++) {
z = z * polynomial<field_value_type>({-arithmetic_sequence[i], field_value_type::one()});
z = z *
polynomial<field_value_type>({-arithmetic_sequence_element(i), field_value_type::one()});
}
return z;
}
Expand Down
47 changes: 26 additions & 21 deletions libs/math/include/nil/crypto3/math/domains/basic_radix2_domain.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,23 @@ namespace nil {
detail::create_fft_cache<FieldType>(this->m, omega.inversed(), fft_cache->second);
}

void inverse_fft_impl(std::vector<value_type> &a) const {
if (a.size() != this->m) {
if (a.size() < this->m) {
a.resize(this->m, value_type::zero());
} else {
throw std::invalid_argument("basic_radix2: expected a.size() == this->m");
}
}

detail::basic_radix2_fft_cached<FieldType>(a, fft_cache->second);

const field_value_type sconst = field_value_type(this->m).inversed();
for (value_type &a_i : a) {
a_i *= sconst;
}
}

public:
typedef FieldType field_type;
using evaluation_domain<FieldType, ValueType>::evaluate_all_lagrange_polynomials;
Expand Down Expand Up @@ -148,51 +165,39 @@ namespace nil {
}

void inverse_fft(std::vector<value_type> &a) override {
if (a.size() != this->m) {
if (a.size() < this->m) {
a.resize(this->m, value_type::zero());
} else {
throw std::invalid_argument("basic_radix2: expected a.size() == this->m");
}
}

detail::basic_radix2_fft_cached<FieldType>(a, fft_cache->second);

const field_value_type sconst = field_value_type(this->m).inversed();
for (value_type &a_i : a) {
a_i *= sconst;
}
inverse_fft_impl(a);
}

std::vector<field_value_type> evaluate_all_lagrange_polynomials(const field_value_type &t) override {
std::vector<field_value_type>
evaluate_all_lagrange_polynomials(const field_value_type &t) const override {
return detail::basic_radix2_evaluate_all_lagrange_polynomials<FieldType>(this->m, t);
}

std::vector<value_type> evaluate_all_lagrange_polynomials(
const typename std::vector<value_type>::const_iterator &t_powers_begin,
const typename std::vector<value_type>::const_iterator &t_powers_end) override {
const typename std::vector<value_type>::const_iterator &t_powers_end) const override {
if (std::size_t(std::distance(t_powers_begin, t_powers_end)) < this->m) {
throw std::invalid_argument(
"basic_radix2: expected std::distance(t_powers_begin, t_powers_end) >= this->m");
}
std::vector<value_type> tmp(t_powers_begin, t_powers_begin + this->m);
this->inverse_fft(tmp);
inverse_fft_impl(tmp);
return tmp;
}

const field_value_type &get_unity_root() override {
const field_value_type &get_unity_root() const override {
return omega;
}

field_value_type get_domain_element(const std::size_t idx) override {
field_value_type get_domain_element(const std::size_t idx) const override {
return omega.pow(idx);
}

field_value_type compute_vanishing_polynomial(const field_value_type &t) override {
field_value_type compute_vanishing_polynomial(const field_value_type &t) const override {
return (t.pow(this->m)) - field_value_type::one();
}

polynomial<field_value_type> get_vanishing_polynomial() override {
polynomial<field_value_type> get_vanishing_polynomial() const override {
polynomial<field_value_type> z(this->m + 1, field_value_type::zero());
z[this->m] = field_value_type::one();
z[0] = -field_value_type::one();
Expand Down
15 changes: 8 additions & 7 deletions libs/math/include/nil/crypto3/math/domains/evaluation_domain.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,12 +69,12 @@ namespace nil {
/**
* Get the unity root.
*/
virtual const field_value_type &get_unity_root() = 0;
virtual const field_value_type &get_unity_root() const = 0;

/**
* Get the idx-th element in S.
*/
virtual field_value_type get_domain_element(const std::size_t idx) = 0;
virtual field_value_type get_domain_element(const std::size_t idx) const = 0;

/**
* Compute the FFT, over the domain S, of the vector a.
Expand Down Expand Up @@ -105,7 +105,8 @@ namespace nil {
* The output is a vector (b_{0},...,b_{m-1})
* where b_{i} is the evaluation of L_{i,S}(z) at z = t.
*/
virtual std::vector<field_value_type> evaluate_all_lagrange_polynomials(const field_value_type &t) = 0;
virtual std::vector<field_value_type>
evaluate_all_lagrange_polynomials(const field_value_type &t) const = 0;

/**
* Evaluate all Lagrange polynomials and the domain vanishing polynomial at t.
Expand All @@ -115,7 +116,7 @@ namespace nil {
*/
virtual std::vector<field_value_type>
evaluate_all_lagrange_polynomials(const field_value_type &t,
field_value_type &vanishing_polynomial_at_t) {
field_value_type &vanishing_polynomial_at_t) const {
std::vector<field_value_type> result = evaluate_all_lagrange_polynomials(t);
vanishing_polynomial_at_t = compute_vanishing_polynomial(t);
return result;
Expand All @@ -132,17 +133,17 @@ namespace nil {
*/
virtual std::vector<value_type> evaluate_all_lagrange_polynomials(
const typename std::vector<value_type>::const_iterator &t_powers_begin,
const typename std::vector<value_type>::const_iterator &t_powers_end) = 0;
const typename std::vector<value_type>::const_iterator &t_powers_end) const = 0;

/**
* Evaluate the vanishing polynomial of S at the field element t.
*/
virtual field_value_type compute_vanishing_polynomial(const field_value_type &t) = 0;
virtual field_value_type compute_vanishing_polynomial(const field_value_type &t) const = 0;

/**
* Build the vanishing polynomial of S.
*/
virtual polynomial<field_value_type> get_vanishing_polynomial() = 0;
virtual polynomial<field_value_type> get_vanishing_polynomial() const = 0;

/**
* Add the coefficients of the vanishing polynomial of S to the coefficients of the polynomial H.
Expand Down
Loading
Loading