1 #include <cstdlib>
2 #include <vector>
3 #include <fstream>
4 #include <ctime>
5 #ifdef _OPENMP
6 #include <omp.h>
7 #endif
8 using namespace std;
9
10 #include "libnormaliz/libnormaliz.h"
11
12 using namespace libnormaliz;
13
14 typedef long long Integer;
15
rand_simplex(size_t dim,long bound)16 Cone<Integer> rand_simplex(size_t dim, long bound) {
17 vector<vector<Integer> > vertices(dim + 1, vector<Integer>(dim));
18 while (true) { // an eternal loop ...
19 for (size_t i = 0; i <= dim; ++i) {
20 for (size_t j = 0; j < dim; ++j)
21 vertices[i][j] = rand() % (bound + 1);
22 }
23 Cone<Integer> Simplex(Type::polytope, vertices);
24 // we must check the rank and normality
25 if (Simplex.getRank() == dim + 1 && Simplex.isDeg1HilbertBasis())
26 return Simplex;
27 }
28 vector<vector<Integer> > dummy_gen(1, vector<Integer>(1, 1)); // to make the compiler happy
29 return Cone<Integer>(Type::cone, dummy_gen);
30 }
31
exists_jump_over(Cone<Integer> & Polytope,const vector<vector<Integer>> & jump_cands)32 bool exists_jump_over(Cone<Integer>& Polytope, const vector<vector<Integer> >& jump_cands) {
33 vector<vector<Integer> > test_polytope = Polytope.getExtremeRays();
34 test_polytope.resize(test_polytope.size() + 1);
35 for (const auto& jump_cand : jump_cands) {
36 test_polytope[test_polytope.size() - 1] = jump_cand;
37 Cone<Integer> TestCone(Type::cone, test_polytope);
38 if (TestCone.getNrDeg1Elements() != Polytope.getNrDeg1Elements() + 1)
39 continue;
40 if (TestCone.isDeg1HilbertBasis())
41 return true;
42 }
43 return false;
44 }
45
lattice_widths(Cone<Integer> & Polytope)46 vector<Integer> lattice_widths(Cone<Integer>& Polytope) {
47 if (!Polytope.isDeg1ExtremeRays()) {
48 cerr << "Cone in lattice_widths is not defined by lattice polytope" << endl;
49 exit(1);
50 }
51 vector<Integer> widths(Polytope.getNrExtremeRays());
52 for (size_t i = 0; i < Polytope.getNrSupportHyperplanes(); ++i) {
53 widths[i] = 0;
54 for (size_t j = 0; j < Polytope.getNrExtremeRays(); ++j) {
55 // v_scalar_product is a useful function from vector_operations.h
56 Integer test = v_scalar_product(Polytope.getSupportHyperplanes()[i], Polytope.getExtremeRays()[j]);
57 if (test > widths[i])
58 widths[i] = test;
59 }
60 }
61 return widths;
62 }
63
main(int argc,char * argv[])64 int main(int argc, char* argv[]) {
65 time_t ticks;
66 srand(time(&ticks));
67 cout << "Seed " << ticks << endl; // we may want to reproduce the run
68
69 size_t polytope_dim = 4;
70 size_t cone_dim = polytope_dim + 1;
71 long bound = 6;
72 vector<Integer> grading(cone_dim, 0); // at some points we need the explicit grading
73 grading[polytope_dim] = 1;
74
75 size_t nr_simplex = 0; // for the progress report
76
77 while (true) {
78 #ifdef _OPENMP
79 omp_set_num_threads(1);
80 #endif
81 Cone<Integer> Candidate = rand_simplex(polytope_dim, bound);
82 nr_simplex++;
83 if (nr_simplex % 1000 == 0)
84 cout << "simplex " << nr_simplex << endl;
85 vector<vector<Integer> > supp_hyps_moved = Candidate.getSupportHyperplanes();
86 for (auto& i : supp_hyps_moved)
87 i[polytope_dim] += 1;
88 Cone<Integer> Candidate1(Type::inequalities, supp_hyps_moved, Type::grading, to_matrix(grading));
89 if (Candidate1.getNrDeg1Elements() > Candidate.getNrDeg1Elements())
90 continue; // there exists a point of height 1
91 cout << "No ht 1 jump"
92 << " #latt " << Candidate.getNrDeg1Elements() << endl;
93 // move the hyperplanes further outward
94 for (auto& i : supp_hyps_moved)
95 i[polytope_dim] += polytope_dim;
96 Cone<Integer> Candidate2(Type::inequalities, supp_hyps_moved, Type::grading, to_matrix(grading));
97 cout << "Testing " << Candidate2.getNrDeg1Elements() << " jump candidates" << endl;
98 // including the lattice points in P
99 if (exists_jump_over(Candidate, Candidate2.getDeg1Elements()))
100 continue;
101 cout << "No ht <= 1+dim jump" << endl;
102 vector<Integer> widths = lattice_widths(Candidate);
103 for (size_t i = 0; i < supp_hyps_moved.size(); ++i)
104 supp_hyps_moved[i][polytope_dim] += -polytope_dim + (widths[i]) * (polytope_dim - 2);
105 vector<vector<mpz_class> > mpz_supp_hyps;
106 convert(mpz_supp_hyps, supp_hyps_moved);
107 vector<mpz_class> mpz_grading = convertTo<vector<mpz_class> >(grading);
108 #ifdef _OPENMP
109 omp_set_num_threads(4);
110 #endif
111 Cone<mpz_class> Candidate3(Type::inequalities, mpz_supp_hyps, Type::grading, to_matrix(mpz_grading));
112 Candidate3.compute(ConeProperty::Deg1Elements, ConeProperty::Approximate);
113 vector<vector<Integer> > jumps_cand; // for conversion from mpz_class
114 convert(jumps_cand, Candidate3.getDeg1Elements());
115 cout << "Testing " << jumps_cand.size() << " jump candidates" << endl;
116 if (exists_jump_over(Candidate, jumps_cand))
117 continue;
118 cout << "Maximal simplex found" << endl;
119 cout << "Vertices" << endl;
120 Candidate.getExtremeRaysMatrix().pretty_print(cout);
121 cout << "Number of lattice points = " << Candidate.getNrDeg1Elements();
122 cout << " Multiplicity = " << Candidate.getMultiplicity() << endl;
123 } // end while
124 } // end main
125