Zp_field_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
 *
 *    Copyright (C) 2024 Inria
 *
 *    Modification(s):
 *      - YYYY/MM Author: Description of the modification
 */

 
#ifndef MATRIX_FIELD_ZP_OPERATOR_H_
#define MATRIX_FIELD_ZP_OPERATOR_H_
 
#include <stdexcept>
#include <utility>
#include <vector>
#include <limits.h>
 
namespace Gudhi {
namespace persistence_fields {

template <typename Unsigned_integer_type = unsigned int,
          class = std::enable_if_t<std::is_unsigned_v<Unsigned_integer_type> > >
class Zp_field_operators
{
 public:
  using Element = Unsigned_integer_type; 
  using Characteristic = Element;        
  template <class T>
  using isSignedInteger = std::enable_if_t<std::is_signed_v<T> >;

  constexpr static const Characteristic nullCharacteristic = 0;

  Zp_field_operators(Characteristic characteristic = nullCharacteristic) : characteristic_(nullCharacteristic)
  {
    if (characteristic != nullCharacteristic) set_characteristic(characteristic);
  }

  void set_characteristic(Characteristic characteristic)
  {
    if (characteristic <= 1) throw std::invalid_argument("Characteristic must be a strictly positive prime number.");
 
    inverse_.resize(characteristic);
    inverse_[0] = 0;
    for (unsigned int i = 1; i < characteristic; ++i) {
      unsigned int inv = 1;
      unsigned int mult = inv * i;
      while ((mult % characteristic) != 1) {
        ++inv;
        if (mult == characteristic) throw std::invalid_argument("Characteristic must be a prime number.");
        mult = inv * i;
      }
      inverse_[i] = inv;
    }
 
    characteristic_ = characteristic;
  }

  const Characteristic& get_characteristic() const { return characteristic_; }

  Element get_value(Element e) const { return e < characteristic_ ? e : e % characteristic_; }

  template <typename Signed_integer_type, class = isSignedInteger<Signed_integer_type> >
  Element get_value(Signed_integer_type e) const
  {
    if (e < -static_cast<Signed_integer_type>(characteristic_)) e = e % characteristic_;
    if (e < 0) return e += characteristic_;
    return e < static_cast<Signed_integer_type>(characteristic_) ? e : e % characteristic_;
  }

  Element add(Element e1, Element e2) const { return _add(get_value(e1), get_value(e2), characteristic_); }

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

  Element subtract(Element e1, Element e2) const { return _subtract(get_value(e1), get_value(e2), characteristic_); }

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

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

  Element multiply(Element e1, Element e2) const { return _multiply(get_value(e1), get_value(e2), characteristic_); }

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

  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); }

  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); }

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

  Element get_inverse(Element e) const { return inverse_[get_value(e)]; }

  std::pair<Element, Characteristic> get_partial_inverse(Element e, Characteristic productOfCharacteristics) const
  {
    return {get_inverse(e), productOfCharacteristics};
  }

  static constexpr Element get_additive_identity() { return 0; }

  static constexpr Element get_multiplicative_identity() { return 1; }

  static constexpr Element get_partial_multiplicative_identity([[maybe_unused]] Characteristic productOfCharacteristics)
  {
    return 1;
  }

  friend void swap(Zp_field_operators& f1, Zp_field_operators& f2) noexcept
  {
    std::swap(f1.characteristic_, f2.characteristic_);
    f1.inverse_.swap(f2.inverse_);
  }
 
 private:
  Characteristic characteristic_; 
  std::vector<Element> inverse_;  
 
  static Element _add(Element e1, Element e2, Characteristic characteristic)
  {
    if (UINT_MAX - e1 < e2) {
      // automatic unsigned integer overflow behaviour will make it work
      e1 += e2;
      e1 -= characteristic;
      return e1;
    }
 
    e1 += e2;
    if (e1 >= characteristic) e1 -= characteristic;
 
    return e1;
  }
 
  static Element _subtract(Element e1, Element e2, Characteristic characteristic)
  {
    if (e1 < e2) {
      e1 += characteristic;
    }
    e1 -= e2;
 
    return e1;
  }
 
  static Element _multiply(Element e1, Element e2, Characteristic characteristic)
  {
    unsigned int a = e1;
    e1 = 0;
    unsigned int temp_b;
 
    while (a != 0) {
      if (a & 1) {
        if (e2 >= characteristic - e1) e1 -= characteristic;
        e1 += e2;
      }
      a >>= 1;
 
      temp_b = e2;
      if (e2 >= characteristic - e2) temp_b -= characteristic;
      e2 += temp_b;
    }
 
    return e1;
  }
};
 
}  // namespace persistence_fields
}  // namespace Gudhi
 
#endif  // MATRIX_FIELD_ZP_OPERATOR_H_