Multi_field_small_operators.h

Go to the documentation of this file.

/*    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, Clément Maria
 *
 *    Copyright (C) 2022 Inria
 *
 *    Modification(s):
 *      - YYYY/MM Author: Description of the modification
 */

 
#ifndef MATRIX_FIELD_MULTI_SMALL_OPERATORS_H_
#define MATRIX_FIELD_MULTI_SMALL_OPERATORS_H_
 
#include <utility>
#include <vector>
#include <climits>
#include <stdexcept>
#include <numeric>
 
#include <gudhi/Debug_utils.h>
 
namespace Gudhi {
namespace persistence_fields {

class Multi_field_operators_with_small_characteristics
{
 public:
  using Element = unsigned int;   
  using Characteristic = Element; 
 
  inline static const Characteristic nullCharacteristic = 0; 

  Multi_field_operators_with_small_characteristics()
      : productOfAllCharacteristics_(nullCharacteristic) /* , multiplicativeID_(1) */
  {}

  Multi_field_operators_with_small_characteristics(int minCharacteristic, int maxCharacteristic)
      : productOfAllCharacteristics_(nullCharacteristic)  //, multiplicativeID_(1)
  {
    set_characteristic(minCharacteristic, maxCharacteristic);
  }

  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<unsigned int>(maximum); ++i) {
      if (_is_prime(static_cast<int>(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 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_;
    //  }
  }

  [[nodiscard]] const Characteristic& get_characteristic() const { return productOfAllCharacteristics_; }

  [[nodiscard]] Element get_value(Element e) const
  {
    return e < productOfAllCharacteristics_ ? e : e % productOfAllCharacteristics_;
  }

  [[nodiscard]] Element add(Element e1, Element e2) const
  {
    return _add(get_value(e1), get_value(e2), productOfAllCharacteristics_);
  }

  void add_inplace(Element& e1, Element e2) const
  {
    e1 = _add(get_value(e1), get_value(e2), productOfAllCharacteristics_);
  }

  [[nodiscard]] Element subtract(Element e1, Element e2) const
  {
    return _subtract(get_value(e1), get_value(e2), productOfAllCharacteristics_);
  }

  void subtract_inplace_front(Element& e1, Element e2) const
  {
    e1 = _subtract(get_value(e1), get_value(e2), productOfAllCharacteristics_);
  }

  void subtract_inplace_back(Element e1, Element& e2) const
  {
    e2 = _subtract(get_value(e1), get_value(e2), productOfAllCharacteristics_);
  }

  [[nodiscard]] Element multiply(Element e1, Element e2) const
  {
    return _multiply(get_value(e1), get_value(e2), productOfAllCharacteristics_);
  }

  void multiply_inplace(Element& e1, Element e2) const
  {
    e1 = _multiply(get_value(e1), get_value(e2), productOfAllCharacteristics_);
  }

  [[nodiscard]] Element multiply_and_add(Element e, Element m, Element a) const { return get_value((e * m) + a); }

  void multiply_and_add_inplace_front(Element& e, Element m, Element a) const { e = get_value((e * m) + a); }

  void multiply_and_add_inplace_back(Element e, Element m, Element& a) const { a = get_value((e * m) + a); }

  [[nodiscard]] Element add_and_multiply(Element e, Element a, Element m) const { return get_value((e + a) * m); }

  void add_and_multiply_inplace_front(Element& e, Element a, Element m) const { e = get_value((e + a) * m); }

  void add_and_multiply_inplace_back(Element e, Element a, Element& m) const { m = get_value((e + a) * m); }

  [[nodiscard]] bool are_equal(Element e1, Element e2) const { return get_value(e1) == get_value(e2); }

  [[nodiscard]] Element get_inverse(const Element& e) const
  {
    return get_partial_inverse(e, productOfAllCharacteristics_).first;
  }

  std::pair<Element, Characteristic> get_partial_inverse(const Element& e,
                                                         const Characteristic& productOfCharacteristics) const
  {
    GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_,
                std::invalid_argument(
                    "The given product is not the product of a subset of the current Multi-field characteristics."));
 
    Characteristic gcd = std::gcd(e, productOfAllCharacteristics_);
 
    if (gcd == productOfCharacteristics) return {0, get_multiplicative_identity()};  // partial inverse is 0
 
    Characteristic QT = productOfCharacteristics / gcd;
 
    const Characteristic 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 get_additive_identity() { return 0; }

  static constexpr Element get_multiplicative_identity() { return 1; }
 
  // static Element get_multiplicative_identity(){ return multiplicativeID_; }
  [[nodiscard]] Element get_partial_multiplicative_identity(const Characteristic& productOfCharacteristics) const
  {
    GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_,
                std::invalid_argument(
                    "The given product is not the product of a subset of the current Multi-field characteristics."));
 
    if (productOfCharacteristics == nullCharacteristic || productOfCharacteristics == productOfAllCharacteristics_) {
      return get_multiplicative_identity();
    }
    Element multIdentity = 0;
    for (unsigned int idx = 0; idx < primes_.size(); ++idx) {
      if ((productOfCharacteristics % primes_[idx]) == 0) {
        multIdentity = _add(multIdentity, partials_[idx], productOfAllCharacteristics_);
      }
    }
    return multIdentity;
  }

  friend void swap(Multi_field_operators_with_small_characteristics& f1,
                   Multi_field_operators_with_small_characteristics& f2) noexcept
  {
    f1.primes_.swap(f2.primes_);
    std::swap(f1.productOfAllCharacteristics_, f2.productOfAllCharacteristics_);
    f1.partials_.swap(f2.partials_);
  }
 
 private:
  std::vector<unsigned int> primes_;           
  Characteristic productOfAllCharacteristics_; 
  std::vector<Characteristic> partials_;       
 
  // static inline constexpr unsigned int multiplicativeID_ = 1;
 
  static Element _add(Element element, Element v, Characteristic 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;
 
    return element;
  }
 
  static Element _subtract(Element element, Element v, Characteristic characteristic)
  {
    if (element < v) {
      element += characteristic;
    }
    element -= v;
 
    return element;
  }
 
  static Element _multiply(Element a, Element b, Characteristic characteristic)
  {
    Element res = 0;
    Element 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;
  }
 
  static constexpr long int _get_inverse(Element element, Characteristic mod)
  {
    // to solve: Ax + My = 1
    Element M = mod;
    Element A = element;
    long int y = 0, x = 1;
    // extended euclidean 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;
  }
 
  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;
  }
};
 
}  // namespace persistence_fields
}  // namespace Gudhi
 
#endif  // MATRIX_FIELD_MULTI_SMALL_OPERATORS_H_