1 /*
2 * NormalizInterface: GAP wrapper for Normaliz
3 * Copyright (C) 2014 Sebastian Gutsche, Max Horn, Christof Söger
4 *
5 * This program is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU General Public License as
7 * published by the Free Software Foundation; either version 2 of
8 * the License, or (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public
16 * License along with this program; if not, write to the Free
17 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
18 * Boston, MA 02110-1301, USA.
19 */
20
21 /*
22 #! @Chapter Functions
23 #! @Section YOU FORGOT TO SET A SECTION
24 */
25
26 #include "src/compiled.h" /* GAP headers */
27
28 #include "libnormaliz/cone.h"
29 #include "libnormaliz/map_operations.h"
30
31 #include <vector>
32
33 #include <csignal>
34 using std::signal;
35
36 typedef void (*sighandler_t)(int);
37
38 // old versions of libnormaliz (before 2.99.1) did not include such a define
39 #if !defined(NMZ_RELEASE) || NMZ_RELEASE < 30504
40 #error Your Normaliz version is to old! Update to 3.5.4 or newer.
41 #endif
42
43 #define FUNC_BEGIN try {
44
45 #define FUNC_END \
46 } \
47 catch (std::exception & e) \
48 { \
49 ErrorQuit(e.what(), 0, 0); \
50 return Fail; \
51 } \
52 catch (const char * a) \
53 { \
54 ErrorQuit(a, 0, 0); \
55 return Fail; \
56 } \
57 catch (const std::string & a) \
58 { \
59 ErrorQuit(a.c_str(), 0, 0); \
60 return Fail; \
61 } \
62 catch (...) \
63 { \
64 ErrorQuit("PANIC: caught general exception", 0, 0); \
65 return Fail; \
66 }
67
68 #define SIGNAL_HANDLER_BEGIN \
69 sighandler_t current_interpreter_sigint_handler = \
70 signal(SIGINT, signal_handler); \
71 try {
72
73 #define SIGNAL_HANDLER_END \
74 } \
75 catch (libnormaliz::InterruptException & e) \
76 { \
77 signal(SIGINT, current_interpreter_sigint_handler); \
78 libnormaliz::nmz_interrupted = false; \
79 ErrorQuit("computation interrupted", 0, 0); \
80 return 0; \
81 } \
82 catch (...) \
83 { \
84 signal(SIGINT, current_interpreter_sigint_handler); \
85 throw; \
86 } \
87 signal(SIGINT, current_interpreter_sigint_handler);
88
89
90 // Paranoia check
91 #ifdef SYS_IS_64_BIT
92 #if GMP_LIMB_BITS != 64
93 #error GAP compiled in 64 bit mode, but GMP limbs are not 64 bit
94 #endif
95 #else
96 #if GMP_LIMB_BITS != 32
97 #error GAP compiled in 32 bit mode, but GMP limbs are not 32 bit
98 #endif
99 #endif
100
101
102 using libnormaliz::Cone;
103 // using libnormaliz::ConeProperty;
104 using libnormaliz::ConeProperties;
105 using libnormaliz::Sublattice_Representation;
106 using libnormaliz::Type::InputType;
107
108 using std::map;
109 using std::string;
110 using std::vector;
111
signal_handler(int signal)112 static void signal_handler(int signal)
113 {
114 libnormaliz::nmz_interrupted = true;
115 }
116
117 static Obj TheTypeNormalizCone;
118
119 static UInt T_NORMALIZ = 0;
120
IS_CONE(Obj o)121 static inline bool IS_CONE(Obj o)
122 {
123 return TNUM_OBJ(o) == T_NORMALIZ;
124 }
125
126 template <typename Integer>
SET_CONE(Obj o,libnormaliz::Cone<Integer> * p)127 static inline void SET_CONE(Obj o, libnormaliz::Cone<Integer> * p)
128 {
129 ADDR_OBJ(o)[0] = reinterpret_cast<Obj>(p);
130 }
131
132 template <typename Integer>
GET_CONE(Obj o)133 static inline libnormaliz::Cone<Integer> * GET_CONE(Obj o)
134 {
135 return reinterpret_cast<libnormaliz::Cone<Integer> *>(ADDR_OBJ(o)[0]);
136 }
137
NewCone(Cone<mpz_class> * C)138 static Obj NewCone(Cone<mpz_class> * C)
139 {
140 Obj o;
141 o = NewBag(T_NORMALIZ, 1 * sizeof(Obj));
142 SET_CONE<mpz_class>(o, C);
143 return o;
144 }
145
NewProxyCone(Cone<mpz_class> * C,Obj parentCone)146 static Obj NewProxyCone(Cone<mpz_class> * C, Obj parentCone)
147 {
148 Obj o;
149 o = NewBag(T_NORMALIZ, 2 * sizeof(Obj));
150 SET_CONE<mpz_class>(o, C);
151 ADDR_OBJ(o)[1] = parentCone;
152 return o;
153 }
154
155 /* Free function */
NormalizFreeFunc(Obj o)156 static void NormalizFreeFunc(Obj o)
157 {
158 if (SIZE_OBJ(o) != 2 * sizeof(Obj)) {
159 delete GET_CONE<mpz_class>(o);
160 }
161 }
162
163 /* Type object function for the object */
NormalizTypeFunc(Obj o)164 static Obj NormalizTypeFunc(Obj o)
165 {
166 return TheTypeNormalizCone;
167 }
168
NormalizCopyFunc(Obj o,Int mut)169 static Obj NormalizCopyFunc(Obj o, Int mut)
170 {
171 // Cone objects are mathematically immutable, so
172 // we don't need to do anything,
173 return o;
174 }
175
NormalizCleanFunc(Obj o)176 static void NormalizCleanFunc(Obj o)
177 {
178 }
179
NormalizIsMutableObjFuncs(Obj o)180 static Int NormalizIsMutableObjFuncs(Obj o)
181 {
182 // Cone objects are mathematically immutable.
183 return 0L;
184 }
185
186 //
187 // conversion of various C++ types to GAP objects
188 //
189
190 template <typename T, typename U>
191 static Obj NmzToGAP(const std::pair<T, U> & in);
192
193 template <typename T>
194 static Obj NmzToGAP(const vector<T> & in);
195
NmzToGAP(const mpz_t x)196 static Obj NmzToGAP(const mpz_t x)
197 {
198 #if GAP_KERNEL_MAJOR_VERSION >= 7
199 return MakeObjInt((const UInt *)x->_mp_d,
200 x->_mp_size); // added in GAP 4.11
201 #else
202 Obj res;
203 Int size = x->_mp_size;
204 int sign;
205 if (size == 0) {
206 return INTOBJ_INT(0);
207 }
208 else if (size < 0) {
209 size = -size;
210 sign = -1;
211 }
212 else {
213 sign = +1;
214 }
215 if (size == 1) {
216 res = ObjInt_UInt(x->_mp_d[0]);
217 if (sign < 0)
218 res = AInvInt(res);
219 }
220 else {
221 size = sizeof(mp_limb_t) * size;
222 if (sign > 0)
223 res = NewBag(T_INTPOS, size);
224 else
225 res = NewBag(T_INTNEG, size);
226 memcpy(ADDR_INT(res), x->_mp_d, size);
227 }
228 return res;
229 #endif
230 }
231
NmzToGAP(unsigned int x)232 static Obj NmzToGAP(unsigned int x) // key_t = unsigned int
233 {
234 return ObjInt_UInt(x);
235 }
236
NmzToGAP(unsigned long x)237 static Obj NmzToGAP(unsigned long x) // size_t = unsigned long
238 {
239 return ObjInt_UInt(x);
240 }
241
NmzToGAP(int x)242 static Obj NmzToGAP(int x)
243 {
244 return ObjInt_Int(x);
245 }
246
NmzToGAP(long x)247 static Obj NmzToGAP(long x)
248 {
249 return ObjInt_Int(x);
250 }
251
NmzToGAP(mpz_class x)252 static Obj NmzToGAP(mpz_class x)
253 {
254 return NmzToGAP(x.get_mpz_t());
255 }
256
NmzToGAP(mpq_class x)257 static Obj NmzToGAP(mpq_class x)
258 {
259 Obj num = NmzToGAP(x.get_num());
260 Obj den = NmzToGAP(x.get_den());
261 return QUO(num, den);
262 }
263
NmzToGAP(double x)264 static Obj NmzToGAP(double x)
265 {
266 return NEW_MACFLOAT(x);
267 }
268
269 template <typename T>
NmzToGAP(const libnormaliz::STANLEYDATA<T> & sd)270 static Obj NmzToGAP(const libnormaliz::STANLEYDATA<T> & sd)
271 {
272 Obj ret = NEW_PLIST(T_PLIST, 2);
273 ASS_LIST(ret, 1, NmzToGAP(sd.key));
274 ASS_LIST(ret, 2, NmzToGAP(sd.offsets));
275 return ret;
276 }
277
278 /* TODO: HSOP
279 * There are two representations for Hilbert series in Normaliz,
280 * standard and HSOP. Currently, only the standard representation is returned.
281 */
NmzToGAP(const libnormaliz::HilbertSeries & HS)282 static Obj NmzToGAP(const libnormaliz::HilbertSeries & HS)
283 {
284 Obj ret;
285 ret = NEW_PLIST(T_PLIST, 3);
286 ASS_LIST(ret, 1, NmzToGAP(HS.getNum()));
287 ASS_LIST(ret, 2, NmzToGAP(libnormaliz::to_vector(HS.getDenom())));
288 ASS_LIST(ret, 3, NmzToGAP(HS.getShift()));
289 return ret;
290 }
291
292 static Obj
NmzHilbertQuasiPolynomialToGAP(const libnormaliz::HilbertSeries & HS)293 NmzHilbertQuasiPolynomialToGAP(const libnormaliz::HilbertSeries & HS)
294 {
295 Obj ret = NmzToGAP(HS.getHilbertQuasiPolynomial());
296 AddList(ret, NmzToGAP(HS.getHilbertQuasiPolynomialDenom()));
297 return ret;
298 }
299
300 //
301 // generic recursive conversion of C++ containers to GAP objects
302 //
303
304 // convert maps to list of [key,value] pairs
305 template <typename K, typename V>
NmzToGAP(const map<K,V> & in)306 static Obj NmzToGAP(const map<K, V> & in)
307 {
308 Obj list = NEW_PLIST(T_PLIST, in.size());
309 for (const auto & it : in) {
310 AddList(list, NmzToGAP(it));
311 }
312 return list;
313 }
314
315 template <typename T>
NmzToGAP(const std::list<T> & in)316 static Obj NmzToGAP(const std::list<T> & in)
317 {
318 Obj list = NEW_PLIST(T_PLIST, in.size());
319 for (const auto & it : in) {
320 AddList(list, NmzToGAP(it));
321 }
322 return list;
323 }
324
325 template <typename T>
NmzToGAP(const vector<T> & in)326 static Obj NmzToGAP(const vector<T> & in)
327 {
328 const size_t n = in.size();
329 Obj list = NEW_PLIST(T_PLIST, n);
330 for (size_t i = 0; i < n; ++i) {
331 ASS_LIST(list, i + 1, NmzToGAP(in[i]));
332 }
333 return list;
334 }
335
336 template <>
NmzToGAP(const vector<bool> & in)337 Obj NmzToGAP(const vector<bool> & in)
338 {
339 const size_t n = in.size();
340 Obj list = NewBag(T_BLIST, SIZE_PLEN_BLIST(n));
341 SET_LEN_BLIST(list, n);
342 for (size_t i = 0; i < n; ++i) {
343 if (in[i])
344 SET_BIT_BLIST(list, i + 1);
345 }
346 return list;
347 }
348
NmzToGAP(const boost::dynamic_bitset<> & in)349 static Obj NmzToGAP(const boost::dynamic_bitset<> & in)
350 {
351 const size_t n = in.size();
352 Obj list = NewBag(T_BLIST, SIZE_PLEN_BLIST(n));
353 SET_LEN_BLIST(list, n);
354 for (size_t i = 0; i < n; ++i) {
355 if (in.test(i))
356 SET_BIT_BLIST(list, i + 1);
357 }
358 return list;
359 }
360
361 template <typename T, typename U>
NmzToGAP(const std::pair<T,U> & in)362 static Obj NmzToGAP(const std::pair<T, U> & in)
363 {
364 Obj pair = NEW_PLIST(T_PLIST, 2);
365 ASS_LIST(pair, 1, NmzToGAP(in.first));
366 ASS_LIST(pair, 2, NmzToGAP(in.second));
367 return pair;
368 }
369
370 // TODO: generate a "proper" MatrixObj, once available
371 template <typename T>
NmzToGAP(const libnormaliz::Matrix<T> & in)372 static Obj NmzToGAP(const libnormaliz::Matrix<T> & in)
373 {
374 return NmzToGAP(in.get_elements());
375 }
376
GAPNumberToNmz(long & out,Obj x)377 bool GAPNumberToNmz(long & out, Obj x)
378 {
379 if (IS_INTOBJ(x)) {
380 out = INT_INTOBJ(x);
381 return true;
382 }
383 else if (TNUM_OBJ(x) == T_INTPOS || TNUM_OBJ(x) == T_INTNEG) {
384 UInt size = SIZE_INT(x);
385 if (size == 1) {
386 out = *ADDR_INT(x);
387 if (out < 0)
388 return false; // overflow
389 if (TNUM_OBJ(x) == T_INTNEG)
390 out = -out;
391 return true;
392 }
393 }
394 return false;
395 }
396
GAPNumberToNmz(mpz_class & out,Obj x)397 bool GAPNumberToNmz(mpz_class & out, Obj x)
398 {
399 if (IS_INTOBJ(x)) {
400 out = INT_INTOBJ(x);
401 return true;
402 }
403 else if (TNUM_OBJ(x) == T_INTPOS || TNUM_OBJ(x) == T_INTNEG) {
404 UInt size = SIZE_INT(x);
405 mpz_ptr m = out.get_mpz_t();
406 mpz_realloc2(m, size * GMP_NUMB_BITS);
407 memcpy(m->_mp_d, ADDR_INT(x), sizeof(mp_limb_t) * size);
408 m->_mp_size = (TNUM_OBJ(x) == T_INTPOS) ? (Int)size : -(Int)size;
409 return true;
410 }
411 return false;
412 }
413
GAPNumberToNmz(mpq_class & out,Obj x)414 bool GAPNumberToNmz(mpq_class & out, Obj x)
415 {
416 if (IS_INTOBJ(x)) {
417 out = INT_INTOBJ(x);
418 return true;
419 }
420 else if (TNUM_OBJ(x) == T_INTPOS || TNUM_OBJ(x) == T_INTNEG) {
421 out.get_den() = 1;
422 return GAPNumberToNmz(out.get_num(), x);
423 }
424 else if (TNUM_OBJ(x) == T_RAT) {
425 return GAPNumberToNmz(out.get_num(), NUM_RAT(x)) &&
426 GAPNumberToNmz(out.get_den(), DEN_RAT(x));
427 }
428 return false;
429 }
430
431 template <typename Number>
GAPVectorToNmz(vector<Number> & out,Obj V)432 static bool GAPVectorToNmz(vector<Number> & out, Obj V)
433 {
434 if (!IS_PLIST(V) || !IS_DENSE_LIST(V))
435 return false;
436 const int n = LEN_PLIST(V);
437 out.resize(n);
438 for (int i = 0; i < n; ++i) {
439 Obj tmp = ELM_PLIST(V, i + 1);
440 if (!GAPNumberToNmz(out[i], tmp))
441 return false;
442 }
443 return true;
444 }
445
446 template <typename Number>
GAPMatrixToNmz(vector<vector<Number>> & out,Obj M)447 static bool GAPMatrixToNmz(vector<vector<Number>> & out, Obj M)
448 {
449 if (!IS_PLIST(M) || !IS_DENSE_LIST(M))
450 return false;
451 const int nr = LEN_PLIST(M);
452 out.resize(nr);
453 for (int i = 0; i < nr; ++i) {
454 bool okay = GAPVectorToNmz(out[i], ELM_PLIST(M, i + 1));
455 if (!okay)
456 return false;
457 }
458 return true;
459 }
460
461 template <typename Integer>
_NmzConeIntern(Obj input_list)462 static Obj _NmzConeIntern(Obj input_list)
463 {
464 bool has_polynomial_input = false;
465 string polynomial;
466 map<InputType, vector<vector<mpq_class>>> input;
467 const int n = LEN_PLIST(input_list);
468 if (n & 1) {
469 throw std::runtime_error(
470 "Input list must have even number of elements");
471 }
472 for (int i = 0; i < n; i += 2) {
473 Obj type = ELM_PLIST(input_list, i + 1);
474 if (!IS_STRING_REP(type)) {
475 throw std::runtime_error(
476 "Element " + std::to_string(i + 1) +
477 " of the input list must be a type string");
478 }
479 string type_str(CSTR_STRING(type));
480 Obj M = ELM_PLIST(input_list, i + 2);
481 if (type_str.compare("polynomial") == 0) {
482 if (!IS_STRING_REP(M)) {
483 throw std::runtime_error(
484 "Element " + std::to_string(i + 2) +
485 " of the input list must be a string");
486 }
487 polynomial = string(CSTR_STRING(M));
488 has_polynomial_input = true;
489 continue;
490 }
491 vector<vector<mpq_class>> Mat;
492 bool okay = GAPMatrixToNmz(Mat, M);
493 if (!okay) {
494 throw std::runtime_error(
495 "Element " + std::to_string(i + 2) +
496 " of the input list must be an integer matrix");
497 }
498
499 input[libnormaliz::to_type(type_str)] = Mat;
500 }
501
502 Cone<Integer> * C = new Cone<Integer>(input);
503 if (has_polynomial_input) {
504 C->setPolynomial(polynomial);
505 }
506 Obj Cone = NewCone(C);
507 return Cone;
508 }
509
Func_NmzCone(Obj self,Obj input_list)510 static Obj Func_NmzCone(Obj self, Obj input_list)
511 {
512 if (!IS_PLIST(input_list) || !IS_DENSE_LIST(input_list))
513 ErrorQuit("Input argument must be a list", 0, 0);
514
515 FUNC_BEGIN
516 return _NmzConeIntern<mpz_class>(input_list);
517 FUNC_END
518 }
519
Func_NmzCompute(Obj self,Obj cone,Obj to_compute)520 static Obj Func_NmzCompute(Obj self, Obj cone, Obj to_compute)
521 {
522 if (!IS_CONE(cone))
523 ErrorQuit("<cone> must be a Normaliz cone", 0, 0);
524 if (!IS_PLIST(to_compute) || !IS_DENSE_LIST(to_compute))
525 ErrorQuit("<props> must be a list of strings", 0, 0);
526
527 FUNC_BEGIN
528 ConeProperties propsToCompute;
529 // we have a list
530 const int n = LEN_PLIST(to_compute);
531
532 for (int i = 0; i < n; ++i) {
533 Obj prop = ELM_PLIST(to_compute, i + 1);
534 if (!IS_STRING_REP(prop)) {
535 throw std::runtime_error(
536 "Element " + std::to_string(i + 1) +
537 " of the input list must be a type string");
538 }
539 string prop_str(CSTR_STRING(prop));
540 propsToCompute.set(libnormaliz::toConeProperty(prop_str));
541 }
542
543 Cone<mpz_class> * C = GET_CONE<mpz_class>(cone);
544
545 ConeProperties notComputed;
546 SIGNAL_HANDLER_BEGIN
547 notComputed = C->compute(propsToCompute);
548 SIGNAL_HANDLER_END
549
550 // Cone.compute returns the not computed properties
551 // we return a bool, true when everything requested was computed
552 return notComputed.none() ? True : False;
553 FUNC_END
554 }
555
556
557 /*
558 #! @Section Use a NmzCone
559 #! @Arguments cone property
560 #! @Returns whether the cone has already computed the given property
561 #! @Description
562 #! See <Ref Func="NmzConeProperty"/> for a list of recognized properties.
563 #!
564 #! @InsertChunk NmzHasConeProperty_example
565 DeclareGlobalFunction("NmzHasConeProperty");
566 */
FuncNmzHasConeProperty(Obj self,Obj cone,Obj prop)567 static Obj FuncNmzHasConeProperty(Obj self, Obj cone, Obj prop)
568 {
569 if (!IS_CONE(cone))
570 ErrorQuit("<cone> must be a Normaliz cone", 0, 0);
571 if (!IS_STRING_REP(prop))
572 ErrorQuit("<prop> must be a string", 0, 0);
573
574 FUNC_BEGIN
575
576 libnormaliz::ConeProperty::Enum p =
577 libnormaliz::toConeProperty(CSTR_STRING(prop));
578
579 Cone<mpz_class> * C = GET_CONE<mpz_class>(cone);
580 return C->isComputed(p) ? True : False;
581
582 FUNC_END
583 }
584
585 /*
586 #! @Section Use a NmzCone
587 #! @Arguments cone
588 #! @Returns a list of strings representing the known (computed) cone
589 #! properties
590 #! @Description
591 #! Given a Normaliz cone object, return a list of all properties already
592 #! computed for the cone.
593 #!
594 #! @InsertChunk NmzKnownConeProperties_example
595 DeclareGlobalFunction("NmzKnownConeProperties");
596 */
FuncNmzKnownConeProperties(Obj self,Obj cone)597 static Obj FuncNmzKnownConeProperties(Obj self, Obj cone)
598 {
599 if (!IS_CONE(cone))
600 ErrorQuit("<cone> must be a Normaliz cone", 0, 0);
601
602 FUNC_BEGIN
603
604 size_t n = 0;
605 Obj M = NEW_PLIST(T_PLIST, libnormaliz::ConeProperty::EnumSize);
606 Cone<mpz_class> * C = GET_CONE<mpz_class>(cone);
607
608 // FIXME: This code could be a lot simpler if there was
609 // a Cone method for reading the value of is_Computed.
610 for (int i = 0; i < libnormaliz::ConeProperty::EnumSize; ++i) {
611 libnormaliz::ConeProperty::Enum p =
612 (libnormaliz::ConeProperty::Enum)i;
613
614 // skip internal control properties
615 if (p == libnormaliz::ConeProperty::ExplicitHilbertSeries ||
616 p == libnormaliz::ConeProperty::NakedDual)
617 continue;
618
619 if (C->isComputed(p)) {
620 string prop_name(libnormaliz::toString(p));
621
622 ASS_LIST(M, ++n, MakeImmString(prop_name.c_str()));
623 if (p == libnormaliz::ConeProperty::HilbertSeries) {
624 const libnormaliz::HilbertSeries & HS = C->getHilbertSeries();
625 HS.computeHilbertQuasiPolynomial();
626 if (HS.isHilbertQuasiPolynomialComputed()) {
627 ASS_LIST(M, ++n, MakeImmString("HilbertQuasiPolynomial"));
628 }
629 }
630 }
631 }
632 return M;
633
634 FUNC_END
635 }
636
637 template <typename Integer>
_NmzConePropertyImpl(Obj cone,Obj prop)638 static Obj _NmzConePropertyImpl(Obj cone, Obj prop)
639 {
640 Cone<Integer> * C = GET_CONE<Integer>(cone);
641
642 libnormaliz::ConeProperty::Enum p =
643 libnormaliz::toConeProperty(CSTR_STRING(prop));
644 ConeProperties notComputed;
645 SIGNAL_HANDLER_BEGIN
646 notComputed = C->compute(ConeProperties(p));
647 SIGNAL_HANDLER_END
648 #if NMZ_RELEASE >= 30700
649 // workaround a bug where computing HilbertQuasiPolynomial after
650 // HilbertSeries was already computed returned notComputed equal to
651 // NoGradingDenom, even though of course the quasi polynomial was
652 // available.
653 notComputed.reset(libnormaliz::ConeProperty::NoGradingDenom);
654 #endif
655 if (notComputed.any()) {
656 ErrorQuit("Failed to compute %s, missing properties: %s",
657 (Int)(libnormaliz::toString(p).c_str()),
658 (Int)(libnormaliz::toString(notComputed).c_str()));
659 return Fail;
660 }
661
662 // workaround bug in certain Normaliz versions, where the output type for
663 // ClassGroups is reported as libnormaliz::OutputType::Vector, but calling
664 // getVectorConeProperty on it produces an error
665 if (p == libnormaliz::ConeProperty::ClassGroup)
666 return NmzToGAP(C->getClassGroup());
667
668 #if NMZ_RELEASE >= 30700 && NMZ_RELEASE < 30703
669 // workaround bug where getMachineIntegerConeProperty does not support
670 // NumberLatticePoints
671 if (p == libnormaliz::ConeProperty::NumberLatticePoints)
672 return NmzToGAP(C->getNumberLatticePoints());
673 #endif
674
675 // workaround: these two properties are marked as having output type
676 // "void"
677 if (p == libnormaliz::ConeProperty::IsTriangulationNested)
678 return C->isTriangulationNested() ? True : False;
679 if (p == libnormaliz::ConeProperty::IsTriangulationPartial)
680 return C->isTriangulationPartial() ? True : False;
681
682 switch (libnormaliz::output_type(p)) {
683 case libnormaliz::OutputType::Matrix:
684 // TODO: switch to getMatrixConePropertyMatrix ?
685 return NmzToGAP(C->getMatrixConeProperty(p));
686
687 case libnormaliz::OutputType::MatrixFloat:
688 return NmzToGAP(C->getFloatMatrixConeProperty(p));
689
690 case libnormaliz::OutputType::Vector:
691 return NmzToGAP(C->getVectorConeProperty(p));
692
693 case libnormaliz::OutputType::Integer:
694 return NmzToGAP(C->getIntegerConeProperty(p));
695
696 case libnormaliz::OutputType::GMPInteger:
697 return NmzToGAP(C->getGMPIntegerConeProperty(p));
698
699 case libnormaliz::OutputType::Rational:
700 return NmzToGAP(C->getRationalConeProperty(p));
701
702 #if NMZ_RELEASE >= 30700
703 case libnormaliz::OutputType::FieldElem:
704 throw "OutputType::FieldElem not yet supported";
705 #endif
706
707 case libnormaliz::OutputType::Float:
708 return NmzToGAP(C->getFloatConeProperty(p));
709
710 case libnormaliz::OutputType::MachineInteger:
711 return NmzToGAP(C->getMachineIntegerConeProperty(p));
712
713 case libnormaliz::OutputType::Bool:
714 return C->getBooleanConeProperty(p) ? True : False;
715
716 case libnormaliz::OutputType::Complex:
717 // more complex data structures are handled below
718 break;
719
720 case libnormaliz::OutputType::Void:
721 // throw "cone property is input-only";
722 return Fail;
723
724 default:
725 throw "unsupported output_type";
726 }
727
728 switch (p) {
729 #if NMZ_RELEASE >= 30703
730 case libnormaliz::ConeProperty::AmbientAutomorphisms:
731 case libnormaliz::ConeProperty::Automorphisms:
732 case libnormaliz::ConeProperty::CombinatorialAutomorphisms:
733 case libnormaliz::ConeProperty::EuclideanAutomorphisms:
734 case libnormaliz::ConeProperty::RationalAutomorphisms:
735 throw "querying automorphisms not yet supported";
736 #endif
737
738 case libnormaliz::ConeProperty::ConeDecomposition:
739 return NmzToGAP(C->getOpenFacets());
740
741 #if NMZ_RELEASE >= 30600
742 case libnormaliz::ConeProperty::EhrhartQuasiPolynomial:
743 #if NMZ_RELEASE >= 30700
744 return NmzHilbertQuasiPolynomialToGAP(C->getEhrhartSeries());
745 #else
746 throw "Extracting EhrhartQuasiPolynomial requires Normaliz >= 3.7.0";
747 #endif
748 #endif
749
750 case libnormaliz::ConeProperty::EhrhartSeries:
751 #if NMZ_RELEASE >= 30700
752 return NmzToGAP(C->getEhrhartSeries());
753 #else
754 throw "Extracting EhrhartSeries requires Normaliz >= 3.7.0";
755 #endif
756
757 #if NMZ_RELEASE >= 30700
758 case libnormaliz::ConeProperty::FaceLattice:
759 return NmzToGAP(C->getFaceLattice());
760
761 case libnormaliz::ConeProperty::FVector:
762 return NmzToGAP(C->getFVector());
763 #endif
764
765 case libnormaliz::ConeProperty::HilbertQuasiPolynomial:
766 return NmzHilbertQuasiPolynomialToGAP(C->getHilbertSeries());
767
768 case libnormaliz::ConeProperty::HilbertSeries:
769 return NmzToGAP(C->getHilbertSeries());
770
771 case libnormaliz::ConeProperty::InclusionExclusionData:
772 return NmzToGAP(C->getInclusionExclusionData());
773
774 case libnormaliz::ConeProperty::IntegerHull:
775 return NewProxyCone(&(C->getIntegerHullCone()), cone);
776
777 case libnormaliz::ConeProperty::ProjectCone:
778 return NewProxyCone(&(C->getProjectCone()), cone);
779
780 case libnormaliz::ConeProperty::Sublattice:
781 return _NmzBasisChangeIntern(C);
782
783 case libnormaliz::ConeProperty::StanleyDec:
784 return NmzToGAP(C->getStanleyDec());
785
786 case libnormaliz::ConeProperty::Triangulation:
787 return NmzToGAP(C->getTriangulation());
788
789 case libnormaliz::ConeProperty::WeightedEhrhartQuasiPolynomial:
790 return NmzHilbertQuasiPolynomialToGAP(
791 C->getWeightedEhrhartSeries().first);
792
793 case libnormaliz::ConeProperty::WeightedEhrhartSeries:
794 return NmzToGAP(C->getWeightedEhrhartSeries());
795
796 default:
797 throw "unsupported cone property " + libnormaliz::toString(p);
798 }
799
800 return Fail;
801 }
802
Func_NmzConeProperty(Obj self,Obj cone,Obj prop)803 static Obj Func_NmzConeProperty(Obj self, Obj cone, Obj prop)
804 {
805 if (!IS_CONE(cone))
806 ErrorQuit("<cone> must be a Normaliz cone", 0, 0);
807 if (!IS_STRING_REP(prop))
808 ErrorQuit("<prop> must be a string", 0, 0);
809
810 FUNC_BEGIN
811 return _NmzConePropertyImpl<mpz_class>(cone, prop);
812 FUNC_END
813 }
814
815
816 /*
817 #! @Section Use a NmzCone
818 #! @Arguments verboseFlag
819 #! @Returns the previous verbosity
820 #! @Description
821 #! Set the global default verbosity state in libnormaliz.
822 #! This will influence all NmzCone created afterwards, but not any existing
823 #! ones.
824 #!
825 #! See also <Ref Func="NmzSetVerbose"/>
826 DeclareGlobalFunction("NmzSetVerboseDefault");
827 */
FuncNmzSetVerboseDefault(Obj self,Obj value)828 static Obj FuncNmzSetVerboseDefault(Obj self, Obj value)
829 {
830 if (value != True && value != False)
831 ErrorQuit("<value> must be a boolean value", 0, 0);
832 FUNC_BEGIN
833 return libnormaliz::setVerboseDefault(value == True) ? True : False;
834 FUNC_END
835 }
836
837
838 /*
839 #! @Arguments cone verboseFlag
840 #! @Returns the previous verbosity
841 #! @Description
842 #! Set the verbosity state for a cone.
843 #!
844 #! See also <Ref Func="NmzSetVerboseDefault"/>
845 DeclareGlobalFunction("NmzSetVerbose");
846 */
FuncNmzSetVerbose(Obj self,Obj cone,Obj value)847 static Obj FuncNmzSetVerbose(Obj self, Obj cone, Obj value)
848 {
849 if (!IS_CONE(cone))
850 ErrorQuit("<cone> must be a Normaliz cone", 0, 0);
851 if (value != True && value != False)
852 ErrorQuit("<value> must be a boolean value", 0, 0);
853 bool old_value;
854
855 FUNC_BEGIN
856 Cone<mpz_class> * C = GET_CONE<mpz_class>(cone);
857 old_value = C->setVerbose(value == True);
858 return old_value ? True : False;
859 FUNC_END
860 }
861
862 /*
863 #! @Section Cone properties
864 */
865
866
867 template <typename Integer>
_NmzBasisChangeIntern(Cone<Integer> * C)868 static Obj _NmzBasisChangeIntern(Cone<Integer> * C)
869 {
870 Sublattice_Representation<Integer> bc;
871 SIGNAL_HANDLER_BEGIN
872 bc = C->getSublattice();
873 SIGNAL_HANDLER_END
874
875 Obj res = NEW_PLIST(T_PLIST, 3);
876 ASS_LIST(res, 1, NmzToGAP(bc.getEmbedding()));
877 ASS_LIST(res, 2, NmzToGAP(bc.getProjection()));
878 ASS_LIST(res, 3, NmzToGAP(bc.getAnnihilator()));
879 // Dim, Rank, Equations and Congruences are already covered by special
880 // functions The index is not always computed and not so relevant
881 return res;
882 }
883
884
Func_NmzVersion(Obj self)885 static Obj Func_NmzVersion(Obj self)
886 {
887 Obj res = NEW_PLIST(T_PLIST, 3);
888 ASS_LIST(res, 1, INTOBJ_INT(NMZ_VERSION_MAJOR));
889 ASS_LIST(res, 2, INTOBJ_INT(NMZ_VERSION_MINOR));
890 ASS_LIST(res, 3, INTOBJ_INT(NMZ_VERSION_PATCH));
891 return res;
892 }
893
894 // GVAR_FUNC in GAP 4.9 and 4.10 triggers warnings when used in C++ code; so
895 // we force it to the definition from GAP 4.11
896 #undef GVAR_FUNC
897 #define GVAR_FUNC(name, nargs, args) \
898 { \
899 #name, nargs, args, (ObjFunc)Func##name, __FILE__ ":" #name \
900 }
901
902 // Table of functions to export
903 static StructGVarFunc GVarFuncs[] = {
904 GVAR_FUNC(_NmzCone, 1, "list"),
905
906 GVAR_FUNC(_NmzCompute, 2, "cone, props"),
907 GVAR_FUNC(NmzSetVerboseDefault, 1, "value"),
908 GVAR_FUNC(NmzSetVerbose, 2, "cone, value"),
909
910 GVAR_FUNC(NmzHasConeProperty, 2, "cone, prop"),
911 GVAR_FUNC(_NmzConeProperty, 2, "cone, prop"),
912 GVAR_FUNC(NmzKnownConeProperties, 1, "cone"),
913
914 GVAR_FUNC(_NmzVersion, 0, ""),
915
916 { 0 } /* Finish with an empty entry */
917
918 };
919
920 // initialise kernel data structures
InitKernel(StructInitInfo * module)921 static Int InitKernel(StructInitInfo * module)
922 {
923 /* init filters and functions */
924 InitHdlrFuncsFromTable(GVarFuncs);
925
926 InitCopyGVar("TheTypeNormalizCone", &TheTypeNormalizCone);
927
928 T_NORMALIZ = RegisterPackageTNUM("NormalizCone", NormalizTypeFunc);
929
930 InitMarkFuncBags(T_NORMALIZ, &MarkAllSubBags);
931 InitFreeFuncBag(T_NORMALIZ, &NormalizFreeFunc);
932
933 CopyObjFuncs[T_NORMALIZ] = &NormalizCopyFunc;
934 CleanObjFuncs[T_NORMALIZ] = &NormalizCleanFunc;
935 IsMutableObjFuncs[T_NORMALIZ] = &NormalizIsMutableObjFuncs;
936
937 /* return success */
938 return 0;
939 }
940
941 // initialise library data structures
InitLibrary(StructInitInfo * module)942 static Int InitLibrary(StructInitInfo * module)
943 {
944 /* init filters and functions */
945 InitGVarFuncsFromTable(GVarFuncs);
946
947 /* return success */
948 return 0;
949 }
950
951 // table of init functions
952 static StructInitInfo module = {
953 #ifdef NORMALIZSTATIC
954 /* type = */ MODULE_STATIC,
955 #else
956 /* type = */ MODULE_DYNAMIC,
957 #endif
958 /* name = */ "Normaliz",
959 /* revision_c = */ 0,
960 /* revision_h = */ 0,
961 /* version = */ 0,
962 /* crc = */ 0,
963 /* initKernel = */ InitKernel,
964 /* initLibrary = */ InitLibrary,
965 /* checkInit = */ 0,
966 /* preSave = */ 0,
967 /* postSave = */ 0,
968 /* postRestore = */ 0
969 };
970
971 #ifndef NORMALIZSTATIC
Init__Dynamic(void)972 extern "C" StructInitInfo * Init__Dynamic(void)
973 {
974 return &module;
975 }
976 #endif
977
Init__normaliz(void)978 extern "C" StructInitInfo * Init__normaliz(void)
979 {
980 return &module;
981 }
982