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