1 /* Copyright (c) 1997-2021
2 Ewgenij Gawrilow, Michael Joswig, and the polymake team
3 Technische Universität Berlin, Germany
4 https://polymake.org
5
6 This program is free software; you can redistribute it and/or modify it
7 under the terms of the GNU General Public License as published by the
8 Free Software Foundation; either version 2, or (at your option) any
9 later version: http://www.gnu.org/licenses/gpl.txt.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15 --------------------------------------------------------------------------------
16 */
17
18 #include "polymake/client.h"
19 #include "polymake/Rational.h"
20 #include "polymake/Bitset.h"
21 #include "polymake/Set.h"
22 #include "polymake/hash_map"
23 #include "polymake/Vector.h"
24 #include "polymake/Matrix.h"
25 #include "polymake/topaz/SimplicialComplex_as_FaceMap.h"
26 #include "polymake/topaz/HomologyComplex.h"
induced_subcomplex(BigObject p_in,const Set<Int> & V_in,OptionSet options)27 #include "polymake/topaz/IntersectionForm.h"
28
29 namespace polymake { namespace topaz {
30
31 namespace {
32
33 typedef Rational coeff_type;
34 typedef Matrix<coeff_type> cup_type;
35 typedef Array<Int> reordering_type;
36
37 bool pivot(Int& k, const Int i, const cup_type& M, const reordering_type& ind)
38 {
39 const Int n = M.rows();
40 k = i;
41 while (k<n && is_zero(M(ind[k],ind[k])))
42 ++k;
43 return (k<n);
44 }
45
46 bool non_zero(Int& k, const Int i, const cup_type& M, const reordering_type& ind)
47 {
48 const Int n = M.rows();
49 k = i;
50 while (k<n && is_zero(M(ind[k],ind[i])))
51 ++k;
52 return (k<n);
53 }
54
55 void signature(cup_type& M, Int& positive, Int& negative)
56 {
57 const Int n = M.rows();
58 reordering_type ind(n,entire(sequence(0,n)));
59
60 Int i = 0, k;
61 positive = negative = 0;
62
63 while (i < n) {
64 if (pivot(k, i, M, ind)) {
65 // eliminate with pivot element on the diagonal
66 if (k!=i) std::swap(ind[i],ind[k]);
67 coeff_type p(M(ind[i],ind[i]));
68
69 for (Int j = i+1; j < n; ++j) {
70 coeff_type c(M(ind[j],ind[i])/p);
71 M[ind[j]] -= c * M[ind[i]];
72 }
73 for (Int j = i+1; j < n; ++j)
74 M(ind[i], ind[j]) = 0;
75
76 if (p>0)
77 ++positive;
78 else
79 ++negative;
80
81 ++i;
82 } else {
83 // diagonalize hyperbolic pairs
84 if (non_zero(k,i,M,ind)) { // i!=k
85 if (k!=i+1) std::swap(ind[i+1],ind[k]);
86
87 Vector<coeff_type> sum(M[ind[i]]+M[ind[i+1]]), diff(M[ind[i]]-M[ind[i+1]]);
88 M[ind[i]]=sum; M.col(ind[i])=sum; M[ind[i+1]]=diff; M.col(ind[i+1])=diff;
89 M(ind[i],ind[i])*=2; M(ind[i+1],ind[i+1])*=-2;
90 } else // zero row/column
91 ++i;
92
93 }
94 }
95 }
96
97 } // end unnamed namespace
98
99 void intersection_form(BigObject p)
100 {
101 typedef CycleGroup<Integer> cycle_type;
102 const Array<cycle_type> Cycles = p.give("CYCLES");
103
104 const Int d = Cycles.size()-1;
105 if (d%4 != 0)
106 throw std::runtime_error("intersection_form: Dimension " + std::to_string(d) + " not divisible by 4");
107
108 const cycle_type::face_list facets(Cycles[d].faces);
109 const cycle_type::coeff_matrix signs(Cycles[d].coeffs);
110
111 if (signs.rows() != 1)
112 throw std::runtime_error("intersection_form: Expected exactly one top level homology class, found " + std::to_string(signs.rows()));
113
114 auto f = entire(facets);
115 auto s = entire(signs[0]);
116 hash_map<Bitset,Integer> SignedFacets;
117 for ( ; !f.at_end() && !s.at_end(); ++f, ++s)
118 SignedFacets[Bitset(*f)]=*s;
119
120 const Array<cycle_type> CoCycles = p.give("COCYCLES");
121 const cycle_type::face_list small_faces(CoCycles[d/2].faces);
122 const cycle_type::coeff_matrix small_cocycles(CoCycles[d/2].coeffs);
123 const Int n = small_cocycles.rows();
124
125 cup_type Cup(n,n);
126
127 Int parity = 0; // until we are not convinced of the converse we assume that the intersection form is even
128
129 for (auto c1 = entire<indexed>(rows(small_cocycles)); !c1.at_end(); ++c1) {
130 for (auto c2 = entire<indexed>(rows(small_cocycles)); !c2.at_end(); ++c2) {
131 Integer cup_product(0);
132 for (auto x = entire(*c1); !x.at_end(); ++x) {
133 const Set<Int> face_x = small_faces[x.index()];
134 const Bitset bit_face_x(face_x);
135 for (auto y = entire(*c2); !y.at_end(); ++y) {
136 const Set<Int> face_y = small_faces[y.index()];
137 Bitset this_union(bit_face_x+Bitset(face_y));
138 if (face_x.back()==face_y.front() && SignedFacets.find(this_union)!=SignedFacets.end())
139 cup_product+=SignedFacets[this_union]*(*x)*(*y);
140 }
141 }
142 Cup(c1.index(), c2.index())=cup_product;
143 if (c1.index() == c2.index() && cup_product.odd())
144 parity=1;
145 }
146 }
147
148 #if POLYMAKE_DEBUG
149 for (Int i = 0; i < n; ++i)
150 for (Int j = 0; j < i; ++j)
151 if (Cup(i,j) != Cup(j,i)) {
152 std::ostringstream err;
153 err << "resulting cup product matrix not symmetric: [" << i << "," << j << "]\n";
154 wrap(err) << Cup;
155 throw std::runtime_error(err.str());
156 }
157 #endif
158
159 IntersectionForm IF;
160 IF.parity=parity;
161 signature(Cup,IF.positive,IF.negative);
162 p.take("INTERSECTION_FORM") << IF;
163 }
164
165 Function4perl(&intersection_form, "intersection_form(SimplicialComplex)");
166
167 } }
168
169 // Local Variables:
170 // mode:C++
171 // c-basic-offset:3
172 // indent-tabs-mode:nil
173 // End:
174