The Original CHomP Software
chains.h
Go to the documentation of this file.
1
3
23
24// Copyright (C) 1997-2020 by Pawel Pilarczyk.
25//
26// This file is part of my research software package. This is free software:
27// you can redistribute it and/or modify it under the terms of the GNU
28// General Public License as published by the Free Software Foundation,
29// either version 3 of the License, or (at your option) any later version.
30//
31// This software is distributed in the hope that it will be useful,
32// but WITHOUT ANY WARRANTY; without even the implied warranty of
33// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
34// GNU General Public License for more details.
35//
36// You should have received a copy of the GNU General Public License
37// along with this software; see the file "license.txt". If not,
38// please, see <https://www.gnu.org/licenses/>.
39
40// Started in 1997. Last revision: July 28, 2023.
41
42
43#ifndef _CHOMP_HOMOLOGY_CHAINS_H_
44#define _CHOMP_HOMOLOGY_CHAINS_H_
45
46#include "chomp/system/config.h"
47#include "chomp/system/textfile.h"
48#include "chomp/struct/autoarray.h"
49
50#include <iostream>
51#include <fstream>
52#include <cstdlib>
53#include <iomanip>
54#include <algorithm>
55
56namespace chomp {
57namespace homology {
58
59
60// templates defined within this file (in this order):
61
62// a chain, that is, a linear combination of generators
63template <class euclidom>
64class chain;
65// a set of matrices
66template <class euclidom>
67class matrices;
68// a sparse matrix whose entries are the given coefficients
69template <class euclidom>
70class mmatrix;
71// a chain complex over a given Euclidean domain
72template <class euclidom>
73class chaincomplex;
74// a chain homomorphism between two chain complexes
75template <class euclidom>
76class chainmap;
77
78// templates of functions which display the homology to the output streams
79template <class euclidom>
80outputstream &show_homology (outputstream &out, const chain<euclidom> &c);
81template <class euclidom>
82std::ostream &show_homology (std::ostream &out, const chain<euclidom> &c);
83
84
85// --------------------------------------------------
86// --------------------- chains ---------------------
87// --------------------------------------------------
88
92template <class euclidom>
93class chain
94{
95public:
97 chain ();
98
100 chain (const chain<euclidom> &c);
101
103 chain<euclidom> &operator = (const chain<euclidom> &c);
104
106 ~chain ();
107
110 int_t size () const;
111
113 bool empty () const;
114
118 euclidom getcoefficient (int_t n = -1) const;
119
122 int_t findnumber (int_t n) const;
123
126 euclidom coef (int_t i) const;
127
129 int_t num (int_t i) const;
130
133 bool contains_non_invertible () const;
134
141 int_t findbest (chain<euclidom> *table = NULL) const;
142
144 chain<euclidom> &add (int_t n, euclidom e);
145
147 chain<euclidom> &remove (int_t n);
148
153 chain<euclidom> &add (const chain<euclidom> &other,
154 euclidom e, int_t number = -1,
155 chain<euclidom> *table = NULL);
156
161 chain<euclidom> &swap (chain<euclidom> &other,
162 int_t number = -1, int_t othernumber = -1,
163 chain<euclidom> *table = NULL);
164
166 chain<euclidom> &take (chain<euclidom> &c);
167
170 chain<euclidom> &multiply (euclidom e, int_t number = -1);
171
174 outputstream &show (outputstream &out,
175 const char *label = NULL) const;
176
179 std::ostream &show (std::ostream &out, const char *label = NULL) const;
180
181private:
183 int_t len;
184
186 int_t *_n;
187
189 euclidom *_e;
190
193 chain<euclidom> &insertpair (int_t i, int_t n, euclidom e);
194
197 chain<euclidom> &removepair (int_t i);
198
200 chain<euclidom> &swapnumbers (int_t number1, int_t number2);
201
202}; /* class chain */
203
204// --------------------------------------------------
205
206template <class euclidom>
207inline chain<euclidom>::chain (): len(0)
208{
209 return;
210} /* chain<euclidom>::chain */
211
212template <class euclidom>
213inline chain<euclidom>::chain (const chain<euclidom> &c)
214{
215 // copy the length of the chain
216 len = c. len;
217
218 // allocate new tables and copy the data
219 if (len)
220 {
221 _n = new int_t [len];
222 _e = new euclidom [len];
223 if (!_n || !_e)
224 throw "Not enough memory to create a chain copy.";
225 }
226 for (int_t i = 0; i < len; ++ i)
227 {
228 _n [i] = c. _n [i];
229 _e [i] = c. _e [i];
230 }
231 return;
232} /* chain<euclidom>::chain */
233
234template <class euclidom>
235inline chain<euclidom> &chain<euclidom>::operator =
236 (const chain<euclidom> &c)
237{
238 // protect against self-assignment
239 if (&c == this)
240 return *this;
241
242 // make sure the lengths of the chains are the same
243 if (len != c.len)
244 {
245 if (len)
246 {
247 delete [] _n;
248 delete [] _e;
249 }
250 len = c.len;
251 if (len)
252 {
253 _n = new int_t [len];
254 _e = new euclidom [len];
255 if (!_n || !_e)
256 throw "Not enough memory to create a chain copy =.";
257 }
258 }
259
260 // copy the data
261 for (int_t i = 0; i < len; ++ i)
262 {
263 _n [i] = c. _n [i];
264 _e [i] = c. _e [i];
265 }
266 return *this;
267} /* chain<euclidom>::operator = */
268
269template <class euclidom>
270inline chain<euclidom>::~chain ()
271{
272 if (len)
273 {
274 delete [] _n;
275 delete [] _e;
276 }
277 return;
278} /* chain<euclidom>::~chain */
279
280template <class euclidom>
281inline int_t chain<euclidom>::size () const
282{
283 return len;
284} /* chain<euclidom>::size */
285
286template <class euclidom>
287inline bool chain<euclidom>::empty () const
288{
289 return !len;
290} /* chain<euclidom>::empty */
291
292template <class euclidom>
293/*inline*/ euclidom chain<euclidom>::getcoefficient (int_t n) const
294{
295 if (n < 0)
296 {
297 if (len > 0)
298 return _e [0];
299 else
300 {
301 euclidom zero;
302 zero = 0;
303 return zero;
304 }
305 }
306
307 int_t i = 0;
308 while ((i < len) && (_n [i] < n))
309 ++ i;
310 if ((i < len) && (_n [i] == n))
311 return _e [i];
312 euclidom zero;
313 zero = 0;
314 return zero;
315} /* chain<euclidom>::getcoefficient */
316
317template <class euclidom>
318inline int_t chain<euclidom>::findnumber (int_t n) const
319{
320 for (int_t i = 0; i < len; ++ i)
321 {
322 if (_n [i] == n)
323 return i;
324 else if (_n [i] > n)
325 return -1;
326 }
327 return -1;
328} /* chain<euclidom>::findnumber */
329
330template <class euclidom>
331inline euclidom chain<euclidom>::coef (int_t i) const
332{
333 if (i >= len)
334 throw "Wrong coefficient requested from a chain.";
335 return _e [i];
336} /* chain<euclidom>::coef */
337
338template <class euclidom>
339inline int_t chain<euclidom>::num (int_t i) const
340{
341 if (i >= len)
342 throw "Wrong number requested from a chain.";
343 return _n [i];
344} /* chain<euclidom>::num */
345
346template <class euclidom>
347inline bool chain<euclidom>::contains_non_invertible () const
348{
349 for (int_t i = 0; i < len; ++ i)
350 {
351 if (_e [i]. delta () > 1)
352 return true;
353 }
354 return false;
355} /* chain<euclidom>::contains_non_invertible */
356
357template <class euclidom>
358inline int_t chain<euclidom>::findbest (chain<euclidom> *table) const
359{
360 // if the chain is so short that the answer is obvious, return it
361 if (!len)
362 return -1;
363 if (len == 1)
364 return 0;
365
366 // find the number which has the smallest delta function value
367 int_t best_delta = _e [0]. delta ();
368 int_t best_i = 0;
369
370 // go through the whole table to find the first element
371 // with the lowest delta possible (1 is the smallest value)
372 for (int_t i = 1; (i < len) && (best_delta > 1); ++ i)
373 {
374 // compute the value of the function delta
375 int_t this_delta = _e [i]. delta ();
376
377 // if this delta is better, remember it
378 if (this_delta < best_delta)
379 {
380 best_delta = this_delta;
381 best_i = i;
382 }
383 }
384
385 // if no further analysis is required, return the result just now
386 if (!table)
387 return best_i;
388
389 // analyze which element has the shortest corresponding chain
390 // and update "best_i"
391 int_t best_length = table [_n [best_i]]. size ();
392 for (int_t i = best_i + 1; i < len; ++ i)
393 {
394 if (_e [i]. delta () != best_delta)
395 continue;
396 int_t this_length = table [_n [i]]. size ();
397 if (best_length > this_length)
398 {
399 best_length = this_length;
400 best_i = i;
401 }
402 }
403
404 return best_i;
405} /* chain<euclidom>::findbest */
406
407// --------------------------------------------------
408
409template <class euclidom>
410inline chain<euclidom> &chain<euclidom>::insertpair
411 (int_t i, int_t n, euclidom e)
412{
413 // remember the addresses of the old tables
414 int_t *old_n = _n;
415 euclidom *old_e = _e;
416
417 // allocate new tables
418 _n = new int_t [len + 1];
419 _e = new euclidom [len + 1];
420 if (!_n || !_e)
421 throw "Cannot add an element to a chain.";
422
423 // increase the length
424 ++ len;
425
426 // copy the old data and insert the new pair
427 for (int_t j = 0; j < i; ++ j)
428 {
429 _n [j] = old_n [j];
430 _e [j] = old_e [j];
431 }
432 _n [i] = n;
433 _e [i] = e;
434 for (int_t j = i + 1; j < len; ++ j)
435 {
436 _n [j] = old_n [j - 1];
437 _e [j] = old_e [j - 1];
438 }
439
440 // release the previous tables if they were allocated
441 if (len > 1)
442 {
443 delete [] old_n;
444 delete [] old_e;
445 }
446
447 return *this;
448} /* chain<euclidom>::insertpair */
449
450template <class euclidom>
451inline chain<euclidom> &chain<euclidom>::removepair (int_t i)
452{
453 // make sure the chain is not empty
454 if (!len)
455 throw "Trying to remove an element from an empty chain.";
456
457 // remember the addresses of the old tables
458 int_t *old_n = _n;
459 euclidom *old_e = _e;
460
461 // allocate new tables if necessary
462 if (len > 1)
463 {
464 _n = new int_t [len - 1];
465 _e = new euclidom [len - 1];
466 if (!_n || !_e)
467 throw "Cannot add an element to a chain.";
468 }
469
470 // decrease the length
471 -- len;
472
473 // copy the data form the previous tables
474 for (int_t j = 0; j < i; ++ j)
475 {
476 _n [j] = old_n [j];
477 _e [j] = old_e [j];
478 }
479 for (int_t j = i; j < len; ++ j)
480 {
481 _n [j] = old_n [j + 1];
482 _e [j] = old_e [j + 1];
483 }
484
485 // release the previous tables
486 delete [] old_n;
487 delete [] old_e;
488
489 return *this;
490} /* chain<euclidom>::removepair */
491
492template <class euclidom>
493inline chain<euclidom> &chain<euclidom>::swapnumbers (int_t number1,
494 int_t number2)
495{
496 // if the numbers are the same, do nothing
497 if (number1 == number2)
498 return *this;
499
500 // force the first number be less than the second number
501 if (number1 > number2)
502 std::swap (number1, number2);
503
504 // find both numbers or the positions they should be at
505 int_t i1 = 0;
506 while ((i1 < len) && (_n [i1] < number1))
507 ++ i1;
508 int_t i2 = i1;
509 while ((i2 < len) && (_n [i2] < number2))
510 ++ i2;
511
512 // if the first number was found...
513 if ((i1 < len) && (_n [i1] == number1))
514 {
515 // if both numbers were found, exchange their coefficients
516 if ((i2 < len) && (_n [i2] == number2))
517 swapelements (_e [i1], _e [i2]);
518 // if only the first was found, move it to the new position
519 else
520 {
521 euclidom temp = _e [i1];
522 for (int_t i = i1 + 1; i < i2; ++ i)
523 {
524 _n [i - 1] = _n [i];
525 _e [i - 1] = _e [i];
526 }
527 _n [i2 - 1] = number2;
528 _e [i2 - 1] = temp;
529 }
530 }
531
532 // otherwise if the second number only was found, move it to its pos.
533 else if ((i2 < len) && (_n [i2] == number2))
534 {
535 euclidom temp = _e [i2];
536 for (int_t i = i2; i > i1; -- i)
537 {
538 _n [i] = _n [i - 1];
539 _e [i] = _e [i - 1];
540 }
541 _n [i1] = number1;
542 _e [i1] = temp;
543 }
544
545 return *this;
546} /* chain<euclidom>::swapnumbers */
547
548// --------------------------------------------------
549
550template <class euclidom>
551inline chain<euclidom> &chain<euclidom>::add (int_t n, euclidom e)
552{
553 // if the coefficient is zero, ignore the pair
554 if (e == 0)
555 return *this;
556
557 // find the position in the table for adding this pair
558 int_t i = 0;
559 while ((i < len) && (_n [i] < n))
560 ++ i;
561
562 // if an element with this identifier was found, add the coefficients
563 if ((i < len) && (_n [i] == n))
564 {
565 // add the coefficient
566 _e [i] += e;
567
568 // if the coefficient became zero, remove this pair
569 if (_e [i] == 0)
570 return removepair (i);
571
572 // otherwise we are done
573 else
574 return *this;
575 }
576
577 // otherwise insert this pair into the chain
578 return insertpair (i, n, e);
579
580} /* chain<euclidom>::add */
581
582template <class euclidom>
583chain<euclidom> &chain<euclidom>::remove (int_t n)
584{
585 // find the element of the chain to be removed
586 int_t i = 0;
587 while ((i < len) && (_n [i] < n))
588 ++ i;
589
590 // if found, then remove it
591 if ((i < len) && (_n [i] == n))
592 return removepair (i);
593
594 return *this;
595} /* chain<euclidom>::remove */
596
597template <class euclidom>
598inline chain<euclidom> &chain<euclidom>::add (const chain<euclidom> &other,
599 euclidom e, int_t number, chain<euclidom> *table)
600{
601 // if the coefficient is zero or the other chain is zero,
602 // then there is nothing to do
603 if ((e == 0) || !other. len)
604 return *this;
605
606 // prepare big tables for the new chain
607 int_t tablen = len + other. len;
608 int_t *bigntab = new int_t [tablen];
609 euclidom *bigetab = new euclidom [tablen];
610 if (!bigntab || !bigetab)
611 throw "Not enough memory to add chains.";
612
613 // prepare the counters of elements of the two input chains
614 // and of the output chain
615 int_t i = 0, j = 0, k = 0;
616
617 // go through both input chains and compute the output chain
618 while ((i < len) || (j < other. len))
619 {
620 // if the first chain has already been processed
621 if (i >= len)
622 {
623 bigntab [k] = other. _n [j];
624 bigetab [k] = e * other. _e [j];
625 ++ j;
626 if (table)
627 {
628 table [bigntab [k]]. add (number,
629 bigetab [k]);
630 }
631 ++ k;
632 }
633 // if an element from the first chain should be taken
634 else if ((j >= other. len) || (_n [i] < other. _n [j]))
635 {
636 bigntab [k] = _n [i];
637 bigetab [k] = _e [i];
638 ++ i;
639 ++ k;
640 }
641 // if an element from the other chain should be taken
642 else if (_n [i] > other. _n [j])
643 {
644 bigntab [k] = other. _n [j];
645 bigetab [k] = e * other. _e [j];
646 ++ j;
647 if (table)
648 {
649 table [bigntab [k]]. add (number,
650 bigetab [k]);
651 }
652 ++ k;
653 }
654 else // if (_n [i] == other. _n [j])
655 {
656 bigntab [k] = _n [i];
657 euclidom addelem = e * other. _e [j];
658 ++ j;
659 bigetab [k] = _e [i] + addelem;
660 ++ i;
661 if (!(bigetab [k] == 0))
662 {
663 if (table)
664 {
665 table [bigntab [k]]. add (number,
666 addelem);
667 }
668 ++ k;
669 }
670 else if (table)
671 {
672 table [bigntab [k]]. remove (number);
673 }
674 }
675 }
676
677 // release the old tables if they are useless now
678 if ((k != len) || (k == tablen))
679 {
680 if (len)
681 {
682 delete [] _n;
683 delete [] _e;
684 }
685 }
686
687 // use the previous tables and release the big table if beneficial
688 if ((k == len) && (k != tablen))
689 {
690 for (int_t i = 0; i < len; ++ i)
691 {
692 _n [i] = bigntab [i];
693 _e [i] = bigetab [i];
694 }
695 delete [] bigntab;
696 delete [] bigetab;
697 return *this;
698 }
699
700 // if the new chain is empty then accomodate for this
701 if (k == 0)
702 {
703 len = 0;
704 delete [] bigntab;
705 delete [] bigetab;
706 return *this;
707 }
708
709 // if the big tables cannot be used, allocate new tables
710 if (k != tablen)
711 {
712 _n = new int_t [k];
713 _e = new euclidom [k];
714 if (!_n || !_e)
715 throw "Cannot shorten a sum of chains.";
716 len = k;
717 for (int_t i = 0; i < len; ++ i)
718 {
719 _n [i] = bigntab [i];
720 _e [i] = bigetab [i];
721 }
722 delete [] bigntab;
723 delete [] bigetab;
724 }
725
726 // otherwise, simply use the big tables
727 else
728 {
729 len = k;
730 _n = bigntab;
731 _e = bigetab;
732 }
733
734 return *this;
735} /* chain<euclidom>::add */
736
737template <class euclidom>
738inline chain<euclidom> &chain<euclidom>::swap (chain<euclidom> &other,
739 int_t number, int_t othernumber, chain<euclidom> *table)
740{
741 // swap the data of the chains
742 swapelements (_n, other. _n);
743 swapelements (_e, other. _e);
744 swapelements (len, other. len);
745
746 if (!table)
747 return *this;
748
749 // change the numbers in every relevant entry of the table
750 int_t i = 0;
751 int_t j = 0;
752 while ((i < len) || (j < other. len))
753 {
754 // determine which table entry should be modified
755 int_t n;
756 if (i >= len)
757 n = other. _n [j ++];
758 else if (j >= other. len)
759 n = _n [i ++];
760 else if (_n [i] < other. _n [j])
761 n = _n [i ++];
762 else if (other. _n [j] < _n [i])
763 n = other. _n [j ++];
764 else
765 {
766 n = _n [i ++];
767 ++ j;
768 // ++ i;
769 // ++ j;
770 // continue;
771 }
772
773 // swap numbers in that table entry
774 table [n]. swapnumbers (othernumber, number);
775 }
776
777 return *this;
778} /* chain<euclidom>::swap */
779
780template <class euclidom>
781inline chain<euclidom> &chain<euclidom>::take (chain<euclidom> &c)
782{
783 // release the current tables if they were allocated
784 if (len)
785 {
786 delete [] _n;
787 delete [] _e;
788 }
789
790 _n = c. _n;
791 _e = c. _e;
792
793 // copy the length and reset the other length
794 len = c. len;
795 c. len = 0;
796
797 return *this;
798} /* chain<euclidom>::take */
799
800template <class euclidom>
801inline chain<euclidom> &chain<euclidom>::multiply (euclidom e, int_t number)
802{
803 // if there is only one element to be multiplied, find it and do it
804 if (number >= 0)
805 {
806 for (int_t i = 0; i < len; ++ i)
807 {
808 if (_n [i] == number)
809 {
810 if (e == 0)
811 removepair (i);
812 else
813 {
814 _e [i] *= e;
815 // if (_e [i] == 0)
816 // removepair (i);
817 }
818 return *this;
819 }
820 }
821 }
822
823 // if the entire chain has to be multiplied by a non-zero number...
824 else if (!(e == 0))
825 {
826 for (int_t i = 0; i < len; ++ i)
827 {
828 _e [i] *= e;
829 // if (_e [i] == 0)
830 // removepair (i);
831 }
832 }
833
834 // otherwise, if the chain has to be made zero, clean it
835 else
836 {
837 if (len)
838 {
839 delete [] _n;
840 delete [] _e;
841 }
842 len = 0;
843 }
844
845 return *this;
846} /* chain<euclidom>::multiply */
847
848template <class euclidom>
849inline outputstream &chain<euclidom>::show (outputstream &out,
850 const char *label) const
851{
852 if (len <= 0)
853 out << "0";
854 for (int_t i = 0; i < len; ++ i)
855 {
856 euclidom e = _e [i];
857 int_t n = _n [i] + 1;
858
859 if (e == 1)
860 out << (i ? " + " : "") <<
861 (label ? label : "") << n;
862 else if (-e == 1)
863 out << (i ? " - " : "-") <<
864 (label ? label : "") << n;
865 else
866 out << (i ? " + " : "") << e << " * " <<
867 (label ? label : "") << n;
868 }
869 return out;
870} /* chain<euclidom>::show */
871
872template <class euclidom>
873inline std::ostream &chain<euclidom>::show (std::ostream &out, const char *label) const
874{
875 outputstream tout (out);
876 show (tout, label);
877 return out;
878} /* chain<euclidom>::show */
879
880// --------------------------------------------------
881
884template <class euclidom>
885inline std::ostream &operator << (std::ostream &out, const chain<euclidom> &c)
886{
887 c. show (out);
888 return out;
889} /* operator << */
890
893template <class euclidom>
894inline std::istream &operator >> (std::istream &in, chain<euclidom> &c)
895{
896 ignorecomments (in);
897 int_t closing = readparenthesis (in);
898
899 ignorecomments (in);
900 while (in. peek () != closing)
901 {
902 // read the coefficient
903 euclidom e (1);
904 in >> e;
905
906 // read the multiplication symbol
907 ignorecomments (in);
908 if (in. peek () != '*')
909 throw "The multiplication sign '*' while reading a chain.";
910 in. get ();
911
912 // read the identifier
913 ignorecomments (in);
914 int_t n;
915 in >> n;
916 -- n;
917
918 // if the coefficient is zero, then this is wrong
919 if (e == 0)
920 throw "A zero coefficient in a chain detected.";
921
922 // add this element to the chain
923 c. add (e, n);
924
925 // prepare for the next pair to read
926 ignorecomments (in);
927
928 // read a comma or a plus between two elements of the chain
929 if ((in. peek () == ',') || (in. peek () == '+'))
930 {
931 in. get ();
932 ignorecomments (in);
933 }
934 }
935
936 if (closing != EOF)
937 in. get ();
938
939 return in;
940} /* operator >> */
941
942
943// --------------------------------------------------
944// ------------------- simplelist -------------------
945// --------------------------------------------------
946
949template <class element>
950class simplelist
951{
952public:
954 simplelist ();
955
957 ~simplelist ();
958
960 void add (element &m);
961
963 void remove (element &m);
964
970 element *take ();
971
972private:
974 simplelist (const simplelist<element> &s)
975 {
976 throw "Copying constructor not implemented "
977 "for a simple list.";
978 return;
979 }
980
982 simplelist<element> &operator = (const simplelist<element> &s)
983 {
984 throw "Operator = not implemented "
985 "for a simple list.";
986 return *this;
987 }
988
990 int num;
991
993 int cur;
994
996 element **elem;
997
998}; /* class simplelist */
999
1000// --------------------------------------------------
1001
1002template <class element>
1003inline simplelist<element>::simplelist (): num (0), cur (0), elem (NULL)
1004{
1005 return;
1006} /* simplelist::simplelist */
1007
1008template <class element>
1009inline simplelist<element>::~simplelist ()
1010{
1011 if (elem)
1012 delete [] elem;
1013 return;
1014} /* simplelist::~simplelist */
1015
1016template <class element>
1017inline void simplelist<element>::add (element &m)
1018{
1019 element **newelem = new element * [num + 1];
1020 for (int i = 0; i < num; ++ i)
1021 newelem [i] = elem [i];
1022 newelem [num ++] = &m;
1023 delete [] elem;
1024 elem = newelem;
1025 return;
1026} /* simplelist::add */
1027
1028template <class element>
1029inline void simplelist<element>::remove (element &m)
1030{
1031 int pos = 0;
1032 while ((pos < num) && (elem [pos] != &m))
1033 ++ pos;
1034 -- num;
1035 while (pos < num)
1036 {
1037 elem [pos] = elem [pos + 1];
1038 ++ pos;
1039 }
1040 return;
1041} /* simplelist::remove */
1042
1043template <class element>
1044inline element *simplelist<element>::take ()
1045{
1046 if (cur >= num)
1047 {
1048 cur = 0;
1049 return NULL;
1050 }
1051 else
1052 {
1053 return elem [cur ++];
1054 }
1055} /* simplelist::take */
1056
1057
1058// --------------------------------------------------
1059// -------------------- mmatrix ---------------------
1060// --------------------------------------------------
1061
1065template <class euclidom>
1066class mmatrix
1067{
1068public:
1070 mmatrix ();
1071
1074 mmatrix (const mmatrix<euclidom> &m);
1075
1078 mmatrix<euclidom> &operator = (const mmatrix<euclidom> &s);
1079
1081 ~mmatrix ();
1082
1085 void define (int_t numrows, int_t numcols);
1086
1088 void identity (int_t size);
1089
1091 void add (int_t row, int_t col, const euclidom &e);
1092
1096 euclidom get (int_t row, int_t col) const;
1097
1099 const chain<euclidom> &getrow (int_t n) const;
1100
1102 const chain<euclidom> &getcol (int_t n) const;
1103
1105 int_t getnrows () const;
1106
1108 int_t getncols () const;
1109
1112 void addrow (int_t dest, int_t source, const euclidom &e);
1113
1116 void addcol (int_t dest, int_t source, const euclidom &e);
1117
1120 void swaprows (int_t i, int_t j);
1121
1124 void swapcols (int_t i, int_t j);
1125
1127 void multiplyrow (int_t n, const euclidom &e);
1128
1130 void multiplycol (int_t n, const euclidom &e);
1131
1137 int_t findrow (int_t req_elements = 1, int_t start = -1) const;
1138
1144 int_t findcol (int_t req_elements = 1, int_t start = -1) const;
1145
1151 int_t reducerow (int_t n, int_t preferred);
1152
1158 int_t reducecol (int_t n, int_t preferred);
1159
1164 void simple_reductions (bool quiet = false);
1165
1168 void simple_form (bool quiet = false);
1169
1180 int_t arrange_towards_SNF (int_t *invertible_count = 0);
1181
1189 void simple_form_to_SNF (bool quiet = false);
1190
1192 void invert (void);
1193
1196 void multiply (const mmatrix<euclidom> &m1,
1197 const mmatrix<euclidom> &m2);
1198
1205 simplelist<mmatrix<euclidom> > dom_dom, dom_img, img_dom,
1206 img_img;
1207
1211 void submatrix (const mmatrix<euclidom> &matr,
1212 const chain<euclidom> &domain,
1213 const chain<euclidom> &range);
1214
1217 outputstream &show_hom_col (outputstream &out, int_t col,
1218 const chain<euclidom> &range,
1219 const char *txt = NULL) const;
1220
1223 std::ostream &show_hom_col (std::ostream &out, int_t col,
1224 const chain<euclidom> &range,
1225 const char *txt = NULL) const;
1226
1228 outputstream &showrowscols (outputstream &out,
1229 chain<euclidom> *table, int_t tablen,
1230 int_t first = 0, int_t howmany = 0,
1231 const char *label = NULL) const;
1232
1234 outputstream &showrows (outputstream &out, int_t first = 0,
1235 int_t howmany = 0, const char *label = "Row ") const;
1236
1238 std::ostream &showrows (std::ostream &out, int_t first = 0,
1239 int_t howmany = 0, const char *label = "Row ") const;
1240
1242 outputstream &showcols (outputstream &out, int_t first = 0,
1243 int_t howmany = 0, const char *label = "Col ") const;
1244
1246 std::ostream &showcols (std::ostream &out, int_t first = 0,
1247 int_t howmany = 0, const char *label = "Col ") const;
1248
1250 outputstream &showmap (outputstream &out,
1251 const char *maplabel = NULL,
1252 const char *xlabel = NULL, const char *ylabel = NULL)
1253 const;
1254
1256 std::ostream &showmap (std::ostream &out, const char *maplabel = NULL,
1257 const char *xlabel = NULL, const char *ylabel = NULL)
1258 const;
1259
1260private:
1262 int_t nrows;
1263
1265 int_t ncols;
1266
1268 int_t allrows;
1269
1271 int_t allcols;
1272
1274 chain<euclidom> *rows;
1275
1277 chain<euclidom> *cols;
1278
1281 int_t findrowcol (int_t req_elements, int_t start, int which) const;
1282
1285 void increase (int_t numrows, int_t numcols);
1286
1288 void increaserows (int_t numrows);
1289
1292 void increasecols (int_t numcols);
1293
1300 void division_SNF_correction (const euclidom &a, int_t pos1,
1301 const euclidom &b, int_t pos2);
1302
1309 static void extendedGCD (const euclidom &a, const euclidom &b,
1310 euclidom &x, euclidom &y);
1311
1321 void mult_left (const int_t setRows [],
1322 const int_t setCols [],
1323 const mmatrix<euclidom> &M,
1324 const mmatrix<euclidom> &invM,
1325 bool update_linked);
1326
1336 void mult_right (const int_t setRows [],
1337 const int_t setCols [],
1338 const mmatrix<euclidom> &M,
1339 const mmatrix<euclidom> &invM,
1340 bool update_linked);
1341
1342}; /* class mmatrix */
1343
1344// --------------------------------------------------
1345
1346template <class euclidom>
1347inline mmatrix<euclidom>::mmatrix (): nrows (0), ncols (0),
1348 allrows (0), allcols (0), rows (NULL), cols (NULL)
1349{
1350 return;
1351} /* mmatrix<euclidom>::mmatrix */
1352
1353template <class euclidom>
1354inline mmatrix<euclidom>::~mmatrix ()
1355{
1356 if (rows)
1357 delete [] rows;
1358 if (cols)
1359 delete [] cols;
1360 return;
1361} /* mmatrix<euclidom>::~mmatrix */
1362
1363template <class euclidom>
1364inline void mmatrix<euclidom>::define (int_t numrows, int_t numcols)
1365{
1366 // verify that no nonzero entry will be thrown away
1367 if ((nrows > numrows) || (ncols > numcols))
1368 throw "Trying to define a matrix smaller than it really is";
1369
1370 // increase the size of the matrix to fit the defined one
1371 increase (numrows, numcols);
1372
1373 // set the number of rows and the number of columns as requested
1374 nrows = numrows;
1375 ncols = numcols;
1376
1377 return;
1378} /* mmatrix<euclidom>::define */
1379
1380template <class euclidom>
1381inline mmatrix<euclidom>::mmatrix (const mmatrix<euclidom> &m)
1382{
1383 nrows = m.nrows;
1384 ncols = m.ncols;
1385 allrows = m.allrows;
1386 allcols = m.allcols;
1387
1388 rows = NULL;
1389 cols = NULL;
1390 if (m. allrows > 0)
1391 {
1392 chain<euclidom> *newrows = new chain<euclidom> [m. allrows];
1393 if (!newrows)
1394 throw "Not enough memory for matrix rows.";
1395 for (int_t i = 0; i < m. allrows; ++ i)
1396 newrows [i] = m. rows [i];
1397 rows = newrows;
1398 }
1399
1400 if (m. allcols > 0)
1401 {
1402 chain<euclidom> *newcols = new chain<euclidom> [m. allcols];
1403 if (!newcols)
1404 throw "Not enough memory for matrix columns.";
1405 for (int_t i = 0; i < m.allcols; ++ i)
1406 newcols [i] = m. cols [i];
1407 cols = newcols;
1408 }
1409} /* mmatrix<euclidom>::mmatrix */
1410
1411template <class euclidom>
1412inline mmatrix<euclidom> &mmatrix<euclidom>::operator =
1413 (const mmatrix<euclidom> &m)
1414{
1415 // first release allocated tables if any
1416 if (rows)
1417 delete [] rows;
1418 if (cols)
1419 delete [] cols;
1420
1421 nrows = m. nrows;
1422 ncols = m. ncols;
1423 allrows = m. allrows;
1424 allcols = m. allcols;
1425
1426 rows = NULL;
1427 cols = NULL;
1428 if (m. allrows > 0)
1429 {
1430 chain<euclidom> *newrows = new chain<euclidom> [m. allrows];
1431 if (!newrows)
1432 throw "Not enough memory for matrix rows.";
1433 for (int_t i = 0; i < m. allrows; ++ i)
1434 newrows [i] = m. rows [i];
1435 rows = newrows;
1436 }
1437
1438 if (m. allcols > 0)
1439 {
1440 chain<euclidom> *newcols = new chain<euclidom> [m. allcols];
1441 if (!newcols)
1442 throw "Not enough memory for matrix columns.";
1443 for (int_t i = 0; i < m.allcols; ++ i)
1444 newcols [i] = m. cols [i];
1445 cols = newcols;
1446 }
1447
1448 return *this;
1449} /* mmatrix<euclidom>::operator = */
1450
1451template <class euclidom>
1452inline void mmatrix<euclidom>::identity (int_t size)
1453{
1454 if (!nrows && !ncols)
1455 increase (size, size);
1456 else if ((size > nrows) || (size > ncols))
1457 size = (nrows < ncols) ? nrows : ncols;
1458 for (int_t i = 0; i < size; ++ i)
1459 {
1460 euclidom one;
1461 one = 1;
1462 add (i, i, one);
1463 }
1464 return;
1465} /* mmatrix<euclidom>::identity */
1466
1467template <class euclidom>
1468inline void mmatrix<euclidom>::add (int_t row, int_t col, const euclidom &e)
1469// A [r] [c] += e;
1470{
1471 if (row < 0)
1472 throw "Incorrect row number.";
1473 if (col < 0)
1474 throw "Incorrect column number.";
1475 if (row >= nrows)
1476 {
1477 if (row >= allrows)
1478 increaserows (row + row / 2 + 1);
1479 nrows = row + 1;
1480 }
1481 if (col >= ncols)
1482 {
1483 if (col >= allcols)
1484 increasecols (col + col / 2 + 1);
1485 ncols = col + 1;
1486 }
1487 if (e == 0)
1488 return;
1489 cols [col]. add (row, e);
1490 rows [row]. add (col, e);
1491 return;
1492} /* mmatrix<euclidom>::add */
1493
1494template <class euclidom>
1495inline euclidom mmatrix<euclidom>::get (int_t row, int_t col) const
1496// return (A [r] [c]);
1497{
1498 if ((row >= nrows) || (col >= ncols))
1499 {
1500 euclidom zero;
1501 zero = 0;
1502 return zero;
1503 }
1504 if (row >= 0)
1505 return rows [row]. getcoefficient (col);
1506 else if (col >= 0)
1507 return cols [col]. getcoefficient (row);
1508 else
1509 {
1510 euclidom zero;
1511 zero = 0;
1512 return zero;
1513 }
1514} /* mmatrix<euclidom>::get */
1515
1516template <class euclidom>
1517inline const chain<euclidom> &mmatrix<euclidom>::getrow (int_t n) const
1518{
1519 if ((n < 0) || (n >= nrows))
1520 throw "Incorrect row number.";
1521 return rows [n];
1522} /* mmatrix<euclidom>::getrow */
1523
1524template <class euclidom>
1525inline const chain<euclidom> &mmatrix<euclidom>::getcol (int_t n) const
1526{
1527 if ((n < 0) || (n >= ncols))
1528 throw "Incorrect column number.";
1529 return cols [n];
1530} /* mmatrix<euclidom>::getcol */
1531
1532template <class euclidom>
1533inline int_t mmatrix<euclidom>::getnrows () const
1534{
1535 return nrows;
1536} /* mmatrix<euclidom>::getnrows */
1537
1538template <class euclidom>
1539inline int_t mmatrix<euclidom>::getncols () const
1540{
1541 return ncols;
1542} /* mmatrix<euclidom>::getncols */
1543
1544template <class euclidom>
1545inline void mmatrix<euclidom>::addrow (int_t dest, int_t source,
1546 const euclidom &e)
1547{
1548 // check if the parameters are not out of range
1549 if ((dest < 0) || (dest >= nrows) || (source < 0) ||
1550 (source >= nrows))
1551 throw "Trying to add rows out of range.";
1552
1553 // add this row
1554 rows [dest]. add (rows [source], e, dest, cols);
1555
1556 // update the other matrices
1557 mmatrix<euclidom> *m;
1558 while ((m = img_img. take ()) != NULL)
1559 if (m -> rows)
1560 m -> rows [dest]. add (m -> rows [source], e,
1561 dest, m -> cols);
1562
1563 while ((m = img_dom. take ()) != NULL)
1564 if (m -> cols)
1565 m -> cols [source]. add (m -> cols [dest], -e,
1566 source, m -> rows);
1567
1568 return;
1569} /* mmatrix<euclidom>::addrow */
1570
1571template <class euclidom>
1572inline void mmatrix<euclidom>::addcol (int_t dest, int_t source,
1573 const euclidom &e)
1574{
1575 // check if the parameters are not out of range
1576 if ((dest < 0) || (dest >= ncols) || (source < 0) ||
1577 (source >= ncols))
1578 throw "Trying to add columns out of range.";
1579
1580 // add this column
1581 cols [dest]. add (cols [source], e, dest, rows);
1582
1583 // update the other matrices
1584 mmatrix<euclidom> *m;
1585 while ((m = dom_dom. take ()) != NULL)
1586 if (m -> cols)
1587 m -> cols [dest]. add (m -> cols [source], e,
1588 dest, m -> rows);
1589
1590 while ((m = dom_img. take ()) != NULL)
1591 if (m -> rows)
1592 m -> rows [source]. add (m -> rows [dest], -e,
1593 source, m -> cols);
1594
1595 return;
1596} /* mmatrix<euclidom>::addcol */
1597
1598template <class euclidom>
1599inline void mmatrix<euclidom>::swaprows (int_t i, int_t j)
1600{
1601 // in the trivial case nothing needs to be done
1602 if (i == j)
1603 return;
1604
1605 // check if the parameters are not out of range
1606 if ((i < 0) || (i >= nrows) || (j < 0) || (j >= nrows))
1607 throw "Trying to swap rows out of range.";
1608
1609 // swap the rows
1610 rows [i]. swap (rows [j], i, j, cols);
1611
1612 // update the other matrices
1613 mmatrix<euclidom> *m;
1614 while ((m = img_img. take ()) != NULL)
1615 if ((m -> rows) && (m -> nrows))
1616 m -> rows [i]. swap (m -> rows [j], i, j, m -> cols);
1617
1618 while ((m = img_dom. take ()) != NULL)
1619 if ((m -> cols) && (m -> ncols))
1620 m -> cols [i]. swap (m -> cols [j], i, j, m -> rows);
1621
1622 return;
1623} /* mmatrix<euclidom>::swaprows */
1624
1625template <class euclidom>
1626inline void mmatrix<euclidom>::swapcols (int_t i, int_t j)
1627{
1628 // in the trivial case nothing needs to be done
1629 if (i == j)
1630 return;
1631
1632 // check if the parameters are not out of range
1633 if ((i < 0) || (i >= ncols) || (j < 0) || (j >= ncols))
1634 throw "Trying to swap cols out of range.";
1635
1636 // swap the columns
1637 cols [i]. swap (cols [j], i, j, rows);
1638
1639 // update the other matrices
1640 mmatrix<euclidom> *m;
1641 while ((m = dom_dom. take ()) != NULL)
1642 if ((m -> cols) && (m -> ncols))
1643 m -> cols [i]. swap (m -> cols [j], i, j, m -> rows);
1644
1645 while ((m = dom_img. take ()) != NULL)
1646 if ((m -> rows) && (m -> nrows))
1647 m -> rows [i]. swap (m -> rows [j], i, j, m -> cols);
1648
1649 return;
1650} /* mmatrix<euclidom>::swapcols */
1651
1652template <class euclidom>
1653inline void mmatrix<euclidom>::multiplyrow (int_t n, const euclidom &e)
1654{
1655 // retrieve the row
1656 chain<euclidom> &therow = rows [n];
1657
1658 // multiply the row
1659 therow. multiply (e);
1660
1661 // multiply the corresponding entries in all the columns
1662 for (int_t i = 0; i < therow. size (); ++ i)
1663 cols [therow. num (i)]. multiply (e, n);
1664
1665 return;
1666} /* mmatrix<euclidom>::multiplyrow */
1667
1668template <class euclidom>
1669inline void mmatrix<euclidom>::multiplycol (int_t n, const euclidom &e)
1670{
1671 // retrieve the row
1672 chain<euclidom> &thecol = cols [n];
1673
1674 // multiply the row
1675 thecol. multiply (e);
1676
1677 // multiply the corresponding entries in all the rows
1678 for (int_t i = 0; i < thecol. size (); ++ i)
1679 rows [thecol. num (i)]. multiply (e, n);
1680
1681 return;
1682} /* mmatrix<euclidom>::multiplycol */
1683
1684template <class euclidom>
1685inline int_t mmatrix<euclidom>::findrowcol (int_t req_elements, int_t start,
1686 int which) const
1687{
1688 // start at the starting point
1689 int_t i = start;
1690 int_t random_i = -1;
1691
1692 // set loop to none
1693 int_t loopcounter = 0;
1694
1695 // if a random start is requested, initialize it and set loop to 1
1696 if (start < 0)
1697 {
1698 i = random_i = std::rand () % (which ? nrows : ncols);
1699 loopcounter = 1;
1700 }
1701
1702 // select some candidates
1703 int_t candidate = -1;
1704 int_t candidates_left = 3;
1705
1706 // if the table has one of its dimensions trivial, return the result
1707 if (which ? !ncols : !nrows)
1708 return (((req_elements > 0) ||
1709 (i >= (which ? nrows : ncols))) ? -1 : i);
1710
1711 // while the current position has not gone beyond the table
1712 while (i < (which ? nrows : ncols))
1713 {
1714 // if there is an appropriate row/column, return its number
1715 int_t l = (which ? rows : cols) [i]. size ();
1716 if ((req_elements >= 0) && (l >= req_elements))
1717 return i;
1718 else if ((req_elements < 0) && !l)
1719 {
1720 if (!candidates_left || !(which ? rows : cols) [i].
1721 contains_non_invertible ())
1722 return i;
1723 else
1724 {
1725 candidate = i;
1726 -- candidates_left;
1727 if (start < 0)
1728 {
1729 random_i = std::rand () %
1730 (which ? nrows : ncols);
1731 i = random_i - 1;
1732 loopcounter = 1;
1733 }
1734 }
1735 }
1736
1737 // otherwise increase the counter and rewind to 0 if needed
1738 if (++ i >= (which ? nrows : ncols))
1739 if (loopcounter --)
1740 i = 0;
1741
1742 // if the searching was started at a random position,
1743 // finish earlier
1744 if ((random_i >= 0) && !loopcounter && (i >= random_i))
1745 break;
1746 }
1747
1748 // if not found, return the recent candidate (or -1 if none)
1749 return candidate;
1750} /* mmatrix<euclidom>::findrowcol */
1751
1752template <class euclidom>
1753inline int_t mmatrix<euclidom>::findrow (int_t req_elements, int_t start) const
1754{
1755 return findrowcol (req_elements, start, 1);
1756} /* mmatrix<euclidom>::findrow */
1757
1758template <class euclidom>
1759inline int_t mmatrix<euclidom>::findcol (int_t req_elements, int_t start) const
1760{
1761 return findrowcol (req_elements, start, 0);
1762} /* mmatrix<euclidom>::findcol */
1763
1764template <class euclidom>
1765inline int_t mmatrix<euclidom>::reducerow (int_t n, int_t preferred)
1766{
1767 if (n >= nrows)
1768 throw "Trying to reduce a row out of range.";
1769
1770 int_t the_other = -1;
1771
1772 // repeat until the row contains at most one nonzero entry
1773 int_t len;
1774 while ((len = rows [n]. size ()) > 1)
1775 {
1776 // copy the row to a local structure
1777 chain<euclidom> local (rows [n]);
1778
1779 // find the best element in this row
1780 int_t best_i = local. findbest (cols);
1781
1782 // find the number of the preferred element in the row
1783 int_t preferred_i = (preferred < 0) ? -1 :
1784 local. findnumber (preferred);
1785
1786 // check if the element in the preferred column
1787 // is equally good as the one already chosen
1788 if ((preferred_i >= 0) &&
1789 (local. coef (preferred_i). delta () ==
1790 local. coef (best_i). delta ()))
1791 best_i = preferred_i;
1792
1793 // remember the column number containing this element
1794 the_other = local. num (best_i);
1795
1796 // for each column
1797 for (int_t i = 0; i < len; ++ i)
1798 {
1799 // if this column is the chosen one, do nothing
1800 if (i == best_i)
1801 continue;
1802
1803 // compute the quotient of two elements
1804 euclidom quotient = local. coef (i) /
1805 local. coef (best_i);
1806
1807 // subtract the chosen column from the other one
1808 addcol (local. num (i), local. num (best_i),
1809 -quotient);
1810 }
1811 }
1812
1813 return the_other;
1814} /* mmatrix<euclidom>::reducerow */
1815
1816template <class euclidom>
1817inline int_t mmatrix<euclidom>::reducecol (int_t n, int_t preferred)
1818{
1819 if (n >= ncols)
1820 throw "Trying to reduce a column out of range.";
1821
1822 int_t the_other = -1;
1823
1824 // repeat until the column contains at most one nonzero entry
1825 int_t len;
1826 while ((len = cols [n]. size ()) > 1)
1827 {
1828 // copy the column to a local structure
1829 chain<euclidom> local (cols [n]);
1830
1831 // find the best element in this column
1832 int_t best_i = local. findbest (rows);
1833
1834 // find the number of the preferred element in the row
1835 int_t preferred_i = (preferred < 0) ? -1 :
1836 local. findnumber (preferred);
1837
1838 // check if the element in the preferred column
1839 // is equally good as the one already chosen
1840 if ((preferred_i >= 0) &&
1841 (local. coef (preferred_i). delta () ==
1842 local. coef (best_i). delta ()))
1843 best_i = preferred_i;
1844
1845 // remember the row number containing this element
1846 the_other = local. num (best_i);
1847
1848 // for each row
1849 for (int_t i = 0; i < len; ++ i)
1850 {
1851 // if this row is the chosen one, do nothing
1852 if (i == best_i)
1853 continue;
1854
1855 // compute the quotient of two elements
1856 euclidom quotient = local. coef (i) /
1857 local. coef (best_i);
1858
1859 // subtract the chosen row from the other one
1860 addrow (local. num (i), local. num (best_i),
1861 -quotient);
1862 }
1863 }
1864
1865 return the_other;
1866} /* mmatrix<euclidom>::reducecol */
1867
1868template <class euclidom>
1869inline outputstream &mmatrix<euclidom>::showrowscols (outputstream &out,
1870 chain<euclidom> *table, int_t tablen, int_t first, int_t howmany,
1871 const char *label) const
1872{
1873 if ((first < 0) || (first >= tablen))
1874 first = 0;
1875 if ((howmany <= 0) || (first + howmany > tablen))
1876 howmany = tablen - first;
1877 for (int_t i = 0; i < howmany; ++ i)
1878 out << (label ? label : "") << (first + i + 1) << ": " <<
1879 table [first + i] << '\n';
1880 out << '\n';
1881 return out;
1882} /* matrix<euclidom>::showrowscols */
1883
1884template <class euclidom>
1885inline outputstream &mmatrix<euclidom>::showrows (outputstream &out,
1886 int_t first, int_t howmany, const char *label) const
1887{
1888 return showrowscols (out, rows, nrows, first, howmany, label);
1889} /* mmatrix<euclidom>::showrows */
1890
1891template <class euclidom>
1892inline std::ostream &mmatrix<euclidom>::showrows (std::ostream &out,
1893 int_t first, int_t howmany, const char *label) const
1894{
1895 outputstream tout (out);
1896 showrows (tout, first, howmany, label);
1897 return out;
1898} /* mmatrix<euclidom>::showrows */
1899
1900template <class euclidom>
1901inline outputstream &mmatrix<euclidom>::showcols (outputstream &out,
1902 int_t first, int_t howmany, const char *label) const
1903{
1904 return showrowscols (out, cols, ncols, first, howmany, label);
1905} /* mmatrix<euclidom>::showcols */
1906
1907template <class euclidom>
1908inline std::ostream &mmatrix<euclidom>::showcols (std::ostream &out,
1909 int_t first, int_t howmany, const char *label) const
1910{
1911 outputstream tout (out);
1912 showcols (tout, first, howmany, label);
1913 return out;
1914} /* mmatrix<euclidom>::showcols */
1915
1916template <class euclidom>
1917inline outputstream &mmatrix<euclidom>::showmap (outputstream &out,
1918 const char *maplabel, const char *xlabel, const char *ylabel) const
1919{
1920 if (ncols <= 0)
1921 {
1922 if (maplabel && (maplabel [0] == '\t'))
1923 out << "\t0\n";
1924 else
1925 out << "0\n";
1926 }
1927 for (int_t i = 0; i < ncols; ++ i)
1928 {
1929 out << (maplabel ? maplabel : "f") << " (" <<
1930 (xlabel ? xlabel : "") << (i + 1) << ") = ";
1931 cols [i]. show (out, ylabel) << '\n';
1932 }
1933 return out;
1934} /* mmatrix<euclidom>::showmap */
1935
1936template <class euclidom>
1937inline std::ostream &mmatrix<euclidom>::showmap (std::ostream &out,
1938 const char *maplabel, const char *xlabel, const char *ylabel) const
1939{
1940 outputstream tout (out);
1941 showmap (tout, maplabel, xlabel, ylabel);
1942 return out;
1943} /* mmatrix<euclidom>::showmap */
1944
1945template <class euclidom>
1946inline void mmatrix<euclidom>::simple_reductions (bool quiet)
1947{
1948 // if the matrix has no rows or no columns,
1949 // simple reductions make no sense
1950 if (!nrows || !ncols)
1951 return;
1952
1953 // prepare counters
1954 long countreduced = 0;
1955 long count = 4 * ((nrows > ncols) ? ncols : nrows);
1956
1957 // prepare quazi-random numbers
1958 int_t nr = static_cast<int_t> (std::rand () % nrows);
1959 if (nr < 0)
1960 nr = 1 - nr;
1961 int_t nr_count = 0;
1962 int_t nr_add = 0;
1963
1964 // try quazi-random reductions
1965 while (count --)
1966 {
1967 // try a simple reduction
1968 if ((rows [nr]. size () == 1) &&
1969 (rows [nr]. coef (0). delta () == 1) &&
1970 (cols [rows [nr]. num (0)]. size () > 1))
1971 {
1972 ++ countreduced;
1973 reducecol (rows [nr]. num (0), -1);
1974 }
1975
1976 // try a new semi-random position
1977 if (nr_count)
1978 -- nr_count;
1979 else
1980 {
1981 nr_add = ((std::rand () >> 2) + 171) % nrows;
1982 if (nr_add < 1)
1983 nr_add = 1;
1984 nr_count = 100;
1985 }
1986 nr += nr_add;
1987 if (nr >= nrows)
1988 nr -= nrows;
1989
1990 // display a counter
1991 if (!quiet && !(count % 373))
1992 scon << std::setw (12) << count <<
1993 "\b\b\b\b\b\b\b\b\b\b\b\b";
1994 }
1995
1996 if (!quiet)
1997 sout << countreduced << " +";
1998
1999 return;
2000} /* mmatrix<euclidom>::simple_reductions */
2001
2002template <class euclidom>
2003inline void mmatrix<euclidom>::simple_form (bool quiet)
2004{
2005 // if the matrix has no rows or columns, it is already in simple form
2006 if (!nrows || !ncols)
2007 return;
2008
2009 // make some random simple reductions before proceeding
2010 simple_reductions (quiet);
2011
2012 // prepare a counter
2013 long count = 0;
2014
2015 // prepare variables for chosen row and column numbers
2016 int_t row, col, prev_row, prev_col;
2017
2018 // find a row or a column with at least two nonzero entries
2019 row = -1;
2020 col = findcol (2);
2021 prev_row = -1, prev_col = -1;
2022 if (col < 0)
2023 row = findrow (2);
2024
2025 // repeat while such a row or a column can be found
2026 while ((row >= 0) || (col >= 0))
2027 {
2028 // reduce rows and columns until a single entry remains
2029 while ((row >= 0) || (col >= 0))
2030 if (row >= 0)
2031 {
2032 col = reducerow (row, prev_col);
2033 prev_row = row;
2034 row = -1;
2035 }
2036 else if (col >= 0)
2037 {
2038 row = reducecol (col, prev_row);
2039 prev_col = col;
2040 col = -1;
2041 }
2042
2043 // update the counter and display it if requested to
2044 ++ count;
2045 if (!quiet && !(count % 373))
2046 scon << std::setw (12) << count <<
2047 "\b\b\b\b\b\b\b\b\b\b\b\b";
2048
2049 // find another row or column with at least 2 nonzero entries
2050 col = findcol (2);
2051 if (col < 0)
2052 row = findrow (2);
2053 }
2054
2055 if (!quiet)
2056 sout << " " << count << " reductions made. ";
2057
2058 return;
2059} /* mmatrix<euclidom>::simple_form */
2060
2061template <class euclidom>
2062inline int_t mmatrix<euclidom>::arrange_towards_SNF (int_t *invertible_count)
2063{
2064 // prepare a counter for nonzero entries of the diagonal
2065 int_t cur = 0;
2066
2067 // move all the invertible nonzero entries to the front
2068 for (int_t n = 0; n < ncols; ++ n)
2069 {
2070 if (cols [n]. empty ())
2071 continue;
2072 if (cols [n]. coef (0). delta () != 1)
2073 continue;
2074 int_t r = cols [n]. num (0);
2075 if (n != cur)
2076 swapcols (n, cur);
2077 if (r != cur)
2078 swaprows (r, cur);
2079 ++ cur;
2080 }
2081
2082 // store the number of invertible entries at the diagonal
2083 if (invertible_count)
2084 *invertible_count = cur;
2085
2086 // move all the remaining nonzero entries to the front
2087 for (int_t n = cur; n < ncols; ++ n)
2088 {
2089 if (cols [n]. empty ())
2090 continue;
2091 int_t r = cols [n]. num (0);
2092 if (n != cur)
2093 swapcols (n, cur);
2094 if (r != cur)
2095 swaprows (r, cur);
2096 ++ cur;
2097 }
2098
2099 return cur;
2100} /* mmatrix<euclidom>::arrange_towards_SNF */
2101
2102template <class euclidom>
2103inline void mmatrix<euclidom>::simple_form_to_SNF (bool quiet)
2104{
2105 // arrange towards SNF
2106 if (!quiet)
2107 sout << "Determining the 'diagonal'... ";
2108 int_t indexBegin = 0;
2109 int_t indexEnd = arrange_towards_SNF (&indexBegin);
2110 if (!quiet)
2111 {
2112 sout << indexBegin << " invertible + " <<
2113 (indexEnd - indexBegin) << " noninvertible coefs.\n";
2114 }
2115
2116 // check the division condition and make corrections where necessary
2117 if (!quiet)
2118 sout << "Correcting the division condition... ";
2119 int_t countCorrections = 0;
2120 bool divisionOK = false;
2121 euclidom zero;
2122 zero = 0;
2123 while (!divisionOK)
2124 {
2125 divisionOK = true;
2126 for (int_t index = indexBegin + 1; index < indexEnd; ++ index)
2127 {
2128 euclidom e1 (this -> get (index - 1, index - 1));
2129 euclidom e2 (this -> get (index, index));
2130 if (e2 % e1 == zero)
2131 continue;
2132 divisionOK = false;
2133 // sout << "\nDEBUG: " << e1 << " does not divide " <<
2134 // e2 << " at position " << (index - 1) << ".\n";
2135 division_SNF_correction (e1, index - 1, e2, index);
2136 ++ countCorrections;
2137 }
2138 }
2139 sout << countCorrections << " corrections.\n";
2140 return;
2141} /* mmatrix<euclidom>::simple_form_to_SNF */
2142
2143template <class euclidom>
2144inline void mmatrix<euclidom>::division_SNF_correction
2145 (const euclidom &a, int_t pos1, const euclidom &b, int_t pos2)
2146{
2147// sout << "=== DEBUG: division_SNF_correction. ===\n";
2148// sout << "DEBUG: Initial matrix:\n" << *this;
2149
2150 // compute the GCD and the Bezout's identity coefficients
2151 euclidom sigma;
2152 euclidom tau;
2153 extendedGCD (a, b, sigma, tau);
2154 euclidom beta (a * sigma + b * tau);
2155 euclidom alpha (a / beta);
2156 euclidom gamma (b / beta);
2157// sout << "DEBUG: a = " << a << ", b = " << b << ", beta = " << beta <<
2158// " = " << a << "*" << sigma << "+" << b << "*" << tau <<
2159// ", a/beta = " << alpha << ", b/beta = " << gamma <<
2160// ".\nsigma = " << sigma << ", tau = " << tau <<
2161// ", alpha = " << alpha << ", gamma = " << gamma << ".\n";
2162
2163 // add the second column to the first
2164 euclidom one;
2165 one = 1;
2166 this -> addcol (pos1, pos2, one);
2167
2168 // multiply the selected columns/rows by a special 2x2 matrix:
2169 // prepare a chain that defines the columns and rows to be affected
2170 int_t setRowsCols [2];
2171 setRowsCols [0] = pos1;
2172 setRowsCols [1] = pos2;
2173
2174 // prepare a matrix to multiply from the left
2175 mmatrix<euclidom> M;
2176 M. define (2, 2);
2177 M. add (0, 0, sigma);
2178 M. add (0, 1, tau);
2179 M. add (1, 0, -gamma);
2180 M. add (1, 1, alpha);
2181// sout << "DEBUG: Matrix M:\n" << M;
2182
2183 // prepare the inverse of this matrix
2184 mmatrix<euclidom> invM;
2185 invM. define (2, 2);
2186 invM. add (0, 0, alpha);
2187 invM. add (0, 1, -tau);
2188 invM. add (1, 0, gamma);
2189 invM. add (1, 1, sigma);
2190// sout << "DEBUG: Inverse of the matrix M:\n" << invM;
2191
2192 // multiply the current matrix from the left by the special matrix
2193 this -> mult_left (setRowsCols, setRowsCols, M, invM, true);
2194
2195 // add the first column to the second one with a special coefficient
2196 // to make this part of the matirx diagonal again
2197 this -> addcol (pos2, pos1, (-tau) * gamma);
2198
2199// sout << "DEBUG: Final matrix:\n" << *this;
2200 return;
2201} /* mmatrix<euclidom>::division_SNF_correction */
2202
2203template <class euclidom>
2204inline void mmatrix<euclidom>::mult_left
2205 (const int_t setRows [],
2206 const int_t setCols [],
2207 const mmatrix<euclidom> &M,
2208 const mmatrix<euclidom> &invM,
2209 bool update_linked)
2210{
2211 euclidom zero;
2212 zero = 0;
2213 euclidom one;
2214 one = 1;
2215
2216 // compute the possibly changed rows of the matrix
2217 int_t size = M. getnrows ();
2218 auto_array<chain<euclidom> > newRows_p (new chain<euclidom> [size]);
2219 chain<euclidom> *newRows = newRows_p. get ();
2220 for (int_t row = 0; row < size; ++ row)
2221 {
2222 // determine the numbers of columns to process
2223 chain<euclidom> affected;
2224 for (int_t col = 0; col < size; ++ col)
2225 {
2226 const chain<euclidom> &rowChain =
2227 this -> getrow (setCols [col]);
2228 int_t rowSize = rowChain. size ();
2229 for (int cur = 0; cur < rowSize; ++ cur)
2230 {
2231 int_t num (rowChain. num (cur));
2232 if (affected. findnumber (num) < 0)
2233 affected. add (num, one);
2234 }
2235 }
2236
2237 // multiply these columns by the current row of M
2238 int_t col_count = affected. size ();
2239 for (int_t col = 0; col < col_count; ++ col)
2240 {
2241 int_t col_nr = affected. num (col);
2242 euclidom e;
2243 e = 0;
2244 for (int_t k = 0; k < size; ++ k)
2245 {
2246 int_t row_k = setCols [k]; // sic!
2247 e += M. get (row, k) *
2248 this -> get (row_k, col_nr);
2249 }
2250 if (!(e == 0))
2251 newRows [row]. add (col_nr, e);
2252 }
2253 }
2254
2255 // replace the previous rows in the matrix with the new ones
2256 for (int_t row = 0; row < size; ++ row)
2257 {
2258 int_t row_nr = setRows [row];
2259 const chain<euclidom> &row_ch (newRows [row]);
2260
2261 // eliminate entries that do not appear in the new row
2262 const chain<euclidom> row_prev (this -> getrow (row_nr));
2263 int_t len_prev = row_prev. size ();
2264 for (int_t i = 0; i < len_prev; ++ i)
2265 {
2266 int_t col_nr (row_prev. num (i));
2267 euclidom e (row_ch. getcoefficient (col_nr));
2268 if (e == zero)
2269 {
2270 euclidom coef (this -> get (row_nr, col_nr));
2271 this -> add (row_nr, col_nr, -coef);
2272 }
2273 }
2274
2275 // set the other entries to the ones in the new row
2276 int_t len = row_ch. size ();
2277 for (int_t i = 0; i < len; ++ i)
2278 {
2279 int_t col_nr (row_ch. num (i));
2280 euclidom e (row_ch. coef (i) -
2281 this -> get (row_nr, col_nr));
2282 if (!(e == zero))
2283 this -> add (row_nr, col_nr, e);
2284 }
2285 }
2286
2287 // update the other matrices if necessary
2288 if (!update_linked)
2289 return;
2290 mmatrix<euclidom> *m;
2291 while ((m = img_img. take ()) != NULL)
2292 if ((m -> rows) && (m -> nrows))
2293 m -> mult_left (setRows, setCols, M, invM, false);
2294 while ((m = img_dom. take ()) != NULL)
2295 if ((m -> cols) && (m -> ncols))
2296 m -> mult_right (setCols, setRows, invM, M, false);
2297
2298 return;
2299} /* mmatrix<euclidom>::mult_left */
2300
2301template <class euclidom>
2302inline void mmatrix<euclidom>::mult_right
2303 (const int_t setRows [],
2304 const int_t setCols [],
2305 const mmatrix<euclidom> &M,
2306 const mmatrix<euclidom> &invM,
2307 bool update_linked)
2308{
2309 euclidom zero;
2310 zero = 0;
2311 euclidom one;
2312 one = 1;
2313
2314 // compute the possibly changed columns of the matrix
2315 int_t size = M. getncols ();
2316 auto_array<chain<euclidom> > newCols_p (new chain<euclidom> [size]);
2317 chain<euclidom> *newCols = newCols_p. get ();
2318 for (int_t col = 0; col < size; ++ col)
2319 {
2320 // determine the numbers of rows to process
2321 chain<euclidom> affected;
2322 for (int_t row = 0; row < size; ++ row)
2323 {
2324 const chain<euclidom> &colChain =
2325 this -> getcol (setRows [row]);
2326 int_t colSize = colChain. size ();
2327 for (int cur = 0; cur < colSize; ++ cur)
2328 {
2329 int_t num (colChain. num (cur));
2330 if (affected. findnumber (num) < 0)
2331 affected. add (num, one);
2332 }
2333 }
2334
2335 // multiply these columns by the current row of M
2336 int_t row_count = affected. size ();
2337 for (int_t row = 0; row < row_count; ++ row)
2338 {
2339 int_t row_nr = affected. num (row);
2340 euclidom e;
2341 e = 0;
2342 for (int_t k = 0; k < size; ++ k)
2343 {
2344 int_t col_k = setRows [k]; // sic!
2345 e += this -> get (row_nr, col_k) *
2346 M. get (k, col);
2347 }
2348 if (!(e == 0))
2349 newCols [col]. add (row_nr, e);
2350 }
2351 }
2352
2353 // replace the previous columns in the matrix with the new ones
2354 for (int_t col = 0; col < size; ++ col)
2355 {
2356 int_t col_nr = setCols [col];
2357 const chain<euclidom> &col_ch (newCols [col]);
2358
2359 // eliminate entries that do not appear in the new column
2360 const chain<euclidom> col_prev (this -> getcol (col_nr));
2361 int_t len_prev = col_prev. size ();
2362 for (int_t i = 0; i < len_prev; ++ i)
2363 {
2364 int_t row_nr (col_prev. num (i));
2365 euclidom e (col_ch. getcoefficient (row_nr));
2366 if (e == zero)
2367 {
2368 euclidom coef (this -> get (row_nr, col_nr));
2369 this -> add (row_nr, col_nr, -coef);
2370 }
2371 }
2372
2373 // set the other entries to the ones in the new column
2374 int_t len = col_ch. size ();
2375 for (int_t i = 0; i < len; ++ i)
2376 {
2377 int_t row_nr (col_ch. num (i));
2378 euclidom e (col_ch. coef (i) -
2379 this -> get (row_nr, col_nr));
2380 if (!(e == zero))
2381 this -> add (row_nr, col_nr, e);
2382 }
2383 }
2384
2385 // update the other matrices if necessary
2386 if (!update_linked)
2387 return;
2388 mmatrix<euclidom> *m;
2389 while ((m = dom_dom. take ()) != NULL)
2390 if ((m -> cols) && (m -> ncols))
2391 m -> mult_left (setRows, setCols, M, invM, false);
2392 while ((m = dom_img. take ()) != NULL)
2393 if ((m -> rows) && (m -> nrows))
2394 m -> mult_right (setCols, setRows, invM, M, false);
2395
2396 return;
2397} /* mmatrix<euclidom>::mult_right */
2398
2399template <class euclidom>
2400inline void mmatrix<euclidom>::extendedGCD
2401 (const euclidom &a, const euclidom &b, euclidom &x, euclidom &y)
2402{
2403 euclidom aa (a), bb (b);
2404 euclidom xx, yy, lastx, lasty;
2405 xx = 0;
2406 yy = 1;
2407 lastx = 1;
2408 lasty = 0;
2409 while (!(bb == 0))
2410 {
2411 euclidom quotient (aa / bb);
2412 euclidom remainder (aa % bb);
2413 aa = bb;
2414 bb = remainder;
2415 euclidom xxx = lastx - quotient * xx;
2416 lastx = xx;
2417 xx = xxx;
2418 euclidom yyy = lasty - quotient * yy;
2419 lasty = yy;
2420 yy = yyy;
2421 }
2422 x = lastx;
2423 y = lasty;
2424 return;
2425} /* mmatrix<euclidom>::extendedGCD */
2426
2427template <class euclidom>
2428inline void mmatrix<euclidom>::invert (void)
2429{
2430 // check if the matrix is square
2431 if (nrows != ncols)
2432 throw "Trying to invert a non-square matrix.";
2433
2434 // if the matrix is trivial, nothing needs to be done
2435 if (!nrows)
2436 return;
2437
2438 // create the identity matrix of the appropriate size
2439 mmatrix<euclidom> m;
2440 m. identity (ncols);
2441
2442 // transform the matrix to the identity
2443 // by row operations (swapping and adding)
2444 // and apply the same operations to the matrix 'm'
2445 for (int_t col = 0; col < ncols; ++ col)
2446 {
2447 // find an invertible element in the column
2448 int_t len = cols [col]. size ();
2449 if (len <= 0)
2450 {
2451 throw "Matrix inverting: Zero column found.";
2452 }
2453 int_t chosen = 0;
2454 while ((chosen < len) &&
2455 ((cols [col]. num (chosen) < col) ||
2456 (cols [col]. coef (chosen). delta () != 1)))
2457 {
2458 ++ chosen;
2459 }
2460 if (chosen >= len)
2461 {
2462 throw "Matrix inverting: "
2463 "No invertible element in a column.";
2464 }
2465
2466 // make the leading entry equal 1 in the chosen row
2467 euclidom invcoef;
2468 invcoef = 1;
2469 invcoef = invcoef / cols [col]. coef (chosen);
2470 m. multiplyrow (cols [col]. num (chosen), invcoef);
2471 multiplyrow (cols [col]. num (chosen), invcoef);
2472
2473 // move the chosen row to the diagonal position
2474 m. swaprows (col, cols [col]. num (chosen));
2475 swaprows (col, cols [col]. num (chosen));
2476
2477 // zero the column below and above the chosen entry
2478 invcoef = -1;
2479 for (int_t i = 0; i < len; ++ i)
2480 {
2481 if (cols [col]. num (i) == col)
2482 continue;
2483 euclidom coef = invcoef * cols [col]. coef (i);
2484 m. addrow (cols [col]. num (i), col, coef);
2485 addrow (cols [col]. num (i), col, coef);
2486 -- i;
2487 -- len;
2488 }
2489 }
2490
2491 // take the matrix 'm' as the result
2492 for (int_t i = 0; i < ncols; ++ i)
2493 {
2494 cols [i]. take (m. cols [i]);
2495 rows [i]. take (m. rows [i]);
2496 }
2497
2498 return;
2499} /* mmatrix<euclidom>::invert */
2500
2501template <class euclidom>
2502inline void mmatrix<euclidom>::multiply (const mmatrix<euclidom> &m1,
2503 const mmatrix<euclidom> &m2)
2504{
2505 if (m1. ncols != m2. nrows)
2506 throw "Trying to multiply matrices of wrong sizes.";
2507 int_t K = m1. ncols;
2508 define (m1. nrows, m2. ncols);
2509 for (int_t i = 0; i < nrows; ++ i)
2510 {
2511 for (int_t j = 0; j < ncols; ++ j)
2512 {
2513 euclidom e;
2514 e = 0;
2515 for (int_t k = 0; k < K; ++ k)
2516 e += m1. get (i, k) * m2. get (k, j);
2517 add (i, j, e);
2518 }
2519 }
2520 return;
2521} /* mmatrix<euclidom>::multiply */
2522
2523template <class euclidom>
2524inline void mmatrix<euclidom>::submatrix (const mmatrix<euclidom> &matr,
2525 const chain<euclidom> &domain, const chain<euclidom> &range)
2526{
2527 for (int_t i = 0; i < domain. size (); ++ i)
2528 {
2529 for (int_t j = 0; j < range. size (); ++ j)
2530 {
2531 euclidom e = matr. get (domain. num (i),
2532 range. num (j));
2533 if (!(e == 0))
2534 add (i, j, e);
2535 }
2536 }
2537
2538 return;
2539} /* mmatrix<euclidom>::submatrix */
2540
2541template <class euclidom>
2542inline outputstream &mmatrix<euclidom>::show_hom_col (outputstream &out,
2543 int_t col, const chain<euclidom> &range, const char *txt) const
2544{
2545 // keep current position in 'range'
2546 int_t j = 0;
2547
2548 // remember if this is the first displayed element
2549 int_t first = 1;
2550
2551 // go through the column of the matrix
2552 for (int_t i = 0; i < cols [col]. size (); ++ i)
2553 {
2554 // find the current element in the range
2555 while ((j < range. size ()) &&
2556 (range. num (j) < cols [col]. num (i)))
2557 {
2558 ++ j;
2559 }
2560
2561 // if this element was found in the range, display it
2562 if ((j < range. size ()) &&
2563 (range. num (j) == cols [col]. num (i)))
2564 {
2565 if (first)
2566 first = 0;
2567 else
2568 out << " + ";
2569 if (-cols [col]. coef (i) == 1)
2570 out << "-";
2571 else if (cols [col]. coef (i) != 1)
2572 out << cols [col]. coef (i) << " * ";
2573 if (txt)
2574 out << txt;
2575 out << (j + 1);
2576 }
2577 }
2578
2579 // if nothing was shown, display 0
2580 if (first)
2581 out << 0;
2582
2583 return out;
2584} /* mmatrix<euclidom>::show_hom_col */
2585
2586template <class euclidom>
2587inline std::ostream &mmatrix<euclidom>::show_hom_col (std::ostream &out, int_t col,
2588 const chain<euclidom> &range, const char *txt) const
2589{
2590 outputstream tout (out);
2591 show_hom_col (tout, col, range, txt);
2592 return out;
2593} /* mmatrix<euclidom>::show_hom_col */
2594
2595// --------------------------------------------------
2596
2597template <class euclidom>
2598inline void mmatrix<euclidom>::increase (int_t numrows, int_t numcols)
2599{
2600 increaserows (numrows);
2601 increasecols (numcols);
2602 return;
2603} /* mmatrix<euclidom>::increase */
2604
2605template <class euclidom>
2606inline void mmatrix<euclidom>::increaserows (int_t numrows)
2607{
2608 if (allrows >= numrows)
2609 return;
2610 chain<euclidom> *newrows = new chain<euclidom> [numrows];
2611 if (!newrows)
2612 throw "Not enough memory for matrix rows.";
2613 for (int_t i = 0; i < nrows; ++ i)
2614 newrows [i]. take (rows [i]);
2615 if (rows)
2616 delete [] rows;
2617 rows = newrows;
2618 allrows = numrows;
2619 return;
2620} /* mmatrix<euclidom>::increaserows */
2621
2622template <class euclidom>
2623inline void mmatrix<euclidom>::increasecols (int_t numcols)
2624{
2625 if (allcols >= numcols)
2626 return;
2627 chain<euclidom> *newcols = new chain<euclidom> [numcols];
2628 if (!newcols)
2629 throw "Not enough memory for matrix columns.";
2630 for (int_t i = 0; i < ncols; ++ i)
2631 newcols [i]. take (cols [i]);
2632 if (cols)
2633 delete [] cols;
2634 cols = newcols;
2635 allcols = numcols;
2636
2637 return;
2638} /* mmatrix<euclidom>::increasecols */
2639
2640// --------------------------------------------------
2641
2644template <class euclidom>
2645inline std::ostream &operator << (std::ostream &out,
2646 const mmatrix<euclidom> &m)
2647{
2648 return m. showcols (out);
2649} /* operator << */
2650
2651
2652// --------------------------------------------------
2653// ------------------ chaincomplex ------------------
2654// --------------------------------------------------
2655
2660template <class euclidom>
2661class chaincomplex
2662{
2663public:
2668 chaincomplex (int d, int trace_generators = 0, int trace_bases = 0);
2669
2671 ~chaincomplex ();
2672
2677 void def_gen (int q, int_t n);
2678
2681 void add (int q, int_t m, int_t n, const euclidom &e);
2682
2685 euclidom get (int q, int_t row, int_t col) const;
2686
2688 const mmatrix<euclidom> &getboundary (int i) const;
2689
2691 int_t getnumgen (int i) const;
2692
2694 int dim () const;
2695
2698 const chain<euclidom> &gethomgen (int q, int_t i) const;
2699
2701 const mmatrix<euclidom> &gethomgen (int q) const;
2702
2704 const mmatrix<euclidom> &getchgbasis (int q) const;
2705
2708 void simple_form (int which, bool quiet);
2709
2714 void simple_form (const int *level, bool quiet);
2715
2721 int simple_homology (chain<euclidom> &result, int which) const;
2722
2731 int simple_homology (chain<euclidom> *&result,
2732 const int *level = NULL) const;
2733
2737 void take_homology (const chain<euclidom> *hom_chain);
2738
2742 outputstream &show_homology (outputstream &out,
2743 const chain<euclidom> *hom,
2744 const int *level = NULL) const;
2745
2749 std::ostream &show_homology (std::ostream &out,
2750 const chain<euclidom> *hom,
2751 const int *level = NULL) const;
2752
2756 outputstream &show_generators (outputstream &out,
2757 const chain<euclidom> &list, int q) const;
2758
2762 std::ostream &show_generators (std::ostream &out,
2763 const chain<euclidom> &list, int q) const;
2764
2766 outputstream &compute_and_show_homology (outputstream &out,
2767 chain<euclidom> *&hom, const int *level = NULL);
2768
2770 std::ostream &compute_and_show_homology (std::ostream &out,
2771 chain<euclidom> *&hom, const int *level = NULL);
2772
2775 friend class chainmap<euclidom>;
2776
2777private:
2779 chaincomplex (const chaincomplex<euclidom> &m)
2780 {throw "Copy constructor not implemented "
2781 "for a chain complex.";}
2782
2784 chaincomplex<euclidom> &operator =
2785 (const chaincomplex<euclidom> &s)
2786 {throw "Operator = not implemented "
2787 "for a chain complex."; return *this;}
2788
2790 int len;
2791
2794 mmatrix<euclidom> *boundary;
2795
2798 mmatrix<euclidom> *generators;
2799
2801 mmatrix<euclidom> *chgbases;
2802
2805 int *generators_initialized;
2806
2809 int *chgbases_initialized;
2810
2813 int_t *numgen;
2814
2815}; /* class chaincomplex */
2816
2817// --------------------------------------------------
2818
2819template <class euclidom>
2820inline chaincomplex<euclidom>::chaincomplex (int d,
2821 int trace_generators, int trace_bases)
2822{
2823 // set the number of tables to be sufficient for 0 to d inclusive
2824 len = (d >= 0) ? (d + 1) : 0;
2825
2826 // create a table of empty matrices
2827 boundary = len ? new mmatrix<euclidom> [len] : NULL;
2828
2829 // create a table of matrices for tracing generators of homology
2830 generators = (trace_generators && len) ?
2831 new mmatrix<euclidom> [len] : NULL;
2832 if (generators)
2833 generators_initialized = new int [len];
2834 else
2835 generators_initialized = NULL;
2836
2837 // create a table of matrices for tacing base changes
2838 chgbases = (trace_bases && len) ?
2839 new mmatrix<euclidom> [len] : NULL;
2840 if (chgbases)
2841 chgbases_initialized = new int [len];
2842 else
2843 chgbases_initialized = NULL;
2844
2845 // create a table of generator numbers
2846 numgen = len ? new int_t [len] : NULL;
2847
2848 // link the boundary matrices to each other,
2849 // as well as to the generators and to the base change matrices
2850 for (int i = 0; i < len; ++ i)
2851 {
2852 numgen [i] = -1;
2853 if (i < len - 1)
2854 boundary [i]. dom_img. add (boundary [i + 1]);
2855 if (i > 0)
2856 boundary [i]. img_dom. add (boundary [i - 1]);
2857 if (generators)
2858 {
2859 boundary [i]. dom_dom. add (generators [i]);
2860 if (i > 0)
2861 {
2862 boundary [i]. img_dom. add
2863 (generators [i - 1]);
2864 }
2865 generators_initialized [i] = 0;
2866 }
2867 if (chgbases)
2868 {
2869 boundary [i]. dom_img. add (chgbases [i]); // ???
2870 if (i > 0)
2871 {
2872 boundary [i]. img_img. add // ???
2873 (chgbases [i - 1]);
2874 }
2875 chgbases_initialized [i] = 0;
2876 }
2877 }
2878
2879 return;
2880} /* chaincomplex<euclidom>::chaincomplex */
2881
2882template <class euclidom>
2883inline chaincomplex<euclidom>::~chaincomplex ()
2884{
2885 if (boundary)
2886 delete [] boundary;
2887 if (generators)
2888 delete [] generators;
2889 if (chgbases)
2890 delete [] chgbases;
2891 if (numgen)
2892 delete [] numgen;
2893 if (generators_initialized)
2894 delete [] generators_initialized;
2895 if (chgbases_initialized)
2896 delete [] chgbases_initialized;
2897 return;
2898} /* chaincomplex<euclidom>::~chaincomplex */
2899
2900template <class euclidom>
2901inline void chaincomplex<euclidom>::def_gen (int q, int_t n)
2902{
2903 if ((q < 0) || (q >= len))
2904 return;
2905
2906 if (numgen [q] < n)
2907 numgen [q] = n;
2908
2909 if (q == 0)
2910 boundary [0]. define (0, numgen [q]);
2911 if ((q > 0) && (numgen [q - 1] >= 0))
2912 boundary [q]. define (numgen [q - 1], numgen [q]);
2913 if ((q < dim ()) && (numgen [q + 1] >= 0))
2914 boundary [q + 1]. define (numgen [q], numgen [q + 1]);
2915
2916 return;
2917} /* chaincomplex<euclidom>::def_gen */
2918
2919template <class euclidom>
2920inline void chaincomplex<euclidom>::add (int q, int_t m, int_t n,
2921 const euclidom &e)
2922{
2923 if ((q <= 0) || (q >= len))
2924 throw "Trying to add a boundary for dimension out of range";
2925
2926 if (numgen [q] <= n)
2927 numgen [q] = n + 1;
2928 if (numgen [q - 1] <= m)
2929 numgen [q - 1] = m + 1;
2930
2931 boundary [q]. add (m, n, e);
2932 return;
2933} /* chaincomplex<euclidom>::add */
2934
2935template <class euclidom>
2936inline euclidom chaincomplex<euclidom>::get (int q, int_t row, int_t col) const
2937{
2938 if ((q < 0) || (q >= len))
2939 throw "Boundary coefficient out of range from chain cplx.";
2940
2941 return boundary [q]. get (row, col);
2942} /* chaincomplex<euclidom>::get */
2943
2944template <class euclidom>
2945inline const mmatrix<euclidom> &chaincomplex<euclidom>::getboundary (int i)
2946 const
2947{
2948 if ((i < 0) || (i >= len))
2949 throw "Boundary matrix out of range from chain complex.";
2950
2951 return boundary [i];
2952} /* chaincomplex<euclidom>::getboundary */
2953
2954template <class euclidom>
2955inline int_t chaincomplex<euclidom>::getnumgen (int i) const
2956{
2957 if ((i < 0) || (i >= len))
2958 // throw "Level for the number of generators out of range.";
2959 return 0;
2960
2961 if (numgen [i] < 0)
2962 return 0;
2963 else
2964 return numgen [i];
2965} /* chaincomplex<euclidom>::getnumgen */
2966
2967template <class euclidom>
2968inline int chaincomplex<euclidom>::dim () const
2969{
2970 return len - 1;
2971} /* chaincomplex<euclidom>::dim */
2972
2973template <class euclidom>
2974inline const chain<euclidom> &chaincomplex<euclidom>::gethomgen (int q,
2975 int_t i) const
2976{
2977 if ((q < 0) || (q >= len))
2978 throw "Level for homology generators out of range.";
2979 if (!generators || !generators_initialized [q])
2980 throw "Trying to get non-existent homology generators.";
2981 return generators [q]. getcol (i);
2982} /* chaincomplex<euclidom>::gethomgen */
2983
2984template <class euclidom>
2985inline const mmatrix<euclidom> &chaincomplex<euclidom>::gethomgen (int q)
2986 const
2987{
2988 if ((q < 0) || (q >= len))
2989 throw "Level for homology generators out of range.";
2990 if (!generators || !generators_initialized [q])
2991 throw "Trying to get non-existent homology generators.";
2992 return generators [q];
2993} /* chaincomplex<euclidom>::gethomgen */
2994
2995template <class euclidom>
2996inline const mmatrix<euclidom> &chaincomplex<euclidom>::getchgbasis (int q)
2997 const
2998{
2999 if ((q < 0) || (q >= len))
3000 throw "Level for basis change matrix out of range.";
3001 if (!chgbases || !chgbases_initialized [q])
3002 throw "Trying to get non-existent basis change matrix.";
3003 return chgbases [q];
3004} /* chaincomplex<euclidom>::getchgbasis */
3005
3006template <class euclidom>
3007inline void chaincomplex<euclidom>::simple_form (int which, bool quiet)
3008{
3009// if ((which < 0) || (which >= len))
3010// throw "Trying to find the simple form of a wrong matrix.";
3011
3012 // if the generator tracing matrices are not initialized, do it now
3013 if (generators)
3014 {
3015 if (!generators_initialized [which])
3016 generators [which]. identity (numgen [which]);
3017 generators_initialized [which] = 1;
3018 if ((which > 0) && (!generators_initialized [which - 1]))
3019 {
3020 generators [which - 1]. identity
3021 (numgen [which - 1]);
3022 generators_initialized [which - 1] = 1;
3023 }
3024 }
3025
3026 // if the base change tracing matrices are not initialized, do it now
3027 if (chgbases)
3028 {
3029 if (!chgbases_initialized [which])
3030 chgbases [which]. identity (numgen [which]);
3031 chgbases_initialized [which] = 1;
3032 if ((which > 0) && (!chgbases_initialized [which - 1]))
3033 {
3034 chgbases [which - 1]. identity
3035 (numgen [which - 1]);
3036 chgbases_initialized [which - 1] = 1;
3037 }
3038 }
3039
3040 // verify that the adjacent matrices have sufficient size
3041 if (which > 0)
3042 {
3043 int_t n = boundary [which]. getnrows ();
3044 mmatrix<euclidom> &other = boundary [which - 1];
3045 if (other. getncols () < n)
3046 other. define (other. getnrows (), n);
3047 }
3048 if (which < len - 1)
3049 {
3050 int_t n = boundary [which]. getncols ();
3051 mmatrix<euclidom> &other = boundary [which + 1];
3052 if (other. getnrows () < n)
3053 other. define (n, other. getncols ());
3054 }
3055
3056 // compute simple form of the desired boundary matrix
3057 if (!quiet && which)
3058 sout << "Reducing D_" << which << ": ";
3059 boundary [which]. simple_form (quiet);
3060 if (!quiet && which)
3061 sout << '\n';
3062
3063 return;
3064} /* chaincomplex<euclidom>::simple_form */
3065
3066template <class euclidom>
3067inline void chaincomplex<euclidom>::simple_form (const int *level,
3068 bool quiet)
3069{
3070 for (int i = len - 1; i >= 0; -- i)
3071 {
3072 if (!level || level [i] || (i && level [i - 1]))
3073 simple_form (i, quiet);
3074 }
3075 return;
3076} /* chaincomplex<euclidom>::simple_form */
3077
3078template <class euclidom>
3079inline int chaincomplex<euclidom>::simple_homology (chain<euclidom> &result,
3080 int which) const
3081{
3082 int_t g = boundary [which]. findcol (-1, 0);
3083 while (g >= 0)
3084 {
3085 euclidom e;
3086 if (which == dim ())
3087 e = 0;
3088 else
3089 e = boundary [which + 1]. get (g, -1);
3090 if (e == 0)
3091 {
3092 e = 1;
3093 result. add (g, e);
3094 }
3095 else if (e. delta () > 1)
3096 result. add (g, e. normalized ());
3097 g = boundary [which]. findcol (-1, g + 1);
3098 }
3099
3100 return which;
3101} /* chaincomplex<euclidom>::simple_homology */
3102
3103template <class euclidom>
3104inline int chaincomplex<euclidom>::simple_homology (chain<euclidom> *&result,
3105 const int *level) const
3106{
3107 // if the chain complex is empty, then just set the result to NULL
3108 if (!len)
3109 {
3110 result = NULL;
3111 return dim ();
3112 }
3113
3114 result = new chain<euclidom> [len];
3115 if (!result)
3116 throw "Not enough memory to create homology chains.";
3117
3118 for (int q = 0; q < len; ++ q)
3119 {
3120 if (!level || level [q])
3121 simple_homology (result [q], q);
3122 }
3123
3124 return dim ();
3125} /* chaincomplex<euclidom>::simple_homology */
3126
3127template <class euclidom>
3128inline void chaincomplex<euclidom>::take_homology
3129 (const chain<euclidom> *hom_chain)
3130{
3131 if (!hom_chain)
3132 return;
3133 for (int q = 0; q < len; ++ q)
3134 def_gen (q, hom_chain [q]. size ());
3135 return;
3136} /* chaincomplex<euclidom>::take_homology */
3137
3138template <class euclidom>
3139inline outputstream &chaincomplex<euclidom>::show_homology
3140 (outputstream &out, const chain<euclidom> *hom, const int *level)
3141 const
3142{
3143 int max_level = len - 1;
3144 while ((max_level >= 0) && !hom [max_level]. size ())
3145 -- max_level;
3146 ++ max_level;
3147 for (int q = 0; q < max_level; ++ q)
3148 {
3149 if (!level || level [q])
3150 {
3151 out << "H_" << q << " = ";
3152 chomp::homology::show_homology (out, hom [q]);
3153 out << '\n';
3154 }
3155 }
3156 return out;
3157} /* chaincomplex<euclidom>::show_homology */
3158
3159template <class euclidom>
3160inline std::ostream &chaincomplex<euclidom>::show_homology (std::ostream &out,
3161 const chain<euclidom> *hom, const int *level) const
3162{
3163 outputstream tout (out);
3164 show_homology (tout, hom, level);
3165 return out;
3166} /* chaincomplex<euclidom>::show_homology */
3167
3168template <class euclidom>
3169inline outputstream &chaincomplex<euclidom>::show_generators
3170 (outputstream &out, const chain<euclidom> &list, int q) const
3171{
3172 if (!generators || (q < 0) || (q >= len))
3173 return out;
3174 for (int_t i = 0; i < list. size (); ++ i)
3175 out << gethomgen (q, list. num (i)) << '\n';
3176 return out;
3177} /* chaincomplex<euclidom>::show_generators */
3178
3179template <class euclidom>
3180inline std::ostream &chaincomplex<euclidom>::show_generators
3181 (std::ostream &out, const chain<euclidom> &list, int q) const
3182{
3183 outputstream tout (out);
3184 show_generators (tout, list, q);
3185 return out;
3186} /* chaincomplex<euclidom>::show_generators */
3187
3188template <class euclidom>
3189inline outputstream &chaincomplex<euclidom>::compute_and_show_homology
3190 (outputstream &out, chain<euclidom> *&hom, const int *level)
3191{
3192 simple_form (level, false);
3193 simple_homology (hom, level);
3194 show_homology (out, hom, level);
3195 return out;
3196} /* chaincomplex<euclidom>::compute_and_show_homology */
3197
3198template <class euclidom>
3199inline std::ostream &chaincomplex<euclidom>::compute_and_show_homology
3200 (std::ostream &out, chain<euclidom> *&hom, const int *level)
3201{
3202 outputstream tout (out);
3203 compute_and_show_homology (tout, hom, level);
3204 return out;
3205} /* chaincomplex<euclidom>::compute_and_show_homology */
3206
3207// --------------------------------------------------
3208
3210template <class euclidom>
3211inline std::ostream &operator << (std::ostream &out,
3212 const chaincomplex<euclidom> &c)
3213{
3214 out << "chain complex" << '\n';
3215 out << "max dimension " << c. dim () << '\n';
3216 out << "dimension 0: " << c. getnumgen (0) << '\n';
3217 for (int i = 1; i <= c. dim (); ++ i)
3218 {
3219 out << "dimension " << i << ": " << c. getnumgen (i) << '\n';
3220 for (int_t j = 0; j < c. getnumgen (i); ++ j)
3221 {
3222 if (c. getboundary (i). getcol (j). size ())
3223 {
3224 out << "\t# " << (j + 1) << " = " <<
3225 c. getboundary (i). getcol (j) <<
3226 '\n';
3227 }
3228 }
3229 }
3230 return out;
3231} /* operator << */
3232
3233
3234// --------------------------------------------------
3235// -------------------- chainmap --------------------
3236// --------------------------------------------------
3237
3241template <class euclidom>
3242class chainmap
3243{
3244public:
3247 chainmap (const chaincomplex<euclidom> &domain,
3248 const chaincomplex<euclidom> &range);
3249
3251 chainmap (const chainmap<euclidom> &c);
3252
3254 chainmap<euclidom> &operator = (const chainmap<euclidom> &c);
3255
3257 ~chainmap ();
3258
3260 int dim () const;
3261
3263 const mmatrix<euclidom> &operator [] (int i) const;
3264
3267 void add (int q, int_t m, int_t n, euclidom e);
3268
3270 void invert (void);
3271
3274 void compose (const chainmap<euclidom> &m1,
3275 const chainmap<euclidom> &m2);
3276
3280 outputstream &show (outputstream &out,
3281 const char *maplabel = "f", const char *xtxt = NULL,
3282 const char *ytxt = NULL, const int *level = NULL) const;
3283
3287 std::ostream &show (std::ostream &out,
3288 const char *maplabel = "f", const char *xtxt = NULL,
3289 const char *ytxt = NULL, const int *level = NULL) const;
3290
3294 void take_homology (const chainmap<euclidom> &m,
3295 const chain<euclidom> *hom_domain,
3296 const chain<euclidom> *hom_range);
3297
3301 outputstream &show_homology (outputstream &out,
3302 const chain<euclidom> *hom_domain,
3303 const chain<euclidom> *hom_range, const int *level = NULL,
3304 const char *xtxt = NULL, const char *ytxt = NULL) const;
3305
3309 std::ostream &show_homology (std::ostream &out,
3310 const chain<euclidom> *hom_domain,
3311 const chain<euclidom> *hom_range,
3312 const int *level = NULL,
3313 const char *xtxt = NULL, const char *ytxt = NULL)
3314 const;
3315
3316private:
3318 int len;
3319
3321 mmatrix<euclidom> *map;
3322
3323}; /* class chainmap */
3324
3325// --------------------------------------------------
3326
3327template <class euclidom>
3328inline chainmap<euclidom>::chainmap (const chaincomplex<euclidom> &domain,
3329 const chaincomplex<euclidom> &range)
3330{
3331 // set the dimension
3332 len = domain. len;
3333 if (range. len < domain. len)
3334 len = range. len;
3335
3336 // allocate new matrices
3337 if (len)
3338 map = new mmatrix<euclidom> [len];
3339 else
3340 map = NULL;
3341
3342 for (int i = 0; i < len; ++ i)
3343 {
3344 // define the size of the matrix (number of rows and columns)
3345 map [i]. define (range. getnumgen (i),
3346 domain. getnumgen (i));
3347
3348 // link the matrices to the ones in the chain complexes
3349 domain. boundary [i]. dom_dom. add (map [i]);
3350 range. boundary [i]. dom_img. add (map [i]);
3351 if (i < domain. len - 1)
3352 domain. boundary [i + 1]. img_dom. add (map [i]);
3353 if (i < range. len - 1)
3354 range. boundary [i + 1]. img_img. add (map [i]);
3355 }
3356
3357 return;
3358} /* chainmap<euclidom>::chainmap */
3359
3360template <class euclidom>
3361inline chainmap<euclidom>::chainmap (const chainmap<euclidom> &c)
3362{
3363 len = c. len;
3364 if (len)
3365 map = new mmatrix<euclidom> [len];
3366 else
3367 map = 0;
3368
3369 for (int i = 0; i < len; ++ i)
3370 map [i] = c. map [i];
3371
3372 return;
3373} /* chainmap<euclidom>::chainmap */
3374
3375template <class euclidom>
3376inline chainmap<euclidom> &chainmap<euclidom>::operator =
3377 (const chainmap<euclidom> &c)
3378{
3379 if (map)
3380 delete [] map;
3381
3382 len = c. len;
3383 if (len)
3384 map = new mmatrix<euclidom> [len];
3385 else
3386 map = 0;
3387
3388 for (int i = 0; i < len; ++ i)
3389 map [i] = c. map [i];
3390
3391 return *this;
3392} /* chainmap<euclidom>::operator = */
3393
3394template <class euclidom>
3395inline chainmap<euclidom>::~chainmap ()
3396{
3397 if (map)
3398 delete [] map;
3399 return;
3400} /* chainmap<euclidom>::~chainmap */
3401
3402template <class euclidom>
3403inline int chainmap<euclidom>::dim () const
3404{
3405 return len - 1;
3406} /* chainmap<euclidom>::dim */
3407
3408template <class euclidom>
3409inline const mmatrix<euclidom> &chainmap<euclidom>::operator [] (int i) const
3410{
3411// if ((i < 0) || (i >= len))
3412// throw "Chain map level out of range.";
3413 return map [i];
3414} /* chainmap<euclidom>::operator [] */
3415
3416template <class euclidom>
3417inline void chainmap<euclidom>::add (int q, int_t m, int_t n, euclidom e)
3418{
3419 map [q]. add (m, n, e);
3420 return;
3421} /* chainmap<euclidom>::add */
3422
3423template <class euclidom>
3424inline void chainmap<euclidom>::take_homology (const chainmap<euclidom> &m,
3425 const chain<euclidom> *hom_domain, const chain<euclidom> *hom_range)
3426{
3427 if (!hom_domain || !hom_range)
3428 return;
3429
3430 for (int q = 0; q < len; ++ q)
3431 {
3432 int_t dlen = hom_domain [q]. size ();
3433 const chain<euclidom> &r = hom_range [q];
3434 int_t rlen = r. size ();
3435 map [q]. define (rlen, dlen);
3436 // go through the homology generators in the domain
3437 for (int_t i = 0; i < dlen; ++ i)
3438 {
3439 // retrieve the real number of the homology generator
3440 int_t x = hom_domain [q]. num (i);
3441
3442 // get the image of this element by the chain map
3443 const chain<euclidom> &img = m [q]. getcol (x);
3444
3445 // transform numbers in this image to hom. generators
3446 int_t j = 0;
3447 for (int_t k = 0; k < img. size (); ++ k)
3448 {
3449 // find the current element in the range
3450 while ((j < rlen) &&
3451 (r. num (j) < img. num (k)))
3452 ++ j;
3453
3454 // if found in the range, add it
3455 if ((j < rlen) &&
3456 (r. num (j) == img. num (k)))
3457 map [q]. add (j, i, img. coef (k));
3458 }
3459 }
3460 }
3461 return;
3462} /* chainmap<euclidom>::take_homology */
3463
3464template <class euclidom>
3465inline void chainmap<euclidom>::invert (void)
3466{
3467 for (int q = 0; q < len; ++ q)
3468 map [q]. invert ();
3469 return;
3470} /* chainmap<euclidom>::invert */
3471
3472template <class euclidom>
3473inline void chainmap<euclidom>::compose (const chainmap<euclidom> &m1,
3474 const chainmap<euclidom> &m2)
3475{
3476 if ((m1. len < len) || (m2. len < len))
3477 throw "Trying to compose chain maps with too few levels.";
3478 for (int q = 0; q < len; ++ q)
3479 map [q]. multiply (m1. map [q], m2. map [q]);
3480 return;
3481} /* chainmap<euclidom>::compose */
3482
3483template <class euclidom>
3484inline outputstream &chainmap<euclidom>::show (outputstream &out,
3485 const char *maplabel, const char *xtxt, const char *ytxt,
3486 const int *level) const
3487{
3488 for (int q = 0; q < len; ++ q)
3489 {
3490 if (level && !level [q])
3491 continue;
3492 out << "Dim " << q << ":";
3493 map [q]. showmap (out, maplabel, xtxt, ytxt);
3494 }
3495 return out;
3496} /* chainmap<euclidom>::show */
3497
3498template <class euclidom>
3499inline std::ostream &chainmap<euclidom>::show (std::ostream &out,
3500 const char *maplabel, const char *xtxt, const char *ytxt,
3501 const int *level) const
3502{
3503 outputstream tout (out);
3504 show (tout, maplabel, xtxt, ytxt, level);
3505 return out;
3506} /* chainmap<euclidom>::show */
3507
3508template <class euclidom>
3509inline outputstream &chainmap<euclidom>::show_homology (outputstream &out,
3510 const chain<euclidom> *hom_domain,
3511 const chain<euclidom> *hom_range, const int *level,
3512 const char *xtxt, const char *ytxt) const
3513{
3514 int max_len = len - 1;
3515 while ((max_len >= 0) && !hom_domain [max_len]. size ())
3516 -- max_len;
3517 ++ max_len;
3518 for (int q = 0; q < max_len; ++ q)
3519 {
3520 if (!level || level [q])
3521 {
3522 out << "Dim " << q << ":";
3523 int hlen = hom_domain [q]. size ();
3524 if (!hlen)
3525 out << "\t0" << '\n';
3526 for (int i = 0; i < hlen; ++ i)
3527 {
3528 out << "\tf (";
3529 if (xtxt)
3530 out << xtxt;
3531 out << (i + 1) << ") = ";
3532 map [q]. show_hom_col (out,
3533 hom_domain [q]. num (i),
3534 hom_range [q], ytxt);
3535 out << '\n';
3536 }
3537 }
3538 }
3539 return out;
3540} /* chainmap<euclidom>::show_homology */
3541
3542template <class euclidom>
3543inline std::ostream &chainmap<euclidom>::show_homology (std::ostream &out,
3544 const chain<euclidom> *hom_domain,
3545 const chain<euclidom> *hom_range, const int *level,
3546 const char *xtxt, const char *ytxt) const
3547{
3548 outputstream tout (out);
3549 show_homology (tout, hom_domain, hom_range, level, xtxt, ytxt);
3550 return out;
3551} /* chainmap<euclidom>::show_homology */
3552
3553// --------------------------------------------------
3554
3556template <class euclidom>
3557inline std::ostream &operator << (std::ostream &out,
3558 const chainmap<euclidom> &m)
3559{
3560 out << "chain map" << '\n';
3561 for (int i = 0; i <= m. dim (); ++ i)
3562 {
3563 out << "dimension " << i << '\n';
3564 for (int_t j = 0; j < m [i]. getncols (); ++ j)
3565 {
3566 if (m [i]. getcol (j). size ())
3567 {
3568 out << "\t# " << (j + 1) << " = " <<
3569 m [i]. getcol (j) << '\n';
3570 }
3571 }
3572 }
3573 return out;
3574} /* operator << */
3575
3576
3577// --------------------------------------------------
3578// -------------- displaying homology ---------------
3579// --------------------------------------------------
3580
3583template <class euclidom>
3584inline outputstream &show_homology (outputstream &out,
3585 const chain<euclidom> &c)
3586{
3587 int_t countfree = 0;
3588 bool writeplus = false;
3589 for (int_t i = 0; i < c. size (); ++ i)
3590 {
3591 if (c. coef (i). delta () == 1)
3592 ++ countfree;
3593 else
3594 {
3595 out << (writeplus ? " + " : "") <<
3596 euclidom::ringsymbol () << "_" <<
3597 c. coef (i);
3598 writeplus = true;
3599 }
3600
3601 if (countfree && ((i == c. size () - 1) ||
3602 (c. coef (i + 1). delta () != 1)))
3603 {
3604 out << (writeplus ? " + " : "") <<
3605 euclidom::ringsymbol ();
3606 if (countfree > 1)
3607 out << "^" << countfree;
3608 countfree = 0;
3609 writeplus = true;
3610 }
3611 }
3612
3613 // if there was nothing to show, then just show zero
3614 if (!c. size ())
3615 out << "0";
3616
3617 return out;
3618} /* show_homology */
3619
3622template <class euclidom>
3623inline std::ostream &show_homology (std::ostream &out,
3624 const chain<euclidom> &c)
3625{
3626 outputstream tout (out);
3627 show_homology (tout, c);
3628 return out;
3629} /* show_homology */
3630
3631
3632} // namespace homology
3633} // namespace chomp
3634
3635#endif // _CHOMP_HOMOLOGY_CHAINS_H_
3636
3638
autoarray.h
An auto_array template that mimics selected behaviors of the std::auto_ptr template,...
chomp::homology::auto_array
An auto_array template that mimics selected behaviors of the std::auto_ptr template,...
Definition: autoarray.h:54
chomp::homology::chain
This class defines objects which represent chains as finite sequences of elements identified by integ...
Definition: chains.h:94
chomp::homology::chain::multiply
chain< euclidom > & multiply(euclidom e, int_t number=-1)
Multiplies one or all the coefficients in the chain by the given number.
Definition: chains.h:801
chomp::homology::chain::findbest
int_t findbest(chain< euclidom > *table=NULL) const
Finds the best element in the chain for reduction, that is, the element with minimal value of delta.
Definition: chains.h:358
chomp::homology::chain::chain
chain()
The default constructor.
Definition: chains.h:207
chomp::homology::chain::remove
chain< euclidom > & remove(int_t n)
Removes an element with the given identifier from the chain.
Definition: chains.h:583
chomp::homology::chain::show
outputstream & show(outputstream &out, const char *label=NULL) const
Shows the chain to the output stream.
Definition: chains.h:849
chomp::homology::chain::findnumber
int_t findnumber(int_t n) const
Find the position of an element with the given identifier.
Definition: chains.h:318
chomp::homology::chain::~chain
~chain()
The destructor.
Definition: chains.h:270
chomp::homology::chain::len
int_t len
The length of the list and the length of the table.
Definition: chains.h:183
chomp::homology::chain::coef
euclidom coef(int_t i) const
Returns the coefficient in front of the i-th element in the chain.
Definition: chains.h:331
chomp::homology::chain::insertpair
chain< euclidom > & insertpair(int_t i, int_t n, euclidom e)
Inserts one chain element at the given position.
Definition: chains.h:411
chomp::homology::chain::operator=
chain< euclidom > & operator=(const chain< euclidom > &c)
The assignment operator.
Definition: chains.h:236
chomp::homology::chain::_n
int_t * _n
Identifiers of the elements of the chain (sorted).
Definition: chains.h:186
chomp::homology::chain::getcoefficient
euclidom getcoefficient(int_t n=-1) const
Finds and returns the coefficient in front of the given element.
Definition: chains.h:293
chomp::homology::chain::size
int_t size() const
Returns the size of the chain, that is, the number of elements with non-zero coefficients.
Definition: chains.h:281
chomp::homology::chain::take
chain< euclidom > & take(chain< euclidom > &c)
Takes data from another chain. Destroys the other chain.
Definition: chains.h:781
chomp::homology::chain::add
chain< euclidom > & add(int_t n, euclidom e)
Adds an element algebraically to the chain.
Definition: chains.h:551
chomp::homology::chain::_e
euclidom * _e
Coefficients of the elements of the chain.
Definition: chains.h:189
chomp::homology::chain::num
int_t num(int_t i) const
Returns the number (identifier) of the i-th element in the chain.
Definition: chains.h:339
chomp::homology::chain::contains_non_invertible
bool contains_non_invertible() const
Determines if the chain contains a non-invertible coefficient.
Definition: chains.h:347
chomp::homology::chain::swap
chain< euclidom > & swap(chain< euclidom > &other, int_t number=-1, int_t othernumber=-1, chain< euclidom > *table=NULL)
Swaps one chain with another.
Definition: chains.h:738
chomp::homology::chain::swapnumbers
chain< euclidom > & swapnumbers(int_t number1, int_t number2)
Swaps two numbers (identifiers) in the chain.
Definition: chains.h:493
chomp::homology::chain::removepair
chain< euclidom > & removepair(int_t i)
Removes one chain element at the given position.
Definition: chains.h:451
chomp::homology::chain::empty
bool empty() const
Returns true if the chain is empty (zero), false otherwise.
Definition: chains.h:287
chomp::homology::chaincomplex
This is an implementation of a chain complex over an arbitrary ring.
Definition: chains.h:2662
chomp::homology::chaincomplex::get
euclidom get(int q, int_t row, int_t col) const
Returns an element from the boundary matrix for the given dimension.
Definition: chains.h:2936
chomp::homology::chaincomplex::generators
mmatrix< euclidom > * generators
Matrices which store actual combinations of generators.
Definition: chains.h:2798
chomp::homology::chaincomplex::chaincomplex
chaincomplex(int d, int trace_generators=0, int trace_bases=0)
The default constructor.
Definition: chains.h:2820
chomp::homology::chaincomplex::getchgbasis
const mmatrix< euclidom > & getchgbasis(int q) const
Returns the base change matrix at level q.
Definition: chains.h:2996
chomp::homology::chaincomplex::~chaincomplex
~chaincomplex()
The destructor.
Definition: chains.h:2883
chomp::homology::chaincomplex::add
void add(int q, int_t m, int_t n, const euclidom &e)
Adds a coefficient to the structure: D_q [m, n] += e.
Definition: chains.h:2920
chomp::homology::chaincomplex::show_homology
outputstream & show_homology(outputstream &out, const chain< euclidom > *hom, const int *level=NULL) const
Writes the homology module of the chain complex to an output stream.
Definition: chains.h:3140
chomp::homology::chaincomplex::compute_and_show_homology
outputstream & compute_and_show_homology(outputstream &out, chain< euclidom > *&hom, const int *level=NULL)
Computes the homology and shows the result.
Definition: chains.h:3190
chomp::homology::chaincomplex::chaincomplex
chaincomplex(const chaincomplex< euclidom > &m)
The copy constructor is not implemented.
Definition: chains.h:2779
chomp::homology::chaincomplex::dim
int dim() const
Returns the dimension of the chain complex.
Definition: chains.h:2968
chomp::homology::chaincomplex::simple_form
void simple_form(int which, bool quiet)
Reduces of the given boundary matrix in the chain complex for the purpose of homology computation.
Definition: chains.h:3007
chomp::homology::chaincomplex::getnumgen
int_t getnumgen(int i) const
Returns the number of generators at the given level.
Definition: chains.h:2955
chomp::homology::chaincomplex::chgbases
mmatrix< euclidom > * chgbases
Matrices which store the change of basis to obtain the SNF.
Definition: chains.h:2801
chomp::homology::chaincomplex::gethomgen
const chain< euclidom > & gethomgen(int q, int_t i) const
Returns the given homology generator at level q.
Definition: chains.h:2974
chomp::homology::chaincomplex::getboundary
const mmatrix< euclidom > & getboundary(int i) const
Returns a reference to the given boundary matrix.
Definition: chains.h:2945
chomp::homology::chaincomplex::numgen
int_t * numgen
The numbers of generators in each dimension, or -1's if not defined yet.
Definition: chains.h:2813
chomp::homology::chaincomplex::take_homology
void take_homology(const chain< euclidom > *hom_chain)
Creates a chain complex containing exactly one generator for each homology generator.
Definition: chains.h:3129
chomp::homology::chaincomplex::simple_homology
int simple_homology(chain< euclidom > &result, int which) const
Computes the given level of homology of the chain complex, provided it has been transformed into the ...
Definition: chains.h:3079
chomp::homology::chaincomplex::show_generators
outputstream & show_generators(outputstream &out, const chain< euclidom > &list, int q) const
Writes the homology generators of the homology module to an output stream.
Definition: chains.h:3170
chomp::homology::chaincomplex::boundary
mmatrix< euclidom > * boundary
The matrices of the boundary homomorphism.
Definition: chains.h:2794
chomp::homology::chaincomplex::chgbases_initialized
int * chgbases_initialized
Have the base change tracing matrices been initialized to the identity (of suitable size each)?
Definition: chains.h:2809
chomp::homology::chaincomplex::generators_initialized
int * generators_initialized
Have the generator tracing matrices been initialized to the identity (of suitable size each)?
Definition: chains.h:2805
chomp::homology::chaincomplex::len
int len
The length of the chain complex (i.e., its dimension + 1).
Definition: chains.h:2790
chomp::homology::chaincomplex::def_gen
void def_gen(int q, int_t n)
Defines the number of generators in the given dimension.
Definition: chains.h:2901
chomp::homology::chainmap
This class defines a chain map between two chain complexes.
Definition: chains.h:3243
chomp::homology::chainmap::dim
int dim() const
Returns the dimension of the chain map.
Definition: chains.h:3403
chomp::homology::chainmap::len
int len
The number of matrices (dimension of the chain map + 1).
Definition: chains.h:3318
chomp::homology::chainmap::add
void add(int q, int_t m, int_t n, euclidom e)
Adds a coefficient to the selected matrix of the map: M_q [m, n] += e.
Definition: chains.h:3417
chomp::homology::chainmap::operator=
chainmap< euclidom > & operator=(const chainmap< euclidom > &c)
The assignment operator.
Definition: chains.h:3377
chomp::homology::chainmap::operator[]
const mmatrix< euclidom > & operator[](int i) const
Returns the matrix of the chain map at the given level.
Definition: chains.h:3409
chomp::homology::chainmap::invert
void invert(void)
Inverts the chain map.
Definition: chains.h:3465
chomp::homology::chainmap::show
outputstream & show(outputstream &out, const char *maplabel="f", const char *xtxt=NULL, const char *ytxt=NULL, const int *level=NULL) const
Writes the chain map to an output stream in the text format using specified labels for the map and el...
Definition: chains.h:3484
chomp::homology::chainmap::map
mmatrix< euclidom > * map
The matrices in each dimension.
Definition: chains.h:3321
chomp::homology::chainmap::chainmap
chainmap(const chaincomplex< euclidom > &domain, const chaincomplex< euclidom > &range)
The default constructor of a chain map between the two given chain complexes.
Definition: chains.h:3328
chomp::homology::chainmap::show_homology
outputstream & show_homology(outputstream &out, const chain< euclidom > *hom_domain, const chain< euclidom > *hom_range, const int *level=NULL, const char *xtxt=NULL, const char *ytxt=NULL) const
Writes to an output stream the map induced in homology.
Definition: chains.h:3509
chomp::homology::chainmap::take_homology
void take_homology(const chainmap< euclidom > &m, const chain< euclidom > *hom_domain, const chain< euclidom > *hom_range)
Creates a chain map that represents the map induced in homology by the chain map between the two give...
Definition: chains.h:3424
chomp::homology::chainmap::compose
void compose(const chainmap< euclidom > &m1, const chainmap< euclidom > &m2)
Composes two given chain maps.
Definition: chains.h:3473
chomp::homology::chainmap::~chainmap
~chainmap()
The destructor.
Definition: chains.h:3395
chomp::homology::mmatrix
A class for representing sparse matrices containing elements of the 'euclidom' type.
Definition: chains.h:1067
chomp::homology::mmatrix::img_img
simplelist< mmatrix< euclidom > > img_img
Definition: chains.h:1206
chomp::homology::mmatrix::increaserows
void increaserows(int_t numrows)
Increases tables to be enough to keep the given number of rows.
Definition: chains.h:2606
chomp::homology::mmatrix::getcol
const chain< euclidom > & getcol(int_t n) const
Returns a reference to the entire column stored as a chain.
Definition: chains.h:1525
chomp::homology::mmatrix::show_hom_col
outputstream & show_hom_col(outputstream &out, int_t col, const chain< euclidom > &range, const char *txt=NULL) const
Writes a column with only these coefficients which exist in 'range', and shows their positions in tha...
Definition: chains.h:2542
chomp::homology::mmatrix::mmatrix
mmatrix()
The default constructor.
Definition: chains.h:1347
chomp::homology::mmatrix::define
void define(int_t numrows, int_t numcols)
Defines the number of rows and columns and increases the internal tables if necessary.
Definition: chains.h:1364
chomp::homology::mmatrix::img_dom
simplelist< mmatrix< euclidom > > img_dom
Definition: chains.h:1205
chomp::homology::mmatrix::getnrows
int_t getnrows() const
Returns the number of rows in the matrix.
Definition: chains.h:1533
chomp::homology::mmatrix::get
euclidom get(int_t row, int_t col) const
Returns the value at the desired location of the matrix.
Definition: chains.h:1495
chomp::homology::mmatrix::allrows
int_t allrows
The number of allocated rows.
Definition: chains.h:1268
chomp::homology::mmatrix::arrange_towards_SNF
int_t arrange_towards_SNF(int_t *invertible_count=0)
Swaps rows and columns to make the simple form closer to SNF.
Definition: chains.h:2062
chomp::homology::mmatrix::addrow
void addrow(int_t dest, int_t source, const euclidom &e)
Adds one row to another with a given coefficient.
Definition: chains.h:1545
chomp::homology::mmatrix::nrows
int_t nrows
The number of rows in the matrix.
Definition: chains.h:1262
chomp::homology::mmatrix::dom_img
simplelist< mmatrix< euclidom > > dom_img
Definition: chains.h:1205
chomp::homology::mmatrix::addcol
void addcol(int_t dest, int_t source, const euclidom &e)
Adds one column to another with a given coefficient.
Definition: chains.h:1572
chomp::homology::mmatrix::showrows
outputstream & showrows(outputstream &out, int_t first=0, int_t howmany=0, const char *label="Row ") const
Writes the matrix to an output stream by its rows.
Definition: chains.h:1885
chomp::homology::mmatrix::~mmatrix
~mmatrix()
The destructor of a matrix.
Definition: chains.h:1354
chomp::homology::mmatrix::invert
void invert(void)
Inverts the matrix. Throws an error message on failure.
Definition: chains.h:2428
chomp::homology::mmatrix::reducerow
int_t reducerow(int_t n, int_t preferred)
Reduces the given row of the matrix and updates its columns.
Definition: chains.h:1765
chomp::homology::mmatrix::simple_form
void simple_form(bool quiet=false)
Transforms the matrix to a simple form (nearly SNF).
Definition: chains.h:2003
chomp::homology::mmatrix::swapcols
void swapcols(int_t i, int_t j)
Swaps two columns of the matrix.
Definition: chains.h:1626
chomp::homology::mmatrix::increasecols
void increasecols(int_t numcols)
Increases tables to be enough to keep the given number of columns.
Definition: chains.h:2623
chomp::homology::mmatrix::rows
chain< euclidom > * rows
The rows of the matrix.
Definition: chains.h:1274
chomp::homology::mmatrix::allcols
int_t allcols
The number of allocated columns.
Definition: chains.h:1271
chomp::homology::mmatrix::mult_right
void mult_right(const int_t setRows[], const int_t setCols[], const mmatrix< euclidom > &M, const mmatrix< euclidom > &invM, bool update_linked)
Multiplies the matrix from the right by an invertible square matrix whose possibly nonzero entries ar...
Definition: chains.h:2303
chomp::homology::mmatrix::reducecol
int_t reducecol(int_t n, int_t preferred)
Reduces the given column of the matrix and updates its rows.
Definition: chains.h:1817
chomp::homology::mmatrix::division_SNF_correction
void division_SNF_correction(const euclidom &a, int_t pos1, const euclidom &b, int_t pos2)
Makes a correction of division in the "diagonal" at the given two positions, so that the entry at pos...
Definition: chains.h:2145
chomp::homology::mmatrix::ncols
int_t ncols
The number of columns in the matrix.
Definition: chains.h:1265
chomp::homology::mmatrix::dom_dom
simplelist< mmatrix< euclidom > > dom_dom
This is a list of matrices to be updated together with the changes to the columns or rows of the curr...
Definition: chains.h:1205
chomp::homology::mmatrix::swaprows
void swaprows(int_t i, int_t j)
Swaps two rows of the matrix.
Definition: chains.h:1599
chomp::homology::mmatrix::multiplycol
void multiplycol(int_t n, const euclidom &e)
Multiplies the column by the given coefficient and updates rows.
Definition: chains.h:1669
chomp::homology::mmatrix::showmap
outputstream & showmap(outputstream &out, const char *maplabel=NULL, const char *xlabel=NULL, const char *ylabel=NULL) const
Writes the matrix as a linear map on generators.
Definition: chains.h:1917
chomp::homology::mmatrix::showrowscols
outputstream & showrowscols(outputstream &out, chain< euclidom > *table, int_t tablen, int_t first=0, int_t howmany=0, const char *label=NULL) const
Writes the matrix to an output stream by its rows or columns.
Definition: chains.h:1869
chomp::homology::mmatrix::extendedGCD
static void extendedGCD(const euclidom &a, const euclidom &b, euclidom &x, euclidom &y)
Extended Euclidean algorithm for the computation of the greatest common divisor, together with the co...
Definition: chains.h:2401
chomp::homology::mmatrix::increase
void increase(int_t numrows, int_t numcols)
Increases tables to be enough to keep the given number of columns and rows.
Definition: chains.h:2598
chomp::homology::mmatrix::add
void add(int_t row, int_t col, const euclidom &e)
Adds a value to the given element of the matrix.
Definition: chains.h:1468
chomp::homology::mmatrix::multiplyrow
void multiplyrow(int_t n, const euclidom &e)
Multiplies the row by the given coefficient and updates columns.
Definition: chains.h:1653
chomp::homology::mmatrix::getrow
const chain< euclidom > & getrow(int_t n) const
Returns a reference to the entire row stored as a chain.
Definition: chains.h:1517
chomp::homology::mmatrix::getncols
int_t getncols() const
Returns the number of columns in the matrix.
Definition: chains.h:1539
chomp::homology::mmatrix::cols
chain< euclidom > * cols
The columns of the matrix.
Definition: chains.h:1277
chomp::homology::mmatrix::showcols
outputstream & showcols(outputstream &out, int_t first=0, int_t howmany=0, const char *label="Col ") const
Writes the matrix to an output stream by its rows.
Definition: chains.h:1901
chomp::homology::mmatrix::identity
void identity(int_t size)
Defines the matrix to be the identity of the given size.
Definition: chains.h:1452
chomp::homology::mmatrix::multiply
void multiply(const mmatrix< euclidom > &m1, const mmatrix< euclidom > &m2)
Computes the product of the two given matrices.
Definition: chains.h:2502
chomp::homology::mmatrix::simple_reductions
void simple_reductions(bool quiet=false)
Runs some random simple reductions.
Definition: chains.h:1946
chomp::homology::mmatrix::mult_left
void mult_left(const int_t setRows[], const int_t setCols[], const mmatrix< euclidom > &M, const mmatrix< euclidom > &invM, bool update_linked)
Multiplies the matrix from the left by an invertible square matrix whose possibly nonzero entries are...
Definition: chains.h:2205
chomp::homology::mmatrix::findrowcol
int_t findrowcol(int_t req_elements, int_t start, int which) const
An internal procedure for both findrow and findcol.
Definition: chains.h:1685
chomp::homology::mmatrix::findrow
int_t findrow(int_t req_elements=1, int_t start=-1) const
Finds a row containing at least the required number of nonzero elements, starting at the given row.
Definition: chains.h:1753
chomp::homology::mmatrix::operator=
mmatrix< euclidom > & operator=(const mmatrix< euclidom > &s)
The assignment operator.
Definition: chains.h:1413
chomp::homology::mmatrix::findcol
int_t findcol(int_t req_elements=1, int_t start=-1) const
Finds a column containing at least the required number of nonzero elements, starting at the given col...
Definition: chains.h:1759
chomp::homology::mmatrix::submatrix
void submatrix(const mmatrix< euclidom > &matr, const chain< euclidom > &domain, const chain< euclidom > &range)
Extracts a submatrix from the given matrix.
Definition: chains.h:2524
chomp::homology::mmatrix::simple_form_to_SNF
void simple_form_to_SNF(bool quiet=false)
Transforms the matrix from the simple form to actual SNF.
Definition: chains.h:2103
chomp::homology::outputstream
This class defines an output stream for replacing the standard 'cout'.
Definition: textfile.h:64
chomp::homology::simplelist
This class defines a simple list of pointers to objects of the given type.
Definition: chains.h:951
chomp::homology::simplelist::simplelist
simplelist()
The default constructor of an empty list.
Definition: chains.h:1003
chomp::homology::simplelist::operator=
simplelist< element > & operator=(const simplelist< element > &s)
The assignment operator is not implemented.
Definition: chains.h:982
chomp::homology::simplelist::take
element * take()
A simple internal iterator of the list.
Definition: chains.h:1044
chomp::homology::simplelist::remove
void remove(element &m)
Remove an element from the list.
Definition: chains.h:1029
chomp::homology::simplelist::add
void add(element &m)
Adds an element to the list.
Definition: chains.h:1017
chomp::homology::simplelist::simplelist
simplelist(const simplelist< element > &s)
The copy constructor is not implemented.
Definition: chains.h:974
chomp::homology::simplelist::cur
int cur
The current element in the table.
Definition: chains.h:993
chomp::homology::simplelist::elem
element ** elem
A table of element pointers.
Definition: chains.h:996
chomp::homology::simplelist::num
int num
The number of element pointers stored in the list.
Definition: chains.h:990
chomp::homology::simplelist::~simplelist
~simplelist()
The destructor.
Definition: chains.h:1009
config.h
This file contains some precompiler definitions which indicate the operating system and/or compiler u...
int_t
int int_t
Index type for indexing arrays, counting cubes, etc.
Definition: config.h:115
chomp::homengin::ringsymbol
const char * ringsymbol()
Definition: algstruct.h:216
chomp::homology::show_homology
outputstream & show_homology(outputstream &out, const chain< euclidom > &c)
Shows a chain as a list of generators of one level of a homology module.
Definition: chains.h:3584
chomp::homology::sout
outputstream sout
A replacement for standard output stream, with optional logging and other features provided by the cl...
chomp::homology::swapelements
void swapelements(type &x, type &y)
A simple template for swapping two elements with the use of a temporary variable of the same type and...
Definition: textfile.h:329
chomp::homology::operator<<
std::ostream & operator<<(std::ostream &out, const bincube< Dim, twoPower > &b)
Definition: bincube.h:907
chomp::homology::readparenthesis
int readparenthesis(std::istream &in)
Reads an opening parenthesis from the input file.
Definition: textfile.h:432
chomp::homology::operator>>
std::istream & operator>>(std::istream &in, bincube< Dim, twoPower > &b)
Definition: bincube.h:914
chomp::homology::scon
outputstream scon
The console output stream to which one should put all the junk that spoils the log file,...
chomp::homology::ignorecomments
void ignorecomments(std::istream &in)
Ignores white characters (spaces, tabulators, CRs and LFs), as well as comments from the input text f...
Definition: textfile.h:385
chomp::multiwork::swap
void swap(mwWorkerData &data1, mwWorkerData &data2)
Definition: mwcoord.h:108
chomp
This namespace contains the entire CHomP library interface.
Definition: bitmaps.h:51
textfile.h
This file contains some useful functions related to the text input/output procedures.
Generated by doxygen 1.9.4