diophantine.cpp
2 #include <symengine/constants.h>
4 #include <symengine/mul.h>
5
6 namespace SymEngine
7 {
8
9 bool order(const DenseMatrix &t, const std::vector<DenseMatrix> &basis,
10  unsigned k)
11 {
12  bool eq = true;
13
14  for (unsigned j = 0; j < t.ncols(); j++) {
15  SYMENGINE_ASSERT(is_a<Integer>(*t.get(0, j)));
16  integer_class t_
17  = down_cast<const Integer &>(*t.get(0, j)).as_integer_class();
18
19  SYMENGINE_ASSERT(is_a<Integer>(*basis[k].get(0, j)));
20  integer_class b_ = down_cast<const Integer &>(*basis[k].get(0, j))
21  .as_integer_class();
22
23  if (t_ < b_) {
24  return false;
25  }
26  if (t_ > b_) {
27  eq = false;
28  }
29  }
30
31  return not eq;
32 }
33
34 bool is_minimum(const DenseMatrix &t, const std::vector<DenseMatrix> &basis,
35  unsigned n)
36 {
37  if (n == 0) {
38  return true;
39  }
40
41  return not order(t, basis, n - 1) and is_minimum(t, basis, n - 1);
42 }
43
44 // Solve the diophantine system Ax = 0 and return a basis set for solutions
45 // Reference:
46 // Evelyne Contejean, Herve Devie. An Efficient Incremental Algorithm for
47 // Solving
48 // Systems of Linear Diophantine Equations. Information and computation,
49 // 113(1):143-172,
50 // August 1994.
51 void homogeneous_lde(std::vector<DenseMatrix> &basis, const DenseMatrix &A)
52 {
53  unsigned p = A.nrows();
54  unsigned q = A.ncols();
55
56  SYMENGINE_ASSERT(p > 0 and q > 1);
57
58  DenseMatrix row_zero(1, q);
59  zeros(row_zero);
60
61  DenseMatrix col_zero(p, 1);
62  zeros(col_zero);
63
65  P.push_back(row_zero);
66
68  for (unsigned j = 0; j < q; j++) {
69  Frozen[0][j] = false;
70  }
71
72  std::vector<bool> F(q, false);
73
74  DenseMatrix t, transpose, product, T;
75  RCP<const Integer> dot;
76
77  product = DenseMatrix(p, 1);
78  transpose = DenseMatrix(q, 1);
79
80  while (P.size() > 0) {
81  auto n = P.size() - 1;
82  t = P[n];
83  P.pop_back();
84
85  t.transpose(transpose);
86  A.mul_matrix(transpose, product);
87
88  if (product.eq(col_zero) and not t.eq(row_zero)) {
89  basis.push_back(t);
90  } else {
91  for (unsigned i = 0; i < q; i++) {
92  F[i] = Frozen[n][i];
93  }
94
95  T = t;
96  for (unsigned i = 0; i < q; i++) {
97  SYMENGINE_ASSERT(is_a<Integer>(*T.get(0, i)));
98  T.set(0, i,
100
101  if (i > 0) {
102  SYMENGINE_ASSERT(is_a<Integer>(*T.get(0, i - 1)));
103  T.set(0, i - 1,
104  down_cast<const Integer &>(*T.get(0, i - 1))
105  .subint(*one));
106  }
107
108  dot = zero;
109  for (unsigned j = 0; j < p; j++) {
110  SYMENGINE_ASSERT(is_a<Integer>(*product.get(j, 0)));
111  RCP<const Integer> p_j0
112  = rcp_static_cast<const Integer>(product.get(j, 0));
113
114  SYMENGINE_ASSERT(is_a<Integer>(*A.get(j, i)));
115  RCP<const Integer> A_ji
116  = rcp_static_cast<const Integer>(A.get(j, i));
117
119  }
120
121  if (F[i] == false
122  and ((dot->is_negative()
123  and is_minimum(T, basis,
124  numeric_cast<unsigned>(basis.size())))
125  or t.eq(row_zero))) {
126  P.push_back(T);
127  n = n + 1;
128
129  for (unsigned j = 0; j < q; j++) {
130  Frozen[n - 1][j] = F[j];
131  }
132
133  F[i] = true;
134  }
135  }
136  }
137  }
138 }
139 } // namespace SymEngine
Classes and functions relating to the binary operation of addition.
Main namespace for SymEngine package.