1 // Copyright (c) 2013-2015 Mario Albert
2 // This file is part of the CoCoALib suite of examples.
3 // You are free to use any part of this example in your own programs.
4
5 #include "CoCoA/library.H"
6 #include <bitset>
7 using namespace std;
8 using namespace CoCoA;
9 using namespace CoCoA::Involutive;
10
11 //----------------------------------------------------------------------
12 const string ShortDescription =
13 "In this file we explain how to using some methods which are using the Pommaret Basis. \n";
14
15 const string LongDescription =
16 "In this file we explain how to using some methods which are using the Pommaret Basis. \n";
17 //----------------------------------------------------------------------
18
output(vector<RingElem> vec)19 void output(vector<RingElem> vec)
20 {
21 for(vector<RingElem>::iterator iter = vec.begin(); iter != vec.end(); ++iter)
22 {
23 cout << *iter << endl;
24 }
25 }
26
27
program()28 void program()
29 {
30 GlobalManager CoCoAFoundations;
31
32 cout << ShortDescription << endl;
33 cout << boolalpha << endl;
34 ring Q = RingQQ();
35 cout <<"///////////////A homogenous ideal with a Pommaret basis//////////////////////////////////////////" << endl;
36 SparsePolyRing polyRing6 = NewPolyRing(Q, SymbolRange("x",0,5));
37 vector<RingElem> x = indets(polyRing6);
38 vector<RingElem> polys;
39
40 polys.push_back( power(x[0],2) + power(x[1],2) + power(x[2],2) );
41 polys.push_back( x[0]*x[3] + x[1]*x[4] + x[2]*x[5]);
42 polys.push_back( 50*power(x[0],2) - 2*x[0]*x[4] + 14*x[0]*x[5] + power(x[1],2) - 14*x[1]*x[2] + 2*x[1]*x[3] + 49*power(x[2],2) - 14*x[2]*x[3] + power(x[3],2) + power(x[4],2) + power(x[5],2) );
43 polys.push_back( 29*power(x[0],2) - 10*x[0]*x[1] - 4*x[0]*x[2] - 4*x[0]*x[4] + 10*x[0]*x[5] + 5*power(x[1],2) - 20*x[1]*x[2] + 4*x[1]*x[3] - 2*x[1]*x[5] + 26*power(x[2],2) - 10*x[2]*x[3] + 2*x[2]*x[4] + power(x[3],2) + power(x[4],2) + power(x[5],2) );
44 polys.push_back( 9*power(x[0],2) - 6*x[0]*x[5] + 9*power(x[2],2) + 6*x[2]*x[3] + power(x[3],2) + power(x[4],2) + power(x[5],2) );
45 polys.push_back( 9*power(x[0],2) + 6*x[0]*x[5] + 9*power(x[2],2) - 6*x[2]*x[3] + power(x[3],2) + power(x[4],2) + power(x[5],2) );
46
47
48 cout << "---------computing the Janet basis/Pommaret basis:------------------------" << endl;
49 Involutive::JBMill jbMill(Involutive::JBMill::Builder().setInput(polys));
50 Involutive::PBMill mill(Involutive::PBMill::Converter().setJBMill(jbMill));
51
52 cout << "Janet basis" << endl;
53 vector<RingElem> basis = mill.myReturnPB();
54 for(std::vector<RingElem>::const_iterator i = basis.begin(); i != basis.end(); ++i)
55 {
56 cout << *i << endl;
57 }
58
59 cout << "-----computing the depth of the ideal" << endl;
60 cout << "depth = " << mill.myDepth() << endl << endl;
61
62 cout << "-----computing the projective dimension" << endl;
63 cout << "projDim = " << mill.myProjDim() << endl;
64 cout << "-----test if I is Cohen Macaulay" << endl;
65 cout << "Ideal is a Cohen-Macaulay ideal = " << mill.IamCohenMacaulay() << endl;
66
67
68 cout << "-----regular sequence" << endl;
69 std::vector<RingElem> regularSeq = mill.myRegSeq();
70 for (std::vector<RingElem>::iterator i = regularSeq.begin(); i != regularSeq.end(); ++i)
71 {
72 if(i != regularSeq.begin())
73 {
74 cout << ", " << *i;
75 }
76 else
77 {
78 cout << *i;
79 }
80 }
81 cout << endl;
82
83
84 cout << "-----maxStronglyIndependentSet" << endl;
85 std::vector<RingElem> stronglySet = mill.myMaxStronglyIndependentSet();
86 for (std::vector<RingElem>::iterator i = stronglySet.begin(); i != stronglySet.end(); ++i)
87 {
88 if(i != stronglySet.begin())
89 {
90 cout << ", " << *i;
91 }
92 else
93 {
94 cout << *i;
95 }
96 }
97 cout << endl;
98
99 cout << "---------regularity" << endl;
100 cout << "regularity = " << mill.myRegularity() << endl;
101
102 cout << "----------satiety" << endl;
103 cout << "satiety = " << mill.mySatiety() << endl;
104
105 cout << "----------saturation" << endl;
106 std::vector<RingElem> satur = mill.mySaturation();
107 long counter = 0;
108 for (std::vector<RingElem>::iterator i = satur.begin(); i != satur.end(); ++i)
109 {
110 ++counter;
111 cout << *i << endl;
112 }
113
114 cout << "----------classes stuff" << endl;
115
116 RingElem clsRingElem = x[0]*x[1];
117
118 cout << "cls of " << LPP(clsRingElem) << " = " << mill.myCls(LPP(clsRingElem)) << endl;
119
120 cout << "minimal class of result = " << mill.myMinCls() << endl;
121
122 cout << "all elements with cls 3 in the current basis" << endl;
123 std::vector<RingElem> clsElems = mill.myElementsWithClass(3);
124 for (std::vector<RingElem>::iterator i = clsElems.begin(); i != clsElems.end(); ++i)
125 {
126 cout << *i << endl;
127 }
128
129 cout << endl << endl;
130
131 cout << "extremalBettiNumbers" << endl;
132
133 map<pair<long, long>, long> bettis = mill.myExtremalBettiNumbers();
134 for (map<pair<long, long>, long>::iterator i = bettis.begin(); i != bettis.end(); ++i)
135 {
136 cout << "(" << i->first.first << ", " << i->first.second << ")" << "=" << i->second << endl;
137 }
138 }
139
140 //----------------------------------------------------------------------
141 // Use main() to handle any uncaught exceptions and warn the user about them.
main()142 int main()
143 {
144 try
145 {
146 program();
147 return 0;
148 }
149 catch (const CoCoA::ErrorInfo& err)
150 {
151 cerr << "***ERROR*** UNCAUGHT CoCoA error";
152 ANNOUNCE(cerr, err);
153 }
154 catch (const std::exception& exc)
155 {
156 cerr << "***ERROR*** UNCAUGHT std::exception: " << exc.what() << endl;
157 }
158 catch(...)
159 {
160 cerr << "***ERROR*** UNCAUGHT UNKNOWN EXCEPTION" << endl;
161 }
162 return 1;
163 }
164