diff --git a/src/Persistence_matrix/concept/FieldOperators.h b/src/Persistence_matrix/concept/FieldOperators.h index fc45d8453b..c8f3678f92 100644 --- a/src/Persistence_matrix/concept/FieldOperators.h +++ b/src/Persistence_matrix/concept/FieldOperators.h @@ -42,7 +42,8 @@ class FieldOperators FieldOperators(characteristic_type characteristic = 0); /** - * @brief Sets the characteristic of the field. + * @brief Sets the characteristic of the field. Can eventually be omitted if the characteristic of the class + * is fixed. * * @param characteristic Prime number corresponding to the desired characteristic of the field. */ @@ -107,8 +108,8 @@ class FieldOperators // * Returns the result in the field. // * // * @param e First element. - // * @param m Second element. - // * @param a Third element. + // * @param a Second element. + // * @param m Third element. // * @return `((e + a) * m) % characteristic`, such that the result is positive. // */ // element_type add_and_multiply(element_type e, element_type a, element_type m) const; diff --git a/src/Persistence_matrix/doc/Intro_field_elements_and_operators.h b/src/Persistence_matrix/doc/Intro_field_elements_and_operators.h index 30db408721..3b515c07de 100644 --- a/src/Persistence_matrix/doc/Intro_field_elements_and_operators.h +++ b/src/Persistence_matrix/doc/Intro_field_elements_and_operators.h @@ -17,10 +17,20 @@ namespace persistence_fields { /** \defgroup persistence_fields Persistence Fields * @{ - * \author Hannah Schreiber + * \author Hannah Schreiber, Clément Maria * - * TODO: - * + * Set of classes allowing addition and multiplication, as well as inverse computation, in \f$ \mathbb{F}_p \f$ fields, + * with \f$ p \f$ some prime number, or in multi-fields as defined in \cite boissonnat:hal-00922572. + * + * There are two types of classes: + * - those defining directly a field element, allowing to use them as any integer: the operators are overwritten such + * that the calculation is done in the field. For example, if \f$ e = 2 \f$ is an instanciation of an + * \f$ \mathbb{F}_3 \f$ element class, then `e + 3` returns `2`, + * - those only defining the operators of a field or multi-field. They represent a collection of methods taking + * one or two integers as input and treating them as elements of the field. For example, if \f$ op \f$ is an + * instanciation of a \f$ \mathbb{F}_3 \f$ operator class, `op.add(2, 3)` returns `2`. + * + * The field operator classes all respect the @ref persistence_matrix::FieldOperators concept. * * \subsection fieldsexamples Examples * diff --git a/src/Persistence_matrix/include/gudhi/Fields/Multi_field.h b/src/Persistence_matrix/include/gudhi/Fields/Multi_field.h index ce69f46480..6d0e3fd4b3 100644 --- a/src/Persistence_matrix/include/gudhi/Fields/Multi_field.h +++ b/src/Persistence_matrix/include/gudhi/Fields/Multi_field.h @@ -1,13 +1,19 @@ /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. - * Author(s): Hannah Schreiber + * Author(s): Hannah Schreiber, Clément Maria * - * Copyright (C) 2022 Inria + * Copyright (C) 2022-24 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification */ +/** + * @file Multi_field.h + * @author Hannah Schreiber + * @brief Contains the @ref Multi_field_element class. + */ + #ifndef MATRIX_FIELD_MULTI_H_ #define MATRIX_FIELD_MULTI_H_ @@ -20,409 +26,482 @@ namespace Gudhi { namespace persistence_fields { /** + * @class Multi_field_element Multi_field.h gudhi/Fields/Multi_field.h * @ingroup persistence_fields * - * @brief TODO: - * + * @brief Class representing an element of a multi-field. + * The characteristics will corresponds to all prime numbers in the interval given as template. + * + * @tparam minimum Interval closed lower bound. + * @tparam maximum Interval closed upper bound. */ -template -class Multi_field_element { -public: - using element_type = mpz_class; - - Multi_field_element(); - Multi_field_element(mpz_class element); - Multi_field_element(const Multi_field_element& toCopy); - Multi_field_element(Multi_field_element&& toMove) noexcept; - - friend void operator+=(Multi_field_element& f1, Multi_field_element const& f2){ - mpz_add(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), f2.element_.get_mpz_t()); - mpz_mod(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - } - friend Multi_field_element operator+(Multi_field_element f1, Multi_field_element const& f2){ - f1 += f2; - return f1; - } - friend void operator+=(Multi_field_element& f, mpz_class const v){ - mpz_add(f.element_.get_mpz_t(), f.element_.get_mpz_t(), v.get_mpz_t()); - mpz_mod(f.element_.get_mpz_t(), f.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - } - friend Multi_field_element operator+(Multi_field_element f, mpz_class const v){ - f += v; - return f; - } - friend mpz_class operator+(mpz_class v, Multi_field_element const& f){ - mpz_class e(v); - mpz_add(e.get_mpz_t(), e.get_mpz_t(), f.element_.get_mpz_t()); - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e; - } - - friend void operator-=(Multi_field_element& f1, Multi_field_element const& f2){ - mpz_sub(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), f2.element_.get_mpz_t()); - mpz_mod(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - } - friend Multi_field_element operator-(Multi_field_element f1, Multi_field_element const& f2){ - f1 -= f2; - return f1; - } - friend void operator-=(Multi_field_element& f, mpz_class const v){ - mpz_sub(f.element_.get_mpz_t(), f.element_.get_mpz_t(), v.get_mpz_t()); - mpz_mod(f.element_.get_mpz_t(), f.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - } - friend Multi_field_element operator-(Multi_field_element f, mpz_class const v){ - f -= v; - return f; - } - friend mpz_class operator-(mpz_class v, Multi_field_element const& f){ - mpz_class e(v); - if (e >= productOfAllCharacteristics_) - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - if (f.element_ > e) - mpz_add(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - mpz_sub(e.get_mpz_t(), e.get_mpz_t(), f.element_.get_mpz_t()); - return e; - } - - friend void operator*=(Multi_field_element& f1, Multi_field_element const& f2){ - mpz_mul(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), f2.element_.get_mpz_t()); - mpz_mod(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - } - friend Multi_field_element operator*(Multi_field_element f1, Multi_field_element const& f2){ - f1 *= f2; - return f1; - } - friend void operator*=(Multi_field_element& f, mpz_class const v){ - mpz_mul(f.element_.get_mpz_t(), f.element_.get_mpz_t(), v.get_mpz_t()); - mpz_mod(f.element_.get_mpz_t(), f.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - } - friend Multi_field_element operator*(Multi_field_element f, mpz_class const v){ - f *= v; - return f; - } - friend mpz_class operator*(mpz_class v, Multi_field_element const& f){ - mpz_class e(v); - mpz_mul(e.get_mpz_t(), e.get_mpz_t(), f.element_.get_mpz_t()); - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e; - } - - friend bool operator==(const Multi_field_element& f1, const Multi_field_element& f2){ - return f1.element_ == f2.element_; - } - friend bool operator==(const mpz_class v, const Multi_field_element& f){ - if (v < productOfAllCharacteristics_) return v == f.element_; - mpz_class e(v); - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e == f.element_; - } - friend bool operator==(const Multi_field_element& f, const mpz_class v){ - if (v < productOfAllCharacteristics_) return v == f.element_; - mpz_class e(v); - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e == f.element_; - } - friend bool operator==(const unsigned int v, const Multi_field_element& f){ - mpz_class e(v); - if (e < productOfAllCharacteristics_) return e == f.element_; - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e == f.element_; - } - friend bool operator==(const Multi_field_element& f, const unsigned int v){ - mpz_class e(v); - if (e < productOfAllCharacteristics_) return e == f.element_; - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e == f.element_; - } - friend bool operator!=(const Multi_field_element& f1, const Multi_field_element& f2){ - return !(f1 == f2); - } - friend bool operator!=(const mpz_class v, const Multi_field_element& f){ - return !(v == f); - } - friend bool operator!=(const Multi_field_element& f, const mpz_class v){ - return !(v == f); - } - friend bool operator!=(const unsigned int v, const Multi_field_element& f){ - return !(v == f); - } - friend bool operator!=(const Multi_field_element& f, const unsigned int v){ - return !(v == f); - } - - Multi_field_element& operator=(Multi_field_element other); - Multi_field_element& operator=(const mpz_class value); - operator unsigned int() const; - operator mpz_class() const; - friend void swap(Multi_field_element& f1, Multi_field_element& f2){ - std::swap(f1.element_, f2.element_); - } - - Multi_field_element get_inverse() const; - std::pair get_partial_inverse(const mpz_class& productOfCharacteristics) const; - - static Multi_field_element get_additive_identity(); - static Multi_field_element get_multiplicative_identity(); - static Multi_field_element get_partial_multiplicative_identity(const mpz_class& productOfCharacteristics); - static mpz_class get_characteristic(); - - mpz_class get_value() const; - - static constexpr bool handles_only_z2(){ - return false; - } - -private: - mpz_class element_; - static inline const std::vector primes_ = [](){ - std::vector res; - - unsigned int curr_prime = minimum; - mpz_t tmp_prime; - mpz_init_set_ui(tmp_prime, minimum); - // test if min_prime is prime - int is_prime = mpz_probab_prime_p(tmp_prime, 25); // probabilistic primality test - - if (is_prime == 0) { // min_prime is composite - mpz_nextprime(tmp_prime, tmp_prime); - curr_prime = mpz_get_ui(tmp_prime); - } - - while (curr_prime <= maximum) { - res.push_back(curr_prime); - mpz_nextprime(tmp_prime, tmp_prime); - curr_prime = mpz_get_ui(tmp_prime); - } - mpz_clear(tmp_prime); - - return res; - }(); - static inline const mpz_class productOfAllCharacteristics_ = [](){ - mpz_class res = 1; - for (const auto p : primes_) { - res *= p; - } - - return res; - }(); - static inline const std::vector partials_ = [](){ - std::vector res; - - if (productOfAllCharacteristics_ == 1) return res; - - for (unsigned int i = 0; i < primes_.size(); ++i){ - unsigned int p = primes_[i]; - res.push_back(productOfAllCharacteristics_ / p); - mpz_powm_ui(res.back().get_mpz_t(), res.back().get_mpz_t(), p - 1, - productOfAllCharacteristics_.get_mpz_t()); - } - - return res; - }(); - //If I understood the paper well, multiplicativeID_ always equals to 1. But in Clement's code, - //multiplicativeID_ is computed (see commented lambda function below). TODO: verify with Clement. - static inline const mpz_class multiplicativeID_ = 1;/*[](){ - mpz_class res = 0; - for (unsigned int i = 0; i < partials_.size(); ++i){ - res = (res + partials_[i]) % productOfAllCharacteristics_; - } - - return res; - }();*/ - - static constexpr bool _is_prime(const int p); +template +class Multi_field_element +{ + public: + using element_type = mpz_class; /**< Type for the elements in the field. */ + + /** + * @brief Default constructor. Sets the element to 0. + */ + Multi_field_element(); + /** + * @brief Constructor setting the element to the given value. + * + * @param element Value of the element. + */ + Multi_field_element(mpz_class element); + /** + * @brief Copy constructor. + * + * @param toCopy Element to copy. + */ + Multi_field_element(const Multi_field_element& toCopy); + /** + * @brief Move constructor. + * + * @param toMove Element to move. + */ + Multi_field_element(Multi_field_element&& toMove) noexcept; + + /** + * @brief operator+= + */ + friend void operator+=(Multi_field_element& f1, Multi_field_element const& f2) { + mpz_add(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), f2.element_.get_mpz_t()); + mpz_mod(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + } + /** + * @brief operator+ + */ + friend Multi_field_element operator+(Multi_field_element f1, Multi_field_element const& f2) { + f1 += f2; + return f1; + } + /** + * @brief operator+= + */ + friend void operator+=(Multi_field_element& f, mpz_class const v) { + mpz_add(f.element_.get_mpz_t(), f.element_.get_mpz_t(), v.get_mpz_t()); + mpz_mod(f.element_.get_mpz_t(), f.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + } + /** + * @brief operator+ + */ + friend Multi_field_element operator+(Multi_field_element f, mpz_class const v) { + f += v; + return f; + } + /** + * @brief operator+ + */ + friend mpz_class operator+(mpz_class v, Multi_field_element const& f) { + mpz_class e(v); + mpz_add(e.get_mpz_t(), e.get_mpz_t(), f.element_.get_mpz_t()); + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e; + } + + /** + * @brief operator-= + */ + friend void operator-=(Multi_field_element& f1, Multi_field_element const& f2) { + mpz_sub(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), f2.element_.get_mpz_t()); + mpz_mod(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + } + /** + * @brief operator- + */ + friend Multi_field_element operator-(Multi_field_element f1, Multi_field_element const& f2) { + f1 -= f2; + return f1; + } + /** + * @brief operator-= + */ + friend void operator-=(Multi_field_element& f, mpz_class const v) { + mpz_sub(f.element_.get_mpz_t(), f.element_.get_mpz_t(), v.get_mpz_t()); + mpz_mod(f.element_.get_mpz_t(), f.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + } + /** + * @brief operator- + */ + friend Multi_field_element operator-(Multi_field_element f, mpz_class const v) { + f -= v; + return f; + } + /** + * @brief operator- + */ + friend mpz_class operator-(mpz_class v, Multi_field_element const& f) { + mpz_class e(v); + if (e >= productOfAllCharacteristics_) + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + if (f.element_ > e) mpz_add(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + mpz_sub(e.get_mpz_t(), e.get_mpz_t(), f.element_.get_mpz_t()); + return e; + } + + /** + * @brief operator*= + */ + friend void operator*=(Multi_field_element& f1, Multi_field_element const& f2) { + mpz_mul(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), f2.element_.get_mpz_t()); + mpz_mod(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + } + /** + * @brief operator* + */ + friend Multi_field_element operator*(Multi_field_element f1, Multi_field_element const& f2) { + f1 *= f2; + return f1; + } + /** + * @brief operator*= + */ + friend void operator*=(Multi_field_element& f, mpz_class const v) { + mpz_mul(f.element_.get_mpz_t(), f.element_.get_mpz_t(), v.get_mpz_t()); + mpz_mod(f.element_.get_mpz_t(), f.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + } + /** + * @brief operator* + */ + friend Multi_field_element operator*(Multi_field_element f, mpz_class const v) { + f *= v; + return f; + } + /** + * @brief operator* + */ + friend mpz_class operator*(mpz_class v, Multi_field_element const& f) { + mpz_class e(v); + mpz_mul(e.get_mpz_t(), e.get_mpz_t(), f.element_.get_mpz_t()); + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e; + } + + /** + * @brief operator== + */ + friend bool operator==(const Multi_field_element& f1, const Multi_field_element& f2) { + return f1.element_ == f2.element_; + } + /** + * @brief operator== + */ + friend bool operator==(const mpz_class v, const Multi_field_element& f) { + if (v < productOfAllCharacteristics_) return v == f.element_; + mpz_class e(v); + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e == f.element_; + } + /** + * @brief operator== + */ + friend bool operator==(const Multi_field_element& f, const mpz_class v) { + if (v < productOfAllCharacteristics_) return v == f.element_; + mpz_class e(v); + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e == f.element_; + } + /** + * @brief operator== + */ + friend bool operator==(const unsigned int v, const Multi_field_element& f) { + mpz_class e(v); + if (e < productOfAllCharacteristics_) return e == f.element_; + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e == f.element_; + } + /** + * @brief operator== + */ + friend bool operator==(const Multi_field_element& f, const unsigned int v) { + mpz_class e(v); + if (e < productOfAllCharacteristics_) return e == f.element_; + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e == f.element_; + } + /** + * @brief operator!= + */ + friend bool operator!=(const Multi_field_element& f1, const Multi_field_element& f2) { return !(f1 == f2); } + /** + * @brief operator!= + */ + friend bool operator!=(const mpz_class v, const Multi_field_element& f) { return !(v == f); } + /** + * @brief operator!= + */ + friend bool operator!=(const Multi_field_element& f, const mpz_class v) { return !(v == f); } + /** + * @brief operator!= + */ + friend bool operator!=(const unsigned int v, const Multi_field_element& f) { return !(v == f); } + /** + * @brief operator!= + */ + friend bool operator!=(const Multi_field_element& f, const unsigned int v) { return !(v == f); } + + /** + * @brief Assign operator. + */ + Multi_field_element& operator=(Multi_field_element other); + /** + * @brief Assign operator. + */ + Multi_field_element& operator=(const mpz_class value); + /** + * @brief Swap operator. + */ + friend void swap(Multi_field_element& f1, Multi_field_element& f2) { std::swap(f1.element_, f2.element_); } + + /** + * @brief Casts the element into an unsigned int. + */ + operator unsigned int() const; + /** + * @brief Casts the element into an mpz_class. + */ + operator mpz_class() const; + + /** + * @brief Returns the inverse of the element in the multi-field, see @cite boissonnat:hal-00922572. + * + * @return The inverse. + */ + Multi_field_element get_inverse() const; + /** + * @brief Returns the inverse of the element with respect to a sub-product of the characteristics in the multi-field, + * see @cite boissonnat:hal-00922572. + * + * @param productOfCharacteristics Sub-product of the characteristics. + * @return Pair of the inverse and the characteristic the inverse corresponds to. + */ + std::pair get_partial_inverse(const mpz_class& productOfCharacteristics) const; + + /** + * @brief Returns the additive identity of a field. + * + * @return The additive identity of a field. + */ + static Multi_field_element get_additive_identity(); + /** + * @brief Returns the multiplicative identity of a field. + * + * @return The multiplicative identity of a field. + */ + static Multi_field_element get_multiplicative_identity(); + /** + * @brief Returns the partial multiplicative identity of the multi-field from the given product. + * See @cite boissonnat:hal-00922572 for more details. + * + * @param productOfCharacteristics Product of the different characteristics to take into account in the multi-field. + * @return The partial multiplicative identity of the multi-field. + */ + static Multi_field_element get_partial_multiplicative_identity(const mpz_class& productOfCharacteristics); + /** + * @brief Returns the product of all characteristics. + * + * @return The product of all characteristics. + */ + static mpz_class get_characteristic(); + + /** + * @brief Returns the value of the element. + * + * @return Value of the element. + */ + mpz_class get_value() const; + + // static constexpr bool handles_only_z2() { return false; } + + private: + mpz_class element_; + static inline const std::vector primes_ = []() { + std::vector res; + + unsigned int curr_prime = minimum; + mpz_t tmp_prime; + mpz_init_set_ui(tmp_prime, minimum); + // test if min_prime is prime + int is_prime = mpz_probab_prime_p(tmp_prime, 25); // probabilistic primality test + + if (is_prime == 0) { // min_prime is composite + mpz_nextprime(tmp_prime, tmp_prime); + curr_prime = mpz_get_ui(tmp_prime); + } + + while (curr_prime <= maximum) { + res.push_back(curr_prime); + mpz_nextprime(tmp_prime, tmp_prime); + curr_prime = mpz_get_ui(tmp_prime); + } + mpz_clear(tmp_prime); + + return res; + }(); + static inline const mpz_class productOfAllCharacteristics_ = []() { + mpz_class res = 1; + for (const auto p : primes_) { + res *= p; + } + + return res; + }(); + static inline const std::vector partials_ = []() { + std::vector res; + + if (productOfAllCharacteristics_ == 1) return res; + + for (unsigned int i = 0; i < primes_.size(); ++i) { + unsigned int p = primes_[i]; + res.push_back(productOfAllCharacteristics_ / p); + mpz_powm_ui(res.back().get_mpz_t(), res.back().get_mpz_t(), p - 1, productOfAllCharacteristics_.get_mpz_t()); + } + + return res; + }(); + // If I understood the paper well, multiplicativeID_ always equals to 1. But in Clement's code, + // multiplicativeID_ is computed (see commented lambda function below). TODO: verify with Clement. + static inline const mpz_class multiplicativeID_ = 1; /*[](){ + mpz_class res = 0; + for (unsigned int i = 0; i < partials_.size(); ++i){ + res = (res + partials_[i]) % productOfAllCharacteristics_; + } + + return res; + }();*/ + + static constexpr bool _is_prime(const int p); }; -template -inline Multi_field_element::Multi_field_element() - : element_(0) -{ - static_assert(maximum >= 2, "Characteristics have to be positive."); - static_assert(minimum <= maximum, "The given interval is not valid."); - static_assert(minimum != maximum || _is_prime(minimum), "The given interval does not contain a prime number."); +template +inline Multi_field_element::Multi_field_element() : element_(0) { + static_assert(maximum >= 2, "Characteristics have to be positive."); + static_assert(minimum <= maximum, "The given interval is not valid."); + static_assert(minimum != maximum || _is_prime(minimum), "The given interval does not contain a prime number."); - if (productOfAllCharacteristics_ == 1) - throw std::runtime_error("The given interval does not contain a prime number."); + if (productOfAllCharacteristics_ == 1) + throw std::runtime_error("The given interval does not contain a prime number."); } -template -inline Multi_field_element::Multi_field_element(mpz_class element) - : element_(element) -{ - static_assert(maximum >= 2, "Characteristics has to be positive."); - static_assert(minimum <= maximum, "The given interval is not valid."); - static_assert(minimum != maximum || _is_prime(minimum), "The given interval does not contain a prime number."); +template +inline Multi_field_element::Multi_field_element(mpz_class element) : element_(element) { + static_assert(maximum >= 2, "Characteristics has to be positive."); + static_assert(minimum <= maximum, "The given interval is not valid."); + static_assert(minimum != maximum || _is_prime(minimum), "The given interval does not contain a prime number."); - if (productOfAllCharacteristics_ == 1) - throw std::runtime_error("The given interval does not contain a prime number."); + if (productOfAllCharacteristics_ == 1) + throw std::runtime_error("The given interval does not contain a prime number."); - mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); } -template -inline Multi_field_element::Multi_field_element(const Multi_field_element &toCopy) - : element_(toCopy.element_) -{} - -template -inline Multi_field_element::Multi_field_element(Multi_field_element &&toMove) noexcept - : element_(std::move(toMove.element_)) -{} - -//template -//inline Multi_field_element &Multi_field_element::operator+=(Multi_field_element const &f) -//{ -// mpz_add(element_.get_mpz_t(), element_.get_mpz_t(), f.element_.get_mpz_t()); -// mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); -// return *this; -//} - -//template -//inline Multi_field_element &Multi_field_element::operator+=(mpz_class const v) -//{ -// mpz_add(element_.get_mpz_t(), element_.get_mpz_t(), v.get_mpz_t()); -// mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); -// return *this; -//} - -//template -//inline Multi_field_element &Multi_field_element::operator-=(Multi_field_element const &f) -//{ -// mpz_sub(element_.get_mpz_t(), element_.get_mpz_t(), f.element_.get_mpz_t()); -// mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); -// return *this; -//} - -//template -//inline Multi_field_element &Multi_field_element::operator-=(mpz_class const v) -//{ -// mpz_sub(element_.get_mpz_t(), element_.get_mpz_t(), v.get_mpz_t()); -// mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); -// return *this; -//} - -//template -//inline Multi_field_element &Multi_field_element::operator*=(Multi_field_element const &f) -//{ -// mpz_mul(element_.get_mpz_t(), element_.get_mpz_t(), f.element_.get_mpz_t()); -// mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); -// return *this; -//} - -//template -//inline Multi_field_element &Multi_field_element::operator*=(mpz_class const v) -//{ -// mpz_mul(element_.get_mpz_t(), element_.get_mpz_t(), v.get_mpz_t()); -// mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); -// return *this; -//} - -template -inline Multi_field_element &Multi_field_element::operator=(Multi_field_element other) -{ - std::swap(element_, other.element_); - return *this; +template +inline Multi_field_element::Multi_field_element(const Multi_field_element& toCopy) + : element_(toCopy.element_) {} + +template +inline Multi_field_element::Multi_field_element( + Multi_field_element&& toMove) noexcept + : element_(std::move(toMove.element_)) {} + +template +inline Multi_field_element& Multi_field_element::operator=( + Multi_field_element other) { + std::swap(element_, other.element_); + return *this; } -template -inline Multi_field_element &Multi_field_element::operator=(mpz_class const value) -{ - mpz_mod(element_.get_mpz_t(), value.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return *this; +template +inline Multi_field_element& Multi_field_element::operator=(mpz_class const value) { + mpz_mod(element_.get_mpz_t(), value.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return *this; } -template -inline Multi_field_element::operator unsigned int() const -{ - return element_.get_ui(); +template +inline Multi_field_element::operator unsigned int() const { + return element_.get_ui(); } -template -inline Multi_field_element::operator mpz_class() const -{ - return element_; +template +inline Multi_field_element::operator mpz_class() const { + return element_; } -template -inline Multi_field_element Multi_field_element::get_inverse() const -{ - return get_partial_inverse(productOfAllCharacteristics_).first; +template +inline Multi_field_element Multi_field_element::get_inverse() const { + return get_partial_inverse(productOfAllCharacteristics_).first; } -template -inline std::pair,mpz_class> Multi_field_element::get_partial_inverse(const mpz_class& productOfCharacteristics) const -{ - mpz_class QR; - mpz_gcd(QR.get_mpz_t(), element_.get_mpz_t(), productOfCharacteristics.get_mpz_t()); // QR <- gcd(x,QS) +template +inline std::pair, mpz_class> +Multi_field_element::get_partial_inverse(const mpz_class& productOfCharacteristics) const { + mpz_class QR; + mpz_gcd(QR.get_mpz_t(), element_.get_mpz_t(), productOfCharacteristics.get_mpz_t()); // QR <- gcd(x,QS) - if (QR == productOfCharacteristics) - return {Multi_field_element(), multiplicativeID_}; // partial inverse is 0 + if (QR == productOfCharacteristics) return {Multi_field_element(), multiplicativeID_}; // partial inverse is 0 - mpz_class QT = productOfCharacteristics / QR; + mpz_class QT = productOfCharacteristics / QR; - mpz_class inv_qt; - mpz_invert(inv_qt.get_mpz_t(), element_.get_mpz_t(), QT.get_mpz_t()); + mpz_class inv_qt; + mpz_invert(inv_qt.get_mpz_t(), element_.get_mpz_t(), QT.get_mpz_t()); - auto res = get_partial_multiplicative_identity(QT); - res *= inv_qt; + auto res = get_partial_multiplicative_identity(QT); + res *= inv_qt; - return {res, QT}; + return {res, QT}; } -template -inline Multi_field_element Multi_field_element::get_additive_identity() -{ - return Multi_field_element(); +template +inline Multi_field_element Multi_field_element::get_additive_identity() { + return Multi_field_element(); } -template -inline Multi_field_element Multi_field_element::get_multiplicative_identity() -{ - return Multi_field_element(multiplicativeID_); +template +inline Multi_field_element Multi_field_element::get_multiplicative_identity() { + return Multi_field_element(multiplicativeID_); } -template -inline Multi_field_element Multi_field_element::get_partial_multiplicative_identity(const mpz_class& productOfCharacteristics) -{ - if (productOfCharacteristics == 0) { - return Multi_field_element(multiplicativeID_); - } - Multi_field_element mult; - for (unsigned int idx = 0; idx < primes_.size(); ++idx) { - if ((productOfCharacteristics % primes_[idx]) == 0) { - mult += partials_[idx]; - } - } - return mult; +template +inline Multi_field_element Multi_field_element::get_partial_multiplicative_identity( + const mpz_class& productOfCharacteristics) { + if (productOfCharacteristics == 0) { + return Multi_field_element(multiplicativeID_); + } + Multi_field_element mult; + for (unsigned int idx = 0; idx < primes_.size(); ++idx) { + if ((productOfCharacteristics % primes_[idx]) == 0) { + mult += partials_[idx]; + } + } + return mult; } -template -inline mpz_class Multi_field_element::get_characteristic() -{ - return productOfAllCharacteristics_; +template +inline mpz_class Multi_field_element::get_characteristic() { + return productOfAllCharacteristics_; } -template -inline mpz_class Multi_field_element::get_value() const -{ - return element_; +template +inline mpz_class Multi_field_element::get_value() const { + return element_; } -template -inline constexpr bool Multi_field_element::_is_prime(const int p) -{ - if (p <= 1) return false; - if (p <= 3) return true; - if (p % 2 == 0 || p % 3 == 0) return false; +template +inline constexpr bool Multi_field_element::_is_prime(const int p) { + if (p <= 1) return false; + if (p <= 3) return true; + if (p % 2 == 0 || p % 3 == 0) return false; - for (long i = 5; i * i <= p; i = i + 6) - if (p % i == 0 || p % (i + 2) == 0) - return false; + for (long i = 5; i * i <= p; i = i + 6) + if (p % i == 0 || p % (i + 2) == 0) return false; - return true; + return true; } -} //namespace persistence_fields -} //namespace Gudhi +} // namespace persistence_fields +} // namespace Gudhi #endif // MATRIX_FIELD_MULTI_H_ diff --git a/src/Persistence_matrix/include/gudhi/Fields/Multi_field_operators.h b/src/Persistence_matrix/include/gudhi/Fields/Multi_field_operators.h index 110825bc97..58a6e9473a 100644 --- a/src/Persistence_matrix/include/gudhi/Fields/Multi_field_operators.h +++ b/src/Persistence_matrix/include/gudhi/Fields/Multi_field_operators.h @@ -1,13 +1,19 @@ /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. - * Author(s): Hannah Schreiber + * Author(s): Hannah Schreiber, Clément Maria * - * Copyright (C) 2022 Inria + * Copyright (C) 2022-24 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification */ +/** + * @file Multi_field_operators.h + * @author Hannah Schreiber + * @brief Contains the @ref Multi_field_operators class. + */ + #ifndef MATRIX_FIELD_MULTI_OPERATORS_H_ #define MATRIX_FIELD_MULTI_OPERATORS_H_ @@ -20,220 +26,324 @@ namespace Gudhi { namespace persistence_fields { /** + * @class Multi_field_operators Multi_field_operators.h gudhi/Fields/Multi_field_operators.h * @ingroup persistence_fields * - * @brief TODO: - * + * @brief Class defining operators for a multi-field with "consecutive" charateristic range. */ -class Multi_field_operators { -public: - using element_type = mpz_class; - using characteristic_type = element_type; - - Multi_field_operators() : productOfAllCharacteristics_(0)/* , multiplicativeID_(1) */ {} - Multi_field_operators(int minCharacteristic, int maxCharacteristic) - : productOfAllCharacteristics_(0)//, multiplicativeID_(1) - { - set_characteristic(minCharacteristic, maxCharacteristic); - } - Multi_field_operators(const Multi_field_operators& toCopy) - : primes_(toCopy.primes_), - productOfAllCharacteristics_(toCopy.productOfAllCharacteristics_), - partials_(toCopy.partials_)/* , - multiplicativeID_(toCopy.multiplicativeID_) */ - {} - Multi_field_operators(Multi_field_operators&& toMove) noexcept - : primes_(std::move(toMove.primes_)), - productOfAllCharacteristics_(std::move(toMove.productOfAllCharacteristics_)), - partials_(std::move(toMove.partials_))/* , - multiplicativeID_(std::move(toMove.multiplicativeID_)) */ - {} - - void set_characteristic(int minimum, int maximum){ - if (maximum < 2) - throw std::invalid_argument("Characteristic must be strictly positive"); - if (minimum > maximum) - throw std::invalid_argument("The given interval is not valid."); - if (minimum == maximum && !_is_prime(minimum)) - throw std::invalid_argument("The given interval does not contain a prime number."); - - unsigned int curr_prime = minimum; - mpz_t tmp_prime; - mpz_init_set_ui(tmp_prime, minimum); - // test if min_prime is prime - int is_prime = mpz_probab_prime_p(tmp_prime, 25); // probabilistic primality test - - if (is_prime == 0) { // min_prime is composite - mpz_nextprime(tmp_prime, tmp_prime); - curr_prime = mpz_get_ui(tmp_prime); - } - - primes_.clear(); - while (curr_prime <= static_cast(maximum)) { - primes_.push_back(curr_prime); - mpz_nextprime(tmp_prime, tmp_prime); - curr_prime = mpz_get_ui(tmp_prime); - } - mpz_clear(tmp_prime); - - if (primes_.empty()) - throw std::invalid_argument("The given interval does not contain a prime number."); - - productOfAllCharacteristics_ = 1; - for (const unsigned int p : primes_) { - productOfAllCharacteristics_ *= p; - } - - partials_.resize(primes_.size()); - for (unsigned int i = 0; i < primes_.size(); ++i) { - unsigned int p = primes_[i]; - partials_[i] = productOfAllCharacteristics_ / p; - mpz_powm_ui(partials_[i].get_mpz_t(), partials_[i].get_mpz_t(), p - 1, - productOfAllCharacteristics_.get_mpz_t()); - } - - //If I understood the paper well, multiplicativeID_ always equals to 1. But in Clement's code, - //multiplicativeID_ is computed (see commented loop below). TODO: verify with Clement. - // for (unsigned int i = 0; i < partials_.size(); ++i){ - // multiplicativeID_ = (multiplicativeID_ + partials_[i]) % productOfAllCharacteristics_; - // } - } - characteristic_type get_characteristic() const{ - return productOfAllCharacteristics_; - } - - //TODO: verify if mpz_class is trivial to copy around or not. - //If yes, remove the unnecessery references? - //If not, do the operations in-place instead? - - element_type get_value(element_type e) const{ - if (e >= productOfAllCharacteristics_ || e < -productOfAllCharacteristics_) - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - if (e < 0) - mpz_add(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e; - } - - //r = e1 + e2 - element_type add(element_type e1, const element_type& e2) const{ - mpz_add(e1.get_mpz_t(), e1.get_mpz_t(), e2.get_mpz_t()); - return _get_value(e1); - } - - //r = e1 - e2 - element_type substract(element_type e1, const element_type& e2) const{ - mpz_sub(e1.get_mpz_t(), e1.get_mpz_t(), e2.get_mpz_t()); - return _get_value(e1); - } - - //r = e1 * e2 - element_type multiply(element_type e1, const element_type& e2) const{ - mpz_mul(e1.get_mpz_t(), e1.get_mpz_t(), e2.get_mpz_t()); - return _get_value(e1); - } - - //r = e * m + a - element_type multiply_and_add(element_type e, const element_type& m, const element_type& a) const{ - mpz_mul(e.get_mpz_t(), e.get_mpz_t(), m.get_mpz_t()); - mpz_add(e.get_mpz_t(), e.get_mpz_t(), a.get_mpz_t()); - return _get_value(e); - } - - //r = (e + a) * m - element_type add_and_multiply(element_type e, const element_type& a, const element_type& m) const{ - mpz_add(e.get_mpz_t(), e.get_mpz_t(), a.get_mpz_t()); - mpz_mul(e.get_mpz_t(), e.get_mpz_t(), m.get_mpz_t()); - return _get_value(e); - } - - bool are_equal(element_type e1, element_type e2) const{ - return _get_value(e1) == _get_value(e2); - } - - element_type get_inverse(const element_type& e) const{ - return get_partial_inverse(e, productOfAllCharacteristics_).first; - } - std::pair get_partial_inverse(const element_type& e, const characteristic_type& productOfCharacteristics) const{ - characteristic_type QR; - mpz_gcd(QR.get_mpz_t(), e.get_mpz_t(), productOfCharacteristics.get_mpz_t()); // QR <- gcd(x,QS) - - if (QR == productOfCharacteristics) - return {0, get_multiplicative_identity()}; // partial inverse is 0 - - characteristic_type QT = productOfCharacteristics / QR; - - characteristic_type inv_qt; - mpz_invert(inv_qt.get_mpz_t(), e.get_mpz_t(), QT.get_mpz_t()); - - auto res = get_partial_multiplicative_identity(QT); - mpz_mul(res.get_mpz_t(), res.get_mpz_t(), inv_qt.get_mpz_t()); - - return {_get_value(res), QT}; - } - - static element_type get_additive_identity(){ return 0; } - static element_type get_multiplicative_identity(){ return 1; } - // static element_type get_multiplicative_identity(){ return multiplicativeID_; } - element_type get_partial_multiplicative_identity(const characteristic_type& productOfCharacteristics) const - { - if (productOfCharacteristics == 0) { - return get_multiplicative_identity(); - } - element_type multIdentity(0); - for (unsigned int idx = 0; idx < primes_.size(); ++idx) { - if ((productOfCharacteristics % primes_[idx]) == 0) { - mpz_add(multIdentity.get_mpz_t(), multIdentity.get_mpz_t(), partials_[idx].get_mpz_t()); - } - } - return _get_value(multIdentity); - } - - static constexpr bool handles_only_z2(){ - return false; - } - - Multi_field_operators& operator=(Multi_field_operators other){ - primes_.swap(other.primes_); - productOfAllCharacteristics_ = other.productOfAllCharacteristics_; - partials_.swap(other.partials_); - - return *this; - } - - friend void swap(Multi_field_operators& f1, Multi_field_operators& f2){ - f1.primes_.swap(f2.primes_); - std::swap(f1.productOfAllCharacteristics_, f2.productOfAllCharacteristics_); - f1.partials_.swap(f2.partials_); - } - -private: - std::vector primes_; - characteristic_type productOfAllCharacteristics_; - std::vector partials_; - // static const element_type multiplicativeID_; - - static constexpr bool _is_prime(const int p){ - if (p <= 1) return false; - if (p <= 3) return true; - if (p % 2 == 0 || p % 3 == 0) return false; - - for (long i = 5; i * i <= p; i = i + 6) - if (p % i == 0 || p % (i + 2) == 0) - return false; - - return true; - } - - element_type& _get_value(element_type& e) const{ - if (e >= productOfAllCharacteristics_ || e < -productOfAllCharacteristics_) - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - if (e < 0) - mpz_add(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e; - } +class Multi_field_operators +{ + public: + using element_type = mpz_class; /**< Type for the elements in the field. */ + using characteristic_type = element_type; /**< Type for the field characteristic. */ + + /** + * @brief Default constructor, sets the product of all characteristics to 0. + */ + Multi_field_operators() : productOfAllCharacteristics_(0) /* , multiplicativeID_(1) */ + {} + /** + * @brief Constructor setting the characteristics to all prime numbers between the two given integers. + * + * @param minCharacteristic Smallest value of a prime. + * @param maxCharacteristic Heighest value of a prime. + */ + Multi_field_operators(int minCharacteristic, int maxCharacteristic) + : productOfAllCharacteristics_(0) //, multiplicativeID_(1) + { + set_characteristic(minCharacteristic, maxCharacteristic); + } + /** + * @brief Copy constructor. + * + * @param toCopy Operators to copy. + */ + Multi_field_operators(const Multi_field_operators& toCopy) + : primes_(toCopy.primes_), + productOfAllCharacteristics_(toCopy.productOfAllCharacteristics_), + partials_(toCopy.partials_) /* , + multiplicativeID_(toCopy.multiplicativeID_) */ + {} + /** + * @brief Move constructor. + * + * @param toMove Operators to move. + */ + Multi_field_operators(Multi_field_operators&& toMove) noexcept + : primes_(std::move(toMove.primes_)), + productOfAllCharacteristics_(std::move(toMove.productOfAllCharacteristics_)), + partials_(std::move(toMove.partials_)) /* , + multiplicativeID_(std::move(toMove.multiplicativeID_)) */ + {} + + /** + * @brief Set the characteristics of the field, which are stored in a single value as a product of all of them. + * The characteristics will be all prime numbers in the given interval. + * + * @param minimum Smallest value of a prime. + * @param maximum Heighest value of a prime. + */ + void set_characteristic(int minimum, int maximum) { + if (maximum < 2) throw std::invalid_argument("Characteristic must be strictly positive"); + if (minimum > maximum) throw std::invalid_argument("The given interval is not valid."); + if (minimum == maximum && !_is_prime(minimum)) + throw std::invalid_argument("The given interval does not contain a prime number."); + + unsigned int curr_prime = minimum; + mpz_t tmp_prime; + mpz_init_set_ui(tmp_prime, minimum); + // test if min_prime is prime + int is_prime = mpz_probab_prime_p(tmp_prime, 25); // probabilistic primality test + + if (is_prime == 0) { // min_prime is composite + mpz_nextprime(tmp_prime, tmp_prime); + curr_prime = mpz_get_ui(tmp_prime); + } + + primes_.clear(); + while (curr_prime <= static_cast(maximum)) { + primes_.push_back(curr_prime); + mpz_nextprime(tmp_prime, tmp_prime); + curr_prime = mpz_get_ui(tmp_prime); + } + mpz_clear(tmp_prime); + + if (primes_.empty()) throw std::invalid_argument("The given interval does not contain a prime number."); + + productOfAllCharacteristics_ = 1; + for (const unsigned int p : primes_) { + productOfAllCharacteristics_ *= p; + } + + partials_.resize(primes_.size()); + for (unsigned int i = 0; i < primes_.size(); ++i) { + unsigned int p = primes_[i]; + partials_[i] = productOfAllCharacteristics_ / p; + mpz_powm_ui(partials_[i].get_mpz_t(), partials_[i].get_mpz_t(), p - 1, productOfAllCharacteristics_.get_mpz_t()); + } + + // If I understood the paper well, multiplicativeID_ always equals to 1. But in Clement's code, + // multiplicativeID_ is computed (see commented loop below). TODO: verify with Clement. + // for (unsigned int i = 0; i < partials_.size(); ++i) { + // multiplicativeID_ = (multiplicativeID_ + partials_[i]) % productOfAllCharacteristics_; + // } + } + /** + * @brief Returns the current characteristics as the product of all of them. + * + * @return The value of the current characteristic. + */ + characteristic_type get_characteristic() const { return productOfAllCharacteristics_; } + + // TODO: verify if mpz_class is trivial to copy around or not. + // If yes, remove the unnecessery references? + // If not, do the operations in-place instead? + + /** + * @brief Returns the value of an element in the field. + * That is the positive value of the integer modulo the current characteristic. + * + * @param e Integer to return the value from. + * @return @p e modulo the current characteristic, such that the result is positive. + */ + element_type get_value(element_type e) const { + if (e >= productOfAllCharacteristics_ || e < -productOfAllCharacteristics_) + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + if (e < 0) mpz_add(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e; + } + + /** + * @brief Returns the sum of two elements in the field. + * + * @param e1 First element. + * @param e2 Second element. + * @return `(e1 + e2) % productOfAllCharacteristics`, such that the result is positive. + */ + element_type add(element_type e1, const element_type& e2) const { + mpz_add(e1.get_mpz_t(), e1.get_mpz_t(), e2.get_mpz_t()); + return _get_value(e1); + } + + /** + * @brief Returns the substraction in the field of the first element by the second element. + * + * @param e1 First element. + * @param e2 Second element. + * @return `(e1 - e2) % productOfAllCharacteristics`, such that the result is positive. + */ + element_type substract(element_type e1, const element_type& e2) const { + mpz_sub(e1.get_mpz_t(), e1.get_mpz_t(), e2.get_mpz_t()); + return _get_value(e1); + } + + /** + * @brief Returns the multiplication of two elements in the field. + * + * @param e1 First element. + * @param e2 Second element. + * @return `(e1 * e2) % productOfAllCharacteristics`, such that the result is positive. + */ + element_type multiply(element_type e1, const element_type& e2) const { + mpz_mul(e1.get_mpz_t(), e1.get_mpz_t(), e2.get_mpz_t()); + return _get_value(e1); + } + + /** + * @brief Multiplies the first element with the second one and adds the third one. Returns the result in the field. + * + * @param e First element. + * @param m Second element. + * @param a Third element. + * @return `(e * m + a) % productOfAllCharacteristics`, such that the result is positive. + */ + element_type multiply_and_add(element_type e, const element_type& m, const element_type& a) const { + mpz_mul(e.get_mpz_t(), e.get_mpz_t(), m.get_mpz_t()); + mpz_add(e.get_mpz_t(), e.get_mpz_t(), a.get_mpz_t()); + return _get_value(e); + } + + /** + * @brief Adds the first element to the second one and multiplies the third one with it. + * Returns the result in the field. + * + * @param e First element. + * @param a Second element. + * @param m Third element. + * @return `((e + a) * m) % productOfAllCharacteristics`, such that the result is positive. + */ + element_type add_and_multiply(element_type e, const element_type& a, const element_type& m) const { + mpz_add(e.get_mpz_t(), e.get_mpz_t(), a.get_mpz_t()); + mpz_mul(e.get_mpz_t(), e.get_mpz_t(), m.get_mpz_t()); + return _get_value(e); + } + + /** + * @brief Returns true if the two given elements are equal in the field, false otherwise. + * + * @param e1 First element to compare. + * @param e2 Second element to compare. + * @return true If `e1 % productOfAllCharacteristics == e2 % productOfAllCharacteristics`. + * @return false Otherwise. + */ + bool are_equal(element_type e1, element_type e2) const { return _get_value(e1) == _get_value(e2); } + + /** + * @brief Returns the inverse of the given element in the sense of @cite boissonnat:hal-00922572 with respect + * to the product of all characteristics. + * + * @param e Element to get the inverse from. + * @return Inverse in the current field. + */ + element_type get_inverse(const element_type& e) const { + return get_partial_inverse(e, productOfAllCharacteristics_).first; + } + /** + * @brief Returns the inverse of the given element in the multi-field corresponding to the given sub-product + * of the product of all characteristics in the multi-field. See @cite boissonnat:hal-00922572 for more details. + * + * @param e Element to get the inverse from. + * @param productOfCharacteristics Product of the different characteristics to take into account in the multi-field. + * @return Pair of the inverse of @p e and the characteristic the inverse is coming from. + */ + std::pair get_partial_inverse( + const element_type& e, const characteristic_type& productOfCharacteristics) const { + characteristic_type QR; + mpz_gcd(QR.get_mpz_t(), e.get_mpz_t(), productOfCharacteristics.get_mpz_t()); // QR <- gcd(x,QS) + + if (QR == productOfCharacteristics) return {0, get_multiplicative_identity()}; // partial inverse is 0 + + characteristic_type QT = productOfCharacteristics / QR; + + characteristic_type inv_qt; + mpz_invert(inv_qt.get_mpz_t(), e.get_mpz_t(), QT.get_mpz_t()); + + auto res = get_partial_multiplicative_identity(QT); + mpz_mul(res.get_mpz_t(), res.get_mpz_t(), inv_qt.get_mpz_t()); + + return {_get_value(res), QT}; + } + + /** + * @brief Returns the additive identity of a field. + * + * @return The additive identity of a field. + */ + static element_type get_additive_identity() { return 0; } + /** + * @brief Returns the multiplicative identity of a field. + * + * @return The multiplicative identity of a field. + */ + static element_type get_multiplicative_identity() { return 1; } + // static element_type get_multiplicative_identity(){ return multiplicativeID_; } + /** + * @brief Returns the partial multiplicative identity of the multi-field from the given product. + * See @cite boissonnat:hal-00922572 for more details. + * + * @param productOfCharacteristics Product of the different characteristics to take into account in the multi-field. + * @return The partial multiplicative identity of the multi-field. + */ + element_type get_partial_multiplicative_identity(const characteristic_type& productOfCharacteristics) const { + if (productOfCharacteristics == 0) { + return get_multiplicative_identity(); + } + element_type multIdentity(0); + for (unsigned int idx = 0; idx < primes_.size(); ++idx) { + if ((productOfCharacteristics % primes_[idx]) == 0) { + mpz_add(multIdentity.get_mpz_t(), multIdentity.get_mpz_t(), partials_[idx].get_mpz_t()); + } + } + return _get_value(multIdentity); + } + + // static constexpr bool handles_only_z2() { return false; } + + /** + * @brief Assign operator. + */ + Multi_field_operators& operator=(Multi_field_operators other) { + primes_.swap(other.primes_); + productOfAllCharacteristics_ = other.productOfAllCharacteristics_; + partials_.swap(other.partials_); + + return *this; + } + /** + * @brief Swap operator. + */ + friend void swap(Multi_field_operators& f1, Multi_field_operators& f2) { + f1.primes_.swap(f2.primes_); + std::swap(f1.productOfAllCharacteristics_, f2.productOfAllCharacteristics_); + f1.partials_.swap(f2.partials_); + } + + private: + std::vector primes_; /**< All characteristics. */ + characteristic_type productOfAllCharacteristics_; /**< Product of all characteristics. */ + std::vector partials_; /**< Partial products of the characteristics. */ + // static const element_type multiplicativeID_; + + static constexpr bool _is_prime(const int p) { + if (p <= 1) return false; + if (p <= 3) return true; + if (p % 2 == 0 || p % 3 == 0) return false; + + for (long i = 5; i * i <= p; i = i + 6) + if (p % i == 0 || p % (i + 2) == 0) return false; + + return true; + } + + element_type& _get_value(element_type& e) const { + if (e >= productOfAllCharacteristics_ || e < -productOfAllCharacteristics_) + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + if (e < 0) mpz_add(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e; + } }; -} //namespace persistence_fields -} //namespace Gudhi +} // namespace persistence_fields +} // namespace Gudhi #endif // MATRIX_FIELD_MULTI_OPERATORS_H_ diff --git a/src/Persistence_matrix/include/gudhi/Fields/Multi_field_shared.h b/src/Persistence_matrix/include/gudhi/Fields/Multi_field_shared.h index 7071a419fd..8a6c772cde 100644 --- a/src/Persistence_matrix/include/gudhi/Fields/Multi_field_shared.h +++ b/src/Persistence_matrix/include/gudhi/Fields/Multi_field_shared.h @@ -1,13 +1,19 @@ /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. - * Author(s): Hannah Schreiber + * Author(s): Hannah Schreiber, Clément Maria * - * Copyright (C) 2022 Inria + * Copyright (C) 2022-24 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification */ +/** + * @file Multi_field_shared.h + * @author Hannah Schreiber + * @brief Contains the @ref Shared_multi_field_element class. + */ + #ifndef MATRIX_FIELD_MULTI_SHARED_H_ #define MATRIX_FIELD_MULTI_SHARED_H_ @@ -20,374 +26,543 @@ namespace Gudhi { namespace persistence_fields { /** + * @class Shared_multi_field_element Multi_field_shared.h gudhi/Fields/Multi_field_shared.h * @ingroup persistence_fields * - * @brief TODO: - * + * @brief Class representing an element of a multi-field. If each instanciation of the class can represent another + * element, they all share the same characteritics. That is if the characteristics are set for one, they will be + * set for all the others. The characteristics can be set before instianciating the elements with the static + * @ref Shared_multi_field_element::initialize method. */ -class Shared_multi_field_element { -public: - using element_type = mpz_class; - - Shared_multi_field_element(); - Shared_multi_field_element(mpz_class element); - Shared_multi_field_element(const Shared_multi_field_element& toCopy); - Shared_multi_field_element(Shared_multi_field_element&& toMove) noexcept; - - static void initialize(unsigned int minimum, unsigned int maximum); - - friend void operator+=(Shared_multi_field_element& f1, Shared_multi_field_element const& f2){ - mpz_add(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), f2.element_.get_mpz_t()); - mpz_mod(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - } - friend Shared_multi_field_element operator+(Shared_multi_field_element f1, Shared_multi_field_element const& f2){ - f1 += f2; - return f1; - } - friend void operator+=(Shared_multi_field_element& f, mpz_class const v){ - mpz_add(f.element_.get_mpz_t(), f.element_.get_mpz_t(), v.get_mpz_t()); - mpz_mod(f.element_.get_mpz_t(), f.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - } - friend Shared_multi_field_element operator+(Shared_multi_field_element f, mpz_class const v){ - f += v; - return f; - } - friend mpz_class operator+(mpz_class v, Shared_multi_field_element const& f){ - mpz_class e(v); - mpz_add(e.get_mpz_t(), e.get_mpz_t(), f.element_.get_mpz_t()); - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e; - } - - friend void operator-=(Shared_multi_field_element& f1, Shared_multi_field_element const& f2){ - mpz_sub(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), f2.element_.get_mpz_t()); - mpz_mod(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - } - friend Shared_multi_field_element operator-(Shared_multi_field_element f1, Shared_multi_field_element const& f2){ - f1 -= f2; - return f1; - } - friend void operator-=(Shared_multi_field_element& f, mpz_class const v){ - mpz_sub(f.element_.get_mpz_t(), f.element_.get_mpz_t(), v.get_mpz_t()); - mpz_mod(f.element_.get_mpz_t(), f.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - } - friend Shared_multi_field_element operator-(Shared_multi_field_element f, mpz_class const v){ - f -= v; - return f; - } - friend mpz_class operator-(mpz_class v, Shared_multi_field_element const& f){ - mpz_class e(v); - if (e >= productOfAllCharacteristics_) - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - if (f.element_ > e) - mpz_add(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - mpz_sub(e.get_mpz_t(), e.get_mpz_t(), f.element_.get_mpz_t()); - return e; - } - - friend void operator*=(Shared_multi_field_element& f1, Shared_multi_field_element const& f2){ - mpz_mul(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), f2.element_.get_mpz_t()); - mpz_mod(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - } - friend Shared_multi_field_element operator*(Shared_multi_field_element f1, Shared_multi_field_element const& f2){ - f1 *= f2; - return f1; - } - friend void operator*=(Shared_multi_field_element& f, mpz_class const v){ - mpz_mul(f.element_.get_mpz_t(), f.element_.get_mpz_t(), v.get_mpz_t()); - mpz_mod(f.element_.get_mpz_t(), f.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - } - friend Shared_multi_field_element operator*(Shared_multi_field_element f, mpz_class const v){ - f *= v; - return f; - } - friend mpz_class operator*(mpz_class v, Shared_multi_field_element const& f){ - mpz_class e(v); - mpz_mul(e.get_mpz_t(), e.get_mpz_t(), f.element_.get_mpz_t()); - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e; - } - - friend bool operator==(const Shared_multi_field_element& f1, const Shared_multi_field_element& f2){ - return f1.element_ == f2.element_; - } - friend bool operator==(const mpz_class v, const Shared_multi_field_element& f){ - if (v < productOfAllCharacteristics_) return v == f.element_; - mpz_class e(v); - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e == f.element_; - } - friend bool operator==(const Shared_multi_field_element& f, const mpz_class v){ - if (v < productOfAllCharacteristics_) return v == f.element_; - mpz_class e(v); - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e == f.element_; - } - friend bool operator==(const unsigned int v, const Shared_multi_field_element& f){ - mpz_class e(v); - if (e < productOfAllCharacteristics_) return e == f.element_; - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e == f.element_; - } - friend bool operator==(const Shared_multi_field_element& f, const unsigned int v){ - mpz_class e(v); - if (e < productOfAllCharacteristics_) return e == f.element_; - mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return e == f.element_; - } - friend bool operator!=(const Shared_multi_field_element& f1, const Shared_multi_field_element& f2){ - return !(f1 == f2); - } - friend bool operator!=(const mpz_class v, const Shared_multi_field_element& f){ - return !(v == f); - } - friend bool operator!=(const Shared_multi_field_element& f, const mpz_class v){ - return !(v == f); - } - friend bool operator!=(const unsigned int v, const Shared_multi_field_element& f){ - return !(v == f); - } - friend bool operator!=(const Shared_multi_field_element& f, const unsigned int v){ - return !(v == f); - } - - Shared_multi_field_element& operator=(Shared_multi_field_element other); - Shared_multi_field_element& operator=(const mpz_class value); - operator unsigned int() const; - operator mpz_class() const; - friend void swap(Shared_multi_field_element& f1, Shared_multi_field_element& f2){ - std::swap(f1.element_, f2.element_); - } - - Shared_multi_field_element get_inverse() const; - std::pair get_partial_inverse(const mpz_class& productOfCharacteristics) const; - - static Shared_multi_field_element get_additive_identity(); - static Shared_multi_field_element get_multiplicative_identity(); - static Shared_multi_field_element get_partial_multiplicative_identity(const mpz_class& productOfCharacteristics); - static mpz_class get_characteristic(); - - mpz_class get_value() const; - - static constexpr bool handles_only_z2(){ - return false; - } - -private: - mpz_class element_; - static inline std::vector primes_; - static inline mpz_class productOfAllCharacteristics_ = 0; - static inline std::vector partials_; - static inline const mpz_class multiplicativeID_ = 1; - - static constexpr bool _is_prime(const int p); +class Shared_multi_field_element +{ + public: + using element_type = mpz_class; /**< Type for the elements in the field. */ + + /** + * @brief Default constructor. Sets the element to 0. + */ + Shared_multi_field_element(); + /** + * @brief Constructor setting the element to the given value. + * + * @param element Value of the element. + */ + Shared_multi_field_element(mpz_class element); + /** + * @brief Copy constructor. + * + * @param toCopy Element to copy. + */ + Shared_multi_field_element(const Shared_multi_field_element& toCopy); + /** + * @brief Move constructor. + * + * @param toMove Element to move. + */ + Shared_multi_field_element(Shared_multi_field_element&& toMove) noexcept; + + /** + * @brief Initialize the multi-field to the characteristics (primes) contained in the given interval. + * Should be called first before constructing the field elements. + * + * @param minimum Lowest value in the interval. + * @param maximum Highest value in the interval. + */ + static void initialize(unsigned int minimum, unsigned int maximum); + + /** + * @brief operator+= + * + * @param f1 Shared_multi_field_element + * @param f2 Shared_multi_field_element + */ + friend void operator+=(Shared_multi_field_element& f1, Shared_multi_field_element const& f2) { + mpz_add(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), f2.element_.get_mpz_t()); + mpz_mod(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + } + /** + * @brief operator+ + * + * @param f1 Shared_multi_field_element + * @param f2 const Shared_multi_field_element & + * @return Shared_multi_field_element + */ + friend Shared_multi_field_element operator+(Shared_multi_field_element f1, Shared_multi_field_element const& f2) { + f1 += f2; + return f1; + } + /** + * @brief operator+= + * + * @param f Shared_multi_field_element& + * @param v const mpz_class + */ + friend void operator+=(Shared_multi_field_element& f, mpz_class const v) { + mpz_add(f.element_.get_mpz_t(), f.element_.get_mpz_t(), v.get_mpz_t()); + mpz_mod(f.element_.get_mpz_t(), f.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + } + /** + * @brief operator+ + * + * @param f Shared_multi_field_element + * @param v const mpz_class + * @return Shared_multi_field_element + */ + friend Shared_multi_field_element operator+(Shared_multi_field_element f, mpz_class const v) { + f += v; + return f; + } + /** + * @brief operator+ + * + * @param v mpz_class + * @param f const Shared_multi_field_element & + * @return mpz_class + */ + friend mpz_class operator+(mpz_class v, Shared_multi_field_element const& f) { + mpz_class e(v); + mpz_add(e.get_mpz_t(), e.get_mpz_t(), f.element_.get_mpz_t()); + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e; + } + + /** + * @brief operator-= + * + * @param f1 Shared_multi_field_element& + * @param f2 const Shared_multi_field_element & + */ + friend void operator-=(Shared_multi_field_element& f1, Shared_multi_field_element const& f2) { + mpz_sub(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), f2.element_.get_mpz_t()); + mpz_mod(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + } + /** + * @brief operator- + * + * @param f1 Shared_multi_field_element + * @param f2 const Shared_multi_field_element & + * @return Shared_multi_field_element + */ + friend Shared_multi_field_element operator-(Shared_multi_field_element f1, Shared_multi_field_element const& f2) { + f1 -= f2; + return f1; + } + /** + * @brief operator-= + * + * @param f Shared_multi_field_element& + * @param v const mpz_class + */ + friend void operator-=(Shared_multi_field_element& f, mpz_class const v) { + mpz_sub(f.element_.get_mpz_t(), f.element_.get_mpz_t(), v.get_mpz_t()); + mpz_mod(f.element_.get_mpz_t(), f.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + } + /** + * @brief operator- + * + * @param f Shared_multi_field_element + * @param v const mpz_class + * @return Shared_multi_field_element + */ + friend Shared_multi_field_element operator-(Shared_multi_field_element f, mpz_class const v) { + f -= v; + return f; + } + /** + * @brief operator- + * + * @param v mpz_class + * @param f const Shared_multi_field_element & + * @return mpz_class + */ + friend mpz_class operator-(mpz_class v, Shared_multi_field_element const& f) { + mpz_class e(v); + if (e >= productOfAllCharacteristics_) + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + if (f.element_ > e) mpz_add(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + mpz_sub(e.get_mpz_t(), e.get_mpz_t(), f.element_.get_mpz_t()); + return e; + } + + /** + * @brief operator*= + * + * @param f1 Shared_multi_field_element& + * @param f2 const Shared_multi_field_element & + */ + friend void operator*=(Shared_multi_field_element& f1, Shared_multi_field_element const& f2) { + mpz_mul(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), f2.element_.get_mpz_t()); + mpz_mod(f1.element_.get_mpz_t(), f1.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + } + /** + * @brief operator* + * + * @param f1 Shared_multi_field_element + * @param f2 const Shared_multi_field_element & + * @return Shared_multi_field_element + */ + friend Shared_multi_field_element operator*(Shared_multi_field_element f1, Shared_multi_field_element const& f2) { + f1 *= f2; + return f1; + } + /** + * @brief operator*= + * + * @param f Shared_multi_field_element& + * @param v const mpz_class + */ + friend void operator*=(Shared_multi_field_element& f, mpz_class const v) { + mpz_mul(f.element_.get_mpz_t(), f.element_.get_mpz_t(), v.get_mpz_t()); + mpz_mod(f.element_.get_mpz_t(), f.element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + } + /** + * @brief operator* + * + * @param f Shared_multi_field_element + * @param v const mpz_class + * @return Shared_multi_field_element + */ + friend Shared_multi_field_element operator*(Shared_multi_field_element f, mpz_class const v) { + f *= v; + return f; + } + /** + * @brief operator* + * + * @param v mpz_class + * @param f const Shared_multi_field_element & + * @return mpz_class + */ + friend mpz_class operator*(mpz_class v, Shared_multi_field_element const& f) { + mpz_class e(v); + mpz_mul(e.get_mpz_t(), e.get_mpz_t(), f.element_.get_mpz_t()); + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e; + } + + /** + * @brief operator== + * + * @param f1 const Shared_multi_field_element& + * @param f2 const Shared_multi_field_element& + * @return bool + */ + friend bool operator==(const Shared_multi_field_element& f1, const Shared_multi_field_element& f2) { + return f1.element_ == f2.element_; + } + /** + * @brief operator== + * + * @param f1 const Shared_multi_field_element& + * @param f2 const Shared_multi_field_element& + * @return bool + */ + friend bool operator==(const mpz_class v, const Shared_multi_field_element& f) { + if (v < productOfAllCharacteristics_) return v == f.element_; + mpz_class e(v); + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e == f.element_; + } + /** + * @brief operator== + * + * @param f const Shared_multi_field_element& + * @param v const mpz_class + * @return bool + */ + friend bool operator==(const Shared_multi_field_element& f, const mpz_class v) { + if (v < productOfAllCharacteristics_) return v == f.element_; + mpz_class e(v); + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e == f.element_; + } + /** + * @brief operator== + * + * @param v const unsigned int + * @param f const Shared_multi_field_element& + * @return bool + */ + friend bool operator==(const unsigned int v, const Shared_multi_field_element& f) { + mpz_class e(v); + if (e < productOfAllCharacteristics_) return e == f.element_; + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e == f.element_; + } + /** + * @brief operator== + * + * @param f const Shared_multi_field_element& + * @param v const unsigned int + * @return bool + */ + friend bool operator==(const Shared_multi_field_element& f, const unsigned int v) { + mpz_class e(v); + if (e < productOfAllCharacteristics_) return e == f.element_; + mpz_mod(e.get_mpz_t(), e.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return e == f.element_; + } + /** + * @brief operator!= + * + * @param f1 const Shared_multi_field_element& + * @param f2 const Shared_multi_field_element& + * @return bool + */ + friend bool operator!=(const Shared_multi_field_element& f1, const Shared_multi_field_element& f2) { + return !(f1 == f2); + } + /** + * @brief operator!= + * + * @param f1 const Shared_multi_field_element& + * @param f2 const Shared_multi_field_element& + * @return bool + */ + friend bool operator!=(const mpz_class v, const Shared_multi_field_element& f) { return !(v == f); } + /** + * @brief operator!= + * + * @param f const Shared_multi_field_element& + * @param v const mpz_class + * @return bool + */ + friend bool operator!=(const Shared_multi_field_element& f, const mpz_class v) { return !(v == f); } + /** + * @brief operator!= + * + * @param v const unsigned int + * @param f const Shared_multi_field_element& + * @return bool + */ + friend bool operator!=(const unsigned int v, const Shared_multi_field_element& f) { return !(v == f); } + /** + * @brief operator!= + * + * @param f const Shared_multi_field_element& + * @param v const unsigned int + * @return bool + */ + friend bool operator!=(const Shared_multi_field_element& f, const unsigned int v) { return !(v == f); } + + /** + * @brief Assign operator. + */ + Shared_multi_field_element& operator=(Shared_multi_field_element other); + /** + * @brief Assign operator. + */ + Shared_multi_field_element& operator=(const mpz_class value); + /** + * @brief Swap operator. + */ + friend void swap(Shared_multi_field_element& f1, Shared_multi_field_element& f2) { + std::swap(f1.element_, f2.element_); + } + + /** + * @brief Casts the element into an unsigned int. + */ + operator unsigned int() const; + /** + * @brief Casts the element into a mpz_class. + */ + operator mpz_class() const; + + /** + * @brief Returns the inverse of the element in the multi-field, see @cite boissonnat:hal-00922572. + * + * @return The inverse. + */ + Shared_multi_field_element get_inverse() const; + /** + * @brief Returns the inverse of the element with respect to a sub-product of the characteristics in the multi-field, + * see @cite boissonnat:hal-00922572. + * + * @param productOfCharacteristics Sub-product of the characteristics. + * @return Pair of the inverse and the characteristic the inverse corresponds to. + */ + std::pair get_partial_inverse(const mpz_class& productOfCharacteristics) const; + + /** + * @brief Returns the additive identity of a field. + * + * @return The additive identity of a field. + */ + static Shared_multi_field_element get_additive_identity(); + /** + * @brief Returns the multiplicative identity of a field. + * + * @return The multiplicative identity of a field. + */ + static Shared_multi_field_element get_multiplicative_identity(); + /** + * @brief Returns the partial multiplicative identity of the multi-field from the given product. + * See @cite boissonnat:hal-00922572 for more details. + * + * @param productOfCharacteristics Product of the different characteristics to take into account in the multi-field. + * @return The partial multiplicative identity of the multi-field. + */ + static Shared_multi_field_element get_partial_multiplicative_identity(const mpz_class& productOfCharacteristics); + /** + * @brief Returns the product of all characteristics. + * + * @return The product of all characteristics. + */ + static mpz_class get_characteristic(); + + /** + * @brief Returns the value of the element. + * + * @return Value of the element. + */ + mpz_class get_value() const; + + // static constexpr bool handles_only_z2() { return false; } + + private: + mpz_class element_; /**< Element. */ + static inline std::vector primes_; /**< All characteristics. */ + static inline mpz_class productOfAllCharacteristics_ = 0; /**< Product of all characteristics. */ + static inline std::vector partials_; /**< Partial products of the characteristics. */ + static inline const mpz_class multiplicativeID_ = 1; /**< Multiplicative identity. */ + + static constexpr bool _is_prime(const int p); }; -inline Shared_multi_field_element::Shared_multi_field_element() - : element_(0) -{} +inline Shared_multi_field_element::Shared_multi_field_element() : element_(0) {} -inline Shared_multi_field_element::Shared_multi_field_element(mpz_class element) - : element_(element) -{ - mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); +inline Shared_multi_field_element::Shared_multi_field_element(mpz_class element) : element_(element) { + mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); } -inline Shared_multi_field_element::Shared_multi_field_element(const Shared_multi_field_element &toCopy) - : element_(toCopy.element_) -{} - -inline Shared_multi_field_element::Shared_multi_field_element(Shared_multi_field_element &&toMove) noexcept - : element_(std::move(toMove.element_)) -{} - -inline void Shared_multi_field_element::initialize(unsigned int minimum, unsigned int maximum) -{ - if (maximum < 2) - throw std::invalid_argument("Characteristic must be strictly positive"); - if (minimum > maximum) - throw std::invalid_argument("The given interval is not valid."); - if (minimum == maximum && !_is_prime(minimum)) - throw std::invalid_argument("The given interval does not contain a prime number."); - - unsigned int curr_prime = minimum; - mpz_t tmp_prime; - mpz_init_set_ui(tmp_prime, minimum); - // test if min_prime is prime - int is_prime = mpz_probab_prime_p(tmp_prime, 25); // probabilistic primality test - - if (is_prime == 0) { // min_prime is composite - mpz_nextprime(tmp_prime, tmp_prime); - curr_prime = mpz_get_ui(tmp_prime); - } - - primes_.clear(); - while (curr_prime <= maximum) { - primes_.push_back(curr_prime); - mpz_nextprime(tmp_prime, tmp_prime); - curr_prime = mpz_get_ui(tmp_prime); - } - mpz_clear(tmp_prime); - - if (primes_.empty()) - throw std::invalid_argument("The given interval does not contain a prime number."); - - productOfAllCharacteristics_ = 1; - for (const unsigned int p : primes_) { - productOfAllCharacteristics_ *= p; - } - - partials_.resize(primes_.size()); - for (unsigned int i = 0; i < primes_.size(); ++i) { - unsigned int p = primes_[i]; - partials_[i] = productOfAllCharacteristics_ / p; - mpz_powm_ui(partials_[i].get_mpz_t(), partials_[i].get_mpz_t(), p - 1, - productOfAllCharacteristics_.get_mpz_t()); - } - - //If I understood the paper well, multiplicativeID_ always equals to 1. But in Clement's code, - //multiplicativeID_ is computed (see commented loop below). TODO: verify with Clement. -// for (unsigned int i = 0; i < partials_.size(); ++i){ -// multiplicativeID_ = (multiplicativeID_ + partials_[i]) % productOfAllCharacteristics_; -// } +inline Shared_multi_field_element::Shared_multi_field_element(const Shared_multi_field_element& toCopy) + : element_(toCopy.element_) {} + +inline Shared_multi_field_element::Shared_multi_field_element(Shared_multi_field_element&& toMove) noexcept + : element_(std::move(toMove.element_)) {} + +inline void Shared_multi_field_element::initialize(unsigned int minimum, unsigned int maximum) { + if (maximum < 2) throw std::invalid_argument("Characteristic must be strictly positive"); + if (minimum > maximum) throw std::invalid_argument("The given interval is not valid."); + if (minimum == maximum && !_is_prime(minimum)) + throw std::invalid_argument("The given interval does not contain a prime number."); + + unsigned int curr_prime = minimum; + mpz_t tmp_prime; + mpz_init_set_ui(tmp_prime, minimum); + // test if min_prime is prime + int is_prime = mpz_probab_prime_p(tmp_prime, 25); // probabilistic primality test + + if (is_prime == 0) { // min_prime is composite + mpz_nextprime(tmp_prime, tmp_prime); + curr_prime = mpz_get_ui(tmp_prime); + } + + primes_.clear(); + while (curr_prime <= maximum) { + primes_.push_back(curr_prime); + mpz_nextprime(tmp_prime, tmp_prime); + curr_prime = mpz_get_ui(tmp_prime); + } + mpz_clear(tmp_prime); + + if (primes_.empty()) throw std::invalid_argument("The given interval does not contain a prime number."); + + productOfAllCharacteristics_ = 1; + for (const unsigned int p : primes_) { + productOfAllCharacteristics_ *= p; + } + + partials_.resize(primes_.size()); + for (unsigned int i = 0; i < primes_.size(); ++i) { + unsigned int p = primes_[i]; + partials_[i] = productOfAllCharacteristics_ / p; + mpz_powm_ui(partials_[i].get_mpz_t(), partials_[i].get_mpz_t(), p - 1, productOfAllCharacteristics_.get_mpz_t()); + } + + // If I understood the paper well, multiplicativeID_ always equals to 1. But in Clement's code, + // multiplicativeID_ is computed (see commented loop below). TODO: verify with Clement. + // for (unsigned int i = 0; i < partials_.size(); ++i){ + // multiplicativeID_ = (multiplicativeID_ + partials_[i]) % productOfAllCharacteristics_; + // } } -//inline Shared_multi_field_element &Shared_multi_field_element::operator+=(Shared_multi_field_element const &f) -//{ -// mpz_add(element_.get_mpz_t(), element_.get_mpz_t(), f.element_.get_mpz_t()); -// mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); -// return *this; -//} - -//inline Shared_multi_field_element &Shared_multi_field_element::operator+=(mpz_class const v) -//{ -// mpz_add(element_.get_mpz_t(), element_.get_mpz_t(), v.get_mpz_t()); -// mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); -// return *this; -//} - -//inline Shared_multi_field_element &Shared_multi_field_element::operator-=(Shared_multi_field_element const &f) -//{ -// mpz_sub(element_.get_mpz_t(), element_.get_mpz_t(), f.element_.get_mpz_t()); -// mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); -// return *this; -//} - -//inline Shared_multi_field_element &Shared_multi_field_element::operator-=(mpz_class const v) -//{ -// mpz_sub(element_.get_mpz_t(), element_.get_mpz_t(), v.get_mpz_t()); -// mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); -// return *this; -//} - -//inline Shared_multi_field_element &Shared_multi_field_element::operator*=(Shared_multi_field_element const &f) -//{ -// mpz_mul(element_.get_mpz_t(), element_.get_mpz_t(), f.element_.get_mpz_t()); -// mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); -// return *this; -//} - -//inline Shared_multi_field_element &Shared_multi_field_element::operator*=(mpz_class const v) -//{ -// mpz_mul(element_.get_mpz_t(), element_.get_mpz_t(), v.get_mpz_t()); -// mpz_mod(element_.get_mpz_t(), element_.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); -// return *this; -//} - -inline Shared_multi_field_element &Shared_multi_field_element::operator=(Shared_multi_field_element other) -{ - std::swap(element_, other.element_); - return *this; +inline Shared_multi_field_element& Shared_multi_field_element::operator=(Shared_multi_field_element other) { + std::swap(element_, other.element_); + return *this; } -inline Shared_multi_field_element &Shared_multi_field_element::operator=(mpz_class const value) -{ - mpz_mod(element_.get_mpz_t(), value.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); - return *this; +inline Shared_multi_field_element& Shared_multi_field_element::operator=(mpz_class const value) { + mpz_mod(element_.get_mpz_t(), value.get_mpz_t(), productOfAllCharacteristics_.get_mpz_t()); + return *this; } -inline Shared_multi_field_element::operator unsigned int() const -{ - return element_.get_ui(); -} +inline Shared_multi_field_element::operator unsigned int() const { return element_.get_ui(); } -inline Shared_multi_field_element::operator mpz_class() const -{ - return element_; -} +inline Shared_multi_field_element::operator mpz_class() const { return element_; } -inline Shared_multi_field_element Shared_multi_field_element::get_inverse() const -{ - return get_partial_inverse(productOfAllCharacteristics_).first; +inline Shared_multi_field_element Shared_multi_field_element::get_inverse() const { + return get_partial_inverse(productOfAllCharacteristics_).first; } -inline std::pair Shared_multi_field_element::get_partial_inverse(const mpz_class& productOfCharacteristics) const -{ - mpz_class QR; - mpz_gcd(QR.get_mpz_t(), element_.get_mpz_t(), productOfCharacteristics.get_mpz_t()); // QR <- gcd(x,QS) +inline std::pair Shared_multi_field_element::get_partial_inverse( + const mpz_class& productOfCharacteristics) const { + mpz_class QR; + mpz_gcd(QR.get_mpz_t(), element_.get_mpz_t(), productOfCharacteristics.get_mpz_t()); // QR <- gcd(x,QS) - if (QR == productOfCharacteristics) - return {Shared_multi_field_element(), multiplicativeID_}; // partial inverse is 0 + if (QR == productOfCharacteristics) return {Shared_multi_field_element(), multiplicativeID_}; // partial inverse is 0 - mpz_class QT = productOfCharacteristics / QR; + mpz_class QT = productOfCharacteristics / QR; - mpz_class inv_qt; - mpz_invert(inv_qt.get_mpz_t(), element_.get_mpz_t(), QT.get_mpz_t()); + mpz_class inv_qt; + mpz_invert(inv_qt.get_mpz_t(), element_.get_mpz_t(), QT.get_mpz_t()); - auto res = get_partial_multiplicative_identity(QT); - res *= inv_qt; + auto res = get_partial_multiplicative_identity(QT); + res *= inv_qt; - return {res, QT}; + return {res, QT}; } -inline Shared_multi_field_element Shared_multi_field_element::get_additive_identity() -{ - return Shared_multi_field_element(); +inline Shared_multi_field_element Shared_multi_field_element::get_additive_identity() { + return Shared_multi_field_element(); } -inline Shared_multi_field_element Shared_multi_field_element::get_multiplicative_identity() -{ - return Shared_multi_field_element(multiplicativeID_); +inline Shared_multi_field_element Shared_multi_field_element::get_multiplicative_identity() { + return Shared_multi_field_element(multiplicativeID_); } -inline Shared_multi_field_element Shared_multi_field_element::get_partial_multiplicative_identity(const mpz_class& productOfCharacteristics) -{ - if (productOfCharacteristics == 0) { - return Shared_multi_field_element(multiplicativeID_); - } - Shared_multi_field_element mult; - for (unsigned int idx = 0; idx < primes_.size(); ++idx) { - if ((productOfCharacteristics % primes_[idx]) == 0) { - mult += partials_[idx]; - } - } - return mult; +inline Shared_multi_field_element Shared_multi_field_element::get_partial_multiplicative_identity( + const mpz_class& productOfCharacteristics) { + if (productOfCharacteristics == 0) { + return Shared_multi_field_element(multiplicativeID_); + } + Shared_multi_field_element mult; + for (unsigned int idx = 0; idx < primes_.size(); ++idx) { + if ((productOfCharacteristics % primes_[idx]) == 0) { + mult += partials_[idx]; + } + } + return mult; } -inline mpz_class Shared_multi_field_element::get_characteristic() -{ - return productOfAllCharacteristics_; -} +inline mpz_class Shared_multi_field_element::get_characteristic() { return productOfAllCharacteristics_; } -inline mpz_class Shared_multi_field_element::get_value() const -{ - return element_; -} +inline mpz_class Shared_multi_field_element::get_value() const { return element_; } -inline constexpr bool Shared_multi_field_element::_is_prime(const int p) -{ - if (p <= 1) return false; - if (p <= 3) return true; - if (p % 2 == 0 || p % 3 == 0) return false; +inline constexpr bool Shared_multi_field_element::_is_prime(const int p) { + if (p <= 1) return false; + if (p <= 3) return true; + if (p % 2 == 0 || p % 3 == 0) return false; - for (long i = 5; i * i <= p; i = i + 6) - if (p % i == 0 || p % (i + 2) == 0) - return false; + for (long i = 5; i * i <= p; i = i + 6) + if (p % i == 0 || p % (i + 2) == 0) return false; - return true; + return true; } -} //namespace persistence_fields -} //namespace Gudhi +} // namespace persistence_fields +} // namespace Gudhi #endif // MATRIX_FIELD_MULTI_SHARED_H_ diff --git a/src/Persistence_matrix/include/gudhi/Fields/Multi_field_small.h b/src/Persistence_matrix/include/gudhi/Fields/Multi_field_small.h index 1cfa3a0089..8d15a2cc93 100644 --- a/src/Persistence_matrix/include/gudhi/Fields/Multi_field_small.h +++ b/src/Persistence_matrix/include/gudhi/Fields/Multi_field_small.h @@ -1,13 +1,19 @@ /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. - * Author(s): Hannah Schreiber + * Author(s): Hannah Schreiber, Clément Maria * - * Copyright (C) 2022 Inria + * Copyright (C) 2022-24 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification */ +/** + * @file Multi_field_small.h + * @author Hannah Schreiber + * @brief Contains the @ref Multi_field_element_with_small_characteristics class. + */ + #ifndef MATRIX_FIELD_MULTI_SMALL_H_ #define MATRIX_FIELD_MULTI_SMALL_H_ @@ -20,483 +26,595 @@ namespace Gudhi { namespace persistence_fields { -//productOfAllCharacteristics_ has to fit in an unsigned int, ie, productOfAllCharacteristics_ < 4294967295 /** + * @class Multi_field_element_with_small_characteristics Multi_field_small.h gudhi/Fields/Multi_field_small.h * @ingroup persistence_fields * - * @brief TODO: + * @brief Class representing an element of a multi-field, such that the product of all characteristics fits into + * an unsigned int. The characteristics will corresponds to all prime numbers in the interval given as template. * + * @tparam minimum Interval closed lower bound. + * @tparam maximum Interval closed upper bound. */ -template +template class Multi_field_element_with_small_characteristics { -public: - using element_type = unsigned int; - template - using isInteger = std::enable_if_t >; - - Multi_field_element_with_small_characteristics(); - Multi_field_element_with_small_characteristics(unsigned int element); - Multi_field_element_with_small_characteristics(int element); - Multi_field_element_with_small_characteristics(const Multi_field_element_with_small_characteristics& toCopy); - Multi_field_element_with_small_characteristics(Multi_field_element_with_small_characteristics&& toMove) noexcept; - - friend void operator+=(Multi_field_element_with_small_characteristics& f1, Multi_field_element_with_small_characteristics const& f2){ - f1.element_ = _add(f1.element_, f2.element_); - } - friend Multi_field_element_with_small_characteristics operator+(Multi_field_element_with_small_characteristics f1, Multi_field_element_with_small_characteristics const& f2){ - f1 += f2; - return f1; - } - friend void operator+=(Multi_field_element_with_small_characteristics& f, unsigned int const v){ - f.element_ = _add(f.element_, v < productOfAllCharacteristics_ ? v : v % productOfAllCharacteristics_); - } - //v is assumed to be positive and will be casted into an unsigned int - template > - friend Multi_field_element_with_small_characteristics operator+(Multi_field_element_with_small_characteristics f, const Integer_type v){ - f += v; - return f; - } - template > - friend Integer_type operator+(Integer_type v, Multi_field_element_with_small_characteristics const& f){ - v += f.element_; - v %= productOfAllCharacteristics_; - return v; - } - - friend void operator-=(Multi_field_element_with_small_characteristics& f1, Multi_field_element_with_small_characteristics const& f2){ - f1.element_ = _substract(f1.element_, f2.element_); - } - friend Multi_field_element_with_small_characteristics operator-(Multi_field_element_with_small_characteristics f1, Multi_field_element_with_small_characteristics const& f2){ - f1 -= f2; - return f1; - } - friend void operator-=(Multi_field_element_with_small_characteristics& f, unsigned int const v){ - f.element_ = _substract(f.element_, v < productOfAllCharacteristics_ ? v : v % productOfAllCharacteristics_); - } - //v is assumed to be positive and will be casted into an unsigned int - template > - friend Multi_field_element_with_small_characteristics operator-(Multi_field_element_with_small_characteristics f, const Integer_type v){ - f -= v; - return f; - } - //v is assumed to be positive - template > - friend Integer_type operator-(Integer_type v, Multi_field_element_with_small_characteristics const& f){ - if (v >= productOfAllCharacteristics_) v %= productOfAllCharacteristics_; - if (f.element_ > v) v += productOfAllCharacteristics_; - v -= f.element_; - return v; - } - - friend void operator*=(Multi_field_element_with_small_characteristics& f1, Multi_field_element_with_small_characteristics const& f2){ - f1.element_ = _multiply(f1.element_, f2.element_); - } - friend Multi_field_element_with_small_characteristics operator*(Multi_field_element_with_small_characteristics f1, Multi_field_element_with_small_characteristics const& f2){ - f1 *= f2; - return f1; - } - friend void operator*=(Multi_field_element_with_small_characteristics& f, unsigned int const v){ - f.element_ = _multiply(f.element_, v < productOfAllCharacteristics_ ? v : v % productOfAllCharacteristics_); - } - //v is assumed to be positive and will be casted into an unsigned int - template > - friend Multi_field_element_with_small_characteristics operator*(Multi_field_element_with_small_characteristics f, const Integer_type v){ - f *= v; - return f; - } - //uses bitwise operations on v, so be carefull with signed integers - template > - friend Integer_type operator*(Integer_type v, Multi_field_element_with_small_characteristics const& f){ - unsigned int b = f.element_; - unsigned int res = 0; - unsigned int temp_b; - - while (v != 0) { - if (v & 1) { - if (b >= productOfAllCharacteristics_ - res) - res -= productOfAllCharacteristics_; - res += b; - } - v >>= 1; - - temp_b = b; - if (b >= productOfAllCharacteristics_ - b) - temp_b -= productOfAllCharacteristics_; - b += temp_b; - } - - return res; - } - - friend bool operator==(const Multi_field_element_with_small_characteristics& f1, const Multi_field_element_with_small_characteristics& f2){ - return f1.element_ == f2.element_; - } - //v is assumed to be positive - template > - friend bool operator==(const Integer_type v, const Multi_field_element_with_small_characteristics& f){ - if (v < productOfAllCharacteristics_) return v == f.element_; - return (v % productOfAllCharacteristics_) == f.element_; - } - //v is assumed to be positive - template > - friend bool operator==(const Multi_field_element_with_small_characteristics& f, const Integer_type v){ - if (v < productOfAllCharacteristics_) return v == f.element_; - return (v % productOfAllCharacteristics_) == f.element_; - } - friend bool operator!=(const Multi_field_element_with_small_characteristics& f1, const Multi_field_element_with_small_characteristics& f2){ - return !(f1 == f2); - } - //v is assumed to be positive - template > - friend bool operator!=(const Integer_type v, const Multi_field_element_with_small_characteristics& f){ - return !(v == f); - } - //v is assumed to be positive - template > - friend bool operator!=(const Multi_field_element_with_small_characteristics& f, const Integer_type v){ - return !(v == f); - } - - Multi_field_element_with_small_characteristics& operator=(Multi_field_element_with_small_characteristics other); - Multi_field_element_with_small_characteristics& operator=(const unsigned int value); - operator unsigned int() const; - friend void swap(Multi_field_element_with_small_characteristics& f1, Multi_field_element_with_small_characteristics& f2){ - std::swap(f1.element_, f2.element_); - } - - Multi_field_element_with_small_characteristics get_inverse() const; - std::pair get_partial_inverse(unsigned int productOfCharacteristics) const; - - static Multi_field_element_with_small_characteristics get_additive_identity(); - static Multi_field_element_with_small_characteristics get_multiplicative_identity(); - static Multi_field_element_with_small_characteristics get_partial_multiplicative_identity(const mpz_class& productOfCharacteristics); - static constexpr unsigned int get_characteristic(); - - unsigned int get_value() const; - - static constexpr bool handles_only_z2(){ - return false; - } - -private: - static constexpr bool _is_prime(const int p); - static constexpr unsigned int _multiply(unsigned int a, unsigned int b); - static constexpr unsigned int _add(unsigned int element, unsigned int v); - static constexpr unsigned int _substract(unsigned int element, unsigned int v); - static constexpr int _get_inverse(unsigned int element, const unsigned int mod); - - unsigned int element_; - static inline const std::vector primes_ = [](){ - std::vector res; - for (unsigned int i = minimum; i <= maximum; ++i){ - if (_is_prime(i)){ - res.push_back(i); - } - } - return res; - }(); - static inline constexpr unsigned int productOfAllCharacteristics_ = [](){ - unsigned int res = 1; - for (unsigned int i = minimum; i <= maximum; ++i){ - if (_is_prime(i)){ - res *= i; - } - } - return res; - }(); - static inline const std::vector partials_ = [](){ - std::vector res; - - if (productOfAllCharacteristics_ == 1) return res; - - for (unsigned int i = 0; i < primes_.size(); ++i){ - unsigned int p = primes_[i]; - unsigned int base = productOfAllCharacteristics_ / p; - unsigned int exp = p - 1; - res.push_back(1); - - while (exp > 0) { - // If exp is odd, multiply with result - if (exp & 1) - res.back() = _multiply(res.back(), base); - // y must be even now - exp = exp >> 1; - base = _multiply(base, base); - } - } - - return res; - }(); - //If I understood the paper well, multiplicativeID_ always equals to 1. But in Clement's code, - //multiplicativeID_ is computed (see commented lambda function below). TODO: verify with Clement. - static inline constexpr unsigned int multiplicativeID_ = 1; /*= [](){ - unsigned int res = 0; - for (unsigned int i = 0; i < partials_.size(); ++i){ - res = (res + partials_[i]) % productOfAllCharacteristics_; - } - - return res; - }();*/ + public: + using element_type = unsigned int; /**< Type for the elements in the field. */ + template + using isInteger = std::enable_if_t >; + + /** + * @brief Default constructor. Sets the element to 0. + */ + Multi_field_element_with_small_characteristics(); + /** + * @brief Constructor setting the element to the given value. + * + * @param element Value of the element. + */ + Multi_field_element_with_small_characteristics(unsigned int element); + /** + * @brief Constructor setting the element to the given value. + * + * @param element Value of the element. + */ + Multi_field_element_with_small_characteristics(int element); + /** + * @brief Copy constructor. + * + * @param toCopy Element to copy. + */ + Multi_field_element_with_small_characteristics(const Multi_field_element_with_small_characteristics& toCopy); + /** + * @brief Move constructor. + * + * @param toMove Element to move. + */ + Multi_field_element_with_small_characteristics(Multi_field_element_with_small_characteristics&& toMove) noexcept; + + /** + * @brief operator+= + */ + friend void operator+=(Multi_field_element_with_small_characteristics& f1, + Multi_field_element_with_small_characteristics const& f2) { + f1.element_ = _add(f1.element_, f2.element_); + } + /** + * @brief operator+ + */ + friend Multi_field_element_with_small_characteristics operator+( + Multi_field_element_with_small_characteristics f1, Multi_field_element_with_small_characteristics const& f2) { + f1 += f2; + return f1; + } + /** + * @brief operator+= + */ + friend void operator+=(Multi_field_element_with_small_characteristics& f, unsigned int const v) { + f.element_ = _add(f.element_, v < productOfAllCharacteristics_ ? v : v % productOfAllCharacteristics_); + } + /** + * @brief operator+ + * + * @warning @p v is assumed to be positive and will be casted into an unsigned int. + */ + template > + friend Multi_field_element_with_small_characteristics operator+(Multi_field_element_with_small_characteristics f, + const Integer_type v) { + f += v; + return f; + } + /** + * @brief operator+ + */ + template > + friend Integer_type operator+(Integer_type v, Multi_field_element_with_small_characteristics const& f) { + v += f.element_; + v %= productOfAllCharacteristics_; + return v; + } + + /** + * @brief operator-= + */ + friend void operator-=(Multi_field_element_with_small_characteristics& f1, + Multi_field_element_with_small_characteristics const& f2) { + f1.element_ = _substract(f1.element_, f2.element_); + } + /** + * @brief operator- + */ + friend Multi_field_element_with_small_characteristics operator-( + Multi_field_element_with_small_characteristics f1, Multi_field_element_with_small_characteristics const& f2) { + f1 -= f2; + return f1; + } + /** + * @brief operator-= + */ + friend void operator-=(Multi_field_element_with_small_characteristics& f, unsigned int const v) { + f.element_ = _substract(f.element_, v < productOfAllCharacteristics_ ? v : v % productOfAllCharacteristics_); + } + /** + * @brief operator- + * + * @warning @p v is assumed to be positive and will be casted into an unsigned int. + */ + template > + friend Multi_field_element_with_small_characteristics operator-(Multi_field_element_with_small_characteristics f, + const Integer_type v) { + f -= v; + return f; + } + /** + * @brief operator- + * + * @warning @p v is assumed to be positive. + */ + template > + friend Integer_type operator-(Integer_type v, Multi_field_element_with_small_characteristics const& f) { + if (v >= productOfAllCharacteristics_) v %= productOfAllCharacteristics_; + if (f.element_ > v) v += productOfAllCharacteristics_; + v -= f.element_; + return v; + } + + /** + * @brief operator*= + */ + friend void operator*=(Multi_field_element_with_small_characteristics& f1, + Multi_field_element_with_small_characteristics const& f2) { + f1.element_ = _multiply(f1.element_, f2.element_); + } + /** + * @brief operator* + */ + friend Multi_field_element_with_small_characteristics operator*( + Multi_field_element_with_small_characteristics f1, Multi_field_element_with_small_characteristics const& f2) { + f1 *= f2; + return f1; + } + /** + * @brief operator*= + */ + friend void operator*=(Multi_field_element_with_small_characteristics& f, unsigned int const v) { + f.element_ = _multiply(f.element_, v < productOfAllCharacteristics_ ? v : v % productOfAllCharacteristics_); + } + /** + * @brief operator* + * + * @warning @p v is assumed to be positive and will be casted into an unsigned int. + */ + template > + friend Multi_field_element_with_small_characteristics operator*(Multi_field_element_with_small_characteristics f, + const Integer_type v) { + f *= v; + return f; + } + /** + * @brief operator* + * + * @warning Uses bitwise operations on @p v, so be carefull with signed integers. + */ + template > + friend Integer_type operator*(Integer_type v, Multi_field_element_with_small_characteristics const& f) { + unsigned int b = f.element_; + unsigned int res = 0; + unsigned int temp_b; + + while (v != 0) { + if (v & 1) { + if (b >= productOfAllCharacteristics_ - res) res -= productOfAllCharacteristics_; + res += b; + } + v >>= 1; + + temp_b = b; + if (b >= productOfAllCharacteristics_ - b) temp_b -= productOfAllCharacteristics_; + b += temp_b; + } + + return res; + } + + /** + * @brief operator== + */ + friend bool operator==(const Multi_field_element_with_small_characteristics& f1, + const Multi_field_element_with_small_characteristics& f2) { + return f1.element_ == f2.element_; + } + /** + * @brief operator== + * + * @warning @p v is assumed to be positive. + */ + template > + friend bool operator==(const Integer_type v, const Multi_field_element_with_small_characteristics& f) { + if (v < productOfAllCharacteristics_) return v == f.element_; + return (v % productOfAllCharacteristics_) == f.element_; + } + /** + * @brief operator== + * + * @warning @p v is assumed to be positive. + */ + template > + friend bool operator==(const Multi_field_element_with_small_characteristics& f, const Integer_type v) { + if (v < productOfAllCharacteristics_) return v == f.element_; + return (v % productOfAllCharacteristics_) == f.element_; + } + /** + * @brief operator!= + */ + friend bool operator!=(const Multi_field_element_with_small_characteristics& f1, + const Multi_field_element_with_small_characteristics& f2) { + return !(f1 == f2); + } + /** + * @brief operator!= + * + * @warning @p v is assumed to be positive. + */ + template > + friend bool operator!=(const Integer_type v, const Multi_field_element_with_small_characteristics& f) { + return !(v == f); + } + /** + * @brief operator!= + * + * @warning @p v is assumed to be positive. + */ + template > + friend bool operator!=(const Multi_field_element_with_small_characteristics& f, const Integer_type v) { + return !(v == f); + } + + /** + * @brief Assign operator. + */ + Multi_field_element_with_small_characteristics& operator=(Multi_field_element_with_small_characteristics other); + /** + * @brief Assign operator. + */ + Multi_field_element_with_small_characteristics& operator=(const unsigned int value); + /** + * @brief Swap operator. + */ + friend void swap(Multi_field_element_with_small_characteristics& f1, + Multi_field_element_with_small_characteristics& f2) { + std::swap(f1.element_, f2.element_); + } + + /** + * @brief Casts the element into an unsigned int. + */ + operator unsigned int() const; + + /** + * @brief Returns the inverse of the element in the multi-field, see @cite boissonnat:hal-00922572. + * + * @return The inverse. + */ + Multi_field_element_with_small_characteristics get_inverse() const; + /** + * @brief Returns the inverse of the element with respect to a sub-product of the characteristics in the multi-field, + * see @cite boissonnat:hal-00922572. + * + * @param productOfCharacteristics Sub-product of the characteristics. + * @return Pair of the inverse and the characteristic the inverse corresponds to. + */ + std::pair get_partial_inverse( + unsigned int productOfCharacteristics) const; + + /** + * @brief Returns the additive identity of a field. + * + * @return The additive identity of a field. + */ + static Multi_field_element_with_small_characteristics get_additive_identity(); + /** + * @brief Returns the multiplicative identity of a field. + * + * @return The multiplicative identity of a field. + */ + static Multi_field_element_with_small_characteristics get_multiplicative_identity(); + /** + * @brief Returns the partial multiplicative identity of the multi-field from the given product. + * See @cite boissonnat:hal-00922572 for more details. + * + * @param productOfCharacteristics Product of the different characteristics to take into account in the multi-field. + * @return The partial multiplicative identity of the multi-field. + */ + static Multi_field_element_with_small_characteristics get_partial_multiplicative_identity( + const mpz_class& productOfCharacteristics); + /** + * @brief Returns the product of all characteristics. + * + * @return The product of all characteristics. + */ + static constexpr unsigned int get_characteristic(); + + /** + * @brief Returns the value of the element. + * + * @return Value of the element. + */ + unsigned int get_value() const; + + // static constexpr bool handles_only_z2() { return false; } + + private: + static constexpr bool _is_prime(const int p); + static constexpr unsigned int _multiply(unsigned int a, unsigned int b); + static constexpr unsigned int _add(unsigned int element, unsigned int v); + static constexpr unsigned int _substract(unsigned int element, unsigned int v); + static constexpr int _get_inverse(unsigned int element, const unsigned int mod); + + unsigned int element_; + static inline const std::vector primes_ = []() { + std::vector res; + for (unsigned int i = minimum; i <= maximum; ++i) { + if (_is_prime(i)) { + res.push_back(i); + } + } + return res; + }(); + static inline constexpr unsigned int productOfAllCharacteristics_ = []() { + unsigned int res = 1; + for (unsigned int i = minimum; i <= maximum; ++i) { + if (_is_prime(i)) { + res *= i; + } + } + return res; + }(); + static inline const std::vector partials_ = []() { + std::vector res; + + if (productOfAllCharacteristics_ == 1) return res; + + for (unsigned int i = 0; i < primes_.size(); ++i) { + unsigned int p = primes_[i]; + unsigned int base = productOfAllCharacteristics_ / p; + unsigned int exp = p - 1; + res.push_back(1); + + while (exp > 0) { + // If exp is odd, multiply with result + if (exp & 1) res.back() = _multiply(res.back(), base); + // y must be even now + exp = exp >> 1; + base = _multiply(base, base); + } + } + + return res; + }(); + // If I understood the paper well, multiplicativeID_ always equals to 1. But in Clement's code, + // multiplicativeID_ is computed (see commented lambda function below). TODO: verify with Clement. + static inline constexpr unsigned int multiplicativeID_ = 1; /*= [](){ + unsigned int res = 0; + for (unsigned int i = 0; i < partials_.size(); ++i){ + res = (res + partials_[i]) % productOfAllCharacteristics_; + } + + return res; + }();*/ }; -template -inline Multi_field_element_with_small_characteristics::Multi_field_element_with_small_characteristics() - : element_(0) -{ - static_assert(maximum >= 2, "Characteristics have to be positive."); - static_assert(minimum <= maximum, "The given interval is not valid."); - static_assert(minimum != maximum || _is_prime(minimum), "The given interval does not contain a prime number."); - static_assert(productOfAllCharacteristics_ != 1, "The given interval does not contain a prime number."); +template +inline Multi_field_element_with_small_characteristics::Multi_field_element_with_small_characteristics() + : element_(0) { + static_assert(maximum >= 2, "Characteristics have to be positive."); + static_assert(minimum <= maximum, "The given interval is not valid."); + static_assert(minimum != maximum || _is_prime(minimum), "The given interval does not contain a prime number."); + static_assert(productOfAllCharacteristics_ != 1, "The given interval does not contain a prime number."); } -template -inline Multi_field_element_with_small_characteristics::Multi_field_element_with_small_characteristics(unsigned int element) - : element_(element % productOfAllCharacteristics_) -{ - static_assert(maximum >= 2, "Characteristics has to be positive."); - static_assert(minimum <= maximum, "The given interval is not valid."); - static_assert(minimum != maximum || _is_prime(minimum), "The given interval does not contain a prime number."); - static_assert(productOfAllCharacteristics_ != 1, "The given interval does not contain a prime number."); +template +inline Multi_field_element_with_small_characteristics::Multi_field_element_with_small_characteristics( + unsigned int element) + : element_(element % productOfAllCharacteristics_) { + static_assert(maximum >= 2, "Characteristics has to be positive."); + static_assert(minimum <= maximum, "The given interval is not valid."); + static_assert(minimum != maximum || _is_prime(minimum), "The given interval does not contain a prime number."); + static_assert(productOfAllCharacteristics_ != 1, "The given interval does not contain a prime number."); } -template -inline Multi_field_element_with_small_characteristics::Multi_field_element_with_small_characteristics(int element) - : element_(element % productOfAllCharacteristics_) -{ - static_assert(maximum >= 2, "Characteristics has to be positive."); - static_assert(minimum <= maximum, "The given interval is not valid."); - static_assert(minimum != maximum || _is_prime(minimum), "The given interval does not contain a prime number."); - static_assert(productOfAllCharacteristics_ != 1, "The given interval does not contain a prime number."); +template +inline Multi_field_element_with_small_characteristics::Multi_field_element_with_small_characteristics( + int element) + : element_(element % productOfAllCharacteristics_) { + static_assert(maximum >= 2, "Characteristics has to be positive."); + static_assert(minimum <= maximum, "The given interval is not valid."); + static_assert(minimum != maximum || _is_prime(minimum), "The given interval does not contain a prime number."); + static_assert(productOfAllCharacteristics_ != 1, "The given interval does not contain a prime number."); } -template -inline Multi_field_element_with_small_characteristics::Multi_field_element_with_small_characteristics(const Multi_field_element_with_small_characteristics &toCopy) - : element_(toCopy.element_) -{} - -template -inline Multi_field_element_with_small_characteristics::Multi_field_element_with_small_characteristics(Multi_field_element_with_small_characteristics &&toMove) noexcept - : element_(std::exchange(toMove.element_, 0)) -{} - -//template -//inline Multi_field_element_with_small_characteristics &Multi_field_element_with_small_characteristics::operator+=(Multi_field_element_with_small_characteristics const &f) -//{ -// _add(f.element_); -// return *this; -//} - -//template -//inline Multi_field_element_with_small_characteristics &Multi_field_element_with_small_characteristics::operator+=(unsigned int const v) -//{ -// _add(v % productOfAllCharacteristics_); -// return *this; -//} - -//template -//inline Multi_field_element_with_small_characteristics &Multi_field_element_with_small_characteristics::operator-=(Multi_field_element_with_small_characteristics const &f) -//{ -// _substract(f.element_); -// return *this; -//} - -//template -//inline Multi_field_element_with_small_characteristics &Multi_field_element_with_small_characteristics::operator-=(unsigned int const v) -//{ -// _substract(v % productOfAllCharacteristics_); -// return *this; -//} - -//template -//inline Multi_field_element_with_small_characteristics &Multi_field_element_with_small_characteristics::operator*=(Multi_field_element_with_small_characteristics const &f) -//{ -// element_ = _multiply(element_, f.element_); -// return *this; -//} - -//template -//inline Multi_field_element_with_small_characteristics &Multi_field_element_with_small_characteristics::operator*=(unsigned int const v) -//{ -// element_ = _multiply(element_, v % productOfAllCharacteristics_); -// return *this; -//} - -template -inline Multi_field_element_with_small_characteristics &Multi_field_element_with_small_characteristics::operator=(Multi_field_element_with_small_characteristics other) -{ - std::swap(element_, other.element_); - return *this; +template +inline Multi_field_element_with_small_characteristics::Multi_field_element_with_small_characteristics( + const Multi_field_element_with_small_characteristics& toCopy) + : element_(toCopy.element_) {} + +template +inline Multi_field_element_with_small_characteristics::Multi_field_element_with_small_characteristics( + Multi_field_element_with_small_characteristics&& toMove) noexcept + : element_(std::exchange(toMove.element_, 0)) {} + +template +inline Multi_field_element_with_small_characteristics& +Multi_field_element_with_small_characteristics::operator=( + Multi_field_element_with_small_characteristics other) { + std::swap(element_, other.element_); + return *this; } -template -inline Multi_field_element_with_small_characteristics &Multi_field_element_with_small_characteristics::operator=(unsigned int const value) -{ - element_ = value % productOfAllCharacteristics_; - return *this; +template +inline Multi_field_element_with_small_characteristics& +Multi_field_element_with_small_characteristics::operator=(unsigned int const value) { + element_ = value % productOfAllCharacteristics_; + return *this; } -template -inline Multi_field_element_with_small_characteristics::operator unsigned int() const -{ - return element_; +template +inline Multi_field_element_with_small_characteristics::operator unsigned int() const { + return element_; } -template -inline Multi_field_element_with_small_characteristics Multi_field_element_with_small_characteristics::get_inverse() const -{ - return get_partial_inverse(productOfAllCharacteristics_).first; +template +inline Multi_field_element_with_small_characteristics +Multi_field_element_with_small_characteristics::get_inverse() const { + return get_partial_inverse(productOfAllCharacteristics_).first; } -template -inline std::pair, unsigned int> Multi_field_element_with_small_characteristics::get_partial_inverse(unsigned int productOfCharacteristics) const -{ - unsigned int gcd = std::gcd(element_, productOfAllCharacteristics_); +template +inline std::pair, unsigned int> +Multi_field_element_with_small_characteristics::get_partial_inverse( + unsigned int productOfCharacteristics) const { + unsigned int gcd = std::gcd(element_, productOfAllCharacteristics_); - if (gcd == productOfCharacteristics) - return {Multi_field_element_with_small_characteristics(), multiplicativeID_}; // partial inverse is 0 + if (gcd == productOfCharacteristics) + return {Multi_field_element_with_small_characteristics(), multiplicativeID_}; // partial inverse is 0 - unsigned int QT = productOfCharacteristics / gcd; + unsigned int QT = productOfCharacteristics / gcd; - const unsigned int inv_qt = _get_inverse(element_, QT); + const unsigned int inv_qt = _get_inverse(element_, QT); - auto res = get_partial_multiplicative_identity(QT); - res *= inv_qt; + auto res = get_partial_multiplicative_identity(QT); + res *= inv_qt; - return {res, QT}; + return {res, QT}; } -template -inline Multi_field_element_with_small_characteristics Multi_field_element_with_small_characteristics::get_additive_identity() -{ - return Multi_field_element_with_small_characteristics(); +template +inline Multi_field_element_with_small_characteristics +Multi_field_element_with_small_characteristics::get_additive_identity() { + return Multi_field_element_with_small_characteristics(); } -template -inline Multi_field_element_with_small_characteristics Multi_field_element_with_small_characteristics::get_multiplicative_identity() -{ - return Multi_field_element_with_small_characteristics(multiplicativeID_); +template +inline Multi_field_element_with_small_characteristics +Multi_field_element_with_small_characteristics::get_multiplicative_identity() { + return Multi_field_element_with_small_characteristics(multiplicativeID_); } -template -inline Multi_field_element_with_small_characteristics Multi_field_element_with_small_characteristics::get_partial_multiplicative_identity(const mpz_class& productOfCharacteristics) -{ - if (productOfCharacteristics == 0) { - return Multi_field_element_with_small_characteristics(multiplicativeID_); - } - Multi_field_element_with_small_characteristics mult; - for (unsigned int idx = 0; idx < primes_.size(); ++idx) { - if ((productOfCharacteristics % primes_[idx]) == 0) { - mult += partials_[idx]; - } - } - return mult; +template +inline Multi_field_element_with_small_characteristics +Multi_field_element_with_small_characteristics::get_partial_multiplicative_identity( + const mpz_class& productOfCharacteristics) { + if (productOfCharacteristics == 0) { + return Multi_field_element_with_small_characteristics(multiplicativeID_); + } + Multi_field_element_with_small_characteristics mult; + for (unsigned int idx = 0; idx < primes_.size(); ++idx) { + if ((productOfCharacteristics % primes_[idx]) == 0) { + mult += partials_[idx]; + } + } + return mult; } -template -inline constexpr unsigned int Multi_field_element_with_small_characteristics::get_characteristic() -{ - return productOfAllCharacteristics_; +template +inline constexpr unsigned int Multi_field_element_with_small_characteristics::get_characteristic() { + return productOfAllCharacteristics_; } -template -inline unsigned int Multi_field_element_with_small_characteristics::get_value() const -{ - return element_; +template +inline unsigned int Multi_field_element_with_small_characteristics::get_value() const { + return element_; } -template -inline constexpr unsigned int Multi_field_element_with_small_characteristics::_add(unsigned int element, unsigned int v) -{ - if (UINT_MAX - element < v) { - //automatic unsigned integer overflow behaviour will make it work - element += v; - element -= productOfAllCharacteristics_; - return element; - } +template +inline constexpr unsigned int Multi_field_element_with_small_characteristics::_add( + unsigned int element, unsigned int v) { + if (UINT_MAX - element < v) { + // automatic unsigned integer overflow behaviour will make it work + element += v; + element -= productOfAllCharacteristics_; + return element; + } - element += v; - if (element >= productOfAllCharacteristics_) element -= productOfAllCharacteristics_; + element += v; + if (element >= productOfAllCharacteristics_) element -= productOfAllCharacteristics_; - return element; + return element; } -template -inline constexpr unsigned int Multi_field_element_with_small_characteristics::_substract(unsigned int element, unsigned int v) -{ - if (element < v){ - element += productOfAllCharacteristics_; - } - element -= v; +template +inline constexpr unsigned int Multi_field_element_with_small_characteristics::_substract( + unsigned int element, unsigned int v) { + if (element < v) { + element += productOfAllCharacteristics_; + } + element -= v; - return element; + return element; } -template -inline constexpr unsigned int Multi_field_element_with_small_characteristics::_multiply(unsigned int a, unsigned int b) -{ - unsigned int res = 0; - unsigned int temp_b = 0; - - if (b < a) std::swap(a, b); - - while (a != 0) { - if (a & 1) { - /* Add b to res, modulo m, without overflow */ - if (b >= productOfAllCharacteristics_ - res) - res -= productOfAllCharacteristics_; - res += b; - } - a >>= 1; - - /* Double b, modulo m */ - temp_b = b; - if (b >= productOfAllCharacteristics_ - b) - temp_b -= productOfAllCharacteristics_; - b += temp_b; - } - return res; +template +inline constexpr unsigned int Multi_field_element_with_small_characteristics::_multiply( + unsigned int a, unsigned int b) { + unsigned int res = 0; + unsigned int temp_b = 0; + + if (b < a) std::swap(a, b); + + while (a != 0) { + if (a & 1) { + /* Add b to res, modulo m, without overflow */ + if (b >= productOfAllCharacteristics_ - res) res -= productOfAllCharacteristics_; + res += b; + } + a >>= 1; + + /* Double b, modulo m */ + temp_b = b; + if (b >= productOfAllCharacteristics_ - b) temp_b -= productOfAllCharacteristics_; + b += temp_b; + } + return res; } -template -inline constexpr int Multi_field_element_with_small_characteristics::_get_inverse(unsigned int element, const unsigned int mod) -{ - //to solve: Ax + My = 1 - int M = mod; - int A = element; - int y = 0, x = 1; - //extended euclidien division - while (A > 1) { - int quotient = A / M; - int temp = M; - - M = A % M, A = temp; - temp = y; - - y = x - quotient * y; - x = temp; - } - - if (x < 0) - x += mod; - - return x; +template +inline constexpr int Multi_field_element_with_small_characteristics::_get_inverse( + unsigned int element, const unsigned int mod) { + // to solve: Ax + My = 1 + int M = mod; + int A = element; + int y = 0, x = 1; + // extended euclidien division + while (A > 1) { + int quotient = A / M; + int temp = M; + + M = A % M, A = temp; + temp = y; + + y = x - quotient * y; + x = temp; + } + + if (x < 0) x += mod; + + return x; } -template -inline constexpr bool Multi_field_element_with_small_characteristics::_is_prime(const int p) -{ - if (p <= 1) return false; - if (p <= 3) return true; - if (p % 2 == 0 || p % 3 == 0) return false; +template +inline constexpr bool Multi_field_element_with_small_characteristics::_is_prime(const int p) { + if (p <= 1) return false; + if (p <= 3) return true; + if (p % 2 == 0 || p % 3 == 0) return false; - for (long i = 5; i * i <= p; i = i + 6) - if (p % i == 0 || p % (i + 2) == 0) - return false; + for (long i = 5; i * i <= p; i = i + 6) + if (p % i == 0 || p % (i + 2) == 0) return false; - return true; + return true; } -} //namespace persistence_fields -} //namespace Gudhi +} // namespace persistence_fields +} // namespace Gudhi #endif // MATRIX_FIELD_MULTI_SMALL_H_ diff --git a/src/Persistence_matrix/include/gudhi/Fields/Multi_field_small_operators.h b/src/Persistence_matrix/include/gudhi/Fields/Multi_field_small_operators.h index 058c641027..d3b9024505 100644 --- a/src/Persistence_matrix/include/gudhi/Fields/Multi_field_small_operators.h +++ b/src/Persistence_matrix/include/gudhi/Fields/Multi_field_small_operators.h @@ -1,13 +1,19 @@ /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. - * Author(s): Hannah Schreiber + * Author(s): Hannah Schreiber, Clément Maria * - * Copyright (C) 2022 Inria + * Copyright (C) 2022-24 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification */ +/** + * @file Multi_field_small_operators.h + * @author Hannah Schreiber + * @brief Contains the @ref Multi_field_operators_with_small_characteristics class. + */ + #ifndef MATRIX_FIELD_MULTI_SMALL_OPERATORS_H_ #define MATRIX_FIELD_MULTI_SMALL_OPERATORS_H_ @@ -20,274 +26,381 @@ namespace Gudhi { namespace persistence_fields { -//productOfAllCharacteristics_ ^ 2 has to fit in an unsigned int /** + * @class Multi_field_operators_with_small_characteristics Multi_field_small_operators.h \ + * gudhi/Fields/Multi_field_small_operators.h * @ingroup persistence_fields * - * @brief TODO: - * + * @brief Class defining operators for a multi-field with "consecutive" charateristic range, such that + * `productOfAllCharacteristics ^ 2` fits into an unsigned int. */ -class Multi_field_operators_with_small_characteristics { -public: - using element_type = unsigned int; - using characteristic_type = element_type; - - Multi_field_operators_with_small_characteristics() : productOfAllCharacteristics_(0)/* , multiplicativeID_(1) */ {} - Multi_field_operators_with_small_characteristics(int minCharacteristic, int maxCharacteristic) - : productOfAllCharacteristics_(0)//, multiplicativeID_(1) - { - set_characteristic(minCharacteristic, maxCharacteristic); - } - Multi_field_operators_with_small_characteristics(const Multi_field_operators_with_small_characteristics& toCopy) - : primes_(toCopy.primes_), - productOfAllCharacteristics_(toCopy.productOfAllCharacteristics_), - partials_(toCopy.partials_)/* , - multiplicativeID_(toCopy.multiplicativeID_) */ - {} - Multi_field_operators_with_small_characteristics(Multi_field_operators_with_small_characteristics&& toMove) noexcept - : primes_(std::move(toMove.primes_)), - productOfAllCharacteristics_(std::move(toMove.productOfAllCharacteristics_)), - partials_(std::move(toMove.partials_))/* , - multiplicativeID_(std::move(toMove.multiplicativeID_)) */ - {} - - void set_characteristic(int minimum, int maximum){ - if (maximum < 2) - throw std::invalid_argument("Characteristic must be strictly positive"); - if (minimum > maximum) - throw std::invalid_argument("The given interval is not valid."); - if (minimum == maximum && !_is_prime(minimum)) - throw std::invalid_argument("The given interval does not contain a prime number."); - - productOfAllCharacteristics_ = 1; - primes_.clear(); - for (unsigned int i = minimum; i <= static_cast(maximum); ++i){ - if (_is_prime(i)){ - primes_.push_back(i); - productOfAllCharacteristics_ *= i; - } - } - - if (primes_.empty()) - throw std::invalid_argument("The given interval does not contain a prime number."); - - partials_.resize(primes_.size()); - for (unsigned int i = 0; i < primes_.size(); ++i){ - unsigned int p = primes_[i]; - characteristic_type base = productOfAllCharacteristics_ / p; - unsigned int exp = p - 1; - partials_[i] = 1; - - while (exp > 0) { - // If exp is odd, multiply with result - if (exp & 1) - partials_[i] = _multiply(partials_[i], base, productOfAllCharacteristics_); - // y must be even now - exp = exp >> 1; // y = y/2 - base = _multiply(base, base, productOfAllCharacteristics_); - } - } - - //If I understood the paper well, multiplicativeID_ always equals to 1. But in Clement's code, - //multiplicativeID_ is computed (see commented loop below). TODO: verify with Clement. - // for (unsigned int i = 0; i < partials_.size(); ++i){ - // multiplicativeID_ = (multiplicativeID_ + partials_[i]) % productOfAllCharacteristics_; - // } - } - characteristic_type get_characteristic() const{ - return productOfAllCharacteristics_; - } - - element_type get_value(element_type e) const{ - return e < productOfAllCharacteristics_ ? e : e % productOfAllCharacteristics_; - } - - //r = e1 + e2 - element_type add(element_type e1, element_type e2) const{ - return _add(get_value(e1), get_value(e2), productOfAllCharacteristics_); - } - - //r = e1 - e2 - element_type substract(element_type e1, element_type e2) const{ - return _substract(get_value(e1), get_value(e2), productOfAllCharacteristics_); - } - - //r = e1 * e2 - element_type multiply(element_type e1, element_type e2) const{ - return _multiply(get_value(e1), get_value(e2), productOfAllCharacteristics_); - } - - //r = e * m + a - //WARNING: not overflow safe. - element_type multiply_and_add(element_type e, element_type m, element_type a) const{ - return get_value(e * m + a); - } - - //r = (e + a) * m - //WARNING: not overflow safe. - element_type add_and_multiply(element_type e, element_type a, element_type m) const{ - return get_value((e + a) * m); - } - - bool are_equal(element_type e1, element_type e2) const{ - return get_value(e1) == get_value(e2); - } - - element_type get_inverse(const element_type& e) const{ - return get_partial_inverse(e, productOfAllCharacteristics_).first; - } - std::pair get_partial_inverse(const element_type& e, const characteristic_type& productOfCharacteristics) const{ - characteristic_type gcd = std::gcd(e, productOfAllCharacteristics_); - - if (gcd == productOfCharacteristics) - return {0, get_multiplicative_identity()}; // partial inverse is 0 - - characteristic_type QT = productOfCharacteristics / gcd; - - const characteristic_type inv_qt = _get_inverse(e, QT); - - auto res = get_partial_multiplicative_identity(QT); - res = _multiply(res, inv_qt, productOfAllCharacteristics_); - - return {res, QT}; - } - - static constexpr element_type get_additive_identity(){ return 0; } - static constexpr element_type get_multiplicative_identity(){ return 1; } - // static element_type get_multiplicative_identity(){ return multiplicativeID_; } - element_type get_partial_multiplicative_identity(const characteristic_type& productOfCharacteristics) const - { - if (productOfCharacteristics == 0) { - return get_multiplicative_identity(); - } - element_type multIdentity = 0; - for (unsigned int idx = 0; idx < primes_.size(); ++idx) { - if ((productOfCharacteristics % primes_[idx]) == 0) { - multIdentity = _add(multIdentity, partials_[idx], productOfAllCharacteristics_); - } - } - return multIdentity; - } - - static constexpr bool handles_only_z2(){ - return false; - } - - Multi_field_operators_with_small_characteristics& operator=(Multi_field_operators_with_small_characteristics other){ - primes_.swap(other.primes_); - productOfAllCharacteristics_ = other.productOfAllCharacteristics_; - partials_.swap(other.partials_); - - return *this; - } - - friend void swap(Multi_field_operators_with_small_characteristics& f1, Multi_field_operators_with_small_characteristics& f2){ - f1.primes_.swap(f2.primes_); - std::swap(f1.productOfAllCharacteristics_, f2.productOfAllCharacteristics_); - f1.partials_.swap(f2.partials_); - } - -private: - std::vector primes_; - characteristic_type productOfAllCharacteristics_; - std::vector partials_; - // static inline constexpr unsigned int multiplicativeID_ = 1; - - static element_type _add(element_type element, element_type v, characteristic_type characteristic); - static element_type _substract(element_type element, element_type v, characteristic_type characteristic); - static element_type _multiply(element_type a, element_type b, characteristic_type characteristic); - static constexpr long int _get_inverse(element_type element, characteristic_type mod); - static constexpr bool _is_prime(const int p); +class Multi_field_operators_with_small_characteristics +{ + public: + using element_type = unsigned int; /**< Type for the elements in the field. */ + using characteristic_type = element_type; /**< Type for the field characteristic. */ + + /** + * @brief Default constructor, sets the product of all characteristics to 0. + */ + Multi_field_operators_with_small_characteristics() : productOfAllCharacteristics_(0) /* , multiplicativeID_(1) */ + {} + /** + * @brief Constructor setting the characteristics to all prime numbers between the two given integers. + * The product of all primes to the square has to fit into an unsigned int. + * + * @param minCharacteristic Smallest value of a prime. + * @param maxCharacteristic Heighest value of a prime. + */ + Multi_field_operators_with_small_characteristics(int minCharacteristic, int maxCharacteristic) + : productOfAllCharacteristics_(0) //, multiplicativeID_(1) + { + set_characteristic(minCharacteristic, maxCharacteristic); + } + /** + * @brief Copy constructor. + * + * @param toCopy Operators to copy. + */ + Multi_field_operators_with_small_characteristics(const Multi_field_operators_with_small_characteristics& toCopy) + : primes_(toCopy.primes_), + productOfAllCharacteristics_(toCopy.productOfAllCharacteristics_), + partials_(toCopy.partials_) /* , + multiplicativeID_(toCopy.multiplicativeID_) */ + {} + /** + * @brief Move constructor. + * + * @param toMove Operators to move. + */ + Multi_field_operators_with_small_characteristics(Multi_field_operators_with_small_characteristics&& toMove) noexcept + : primes_(std::move(toMove.primes_)), + productOfAllCharacteristics_(std::move(toMove.productOfAllCharacteristics_)), + partials_(std::move(toMove.partials_)) /* , + multiplicativeID_(std::move(toMove.multiplicativeID_)) */ + {} + + /** + * @brief Set the characteristics of the field, which are stored in a single value as a product of all of them. + * The characteristics will be all prime numbers in the given interval. + * The product of all primes to the square has to fit into an unsigned int. + * + * @param minimum Smallest value of a prime. + * @param maximum Heighest value of a prime. + */ + void set_characteristic(int minimum, int maximum) { + if (maximum < 2) throw std::invalid_argument("Characteristic must be strictly positive"); + if (minimum > maximum) throw std::invalid_argument("The given interval is not valid."); + if (minimum == maximum && !_is_prime(minimum)) + throw std::invalid_argument("The given interval does not contain a prime number."); + + productOfAllCharacteristics_ = 1; + primes_.clear(); + for (unsigned int i = minimum; i <= static_cast(maximum); ++i) { + if (_is_prime(i)) { + primes_.push_back(i); + productOfAllCharacteristics_ *= i; + } + } + + if (primes_.empty()) throw std::invalid_argument("The given interval does not contain a prime number."); + + partials_.resize(primes_.size()); + for (unsigned int i = 0; i < primes_.size(); ++i) { + unsigned int p = primes_[i]; + characteristic_type base = productOfAllCharacteristics_ / p; + unsigned int exp = p - 1; + partials_[i] = 1; + + while (exp > 0) { + // If exp is odd, multiply with result + if (exp & 1) partials_[i] = _multiply(partials_[i], base, productOfAllCharacteristics_); + // y must be even now + exp = exp >> 1; // y = y/2 + base = _multiply(base, base, productOfAllCharacteristics_); + } + } + + // If I understood the paper well, multiplicativeID_ always equals to 1. But in Clement's code, + // multiplicativeID_ is computed (see commented loop below). TODO: verify with Clement. + // for (unsigned int i = 0; i < partials_.size(); ++i){ + // multiplicativeID_ = (multiplicativeID_ + partials_[i]) % productOfAllCharacteristics_; + // } + } + /** + * @brief Returns the current characteristics as the product of all of them. + * + * @return The value of the current characteristic. + */ + characteristic_type get_characteristic() const { return productOfAllCharacteristics_; } + + /** + * @brief Returns the value of an element in the field. + * That is the positive value of the integer modulo the current characteristic. + * + * @param e Integer to return the value from. + * @return @p e modulo the current characteristic, such that the result is positive. + */ + element_type get_value(element_type e) const { + return e < productOfAllCharacteristics_ ? e : e % productOfAllCharacteristics_; + } + + /** + * @brief Returns the sum of two elements in the field. + * + * @param e1 First element. + * @param e2 Second element. + * @return `(e1 + e2) % productOfAllCharacteristics`, such that the result is positive. + */ + element_type add(element_type e1, element_type e2) const { + return _add(get_value(e1), get_value(e2), productOfAllCharacteristics_); + } + + /** + * @brief Returns the substraction in the field of the first element by the second element. + * + * @param e1 First element. + * @param e2 Second element. + * @return `(e1 - e2) % productOfAllCharacteristics`, such that the result is positive. + */ + element_type substract(element_type e1, element_type e2) const { + return _substract(get_value(e1), get_value(e2), productOfAllCharacteristics_); + } + + /** + * @brief Returns the multiplication of two elements in the field. + * + * @param e1 First element. + * @param e2 Second element. + * @return `(e1 * e2) % productOfAllCharacteristics`, such that the result is positive. + */ + element_type multiply(element_type e1, element_type e2) const { + return _multiply(get_value(e1), get_value(e2), productOfAllCharacteristics_); + } + + /** + * @brief Multiplies the first element with the second one and adds the third one. Returns the result in the field. + * + * @warning Not overflow safe. + * + * @param e First element. + * @param m Second element. + * @param a Third element. + * @return `(e * m + a) % productOfAllCharacteristics`, such that the result is positive. + */ + element_type multiply_and_add(element_type e, element_type m, element_type a) const { return get_value(e * m + a); } + + /** + * @brief Adds the first element to the second one and multiplies the third one with it. + * Returns the result in the field. + * + * @warning Not overflow safe. + * + * @param e First element. + * @param a Second element. + * @param m Third element. + * @return `((e + a) * m) % productOfAllCharacteristics`, such that the result is positive. + */ + element_type add_and_multiply(element_type e, element_type a, element_type m) const { return get_value((e + a) * m); } + + /** + * @brief Returns true if the two given elements are equal in the field, false otherwise. + * + * @param e1 First element to compare. + * @param e2 Second element to compare. + * @return true If `e1 % productOfAllCharacteristics == e2 % productOfAllCharacteristics`. + * @return false Otherwise. + */ + bool are_equal(element_type e1, element_type e2) const { return get_value(e1) == get_value(e2); } + + /** + * @brief Returns the inverse of the given element in the sense of @cite boissonnat:hal-00922572 with respect + * to the product of all characteristics. + * + * @param e Element to get the inverse from. + * @return Inverse in the current field. + */ + element_type get_inverse(const element_type& e) const { + return get_partial_inverse(e, productOfAllCharacteristics_).first; + } + /** + * @brief Returns the inverse of the given element in the multi-field corresponding to the given sub-product + * of the product of all characteristics in the multi-field. See @cite boissonnat:hal-00922572 for more details. + * + * @param e Element to get the inverse from. + * @param productOfCharacteristics Product of the different characteristics to take into account in the multi-field. + * @return Pair of the inverse of @p e and the characteristic the inverse is coming from. + */ + std::pair get_partial_inverse( + const element_type& e, const characteristic_type& productOfCharacteristics) const { + characteristic_type gcd = std::gcd(e, productOfAllCharacteristics_); + + if (gcd == productOfCharacteristics) return {0, get_multiplicative_identity()}; // partial inverse is 0 + + characteristic_type QT = productOfCharacteristics / gcd; + + const characteristic_type inv_qt = _get_inverse(e, QT); + + auto res = get_partial_multiplicative_identity(QT); + res = _multiply(res, inv_qt, productOfAllCharacteristics_); + + return {res, QT}; + } + + /** + * @brief Returns the additive identity of a field. + * + * @return The additive identity of a field. + */ + static constexpr element_type get_additive_identity() { return 0; } + /** + * @brief Returns the multiplicative identity of a field. + * + * @return The multiplicative identity of a field. + */ + static constexpr element_type get_multiplicative_identity() { return 1; } + // static element_type get_multiplicative_identity(){ return multiplicativeID_; } + /** + * @brief Returns the partial multiplicative identity of the multi-field from the given product. + * See @cite boissonnat:hal-00922572 for more details. + * + * @param productOfCharacteristics Product of the different characteristics to take into account in the multi-field. + * @return The partial multiplicative identity of the multi-field. + */ + element_type get_partial_multiplicative_identity(const characteristic_type& productOfCharacteristics) const { + if (productOfCharacteristics == 0) { + return get_multiplicative_identity(); + } + element_type multIdentity = 0; + for (unsigned int idx = 0; idx < primes_.size(); ++idx) { + if ((productOfCharacteristics % primes_[idx]) == 0) { + multIdentity = _add(multIdentity, partials_[idx], productOfAllCharacteristics_); + } + } + return multIdentity; + } + + // static constexpr bool handles_only_z2() { return false; } + + /** + * @brief Assign operator. + */ + Multi_field_operators_with_small_characteristics& operator=(Multi_field_operators_with_small_characteristics other) { + primes_.swap(other.primes_); + productOfAllCharacteristics_ = other.productOfAllCharacteristics_; + partials_.swap(other.partials_); + + return *this; + } + /** + * @brief Swap operator. + */ + friend void swap(Multi_field_operators_with_small_characteristics& f1, + Multi_field_operators_with_small_characteristics& f2) { + f1.primes_.swap(f2.primes_); + std::swap(f1.productOfAllCharacteristics_, f2.productOfAllCharacteristics_); + f1.partials_.swap(f2.partials_); + } + + private: + std::vector primes_; /**< All characteristics. */ + characteristic_type productOfAllCharacteristics_; /**< Product of all characteristics. */ + std::vector partials_; /**< Partial products of the characteristics. */ + // static inline constexpr unsigned int multiplicativeID_ = 1; + + static element_type _add(element_type element, element_type v, characteristic_type characteristic); + static element_type _substract(element_type element, element_type v, characteristic_type characteristic); + static element_type _multiply(element_type a, element_type b, characteristic_type characteristic); + static constexpr long int _get_inverse(element_type element, characteristic_type mod); + static constexpr bool _is_prime(const int p); }; -inline Multi_field_operators_with_small_characteristics::element_type Multi_field_operators_with_small_characteristics::_add(element_type element, element_type v, characteristic_type characteristic) -{ - if (UINT_MAX - element < v) { - //automatic unsigned integer overflow behaviour will make it work - element += v; - element -= characteristic; - return element; - } +inline Multi_field_operators_with_small_characteristics::element_type +Multi_field_operators_with_small_characteristics::_add(element_type element, element_type v, + characteristic_type characteristic) { + if (UINT_MAX - element < v) { + // automatic unsigned integer overflow behaviour will make it work + element += v; + element -= characteristic; + return element; + } - element += v; - if (element >= characteristic) element -= characteristic; + element += v; + if (element >= characteristic) element -= characteristic; - return element; + return element; } -inline Multi_field_operators_with_small_characteristics::element_type Multi_field_operators_with_small_characteristics::_substract(element_type element, element_type v, characteristic_type characteristic) -{ - if (element < v){ - element += characteristic; - } - element -= v; +inline Multi_field_operators_with_small_characteristics::element_type +Multi_field_operators_with_small_characteristics::_substract(element_type element, element_type v, + characteristic_type characteristic) { + if (element < v) { + element += characteristic; + } + element -= v; - return element; + return element; } -inline Multi_field_operators_with_small_characteristics::element_type Multi_field_operators_with_small_characteristics::_multiply(element_type a, element_type b, characteristic_type characteristic) -{ - element_type res = 0; - element_type temp_b = 0; - - if (b < a) std::swap(a, b); - - while (a != 0) { - if (a & 1) { - /* Add b to res, modulo m, without overflow */ - if (b >= characteristic - res) - res -= characteristic; - res += b; - } - a >>= 1; - - /* Double b, modulo m */ - temp_b = b; - if (b >= characteristic - b) - temp_b -= characteristic; - b += temp_b; - } - return res; +inline Multi_field_operators_with_small_characteristics::element_type +Multi_field_operators_with_small_characteristics::_multiply(element_type a, element_type b, + characteristic_type characteristic) { + element_type res = 0; + element_type temp_b = 0; + + if (b < a) std::swap(a, b); + + while (a != 0) { + if (a & 1) { + /* Add b to res, modulo m, without overflow */ + if (b >= characteristic - res) res -= characteristic; + res += b; + } + a >>= 1; + + /* Double b, modulo m */ + temp_b = b; + if (b >= characteristic - b) temp_b -= characteristic; + b += temp_b; + } + return res; } -inline constexpr long int Multi_field_operators_with_small_characteristics::_get_inverse(element_type element, characteristic_type mod) -{ - //to solve: Ax + My = 1 - element_type M = mod; - element_type A = element; - long int y = 0, x = 1; - //extended euclidien division - while (A > 1) { - int quotient = A / M; - int temp = M; - - M = A % M, A = temp; - temp = y; - - y = x - quotient * y; - x = temp; - } - - if (x < 0) - x += mod; - - return x; +inline constexpr long int Multi_field_operators_with_small_characteristics::_get_inverse(element_type element, + characteristic_type mod) { + // to solve: Ax + My = 1 + element_type M = mod; + element_type A = element; + long int y = 0, x = 1; + // extended euclidien division + while (A > 1) { + int quotient = A / M; + int temp = M; + + M = A % M, A = temp; + temp = y; + + y = x - quotient * y; + x = temp; + } + + if (x < 0) x += mod; + + return x; } -inline constexpr bool Multi_field_operators_with_small_characteristics::_is_prime(const int p) -{ - if (p <= 1) return false; - if (p <= 3) return true; - if (p % 2 == 0 || p % 3 == 0) return false; +inline constexpr bool Multi_field_operators_with_small_characteristics::_is_prime(const int p) { + if (p <= 1) return false; + if (p <= 3) return true; + if (p % 2 == 0 || p % 3 == 0) return false; - for (long i = 5; i * i <= p; i = i + 6) - if (p % i == 0 || p % (i + 2) == 0) - return false; + for (long i = 5; i * i <= p; i = i + 6) + if (p % i == 0 || p % (i + 2) == 0) return false; - return true; + return true; } -} //namespace persistence_fields -} //namespace Gudhi +} // namespace persistence_fields +} // namespace Gudhi #endif // MATRIX_FIELD_MULTI_SMALL_OPERATORS_H_ diff --git a/src/Persistence_matrix/include/gudhi/Fields/Multi_field_small_shared.h b/src/Persistence_matrix/include/gudhi/Fields/Multi_field_small_shared.h index 4a8da63449..d769ef4061 100644 --- a/src/Persistence_matrix/include/gudhi/Fields/Multi_field_small_shared.h +++ b/src/Persistence_matrix/include/gudhi/Fields/Multi_field_small_shared.h @@ -1,13 +1,19 @@ /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. - * Author(s): Hannah Schreiber + * Author(s): Hannah Schreiber, Clément Maria * - * Copyright (C) 2022 Inria + * Copyright (C) 2022-24 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification */ +/** + * @file Multi_field_small_shared.h + * @author Hannah Schreiber + * @brief Contains the @ref Shared_multi_field_element_with_small_characteristics class. + */ + #ifndef MATRIX_FIELD_MULTI_SMALL_SHARED_H_ #define MATRIX_FIELD_MULTI_SMALL_SHARED_H_ @@ -20,441 +26,610 @@ namespace Gudhi { namespace persistence_fields { -//productOfAllCharacteristics_ ^ 2 has to fit in an unsigned int /** + * @class Shared_multi_field_element_with_small_characteristics Multi_field_small_shared.h \ + * gudhi/Fields/Multi_field_small_shared.h * @ingroup persistence_fields * - * @brief TODO: - * + * @brief Class representing an element of a multi-field, such that `productOfAllCharacteristics ^ 2` fits into + * an unsigned int. If each instanciation of the class can represent another element, they all share the same + * characteritics. That is if the characteristics are set for one, they will be set for all the others. + * The characteristics can be set before instianciating the elements with the static + * @ref Shared_multi_field_element::initialize method. */ -class Shared_multi_field_element_with_small_characteristics { -public: - using element_type = unsigned int; - template - using isInteger = std::enable_if_t >; - - Shared_multi_field_element_with_small_characteristics(); - Shared_multi_field_element_with_small_characteristics(unsigned int element); - Shared_multi_field_element_with_small_characteristics(int element); - Shared_multi_field_element_with_small_characteristics(const Shared_multi_field_element_with_small_characteristics& toCopy); - Shared_multi_field_element_with_small_characteristics(Shared_multi_field_element_with_small_characteristics&& toMove) noexcept; - - static void initialize(unsigned int minimum, unsigned int maximum); - - friend void operator+=(Shared_multi_field_element_with_small_characteristics& f1, Shared_multi_field_element_with_small_characteristics const& f2){ - f1.element_ = _add(f1.element_, f2.element_); - } - friend Shared_multi_field_element_with_small_characteristics operator+(Shared_multi_field_element_with_small_characteristics f1, Shared_multi_field_element_with_small_characteristics const& f2){ - f1 += f2; - return f1; - } - friend void operator+=(Shared_multi_field_element_with_small_characteristics& f, unsigned int const v){ - f.element_ = _add(f.element_, v < productOfAllCharacteristics_ ? v : v % productOfAllCharacteristics_); - } - //v is assumed to be positive and will be casted into an unsigned int - template > - friend Shared_multi_field_element_with_small_characteristics operator+(Shared_multi_field_element_with_small_characteristics f, const Integer_type v){ - f += v; - return f; - } - template > - friend Integer_type operator+(Integer_type v, Shared_multi_field_element_with_small_characteristics const& f){ - v += f.element_; - v %= productOfAllCharacteristics_; - return v; - } - - friend void operator-=(Shared_multi_field_element_with_small_characteristics& f1, Shared_multi_field_element_with_small_characteristics const& f2){ - f1.element_ = _substract(f1.element_, f2.element_); - } - friend Shared_multi_field_element_with_small_characteristics operator-(Shared_multi_field_element_with_small_characteristics f1, Shared_multi_field_element_with_small_characteristics const& f2){ - f1 -= f2; - return f1; - } - friend void operator-=(Shared_multi_field_element_with_small_characteristics& f, unsigned int const v){ - f.element_ = _substract(f.element_, v < productOfAllCharacteristics_ ? v : v % productOfAllCharacteristics_); - } - //v is assumed to be positive and will be casted into an unsigned int - template > - friend Shared_multi_field_element_with_small_characteristics operator-(Shared_multi_field_element_with_small_characteristics f, const Integer_type v){ - f -= v; - return f; - } - //v is assumed to be positive - template > - friend Integer_type operator-(Integer_type v, Shared_multi_field_element_with_small_characteristics const& f){ - if (v >= productOfAllCharacteristics_) v %= productOfAllCharacteristics_; - if (f.element_ > v) v += productOfAllCharacteristics_; - v -= f.element_; - return v; - } - - friend void operator*=(Shared_multi_field_element_with_small_characteristics& f1, Shared_multi_field_element_with_small_characteristics const& f2){ - f1.element_ = _multiply(f1.element_, f2.element_); - } - friend Shared_multi_field_element_with_small_characteristics operator*(Shared_multi_field_element_with_small_characteristics f1, Shared_multi_field_element_with_small_characteristics const& f2){ - f1 *= f2; - return f1; - } - friend void operator*=(Shared_multi_field_element_with_small_characteristics& f, unsigned int const v){ - f.element_ = _multiply(f.element_, v < productOfAllCharacteristics_ ? v : v % productOfAllCharacteristics_); - } - //v is assumed to be positive and will be casted into an unsigned int - template > - friend Shared_multi_field_element_with_small_characteristics operator*(Shared_multi_field_element_with_small_characteristics f, const Integer_type v){ - f *= v; - return f; - } - //uses bitwise operations on v, so be carefull with signed integers - template > - friend Integer_type operator*(Integer_type v, Shared_multi_field_element_with_small_characteristics const& f){ - unsigned int b = f.element_; - unsigned int res = 0; - unsigned int temp_b; - - while (v != 0) { - if (v & 1) { - if (b >= productOfAllCharacteristics_ - res) - res -= productOfAllCharacteristics_; - res += b; - } - v >>= 1; - - temp_b = b; - if (b >= productOfAllCharacteristics_ - b) - temp_b -= productOfAllCharacteristics_; - b += temp_b; - } - - return res; - } - - friend bool operator==(const Shared_multi_field_element_with_small_characteristics& f1, const Shared_multi_field_element_with_small_characteristics& f2){ - return f1.element_ == f2.element_; - } - //v is assumed to be positive - template > - friend bool operator==(const Integer_type v, const Shared_multi_field_element_with_small_characteristics& f){ - if (v < productOfAllCharacteristics_) return v == f.element_; - return (v % productOfAllCharacteristics_) == f.element_; - } - //v is assumed to be positive - template > - friend bool operator==(const Shared_multi_field_element_with_small_characteristics& f, const Integer_type v){ - if (v < productOfAllCharacteristics_) return v == f.element_; - return (v % productOfAllCharacteristics_) == f.element_; - } - friend bool operator!=(const Shared_multi_field_element_with_small_characteristics& f1, const Shared_multi_field_element_with_small_characteristics& f2){ - return !(f1 == f2); - } - //v is assumed to be positive - template > - friend bool operator!=(const Integer_type v, const Shared_multi_field_element_with_small_characteristics& f){ - return !(v == f); - } - //v is assumed to be positive - template > - friend bool operator!=(const Shared_multi_field_element_with_small_characteristics& f, const Integer_type v){ - return !(v == f); - } - - Shared_multi_field_element_with_small_characteristics& operator=(Shared_multi_field_element_with_small_characteristics other); - Shared_multi_field_element_with_small_characteristics& operator=(const unsigned int value); - operator unsigned int() const; - friend void swap(Shared_multi_field_element_with_small_characteristics& f1, Shared_multi_field_element_with_small_characteristics& f2){ - std::swap(f1.element_, f2.element_); - } - - Shared_multi_field_element_with_small_characteristics get_inverse() const; - std::pair get_partial_inverse(unsigned int productOfCharacteristics) const; - - static Shared_multi_field_element_with_small_characteristics get_additive_identity(); - static Shared_multi_field_element_with_small_characteristics get_multiplicative_identity(); - static Shared_multi_field_element_with_small_characteristics get_partial_multiplicative_identity(const mpz_class& productOfCharacteristics); - static unsigned int get_characteristic(); - - unsigned int get_value() const; - - static constexpr bool handles_only_z2(){ - return false; - } - -private: - static constexpr bool _is_prime(const int p); - static unsigned int _multiply(unsigned int a, unsigned int b); - static unsigned int _add(unsigned int element, unsigned int v); - static unsigned int _substract(unsigned int element, unsigned int v); - static constexpr int _get_inverse(unsigned int element, const unsigned int mod); - - unsigned int element_; - static inline std::vector primes_; - static inline unsigned int productOfAllCharacteristics_; - static inline std::vector partials_; - static inline constexpr unsigned int multiplicativeID_ = 1; +class Shared_multi_field_element_with_small_characteristics +{ + public: + using element_type = unsigned int; /**< Type for the elements in the field. */ + template + using isInteger = std::enable_if_t >; + + /** + * @brief Default constructor. Sets the element to 0. + */ + Shared_multi_field_element_with_small_characteristics(); + /** + * @brief Constructor setting the element to the given value. + * + * @param element Value of the element. + */ + Shared_multi_field_element_with_small_characteristics(unsigned int element); + /** + * @brief Constructor setting the element to the given value. + * + * @param element Value of the element. + */ + Shared_multi_field_element_with_small_characteristics(int element); + /** + * @brief Copy constructor. + * + * @param toCopy Element to copy. + */ + Shared_multi_field_element_with_small_characteristics( + const Shared_multi_field_element_with_small_characteristics& toCopy); + /** + * @brief Move constructor. + * + * @param toMove Element to move. + */ + Shared_multi_field_element_with_small_characteristics( + Shared_multi_field_element_with_small_characteristics&& toMove) noexcept; + + /** + * @brief Initialize the multi-field to the characteristics (primes) contained in the given interval. + * Should be called first before constructing the field elements. + * The characteristics must be small enough such that `productOfAllCharacteristics ^ 2` fits into an unsigned int. + * + * @param minimum Lowest value in the interval. + * @param maximum Highest value in the interval. + */ + static void initialize(unsigned int minimum, unsigned int maximum); + + /** + * @brief operator+= + */ + friend void operator+=(Shared_multi_field_element_with_small_characteristics& f1, + Shared_multi_field_element_with_small_characteristics const& f2) { + f1.element_ = _add(f1.element_, f2.element_); + } + /** + * @brief operator+ + */ + friend Shared_multi_field_element_with_small_characteristics operator+( + Shared_multi_field_element_with_small_characteristics f1, + Shared_multi_field_element_with_small_characteristics const& f2) { + f1 += f2; + return f1; + } + /** + * @brief operator+= + */ + friend void operator+=(Shared_multi_field_element_with_small_characteristics& f, unsigned int const v) { + f.element_ = _add(f.element_, v < productOfAllCharacteristics_ ? v : v % productOfAllCharacteristics_); + } + /** + * @brief operator+ + * + * @warning @p v is assumed to be positive and will be casted into an unsigned int. + */ + template > + friend Shared_multi_field_element_with_small_characteristics operator+( + Shared_multi_field_element_with_small_characteristics f, const Integer_type v) { + f += v; + return f; + } + /** + * @brief operator+ + */ + template > + friend Integer_type operator+(Integer_type v, Shared_multi_field_element_with_small_characteristics const& f) { + v += f.element_; + v %= productOfAllCharacteristics_; + return v; + } + + /** + * @brief operator-= + */ + friend void operator-=(Shared_multi_field_element_with_small_characteristics& f1, + Shared_multi_field_element_with_small_characteristics const& f2) { + f1.element_ = _substract(f1.element_, f2.element_); + } + /** + * @brief operator- + */ + friend Shared_multi_field_element_with_small_characteristics operator-( + Shared_multi_field_element_with_small_characteristics f1, + Shared_multi_field_element_with_small_characteristics const& f2) { + f1 -= f2; + return f1; + } + /** + * @brief operator-= + */ + friend void operator-=(Shared_multi_field_element_with_small_characteristics& f, unsigned int const v) { + f.element_ = _substract(f.element_, v < productOfAllCharacteristics_ ? v : v % productOfAllCharacteristics_); + } + /** + * @brief operator- + * + * @warning @p v is assumed to be positive and will be casted into an unsigned int. + */ + template > + friend Shared_multi_field_element_with_small_characteristics operator-( + Shared_multi_field_element_with_small_characteristics f, const Integer_type v) { + f -= v; + return f; + } + /** + * @brief operator- + * + * @warning @p v is assumed to be positive. + */ + template > + friend Integer_type operator-(Integer_type v, Shared_multi_field_element_with_small_characteristics const& f) { + if (v >= productOfAllCharacteristics_) v %= productOfAllCharacteristics_; + if (f.element_ > v) v += productOfAllCharacteristics_; + v -= f.element_; + return v; + } + + /** + * @brief operator*= + */ + friend void operator*=(Shared_multi_field_element_with_small_characteristics& f1, + Shared_multi_field_element_with_small_characteristics const& f2) { + f1.element_ = _multiply(f1.element_, f2.element_); + } + /** + * @brief operator* + */ + friend Shared_multi_field_element_with_small_characteristics operator*( + Shared_multi_field_element_with_small_characteristics f1, + Shared_multi_field_element_with_small_characteristics const& f2) { + f1 *= f2; + return f1; + } + /** + * @brief operator*= + */ + friend void operator*=(Shared_multi_field_element_with_small_characteristics& f, unsigned int const v) { + f.element_ = _multiply(f.element_, v < productOfAllCharacteristics_ ? v : v % productOfAllCharacteristics_); + } + /** + * @brief operator* + * + * @warning @p v is assumed to be positive and will be casted into an unsigned int. + */ + template > + friend Shared_multi_field_element_with_small_characteristics operator*( + Shared_multi_field_element_with_small_characteristics f, const Integer_type v) { + f *= v; + return f; + } + /** + * @brief operator* + * + * @warning Uses bitwise operations on @p v, so be carefull with signed integers. + */ + template > + friend Integer_type operator*(Integer_type v, Shared_multi_field_element_with_small_characteristics const& f) { + unsigned int b = f.element_; + unsigned int res = 0; + unsigned int temp_b; + + while (v != 0) { + if (v & 1) { + if (b >= productOfAllCharacteristics_ - res) res -= productOfAllCharacteristics_; + res += b; + } + v >>= 1; + + temp_b = b; + if (b >= productOfAllCharacteristics_ - b) temp_b -= productOfAllCharacteristics_; + b += temp_b; + } + + return res; + } + + /** + * @brief operator== + */ + friend bool operator==(const Shared_multi_field_element_with_small_characteristics& f1, + const Shared_multi_field_element_with_small_characteristics& f2) { + return f1.element_ == f2.element_; + } + /** + * @brief operator== + * + * @warning @p v is assumed to be positive. + */ + template > + friend bool operator==(const Integer_type v, const Shared_multi_field_element_with_small_characteristics& f) { + if (v < productOfAllCharacteristics_) return v == f.element_; + return (v % productOfAllCharacteristics_) == f.element_; + } + /** + * @brief operator== + * + * @warning @p v is assumed to be positive. + */ + template > + friend bool operator==(const Shared_multi_field_element_with_small_characteristics& f, const Integer_type v) { + if (v < productOfAllCharacteristics_) return v == f.element_; + return (v % productOfAllCharacteristics_) == f.element_; + } + /** + * @brief operator!= + */ + friend bool operator!=(const Shared_multi_field_element_with_small_characteristics& f1, + const Shared_multi_field_element_with_small_characteristics& f2) { + return !(f1 == f2); + } + /** + * @brief operator!= + * + * @warning @p v is assumed to be positive. + */ + template > + friend bool operator!=(const Integer_type v, const Shared_multi_field_element_with_small_characteristics& f) { + return !(v == f); + } + /** + * @brief operator!= + * + * @warning @p v is assumed to be positive. + */ + template > + friend bool operator!=(const Shared_multi_field_element_with_small_characteristics& f, const Integer_type v) { + return !(v == f); + } + + /** + * @brief Assign operator. + */ + Shared_multi_field_element_with_small_characteristics& operator=( + Shared_multi_field_element_with_small_characteristics other); + /** + * @brief Assign operator. + */ + Shared_multi_field_element_with_small_characteristics& operator=(const unsigned int value); + /** + * @brief Swap operator. + */ + friend void swap(Shared_multi_field_element_with_small_characteristics& f1, + Shared_multi_field_element_with_small_characteristics& f2) { + std::swap(f1.element_, f2.element_); + } + + /** + * @brief Casts the element into an unsigned int. + */ + operator unsigned int() const; + + /** + * @brief Returns the inverse of the element in the multi-field, see @cite boissonnat:hal-00922572. + * + * @return The inverse. + */ + Shared_multi_field_element_with_small_characteristics get_inverse() const; + /** + * @brief Returns the inverse of the element with respect to a sub-product of the characteristics in the multi-field, + * see @cite boissonnat:hal-00922572. + * + * @param productOfCharacteristics Sub-product of the characteristics. + * @return Pair of the inverse and the characteristic the inverse corresponds to. + */ + std::pair get_partial_inverse( + unsigned int productOfCharacteristics) const; + + /** + * @brief Returns the additive identity of a field. + * + * @return The additive identity of a field. + */ + static Shared_multi_field_element_with_small_characteristics get_additive_identity(); + /** + * @brief Returns the multiplicative identity of a field. + * + * @return The multiplicative identity of a field. + */ + static Shared_multi_field_element_with_small_characteristics get_multiplicative_identity(); + /** + * @brief Returns the partial multiplicative identity of the multi-field from the given product. + * See @cite boissonnat:hal-00922572 for more details. + * + * @param productOfCharacteristics Product of the different characteristics to take into account in the multi-field. + * @return The partial multiplicative identity of the multi-field. + */ + static Shared_multi_field_element_with_small_characteristics get_partial_multiplicative_identity( + const mpz_class& productOfCharacteristics); + /** + * @brief Returns the product of all characteristics. + * + * @return The product of all characteristics. + */ + static unsigned int get_characteristic(); + + /** + * @brief Returns the value of the element. + * + * @return Value of the element. + */ + unsigned int get_value() const; + + // static constexpr bool handles_only_z2() { return false; } + + private: + static constexpr bool _is_prime(const int p); + static unsigned int _multiply(unsigned int a, unsigned int b); + static unsigned int _add(unsigned int element, unsigned int v); + static unsigned int _substract(unsigned int element, unsigned int v); + static constexpr int _get_inverse(unsigned int element, const unsigned int mod); + + unsigned int element_; /**< Element. */ + static inline std::vector primes_; /**< All characteristics. */ + static inline unsigned int productOfAllCharacteristics_; /**< Product of all characteristics. */ + static inline std::vector partials_; /**< Partial products of the characteristics. */ + static inline constexpr unsigned int multiplicativeID_ = 1; /**< Multiplicative identity. */ }; inline Shared_multi_field_element_with_small_characteristics::Shared_multi_field_element_with_small_characteristics() - : element_(0) -{} - -inline Shared_multi_field_element_with_small_characteristics::Shared_multi_field_element_with_small_characteristics(unsigned int element) - : element_(element % productOfAllCharacteristics_) -{} - -inline Shared_multi_field_element_with_small_characteristics::Shared_multi_field_element_with_small_characteristics(int element) - : element_(element % productOfAllCharacteristics_) -{} - -inline Shared_multi_field_element_with_small_characteristics::Shared_multi_field_element_with_small_characteristics(const Shared_multi_field_element_with_small_characteristics &toCopy) - : element_(toCopy.element_) -{} - -inline Shared_multi_field_element_with_small_characteristics::Shared_multi_field_element_with_small_characteristics(Shared_multi_field_element_with_small_characteristics &&toMove) noexcept - : element_(std::exchange(toMove.element_, 0)) -{} - -inline void Shared_multi_field_element_with_small_characteristics::initialize(unsigned int minimum, unsigned int maximum) -{ - if (maximum < 2) - throw std::invalid_argument("Characteristic must be strictly positive"); - if (minimum > maximum) - throw std::invalid_argument("The given interval is not valid."); - if (minimum == maximum && !_is_prime(minimum)) - throw std::invalid_argument("The given interval does not contain a prime number."); - - productOfAllCharacteristics_ = 1; - primes_.clear(); - for (unsigned int i = minimum; i <= maximum; ++i){ - if (_is_prime(i)){ - primes_.push_back(i); - productOfAllCharacteristics_ *= i; - } - } - - if (primes_.empty()) - throw std::invalid_argument("The given interval does not contain a prime number."); - - partials_.resize(primes_.size()); - for (unsigned int i = 0; i < primes_.size(); ++i){ - unsigned int p = primes_[i]; - unsigned int base = productOfAllCharacteristics_ / p; - unsigned int exp = p - 1; - partials_[i] = 1; - - while (exp > 0) { - // If exp is odd, multiply with result - if (exp & 1) - partials_[i] = _multiply(partials_[i], base); - // y must be even now - exp = exp >> 1; // y = y/2 - base = _multiply(base, base); - } - } - - //If I understood the paper well, multiplicativeID_ always equals to 1. But in Clement's code, - //multiplicativeID_ is computed (see commented loop below). TODO: verify with Clement. -// for (unsigned int i = 0; i < partials_.size(); ++i){ -// multiplicativeID_ = (multiplicativeID_ + partials_[i]) % productOfAllCharacteristics_; -// } + : element_(0) {} + +inline Shared_multi_field_element_with_small_characteristics::Shared_multi_field_element_with_small_characteristics( + unsigned int element) + : element_(element % productOfAllCharacteristics_) {} + +inline Shared_multi_field_element_with_small_characteristics::Shared_multi_field_element_with_small_characteristics( + int element) + : element_(element % productOfAllCharacteristics_) {} + +inline Shared_multi_field_element_with_small_characteristics::Shared_multi_field_element_with_small_characteristics( + const Shared_multi_field_element_with_small_characteristics& toCopy) + : element_(toCopy.element_) {} + +inline Shared_multi_field_element_with_small_characteristics::Shared_multi_field_element_with_small_characteristics( + Shared_multi_field_element_with_small_characteristics&& toMove) noexcept + : element_(std::exchange(toMove.element_, 0)) {} + +inline void Shared_multi_field_element_with_small_characteristics::initialize(unsigned int minimum, + unsigned int maximum) { + if (maximum < 2) throw std::invalid_argument("Characteristic must be strictly positive"); + if (minimum > maximum) throw std::invalid_argument("The given interval is not valid."); + if (minimum == maximum && !_is_prime(minimum)) + throw std::invalid_argument("The given interval does not contain a prime number."); + + productOfAllCharacteristics_ = 1; + primes_.clear(); + for (unsigned int i = minimum; i <= maximum; ++i) { + if (_is_prime(i)) { + primes_.push_back(i); + productOfAllCharacteristics_ *= i; + } + } + + if (primes_.empty()) throw std::invalid_argument("The given interval does not contain a prime number."); + + partials_.resize(primes_.size()); + for (unsigned int i = 0; i < primes_.size(); ++i) { + unsigned int p = primes_[i]; + unsigned int base = productOfAllCharacteristics_ / p; + unsigned int exp = p - 1; + partials_[i] = 1; + + while (exp > 0) { + // If exp is odd, multiply with result + if (exp & 1) partials_[i] = _multiply(partials_[i], base); + // y must be even now + exp = exp >> 1; // y = y/2 + base = _multiply(base, base); + } + } + + // If I understood the paper well, multiplicativeID_ always equals to 1. But in Clement's code, + // multiplicativeID_ is computed (see commented loop below). TODO: verify with Clement. + // for (unsigned int i = 0; i < partials_.size(); ++i){ + // multiplicativeID_ = (multiplicativeID_ + partials_[i]) % productOfAllCharacteristics_; + // } } -//inline Shared_multi_field_element_with_small_characteristics &Shared_multi_field_element_with_small_characteristics::operator+=(Shared_multi_field_element_with_small_characteristics const &f) +// inline Shared_multi_field_element_with_small_characteristics +// &Shared_multi_field_element_with_small_characteristics::operator+=(Shared_multi_field_element_with_small_characteristics +// const &f) //{ // _add(f.element_); // return *this; -//} +// } -//inline Shared_multi_field_element_with_small_characteristics &Shared_multi_field_element_with_small_characteristics::operator+=(unsigned int const v) +// inline Shared_multi_field_element_with_small_characteristics +// &Shared_multi_field_element_with_small_characteristics::operator+=(unsigned int const v) //{ // _add(v % productOfAllCharacteristics_); // return *this; -//} +// } -//inline Shared_multi_field_element_with_small_characteristics &Shared_multi_field_element_with_small_characteristics::operator-=(Shared_multi_field_element_with_small_characteristics const &f) +// inline Shared_multi_field_element_with_small_characteristics +// &Shared_multi_field_element_with_small_characteristics::operator-=(Shared_multi_field_element_with_small_characteristics +// const &f) //{ // _substract(f.element_); // return *this; -//} +// } -//inline Shared_multi_field_element_with_small_characteristics &Shared_multi_field_element_with_small_characteristics::operator-=(unsigned int const v) +// inline Shared_multi_field_element_with_small_characteristics +// &Shared_multi_field_element_with_small_characteristics::operator-=(unsigned int const v) //{ // _substract(v % productOfAllCharacteristics_); // return *this; -//} +// } -//inline Shared_multi_field_element_with_small_characteristics &Shared_multi_field_element_with_small_characteristics::operator*=(Shared_multi_field_element_with_small_characteristics const &f) +// inline Shared_multi_field_element_with_small_characteristics +// &Shared_multi_field_element_with_small_characteristics::operator*=(Shared_multi_field_element_with_small_characteristics +// const &f) //{ // element_ = _multiply(element_, f.element_); // return *this; -//} +// } -//inline Shared_multi_field_element_with_small_characteristics &Shared_multi_field_element_with_small_characteristics::operator*=(unsigned int const v) +// inline Shared_multi_field_element_with_small_characteristics +// &Shared_multi_field_element_with_small_characteristics::operator*=(unsigned int const v) //{ // element_ = _multiply(element_, v % productOfAllCharacteristics_); // return *this; -//} +// } -inline Shared_multi_field_element_with_small_characteristics &Shared_multi_field_element_with_small_characteristics::operator=(Shared_multi_field_element_with_small_characteristics other) -{ - std::swap(element_, other.element_); - return *this; +inline Shared_multi_field_element_with_small_characteristics& +Shared_multi_field_element_with_small_characteristics::operator=( + Shared_multi_field_element_with_small_characteristics other) { + std::swap(element_, other.element_); + return *this; } -inline Shared_multi_field_element_with_small_characteristics &Shared_multi_field_element_with_small_characteristics::operator=(unsigned int const value) -{ - element_ = value % productOfAllCharacteristics_; - return *this; +inline Shared_multi_field_element_with_small_characteristics& +Shared_multi_field_element_with_small_characteristics::operator=(unsigned int const value) { + element_ = value % productOfAllCharacteristics_; + return *this; } -inline Shared_multi_field_element_with_small_characteristics::operator unsigned int() const -{ - return element_; -} +inline Shared_multi_field_element_with_small_characteristics::operator unsigned int() const { return element_; } -inline Shared_multi_field_element_with_small_characteristics Shared_multi_field_element_with_small_characteristics::get_inverse() const -{ - return get_partial_inverse(productOfAllCharacteristics_).first; +inline Shared_multi_field_element_with_small_characteristics +Shared_multi_field_element_with_small_characteristics::get_inverse() const { + return get_partial_inverse(productOfAllCharacteristics_).first; } -inline std::pair Shared_multi_field_element_with_small_characteristics::get_partial_inverse(unsigned int productOfCharacteristics) const -{ - unsigned int gcd = std::gcd(element_, productOfAllCharacteristics_); +inline std::pair +Shared_multi_field_element_with_small_characteristics::get_partial_inverse( + unsigned int productOfCharacteristics) const { + unsigned int gcd = std::gcd(element_, productOfAllCharacteristics_); - if (gcd == productOfCharacteristics) - return {Shared_multi_field_element_with_small_characteristics(), multiplicativeID_}; // partial inverse is 0 + if (gcd == productOfCharacteristics) + return {Shared_multi_field_element_with_small_characteristics(), multiplicativeID_}; // partial inverse is 0 - unsigned int QT = productOfCharacteristics / gcd; + unsigned int QT = productOfCharacteristics / gcd; - const unsigned int inv_qt = _get_inverse(element_, QT); + const unsigned int inv_qt = _get_inverse(element_, QT); - auto res = get_partial_multiplicative_identity(QT); - res *= inv_qt; + auto res = get_partial_multiplicative_identity(QT); + res *= inv_qt; - return {res, QT}; + return {res, QT}; } -inline Shared_multi_field_element_with_small_characteristics Shared_multi_field_element_with_small_characteristics::get_additive_identity() -{ - return Shared_multi_field_element_with_small_characteristics(); +inline Shared_multi_field_element_with_small_characteristics +Shared_multi_field_element_with_small_characteristics::get_additive_identity() { + return Shared_multi_field_element_with_small_characteristics(); } -inline Shared_multi_field_element_with_small_characteristics Shared_multi_field_element_with_small_characteristics::get_multiplicative_identity() -{ - return Shared_multi_field_element_with_small_characteristics(multiplicativeID_); +inline Shared_multi_field_element_with_small_characteristics +Shared_multi_field_element_with_small_characteristics::get_multiplicative_identity() { + return Shared_multi_field_element_with_small_characteristics(multiplicativeID_); } -inline Shared_multi_field_element_with_small_characteristics Shared_multi_field_element_with_small_characteristics::get_partial_multiplicative_identity(const mpz_class& productOfCharacteristics) -{ - if (productOfCharacteristics == 0) { - return Shared_multi_field_element_with_small_characteristics(multiplicativeID_); - } - Shared_multi_field_element_with_small_characteristics mult; - for (unsigned int idx = 0; idx < primes_.size(); ++idx) { - if ((productOfCharacteristics % primes_[idx]) == 0) { - mult += partials_[idx]; - } - } - return mult; +inline Shared_multi_field_element_with_small_characteristics +Shared_multi_field_element_with_small_characteristics::get_partial_multiplicative_identity( + const mpz_class& productOfCharacteristics) { + if (productOfCharacteristics == 0) { + return Shared_multi_field_element_with_small_characteristics(multiplicativeID_); + } + Shared_multi_field_element_with_small_characteristics mult; + for (unsigned int idx = 0; idx < primes_.size(); ++idx) { + if ((productOfCharacteristics % primes_[idx]) == 0) { + mult += partials_[idx]; + } + } + return mult; } -inline unsigned int Shared_multi_field_element_with_small_characteristics::get_characteristic() -{ - return productOfAllCharacteristics_; +inline unsigned int Shared_multi_field_element_with_small_characteristics::get_characteristic() { + return productOfAllCharacteristics_; } -inline unsigned int Shared_multi_field_element_with_small_characteristics::get_value() const -{ - return element_; -} +inline unsigned int Shared_multi_field_element_with_small_characteristics::get_value() const { return element_; } -inline unsigned int Shared_multi_field_element_with_small_characteristics::_add(unsigned int element, unsigned int v) -{ - if (UINT_MAX - element < v) { - //automatic unsigned integer overflow behaviour will make it work - element += v; - element -= productOfAllCharacteristics_; - return element; - } +inline unsigned int Shared_multi_field_element_with_small_characteristics::_add(unsigned int element, unsigned int v) { + if (UINT_MAX - element < v) { + // automatic unsigned integer overflow behaviour will make it work + element += v; + element -= productOfAllCharacteristics_; + return element; + } - element += v; - if (element >= productOfAllCharacteristics_) element -= productOfAllCharacteristics_; + element += v; + if (element >= productOfAllCharacteristics_) element -= productOfAllCharacteristics_; - return element; + return element; } -inline unsigned int Shared_multi_field_element_with_small_characteristics::_substract(unsigned int element, unsigned int v) -{ - if (element < v){ - element += productOfAllCharacteristics_; - } - element -= v; +inline unsigned int Shared_multi_field_element_with_small_characteristics::_substract(unsigned int element, + unsigned int v) { + if (element < v) { + element += productOfAllCharacteristics_; + } + element -= v; - return element; + return element; } -inline unsigned int Shared_multi_field_element_with_small_characteristics::_multiply(unsigned int a, unsigned int b) -{ - unsigned int res = 0; - unsigned int temp_b = 0; - - if (b < a) std::swap(a, b); - - while (a != 0) { - if (a & 1) { - /* Add b to res, modulo m, without overflow */ - if (b >= productOfAllCharacteristics_ - res) - res -= productOfAllCharacteristics_; - res += b; - } - a >>= 1; - - /* Double b, modulo m */ - temp_b = b; - if (b >= productOfAllCharacteristics_ - b) - temp_b -= productOfAllCharacteristics_; - b += temp_b; - } - return res; +inline unsigned int Shared_multi_field_element_with_small_characteristics::_multiply(unsigned int a, unsigned int b) { + unsigned int res = 0; + unsigned int temp_b = 0; + + if (b < a) std::swap(a, b); + + while (a != 0) { + if (a & 1) { + /* Add b to res, modulo m, without overflow */ + if (b >= productOfAllCharacteristics_ - res) res -= productOfAllCharacteristics_; + res += b; + } + a >>= 1; + + /* Double b, modulo m */ + temp_b = b; + if (b >= productOfAllCharacteristics_ - b) temp_b -= productOfAllCharacteristics_; + b += temp_b; + } + return res; } -inline constexpr int Shared_multi_field_element_with_small_characteristics::_get_inverse(unsigned int element, const unsigned int mod) -{ - //to solve: Ax + My = 1 - int M = mod; - int A = element; - int y = 0, x = 1; - //extended euclidien division - while (A > 1) { - int quotient = A / M; - int temp = M; - - M = A % M, A = temp; - temp = y; - - y = x - quotient * y; - x = temp; - } - - if (x < 0) - x += mod; - - return x; +inline constexpr int Shared_multi_field_element_with_small_characteristics::_get_inverse(unsigned int element, + const unsigned int mod) { + // to solve: Ax + My = 1 + int M = mod; + int A = element; + int y = 0, x = 1; + // extended euclidien division + while (A > 1) { + int quotient = A / M; + int temp = M; + + M = A % M, A = temp; + temp = y; + + y = x - quotient * y; + x = temp; + } + + if (x < 0) x += mod; + + return x; } -inline constexpr bool Shared_multi_field_element_with_small_characteristics::_is_prime(const int p) -{ - if (p <= 1) return false; - if (p <= 3) return true; - if (p % 2 == 0 || p % 3 == 0) return false; +inline constexpr bool Shared_multi_field_element_with_small_characteristics::_is_prime(const int p) { + if (p <= 1) return false; + if (p <= 3) return true; + if (p % 2 == 0 || p % 3 == 0) return false; - for (long i = 5; i * i <= p; i = i + 6) - if (p % i == 0 || p % (i + 2) == 0) - return false; + for (long i = 5; i * i <= p; i = i + 6) + if (p % i == 0 || p % (i + 2) == 0) return false; - return true; + return true; } -} //namespace persistence_fields -} //namespace Gudhi +} // namespace persistence_fields +} // namespace Gudhi #endif // MATRIX_FIELD_MULTI_SMALL_SHARED_H_ diff --git a/src/Persistence_matrix/include/gudhi/Fields/Z2_field.h b/src/Persistence_matrix/include/gudhi/Fields/Z2_field.h index b65d7f0364..40865cf36f 100644 --- a/src/Persistence_matrix/include/gudhi/Fields/Z2_field.h +++ b/src/Persistence_matrix/include/gudhi/Fields/Z2_field.h @@ -8,6 +8,12 @@ * - YYYY/MM Author: Description of the modification */ +/** + * @file Z2_field.h + * @author Hannah Schreiber + * @brief Contains the @ref Z2_field_element class. + */ + #ifndef MATRIX_FIELD_Z2_H_ #define MATRIX_FIELD_Z2_H_ @@ -17,208 +23,327 @@ namespace Gudhi { namespace persistence_fields { /** + * @class Z2_field_element Z2_field.h gudhi/Fields/Z2_field.h * @ingroup persistence_fields * - * @brief TODO: - * + * @brief Class representing an element of the \f$ \mathbb{F}_2 \f$ field. */ -class Z2_field_element { -public: - using element_type = unsigned int; - template - using isInteger = std::enable_if_t >; - - Z2_field_element(); - Z2_field_element(unsigned int element); - Z2_field_element(int element); - Z2_field_element(const Z2_field_element& toCopy); - Z2_field_element(Z2_field_element&& toMove) noexcept; - - friend void operator+=(Z2_field_element& f1, Z2_field_element const& f2){ - f1.element_ = (f1.element_ != f2.element_); - } - friend Z2_field_element operator+(Z2_field_element f1, Z2_field_element const& f2){ - f1 += f2; - return f1; - } - friend void operator+=(Z2_field_element& f, unsigned int const v){ - f.element_ = (f.element_ != (v % 2)); - } - //v will be casted into an unsigned int - template > - friend Z2_field_element operator+(Z2_field_element f, const Integer_type v){ - f += v; - return f; - } - template > - friend Integer_type operator+(const Integer_type v, Z2_field_element const& f){ - return f.element_ != (v % 2); - } - - friend void operator-=(Z2_field_element& f1, Z2_field_element const& f2){ - f1.element_ = (f1.element_ != f2.element_); - } - friend Z2_field_element operator-(Z2_field_element f1, Z2_field_element const& f2){ - f1 -= f2; - return f1; - } - friend void operator-=(Z2_field_element& f, unsigned int const v){ - f.element_ = (f.element_ != (v % 2)); - } - //v will be casted into an unsigned int - template > - friend Z2_field_element operator-(Z2_field_element f, const Integer_type v){ - f -= v; - return f; - } - template > - friend Integer_type operator-(const Integer_type v, Z2_field_element const& f){ - return f.element_ != (v % 2); - } - - friend void operator*=(Z2_field_element& f1, Z2_field_element const& f2){ - f1.element_ = (f1.element_ && f2.element_); - } - friend Z2_field_element operator*(Z2_field_element f1, Z2_field_element const& f2){ - f1 *= f2; - return f1; - } - friend void operator*=(Z2_field_element& f, unsigned int const v){ - f.element_ = (f.element_ && (v % 2)); - } - //v will be casted into an unsigned int - template > - friend Z2_field_element operator*(Z2_field_element f, const Integer_type v){ - f *= v; - return f; - } - template > - friend Integer_type operator*(const Integer_type v, Z2_field_element const& f){ - return f.element_ && (v % 2); - } - - friend bool operator==(const Z2_field_element& f1, const Z2_field_element& f2){ - return f1.element_ == f2.element_; - } - template > - friend bool operator==(const Integer_type v, const Z2_field_element& f){ - return (v % 2) == f.element_; - } - template > - friend bool operator==(const Z2_field_element& f, const Integer_type v){ - return (v % 2) == f.element_; - } - friend bool operator!=(const Z2_field_element& f1, const Z2_field_element& f2){ - return !(f1 == f2); - } - template > - friend bool operator!=(const Integer_type v, const Z2_field_element& f){ - return !(v == f); - } - template > - friend bool operator!=(const Z2_field_element& f, const Integer_type v){ - return !(v == f); - } - - Z2_field_element& operator=(Z2_field_element other); - Z2_field_element& operator=(const unsigned int value); - operator unsigned int() const; - friend void swap(Z2_field_element& f1, Z2_field_element& f2){ - std::swap(f1.element_, f2.element_); - } - - Z2_field_element get_inverse() const; - std::pair get_partial_inverse(unsigned int productOfCharacteristics) const; - - static Z2_field_element get_additive_identity(); - static Z2_field_element get_multiplicative_identity(); - static Z2_field_element get_partial_multiplicative_identity([[maybe_unused]] unsigned int productOfCharacteristics); - static constexpr unsigned int get_characteristic(); - - unsigned int get_value() const; - - static constexpr bool handles_only_z2(){ - return true; - } - -private: - bool element_; +class Z2_field_element +{ + public: + using element_type = unsigned int; /**< Type for the elements in the field. */ + template + using isInteger = std::enable_if_t >; + + /** + * @brief Default constructor. + */ + Z2_field_element(); + /** + * @brief Constructor setting the element to the given value. + * + * @param element Value of the element. + */ + Z2_field_element(unsigned int element); + /** + * @brief Constructor setting the element to the given value. + * + * @param element Value of the element. + */ + Z2_field_element(int element); + /** + * @brief Copy constructor. + * + * @param toCopy Element to copy. + */ + Z2_field_element(const Z2_field_element& toCopy); + /** + * @brief Move constructor. + * + * @param toMove Element to move. + */ + Z2_field_element(Z2_field_element&& toMove) noexcept; + + /** + * @brief operator+= + */ + friend void operator+=(Z2_field_element& f1, Z2_field_element const& f2) { + f1.element_ = (f1.element_ != f2.element_); + } + /** + * @brief operator+ + */ + friend Z2_field_element operator+(Z2_field_element f1, Z2_field_element const& f2) { + f1 += f2; + return f1; + } + /** + * @brief operator+= + */ + friend void operator+=(Z2_field_element& f, unsigned int const v) { f.element_ = (f.element_ != (v % 2)); } + /** + * @brief operator+ + * + * @warning @p v will be casted into an unsigned int. + * + * @tparam Integer_type A native integer type: int, unsigned int, long int, bool, etc. + */ + template > + friend Z2_field_element operator+(Z2_field_element f, const Integer_type v) { + f += v; + return f; + } + /** + * @brief operator+ + * + * @tparam Integer_type A native integer type: int, unsigned int, long int, bool, etc. + */ + template > + friend Integer_type operator+(const Integer_type v, Z2_field_element const& f) { + return f.element_ != (v % 2); + } + + /** + * @brief operator-= + */ + friend void operator-=(Z2_field_element& f1, Z2_field_element const& f2) { + f1.element_ = (f1.element_ != f2.element_); + } + /** + * @brief operator- + */ + friend Z2_field_element operator-(Z2_field_element f1, Z2_field_element const& f2) { + f1 -= f2; + return f1; + } + /** + * @brief operator-= + */ + friend void operator-=(Z2_field_element& f, unsigned int const v) { f.element_ = (f.element_ != (v % 2)); } + /** + * @brief operator- + * + * @warning @p v will be casted into an unsigned int. + * + * @tparam Integer_type A native integer type: int, unsigned int, long int, bool, etc. + */ + template > + friend Z2_field_element operator-(Z2_field_element f, const Integer_type v) { + f -= v; + return f; + } + /** + * @brief operator- + * + * @tparam Integer_type A native integer type: int, unsigned int, long int, bool, etc. + */ + template > + friend Integer_type operator-(const Integer_type v, Z2_field_element const& f) { + return f.element_ != (v % 2); + } + + /** + * @brief operator*= + */ + friend void operator*=(Z2_field_element& f1, Z2_field_element const& f2) { + f1.element_ = (f1.element_ && f2.element_); + } + /** + * @brief operator* + */ + friend Z2_field_element operator*(Z2_field_element f1, Z2_field_element const& f2) { + f1 *= f2; + return f1; + } + /** + * @brief operator*= + */ + friend void operator*=(Z2_field_element& f, unsigned int const v) { f.element_ = (f.element_ && (v % 2)); } + // v will be casted into an unsigned int + /** + * @brief operator* + * + * @warning @p v will be casted into an unsigned int. + * + * @tparam Integer_type A native integer type: int, unsigned int, long int, bool, etc. + */ + template > + friend Z2_field_element operator*(Z2_field_element f, const Integer_type v) { + f *= v; + return f; + } + /** + * @brief operator* + * + * @tparam Integer_type A native integer type: int, unsigned int, long int, bool, etc. + */ + template > + friend Integer_type operator*(const Integer_type v, Z2_field_element const& f) { + return f.element_ && (v % 2); + } + + /** + * @brief operator== + */ + friend bool operator==(const Z2_field_element& f1, const Z2_field_element& f2) { return f1.element_ == f2.element_; } + /** + * @brief operator== + * + * @tparam Integer_type A native integer type: int, unsigned int, long int, bool, etc. + */ + template > + friend bool operator==(const Integer_type v, const Z2_field_element& f) { + return (v % 2) == f.element_; + } + /** + * @brief operator== + * + * @tparam Integer_type A native integer type: int, unsigned int, long int, bool, etc. + */ + template > + friend bool operator==(const Z2_field_element& f, const Integer_type v) { + return (v % 2) == f.element_; + } + /** + * @brief operator!= + */ + friend bool operator!=(const Z2_field_element& f1, const Z2_field_element& f2) { return !(f1 == f2); } + /** + * @brief operator!= + * + * @tparam Integer_type A native integer type: int, unsigned int, long int, bool, etc. + */ + template > + friend bool operator!=(const Integer_type v, const Z2_field_element& f) { + return !(v == f); + } + /** + * @brief operator!= + * + * @tparam Integer_type A native integer type: int, unsigned int, long int, bool, etc. + */ + template > + friend bool operator!=(const Z2_field_element& f, const Integer_type v) { + return !(v == f); + } + + /** + * @brief Assign operator. + */ + Z2_field_element& operator=(Z2_field_element other); + /** + * @brief Assign operator. + */ + Z2_field_element& operator=(const unsigned int value); + /** + * @brief Swap operator. + */ + friend void swap(Z2_field_element& f1, Z2_field_element& f2) { std::swap(f1.element_, f2.element_); } + + /** + * @brief Casts the element into an unsigned int. + */ + operator unsigned int() const; + + /** + * @brief Returns the inverse of the element. + * + * @return The inverse. + */ + Z2_field_element get_inverse() const; + /** + * @brief For interface purposes with multi-fields. Returns the inverse together with the second argument. + * + * @param productOfCharacteristics Some value. + * @return Pair whose first element is the inverse of @p e and the second element is @p productOfCharacteristics. + */ + std::pair get_partial_inverse(unsigned int productOfCharacteristics) const; + + /** + * @brief Returns the additive identity of the field. + * + * @return false. + */ + static Z2_field_element get_additive_identity(); + /** + * @brief Returns the multiplicative identity of the field. + * + * @return true. + */ + static Z2_field_element get_multiplicative_identity(); + /** + * @brief For interface purposes with multi-fields. Returns the multiplicative identity of the field. + * + * @param productOfCharacteristics Some value. + * @return true. + */ + static Z2_field_element get_partial_multiplicative_identity([[maybe_unused]] unsigned int productOfCharacteristics); + /** + * @brief Returns the characteristic of the field, that is `2`. + * + * @return 2. + */ + static constexpr unsigned int get_characteristic(); + + /** + * @brief Returns the value of the element. + * + * @return Value of the element. + */ + unsigned int get_value() const; + + // static constexpr bool handles_only_z2() { return true; } + + private: + bool element_; }; -inline Z2_field_element::Z2_field_element() - : element_(false) -{} +inline Z2_field_element::Z2_field_element() : element_(false) {} -inline Z2_field_element::Z2_field_element(unsigned int element) - : element_(element & 1U) -{} +inline Z2_field_element::Z2_field_element(unsigned int element) : element_(element & 1U) {} -inline Z2_field_element::Z2_field_element(int element) - : element_(element & 1U) -{} +inline Z2_field_element::Z2_field_element(int element) : element_(element & 1U) {} -inline Z2_field_element::Z2_field_element(const Z2_field_element &toCopy) - : element_(toCopy.element_) -{} +inline Z2_field_element::Z2_field_element(const Z2_field_element& toCopy) : element_(toCopy.element_) {} -inline Z2_field_element::Z2_field_element(Z2_field_element &&toMove) noexcept - : element_(std::exchange(toMove.element_, 0)) -{} +inline Z2_field_element::Z2_field_element(Z2_field_element&& toMove) noexcept + : element_(std::exchange(toMove.element_, 0)) {} -inline Z2_field_element &Z2_field_element::operator=(Z2_field_element other) -{ - std::swap(element_, other.element_); - return *this; +inline Z2_field_element& Z2_field_element::operator=(Z2_field_element other) { + std::swap(element_, other.element_); + return *this; } -inline Z2_field_element &Z2_field_element::operator=(unsigned int const value) -{ - element_ = value & 1U; - return *this; +inline Z2_field_element& Z2_field_element::operator=(unsigned int const value) { + element_ = value & 1U; + return *this; } -inline Z2_field_element::operator unsigned int() const -{ - return element_; -} +inline Z2_field_element::operator unsigned int() const { return element_; } -inline Z2_field_element Z2_field_element::get_inverse() const -{ - return element_ ? Z2_field_element(1) : Z2_field_element(); +inline Z2_field_element Z2_field_element::get_inverse() const { + return element_ ? Z2_field_element(1) : Z2_field_element(); } -inline std::pair -Z2_field_element::get_partial_inverse(unsigned int productOfCharacteristics) const -{ - return {get_inverse(), productOfCharacteristics}; +inline std::pair Z2_field_element::get_partial_inverse( + unsigned int productOfCharacteristics) const { + return {get_inverse(), productOfCharacteristics}; } -inline Z2_field_element Z2_field_element::get_additive_identity() -{ - return Z2_field_element(); -} +inline Z2_field_element Z2_field_element::get_additive_identity() { return Z2_field_element(); } -inline Z2_field_element Z2_field_element::get_multiplicative_identity() -{ - return Z2_field_element(1); -} +inline Z2_field_element Z2_field_element::get_multiplicative_identity() { return Z2_field_element(1); } -inline Z2_field_element Z2_field_element::get_partial_multiplicative_identity([[maybe_unused]] unsigned int productOfCharacteristics) -{ - return Z2_field_element(1); +inline Z2_field_element Z2_field_element::get_partial_multiplicative_identity( + [[maybe_unused]] unsigned int productOfCharacteristics) { + return Z2_field_element(1); } -inline constexpr unsigned int Z2_field_element::get_characteristic() -{ - return 2; -} +inline constexpr unsigned int Z2_field_element::get_characteristic() { return 2; } -inline unsigned int Z2_field_element::get_value() const -{ - return element_; -} +inline unsigned int Z2_field_element::get_value() const { return element_; } -} //namespace persistence_fields -} //namespace Gudhi +} // namespace persistence_fields +} // namespace Gudhi #endif // MATRIX_FIELD_Z2_H_ diff --git a/src/Persistence_matrix/include/gudhi/Fields/Z2_field_operators.h b/src/Persistence_matrix/include/gudhi/Fields/Z2_field_operators.h index c47b86f6ad..68551dbdab 100644 --- a/src/Persistence_matrix/include/gudhi/Fields/Z2_field_operators.h +++ b/src/Persistence_matrix/include/gudhi/Fields/Z2_field_operators.h @@ -8,6 +8,12 @@ * - YYYY/MM Author: Description of the modification */ +/** + * @file Z2_field_operators.h + * @author Hannah Schreiber + * @brief Contains the @ref Z2_field_operators class. + */ + #ifndef MATRIX_FIELD_Z2_OPERATORS_H_ #define MATRIX_FIELD_Z2_OPERATORS_H_ @@ -17,117 +23,212 @@ namespace Gudhi { namespace persistence_fields { /** + * @class Z2_field_operators Z2_field_operators.h gudhi/Fields/Z2_field_operators.h * @ingroup persistence_fields * - * @brief TODO: - * + * @brief Class defining operators for the \f$ \mathbb{F}_2 \f$ field. */ -class Z2_field_operators { -public: - using element_type = bool; - using characteristic_type = unsigned int; - template - using isUnsignedInteger = std::enable_if_t >; - template - using isInteger = std::enable_if_t >; - - Z2_field_operators(){}; - - static constexpr characteristic_type get_characteristic(){ - return 2; - } - - template > - static element_type get_value(Integer_type e){ - if constexpr (std::is_same_v){ - return e; - } else { - return e < 2 && e >= 0 ? e : e % 2; //returns bool, so %2 won't be negative and is optimized to & - } - } - - //r = e1 + e2 - template > - static element_type add(Unsigned_integer_type e1, Unsigned_integer_type e2){ - if constexpr (std::is_same_v){ - return e1 != e2; - } else { - return get_value(e1) != get_value(e2); - } - } - - //r = e1 - e2 - template > - static element_type substract(Unsigned_integer_type e1, Unsigned_integer_type e2){ - if constexpr (std::is_same_v){ - return e1 != e2; - } else { - return get_value(e1) != get_value(e2); - } - } - - //r = e1 * e2 - template > - static element_type multiply(Unsigned_integer_type e1, Unsigned_integer_type e2){ - if constexpr (std::is_same_v){ - return e1 && e2; - } else { - return get_value(e1) ? get_value(e2) : false; - } - } - - //r = e * m + a - template > - static element_type multiply_and_add(Unsigned_integer_type e, Unsigned_integer_type m, Unsigned_integer_type a){ - if constexpr (std::is_same_v){ - return (e && m) != a; - } else { - return multiply(e, m) != get_value(a); - } - } - - //r = (e + a) * m - template > - static element_type add_and_multiply(Unsigned_integer_type e, Unsigned_integer_type a, Unsigned_integer_type m){ - if constexpr (std::is_same_v){ - return (e != a) && m; - } else { - return add(e, a) ? get_value(m) : false; - } - } - - template > - static bool are_equal(Unsigned_integer_type e1, Unsigned_integer_type e2){ - if constexpr (std::is_same_v){ - return e1 == e2; - } else { - return get_value(e1) == get_value(e2); - } - } - - template > - static element_type get_inverse(Unsigned_integer_type e){ - if constexpr (std::is_same_v){ - return e; - } else { - return get_value(e); - } - } - template > - static std::pair get_partial_inverse(Unsigned_integer_type e, characteristic_type productOfCharacteristics){ - return {get_inverse(e), productOfCharacteristics}; - } - - static constexpr element_type get_additive_identity(){ return false; } - static constexpr element_type get_multiplicative_identity(){ return true; } - static constexpr element_type get_partial_multiplicative_identity([[maybe_unused]] characteristic_type productOfCharacteristics){ return true; } - - static constexpr bool handles_only_z2(){ - return true; - } +class Z2_field_operators +{ + public: + using element_type = bool; /**< Type for the elements in the field. */ + using characteristic_type = unsigned int; /**< Type for the field characteristic. */ + template + using isUnsignedInteger = std::enable_if_t >; + template + using isInteger = std::enable_if_t >; + + /** + * @brief Default constructor. + */ + Z2_field_operators(){}; + + /** + * @brief Returns the characteristic of the field, that is `2`. + * + * @return 2. + */ + static constexpr characteristic_type get_characteristic() { return 2; } + + /** + * @brief Returns the value of an integer in the field. + * That is the positive value of the integer modulo the current characteristic. + * + * @tparam Integer_type A native integer type: int, unsigned int, long int, bool, etc. + * @param e Integer to return the value from. + * @return A boolean representing `e % 2`. + */ + template > + static element_type get_value(Integer_type e) { + if constexpr (std::is_same_v) { + return e; + } else { + return e < 2 && e >= 0 ? e : e % 2; // returns bool, so %2 won't be negative and is optimized to & + } + } + + /** + * @brief Returns the sum of two elements in the field. + * + * @tparam Unsigned_integer_type A native unsigned integer type: unsigned int, bool, etc. + * @param e1 First element. + * @param e2 Second element. + * @return `(e1 + e2) % 2` as a boolean. + */ + template > + static element_type add(Unsigned_integer_type e1, Unsigned_integer_type e2) { + if constexpr (std::is_same_v) { + return e1 != e2; + } else { + return get_value(e1) != get_value(e2); + } + } + + /** + * @brief Returns the substraction in the field of the first element by the second element. + * + * @tparam Unsigned_integer_type A native unsigned integer type: unsigned int, bool, etc. + * @param e1 First element. + * @param e2 Second element. + * @return `(e1 - e2) % 2` as a boolean. + */ + template > + static element_type substract(Unsigned_integer_type e1, Unsigned_integer_type e2) { + if constexpr (std::is_same_v) { + return e1 != e2; + } else { + return get_value(e1) != get_value(e2); + } + } + + /** + * @brief Returns the multiplication of two elements in the field. + * + * @tparam Unsigned_integer_type A native unsigned integer type: unsigned int, bool, etc. + * @param e1 First element. + * @param e2 Second element. + * @return `(e1 * e2) % 2` as a boolean. + */ + template > + static element_type multiply(Unsigned_integer_type e1, Unsigned_integer_type e2) { + if constexpr (std::is_same_v) { + return e1 && e2; + } else { + return get_value(e1) ? get_value(e2) : false; + } + } + + /** + * @brief Multiplies the first element with the second one and adds the third one. Returns the result in the field. + * + * @tparam Unsigned_integer_type A native unsigned integer type: unsigned int, bool, etc. + * @param e First element. + * @param m Second element. + * @param a Third element. + * @return `(e * m + a) % 2` as a boolean. + */ + template > + static element_type multiply_and_add(Unsigned_integer_type e, Unsigned_integer_type m, Unsigned_integer_type a) { + if constexpr (std::is_same_v) { + return (e && m) != a; + } else { + return multiply(e, m) != get_value(a); + } + } + + /** + * @brief Adds the first element to the second one and multiplies the third one with it. + * Returns the result in the field. + * + * @tparam Unsigned_integer_type A native unsigned integer type: unsigned int, bool, etc. + * @param e First element. + * @param a Second element. + * @param m Third element. + * @return `((e + a) * m) % 2` as a boolean. + */ + template > + static element_type add_and_multiply(Unsigned_integer_type e, Unsigned_integer_type a, Unsigned_integer_type m) { + if constexpr (std::is_same_v) { + return (e != a) && m; + } else { + return add(e, a) ? get_value(m) : false; + } + } + + /** + * @brief Returns true if the two given elements are equal in the field, false otherwise. + * + * @tparam Unsigned_integer_type A native unsigned integer type: unsigned int, bool, etc. + * @param e1 First element to compare. + * @param e2 Second element to compare. + * @return true If `e1 % 2 == e2 % 2`. + * @return false Otherwise. + */ + template > + static bool are_equal(Unsigned_integer_type e1, Unsigned_integer_type e2) { + if constexpr (std::is_same_v) { + return e1 == e2; + } else { + return get_value(e1) == get_value(e2); + } + } + + /** + * @brief Returns the inverse of the given element in the field. + * + * @tparam Unsigned_integer_type A native unsigned integer type: unsigned int, bool, etc. + * @param e Element to get the inverse from. + * @return Inverse in the current field of `e % 2`. + */ + template > + static element_type get_inverse(Unsigned_integer_type e) { + if constexpr (std::is_same_v) { + return e; + } else { + return get_value(e); + } + } + /** + * @brief For interface purposes with multi-fields. Returns the inverse together with the second argument. + * + * @tparam Unsigned_integer_type A native unsigned integer type: unsigned int, bool, etc. + * @param e Element to get the inverse from. + * @param productOfCharacteristics Some value. + * @return Pair whose first element is the inverse of @p e and the second element is @p productOfCharacteristics. + */ + template > + static std::pair get_partial_inverse( + Unsigned_integer_type e, characteristic_type productOfCharacteristics) { + return {get_inverse(e), productOfCharacteristics}; + } + + /** + * @brief Returns the additive identity of the field. + * + * @return false. + */ + static constexpr element_type get_additive_identity() { return false; } + /** + * @brief Returns the multiplicative identity of the field. + * + * @return true. + */ + static constexpr element_type get_multiplicative_identity() { return true; } + /** + * @brief For interface purposes with multi-fields. Returns the multiplicative identity of the field. + * + * @param productOfCharacteristics Some value. + * @return true. + */ + static constexpr element_type get_partial_multiplicative_identity( + [[maybe_unused]] characteristic_type productOfCharacteristics) { + return true; + } + + // static constexpr bool handles_only_z2() { return true; } }; -} //namespace persistence_fields -} //namespace Gudhi +} // namespace persistence_fields +} // namespace Gudhi #endif // MATRIX_FIELD_Z2_OPERATORS_H_ diff --git a/src/Persistence_matrix/include/gudhi/Fields/Zp_field_operators.h b/src/Persistence_matrix/include/gudhi/Fields/Zp_field_operators.h index 701f49c087..dc7a0deb0e 100644 --- a/src/Persistence_matrix/include/gudhi/Fields/Zp_field_operators.h +++ b/src/Persistence_matrix/include/gudhi/Fields/Zp_field_operators.h @@ -17,6 +17,7 @@ #include namespace Gudhi { +/// Field namespace namespace persistence_fields { /** diff --git a/src/Persistence_matrix/include/gudhi/Persistence_matrix/base_matrix_with_column_compression.h b/src/Persistence_matrix/include/gudhi/Persistence_matrix/base_matrix_with_column_compression.h index 2225335dcb..c46ec645a1 100644 --- a/src/Persistence_matrix/include/gudhi/Persistence_matrix/base_matrix_with_column_compression.h +++ b/src/Persistence_matrix/include/gudhi/Persistence_matrix/base_matrix_with_column_compression.h @@ -8,7 +8,7 @@ * - YYYY/MM Author: Description of the modification */ - /** +/** * @file base_matrix_with_column_compression.h * @author Hannah Schreiber * @brief Contains the @ref Base_matrix_with_column_compression class. @@ -325,7 +325,7 @@ class Base_matrix_with_column_compression : protected Master_matrix::Matrix_row_ */ friend void swap(Base_matrix_with_column_compression& matrix1, Base_matrix_with_column_compression& matrix2) { matrix1.columnToRep_.swap(matrix2.columnToRep_); - swap(matrix1.columnClasses_, matrix2.columnClasses_); + std::swap(matrix1.columnClasses_, matrix2.columnClasses_); matrix1.repToColumn_.swap(matrix2.repToColumn_); // be careful when columnPool_ becomes not static std::swap(matrix1.nextColumnIndex_, matrix2.nextColumnIndex_); std::swap(matrix1.operators_, matrix2.operators_); diff --git a/src/Persistence_matrix/include/gudhi/Persistence_matrix/overlay_ididx_to_matidx.h b/src/Persistence_matrix/include/gudhi/Persistence_matrix/overlay_ididx_to_matidx.h index 0d5740888a..8298c7e65a 100644 --- a/src/Persistence_matrix/include/gudhi/Persistence_matrix/overlay_ididx_to_matidx.h +++ b/src/Persistence_matrix/include/gudhi/Persistence_matrix/overlay_ididx_to_matidx.h @@ -9,7 +9,7 @@ */ /** - * @file overlay_@ref IDIdx_to_@ref MatIdx.h + * @file overlay_ididx_to_matidx.h * @author Hannah Schreiber * @brief Contains the @ref Id_to_index_overlay class. */ diff --git a/src/Persistence_matrix/include/gudhi/Persistence_matrix/overlay_posidx_to_matidx.h b/src/Persistence_matrix/include/gudhi/Persistence_matrix/overlay_posidx_to_matidx.h index a29f374b91..4a3085ea8c 100644 --- a/src/Persistence_matrix/include/gudhi/Persistence_matrix/overlay_posidx_to_matidx.h +++ b/src/Persistence_matrix/include/gudhi/Persistence_matrix/overlay_posidx_to_matidx.h @@ -9,7 +9,7 @@ */ /** - * @file overlay_@ref PosIdx_to_@ref MatIdx.h + * @file overlay_posidx_to_matidx.h * @author Hannah Schreiber * @brief Contains the @ref Position_to_index_overlay class. */