Loading...
Searching...
No Matches
Multi_field_small_operators.h
Go to the documentation of this file.
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.
3 * Author(s): Hannah Schreiber, Clément Maria
4 *
5 * Copyright (C) 2022-24 Inria
6 *
7 * Modification(s):
8 * - YYYY/MM Author: Description of the modification
9 */
10
17#ifndef MATRIX_FIELD_MULTI_SMALL_OPERATORS_H_
18#define MATRIX_FIELD_MULTI_SMALL_OPERATORS_H_
19
20#include <utility>
21#include <vector>
22#include <limits.h>
23#include <stdexcept>
24#include <numeric>
25
26namespace Gudhi {
27namespace persistence_fields {
28
38{
39 public:
40 using element_type = unsigned int;
46 Multi_field_operators_with_small_characteristics() : productOfAllCharacteristics_(0) /* , multiplicativeID_(1) */
47 {}
55 Multi_field_operators_with_small_characteristics(int minCharacteristic, int maxCharacteristic)
56 : productOfAllCharacteristics_(0) //, multiplicativeID_(1)
57 {
58 set_characteristic(minCharacteristic, maxCharacteristic);
59 }
66 : primes_(toCopy.primes_),
67 productOfAllCharacteristics_(toCopy.productOfAllCharacteristics_),
68 partials_(toCopy.partials_) /* ,
69 multiplicativeID_(toCopy.multiplicativeID_) */
70 {}
77 : primes_(std::move(toMove.primes_)),
78 productOfAllCharacteristics_(std::move(toMove.productOfAllCharacteristics_)),
79 partials_(std::move(toMove.partials_)) /* ,
80 multiplicativeID_(std::move(toMove.multiplicativeID_)) */
81 {}
82
91 void set_characteristic(int minimum, int maximum) {
92 if (maximum < 2) throw std::invalid_argument("Characteristic must be strictly positive");
93 if (minimum > maximum) throw std::invalid_argument("The given interval is not valid.");
94 if (minimum == maximum && !_is_prime(minimum))
95 throw std::invalid_argument("The given interval does not contain a prime number.");
96
97 productOfAllCharacteristics_ = 1;
98 primes_.clear();
99 for (unsigned int i = minimum; i <= static_cast<unsigned int>(maximum); ++i) {
100 if (_is_prime(i)) {
101 primes_.push_back(i);
102 productOfAllCharacteristics_ *= i;
103 }
104 }
105
106 if (primes_.empty()) throw std::invalid_argument("The given interval does not contain a prime number.");
107
108 partials_.resize(primes_.size());
109 for (unsigned int i = 0; i < primes_.size(); ++i) {
110 unsigned int p = primes_[i];
111 characteristic_type base = productOfAllCharacteristics_ / p;
112 unsigned int exp = p - 1;
113 partials_[i] = 1;
114
115 while (exp > 0) {
116 // If exp is odd, multiply with result
117 if (exp & 1) partials_[i] = _multiply(partials_[i], base, productOfAllCharacteristics_);
118 // y must be even now
119 exp = exp >> 1; // y = y/2
120 base = _multiply(base, base, productOfAllCharacteristics_);
121 }
122 }
123
124 // If I understood the paper well, multiplicativeID_ always equals to 1. But in Clement's code,
125 // multiplicativeID_ is computed (see commented loop below). TODO: verify with Clement.
126 // for (unsigned int i = 0; i < partials_.size(); ++i){
127 // multiplicativeID_ = (multiplicativeID_ + partials_[i]) % productOfAllCharacteristics_;
128 // }
129 }
135 const characteristic_type& get_characteristic() const { return productOfAllCharacteristics_; }
136
145 return e < productOfAllCharacteristics_ ? e : e % productOfAllCharacteristics_;
146 }
147
156 return _add(get_value(e1), get_value(e2), productOfAllCharacteristics_);
157 }
158
167 e1 = _add(get_value(e1), get_value(e2), productOfAllCharacteristics_);
168 }
169
178 return _substract(get_value(e1), get_value(e2), productOfAllCharacteristics_);
179 }
180
189 e1 = _substract(get_value(e1), get_value(e2), productOfAllCharacteristics_);
190 }
199 e2 = _substract(get_value(e1), get_value(e2), productOfAllCharacteristics_);
200 }
201
210 return _multiply(get_value(e1), get_value(e2), productOfAllCharacteristics_);
211 }
212
221 e1 = _multiply(get_value(e1), get_value(e2), productOfAllCharacteristics_);
222 }
223
235
248 e = get_value(e * m + a);
249 }
262 a = get_value(e * m + a);
263 }
264
277
290 e = get_value((e + a) * m);
291 }
304 m = get_value((e + a) * m);
305 }
306
315 bool are_equal(element_type e1, element_type e2) const { return get_value(e1) == get_value(e2); }
316
325 return get_partial_inverse(e, productOfAllCharacteristics_).first;
326 }
335 std::pair<element_type, characteristic_type> get_partial_inverse(
336 const element_type& e, const characteristic_type& productOfCharacteristics) const {
337 characteristic_type gcd = std::gcd(e, productOfAllCharacteristics_);
338
339 if (gcd == productOfCharacteristics) return {0, get_multiplicative_identity()}; // partial inverse is 0
340
341 characteristic_type QT = productOfCharacteristics / gcd;
342
343 const characteristic_type inv_qt = _get_inverse(e, QT);
344
346 res = _multiply(res, inv_qt, productOfAllCharacteristics_);
347
348 return {res, QT};
349 }
350
356 static constexpr element_type get_additive_identity() { return 0; }
362 static constexpr element_type get_multiplicative_identity() { return 1; }
363 // static element_type get_multiplicative_identity(){ return multiplicativeID_; }
372 if (productOfCharacteristics == 0) {
374 }
375 element_type multIdentity = 0;
376 for (unsigned int idx = 0; idx < primes_.size(); ++idx) {
377 if ((productOfCharacteristics % primes_[idx]) == 0) {
378 multIdentity = _add(multIdentity, partials_[idx], productOfAllCharacteristics_);
379 }
380 }
381 return multIdentity;
382 }
383
384 // static constexpr bool handles_only_z2() { return false; }
385
390 primes_.swap(other.primes_);
391 productOfAllCharacteristics_ = other.productOfAllCharacteristics_;
392 partials_.swap(other.partials_);
393
394 return *this;
395 }
401 f1.primes_.swap(f2.primes_);
402 std::swap(f1.productOfAllCharacteristics_, f2.productOfAllCharacteristics_);
403 f1.partials_.swap(f2.partials_);
404 }
405
406 private:
407 std::vector<unsigned int> primes_;
408 characteristic_type productOfAllCharacteristics_;
409 std::vector<characteristic_type> partials_;
410 // static inline constexpr unsigned int multiplicativeID_ = 1;
411
412 static element_type _add(element_type element, element_type v, characteristic_type characteristic);
413 static element_type _substract(element_type element, element_type v, characteristic_type characteristic);
414 static element_type _multiply(element_type a, element_type b, characteristic_type characteristic);
415 static constexpr long int _get_inverse(element_type element, characteristic_type mod);
416 static constexpr bool _is_prime(const int p);
417};
418
420Multi_field_operators_with_small_characteristics::_add(element_type element, element_type v,
421 characteristic_type characteristic) {
422 if (UINT_MAX - element < v) {
423 // automatic unsigned integer overflow behaviour will make it work
424 element += v;
425 element -= characteristic;
426 return element;
427 }
428
429 element += v;
430 if (element >= characteristic) element -= characteristic;
431
432 return element;
433}
434
436Multi_field_operators_with_small_characteristics::_substract(element_type element, element_type v,
437 characteristic_type characteristic) {
438 if (element < v) {
439 element += characteristic;
440 }
441 element -= v;
442
443 return element;
444}
445
447Multi_field_operators_with_small_characteristics::_multiply(element_type a, element_type b,
448 characteristic_type characteristic) {
449 element_type res = 0;
450 element_type temp_b = 0;
451
452 if (b < a) std::swap(a, b);
453
454 while (a != 0) {
455 if (a & 1) {
456 /* Add b to res, modulo m, without overflow */
457 if (b >= characteristic - res) res -= characteristic;
458 res += b;
459 }
460 a >>= 1;
461
462 /* Double b, modulo m */
463 temp_b = b;
464 if (b >= characteristic - b) temp_b -= characteristic;
465 b += temp_b;
466 }
467 return res;
468}
469
470inline constexpr long int Multi_field_operators_with_small_characteristics::_get_inverse(element_type element,
471 characteristic_type mod) {
472 // to solve: Ax + My = 1
473 element_type M = mod;
474 element_type A = element;
475 long int y = 0, x = 1;
476 // extended euclidien division
477 while (A > 1) {
478 int quotient = A / M;
479 int temp = M;
480
481 M = A % M, A = temp;
482 temp = y;
483
484 y = x - quotient * y;
485 x = temp;
486 }
487
488 if (x < 0) x += mod;
489
490 return x;
491}
492
493inline constexpr bool Multi_field_operators_with_small_characteristics::_is_prime(const int p) {
494 if (p <= 1) return false;
495 if (p <= 3) return true;
496 if (p % 2 == 0 || p % 3 == 0) return false;
497
498 for (long i = 5; i * i <= p; i = i + 6)
499 if (p % i == 0 || p % (i + 2) == 0) return false;
500
501 return true;
502}
503
504} // namespace persistence_fields
505} // namespace Gudhi
506
507#endif // MATRIX_FIELD_MULTI_SMALL_OPERATORS_H_
Class defining operators for a multi-field with "consecutive" charateristic range,...
Definition Multi_field_small_operators.h:38
static constexpr element_type get_multiplicative_identity()
Returns the multiplicative identity of a field.
Definition Multi_field_small_operators.h:362
element_type get_value(element_type e) const
Returns the value of an element in the field. That is the positive value of the integer modulo the cu...
Definition Multi_field_small_operators.h:144
element_type add(element_type e1, element_type e2) const
Returns the sum of two elements in the field.
Definition Multi_field_small_operators.h:155
element_type substract(element_type e1, element_type e2) const
Returns the substraction in the field of the first element by the second element.
Definition Multi_field_small_operators.h:177
element_type multiply_and_add(element_type e, element_type m, element_type a) const
Multiplies the first element with the second one and adds the third one. Returns the result in the fi...
Definition Multi_field_small_operators.h:234
unsigned int element_type
Definition Multi_field_small_operators.h:40
void add_and_multiply_inplace_front(element_type &e, element_type a, element_type m) const
Adds the first element to the second one and multiplies the third one with it, that is ((e + a) * m) ...
Definition Multi_field_small_operators.h:289
element_type characteristic_type
Definition Multi_field_small_operators.h:41
Multi_field_operators_with_small_characteristics(const Multi_field_operators_with_small_characteristics &toCopy)
Copy constructor.
Definition Multi_field_small_operators.h:65
void add_inplace(element_type &e1, element_type e2) const
Stores in the first element the sum of two given elements in the field, that is (e1 + e2) % productOf...
Definition Multi_field_small_operators.h:166
void multiply_inplace(element_type &e1, element_type e2) const
Stores in the first element the multiplication of two given elements in the field,...
Definition Multi_field_small_operators.h:220
void multiply_and_add_inplace_front(element_type &e, element_type m, element_type a) const
Multiplies the first element with the second one and adds the third one, that is (e * m + a) % produc...
Definition Multi_field_small_operators.h:247
Multi_field_operators_with_small_characteristics(int minCharacteristic, int maxCharacteristic)
Constructor setting the characteristics to all prime numbers between the two given integers....
Definition Multi_field_small_operators.h:55
void substract_inplace_back(element_type e1, element_type &e2) const
Stores in the second element the substraction in the field of the first element by the second element...
Definition Multi_field_small_operators.h:198
element_type multiply(element_type e1, element_type e2) const
Returns the multiplication of two elements in the field.
Definition Multi_field_small_operators.h:209
void multiply_and_add_inplace_back(element_type e, element_type m, element_type &a) const
Multiplies the first element with the second one and adds the third one, that is (e * m + a) % produc...
Definition Multi_field_small_operators.h:261
Multi_field_operators_with_small_characteristics(Multi_field_operators_with_small_characteristics &&toMove) noexcept
Move constructor.
Definition Multi_field_small_operators.h:76
void set_characteristic(int minimum, int maximum)
Set the characteristics of the field, which are stored in a single value as a product of all of them....
Definition Multi_field_small_operators.h:91
element_type add_and_multiply(element_type e, element_type a, element_type m) const
Adds the first element to the second one and multiplies the third one with it. Returns the result in ...
Definition Multi_field_small_operators.h:276
const characteristic_type & get_characteristic() const
Returns the current characteristics as the product of all of them.
Definition Multi_field_small_operators.h:135
Multi_field_operators_with_small_characteristics & operator=(Multi_field_operators_with_small_characteristics other)
Assign operator.
Definition Multi_field_small_operators.h:389
element_type get_inverse(const element_type &e) const
Returns the inverse of the given element in the sense of with respect to the product of all characte...
Definition Multi_field_small_operators.h:324
void substract_inplace_front(element_type &e1, element_type e2) const
Stores in the first element the substraction in the field of the first element by the second element,...
Definition Multi_field_small_operators.h:188
bool are_equal(element_type e1, element_type e2) const
Returns true if the two given elements are equal in the field, false otherwise.
Definition Multi_field_small_operators.h:315
friend void swap(Multi_field_operators_with_small_characteristics &f1, Multi_field_operators_with_small_characteristics &f2)
Swap operator.
Definition Multi_field_small_operators.h:399
element_type get_partial_multiplicative_identity(const characteristic_type &productOfCharacteristics) const
Returns the partial multiplicative identity of the multi-field from the given product....
Definition Multi_field_small_operators.h:371
Multi_field_operators_with_small_characteristics()
Default constructor, sets the product of all characteristics to 0.
Definition Multi_field_small_operators.h:46
static constexpr element_type get_additive_identity()
Returns the additive identity of a field.
Definition Multi_field_small_operators.h:356
std::pair< element_type, characteristic_type > get_partial_inverse(const element_type &e, const characteristic_type &productOfCharacteristics) const
Returns the inverse of the given element in the multi-field corresponding to the given sub-product of...
Definition Multi_field_small_operators.h:335
void add_and_multiply_inplace_back(element_type e, element_type a, element_type &m) const
Adds the first element to the second one and multiplies the third one with it, that is ((e + a) * m) ...
Definition Multi_field_small_operators.h:303
Gudhi namespace.
Definition SimplicialComplexForAlpha.h:14