Multi_field.h
1 /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2  * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
4  *
5  * Copyright (C) 2014 Inria
6  *
7  * Modification(s):
8  * - YYYY/MM Author: Description of the modification
9  */
10
11 #ifndef PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_
12 #define PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_
13
14 #include <gmpxx.h>
15
16 #include <vector>
17 #include <utility>
18
19 namespace Gudhi {
20
21 namespace persistent_cohomology {
22
31 class Multi_field {
32  public:
33  typedef mpz_class Element;
34
35  Multi_field()
36  : prod_characteristics_(0),
37  mult_id_all(0),
39  }
40
41  /* Initialize the multi-field. The generation of prime numbers might fail with
42  * a very small probability.*/
43  void init(int min_prime, int max_prime) {
44  if (max_prime < 2) {
45  std::cerr << "There is no prime less than " << max_prime << std::endl;
46  }
47  if (min_prime > max_prime) {
48  std::cerr << "No prime in [" << min_prime << ":" << max_prime << "]"
49  << std::endl;
50  }
51  // fill the list of prime numbers
52  int curr_prime = min_prime;
53  mpz_t tmp_prime;
54  mpz_init_set_ui(tmp_prime, min_prime);
55  // test if min_prime is prime
56  int is_prime = mpz_probab_prime_p(tmp_prime, 25); // probabilistic primality test
57
58  if (is_prime == 0) { // min_prime is composite
59  mpz_nextprime(tmp_prime, tmp_prime);
60  curr_prime = mpz_get_ui(tmp_prime);
61  }
62
63  while (curr_prime <= max_prime) {
64  primes_.push_back(curr_prime);
65  mpz_nextprime(tmp_prime, tmp_prime);
66  curr_prime = mpz_get_ui(tmp_prime);
67  }
68  mpz_clear(tmp_prime);
69  // set m to primorial(bound_prime)
70  prod_characteristics_ = 1;
71  for (auto p : primes_) {
72  prod_characteristics_ *= p;
73  }
74
75  // Uvect_
76  Element Ui;
77  Element tmp_elem;
78  for (auto p : primes_) {
79  assert(p > 0); // division by zero + non negative values
80  tmp_elem = prod_characteristics_ / p;
81  // Element tmp_elem_bis = 10;
82  mpz_powm_ui(tmp_elem.get_mpz_t(), tmp_elem.get_mpz_t(), p - 1,
83  prod_characteristics_.get_mpz_t());
84  Uvect_.push_back(tmp_elem);
85  }
86  mult_id_all = 0;
87  for (auto uvect : Uvect_) {
88  assert(prod_characteristics_ > 0); // division by zero + non negative values
89  mult_id_all = (mult_id_all + uvect) % prod_characteristics_;
90  }
91  }
92
94  const Element& additive_identity() const {
96  }
98  const Element& multiplicative_identity() const {
99  return mult_id_all;
100  } // 1 everywhere
101
102  Element multiplicative_identity(Element Q) {
103  if (Q == prod_characteristics_) {
104  return multiplicative_identity();
105  }
106
107  assert(prod_characteristics_ > 0); // division by zero + non negative values
108  Element mult_id = 0;
109  for (unsigned int idx = 0; idx < primes_.size(); ++idx) {
110  assert(primes_[idx] > 0); // division by zero + non negative values
111  if ((Q % primes_[idx]) == 0) {
112  mult_id = (mult_id + Uvect_[idx]) % prod_characteristics_;
113  }
114  }
115  return mult_id;
116  }
117
119  Element times(const Element& y, const Element& w) {
120  return plus_times_equal(0, y, w);
121  }
122
123  Element plus_equal(const Element& x, const Element& y) {
124  return plus_times_equal(x, y, (Element)1);
125  }
126
128  const Element& characteristic() const {
129  return prod_characteristics_;
130  }
131
133  std::pair<Element, Element> inverse(Element x, Element QS) {
134  Element QR;
135  mpz_gcd(QR.get_mpz_t(), x.get_mpz_t(), QS.get_mpz_t()); // QR <- gcd(x,QS)
136  if (QR == QS)
137  return std::pair<Element, Element>(additive_identity(), multiplicative_identity()); // partial inverse is 0
138  Element QT = QS / QR;
139  Element inv_qt;
140  mpz_invert(inv_qt.get_mpz_t(), x.get_mpz_t(), QT.get_mpz_t());
141
142  assert(prod_characteristics_ > 0); // division by zero + non negative values
143  return { (inv_qt * multiplicative_identity(QT)) % prod_characteristics_, QT };
144  }
146  Element times_minus(const Element& x, const Element& y) {
147  assert(prod_characteristics_ > 0); // division by zero + non negative values
148  /* This assumes that (x*y)%pc cannot be zero, but Field_Zp has specific code for the 0 case ??? */
149  return prod_characteristics_ - ((x * y) % prod_characteristics_);
150  }
151
153  Element plus_times_equal(const Element& x, const Element& y, const Element& w) {
154  assert(prod_characteristics_ > 0); // division by zero + non negative values
155  Element result = (x + w * y) % prod_characteristics_;
156  if (result < 0)
157  result += prod_characteristics_;
158  return result;
159  }
160
161  Element prod_characteristics_; // product of characteristics of the fields
162  // represented by the multi-field class
163  std::vector<int> primes_; // all the characteristics of the fields
164  std::vector<Element> Uvect_;
165  Element mult_id_all;
167 };
168
169 } // namespace persistent_cohomology
170
171 } // namespace Gudhi
172
173 #endif // PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_
Structure representing coefficients in a set of finite fields simultaneously using the chinese remain...
Definition: Multi_field.h:31
Element times(const Element &y, const Element &w)
Definition: Multi_field.h:119
const Element & characteristic() const
Returns the characteristic of the field.
Definition: Multi_field.h:128
const Element & multiplicative_identity() const
Returns the multiplicative identity of the field.
Definition: Multi_field.h:98
Element plus_times_equal(const Element &x, const Element &y, const Element &w)
Definition: Multi_field.h:153