1 | ////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="version hdepth.lib 4.0.0.0 Jun_2013"; // |
---|
3 | category="Commutative Algebra"; |
---|
4 | info=" |
---|
5 | LIBRARY: hdepth.lib Procedures for computing hdepth_1 |
---|
6 | AUTHORS: Popescu, A., popescu@mathematik.uni-kl.de |
---|
7 | |
---|
8 | SEE ALSO: 'An algorithm to compute the Hilbert depth', Adrian Popescu, arxiv/AC/1307.6084 |
---|
9 | |
---|
10 | KEYWORDS: hdepth, library |
---|
11 | |
---|
12 | |
---|
13 | PROCEDURES: |
---|
14 | hdepth(M [,debug]); hdepth_1 computation of a module M (wrt Z-grading) |
---|
15 | hdepth_p(g, d, debug) the minimum number t <= d s.t. 1/g^t is positive |
---|
16 | "; |
---|
17 | |
---|
18 | /////////////////////////////////////////////////////////////////////////////////// |
---|
19 | static proc myinverse(poly p, int bound) |
---|
20 | "USAGE: myinverse(p,bound), p polynomial in one variable with p(0) nonzero, bound a nonnegative integer |
---|
21 | RETURN: poly, the inverse of the poly p in the power series ring till order bound |
---|
22 | " |
---|
23 | { |
---|
24 | if(bound<=1) |
---|
25 | { |
---|
26 | ERROR("My inverse : negative bound in the inverse"); |
---|
27 | } |
---|
28 | if(p == 0) |
---|
29 | { |
---|
30 | ERROR("My inverse : p is 0"); |
---|
31 | } |
---|
32 | poly original; |
---|
33 | original = p; |
---|
34 | if(leadcoef(p) == 0) |
---|
35 | { |
---|
36 | ERROR("My inverse : the power series is not a unit."); |
---|
37 | } |
---|
38 | poly q = 1/p[1]; |
---|
39 | poly res = q; |
---|
40 | p = q * (p[1] - jet(p,bound)); |
---|
41 | poly s = p; |
---|
42 | while(p != 0) |
---|
43 | { |
---|
44 | res = res + q * p; |
---|
45 | p = jet(p*s,bound); |
---|
46 | } |
---|
47 | //TEST |
---|
48 | if(jet(original*res,bound) != poly(1)) |
---|
49 | { |
---|
50 | ERROR("Myinverse does not work properly."); |
---|
51 | } |
---|
52 | return(res); |
---|
53 | } |
---|
54 | |
---|
55 | /////////////////////////////////////////////////////////////////////////////////// |
---|
56 | static proc hilbconstruct(intvec v) |
---|
57 | "USAGE: hilbconstruct(v), v is the result of hilb(M,2) |
---|
58 | RETURN: poly, the Hilbert Series of M |
---|
59 | AASUME: the ring when called is R = 0,t,ds; |
---|
60 | " |
---|
61 | { |
---|
62 | poly f; |
---|
63 | int i; |
---|
64 | for(i=0;i<size(v)-1;i++) |
---|
65 | { |
---|
66 | f=f+v[i+1]*t^i; |
---|
67 | } |
---|
68 | return(f); |
---|
69 | } |
---|
70 | /////////////////////////////////////////////////////////////////////////////////// |
---|
71 | static proc positiv(poly f) |
---|
72 | "USAGE: positiv(f), f is a polynomial |
---|
73 | RETURN: int, 1 if all the coefficients of f are positive, 0 else |
---|
74 | " |
---|
75 | { |
---|
76 | int pos=1; |
---|
77 | while( (f!=0) && (pos==1) ) |
---|
78 | { |
---|
79 | if(leadcoef(f)<0) |
---|
80 | { |
---|
81 | pos=0; |
---|
82 | } |
---|
83 | f=f-lead(f); |
---|
84 | } |
---|
85 | return(pos); |
---|
86 | } |
---|
87 | /////////////////////////////////////////////////////////////////////////////////// |
---|
88 | static proc sumcoef(poly f) |
---|
89 | "USAGE: sumcoef(f), f is a polynomial |
---|
90 | RETURN: number, the sum of the coefficients |
---|
91 | " |
---|
92 | { |
---|
93 | number c; |
---|
94 | while(f!=0) |
---|
95 | { |
---|
96 | c = c+leadcoef(f); |
---|
97 | f=f-lead(f); |
---|
98 | } |
---|
99 | return(int(c)); |
---|
100 | } |
---|
101 | /////////////////////////////////////////////////////////////////////////////////// |
---|
102 | proc hdepth_p(poly g, int d, int debug) |
---|
103 | "USAGE: hdepth_p(g,d,debug), g is the Hilbert Series of a module M and d is the dimension of M, for debug = 0 the steps will be printed. |
---|
104 | RETURN: int, the minimum number t <= d s.t. 1/g^t is positive |
---|
105 | " |
---|
106 | { |
---|
107 | int dd = d; |
---|
108 | if(debug == 0) |
---|
109 | {"G(t)=",g;} |
---|
110 | if(positiv(g)==1) |
---|
111 | { |
---|
112 | if(debug == 0) |
---|
113 | {return("hdepth =",dd);} |
---|
114 | else |
---|
115 | {return(dd);} |
---|
116 | } |
---|
117 | poly f=g; |
---|
118 | number ag; |
---|
119 | int c1; |
---|
120 | int bound; |
---|
121 | bound = deg(g); |
---|
122 | while(dd >= 0) |
---|
123 | { |
---|
124 | dd = dd-1; |
---|
125 | f = jet( g*myinverse( (1-t)^(d-dd),2*bound ) , bound ); |
---|
126 | if(positiv(f) == 1) |
---|
127 | { |
---|
128 | if(debug == 0) |
---|
129 | { |
---|
130 | "G(t)/(1-t)^",d-dd,"=",f,"+..."; |
---|
131 | return("hdepth =",dd); |
---|
132 | } |
---|
133 | else |
---|
134 | {return(d);} |
---|
135 | } |
---|
136 | c1=sumcoef(f); |
---|
137 | if(c1<0) |
---|
138 | { |
---|
139 | while(c1<0) |
---|
140 | { |
---|
141 | bound=bound+1; |
---|
142 | f=jet( g*myinverse( (1-t)^(d-dd),2*bound ) , bound ); |
---|
143 | c1=sumcoef(f); |
---|
144 | } |
---|
145 | } |
---|
146 | if(debug == 0) |
---|
147 | {"G(t)/(1-t)^",d-dd,"=",f,"+...";} |
---|
148 | } |
---|
149 | ERROR("g was not a Hilbert Series since the coefficient sum is not > 0"); |
---|
150 | } |
---|
151 | example |
---|
152 | { |
---|
153 | "EXAMPLE:";echo=2; |
---|
154 | ring R = 0,t,ds; |
---|
155 | poly f = 2-3t-2t2+2t3+4t4; |
---|
156 | hdepth_p(f,5,0); |
---|
157 | hdepth_p(f,5,1); |
---|
158 | } |
---|
159 | |
---|
160 | /////////////////////////////////////////////////////////////////////////////// |
---|
161 | proc hdepth(module M, list #) |
---|
162 | "USAGE: hdepth(M [,debug]); M is a module, if one want to print the steps debug = 0 |
---|
163 | RETURN: int |
---|
164 | PURPOSE: compute the hdepth_1 of a module M |
---|
165 | EXAMPLE: example hdepth; shows examples |
---|
166 | " |
---|
167 | { |
---|
168 | int debug; |
---|
169 | if(size(#)>0) |
---|
170 | { |
---|
171 | if(typeof(#[1])=="int") |
---|
172 | {debug = #[1];} |
---|
173 | } |
---|
174 | else |
---|
175 | {debug = 1;} |
---|
176 | M = std(M); |
---|
177 | int d=nvars(basering)-dim(M); |
---|
178 | intvec v=hilb(M,2); |
---|
179 | ring R = 0,t,ds; |
---|
180 | poly hp=hilbconstruct(v); |
---|
181 | if(debug == 0) |
---|
182 | {"dim =",d;} |
---|
183 | return(hdepth_p(hp,d,debug)); |
---|
184 | } |
---|
185 | example |
---|
186 | { |
---|
187 | "EXAMPLE:";echo=2; |
---|
188 | ring R = 0,(x(1..10)),dp; |
---|
189 | ideal i=maxideal(1); |
---|
190 | module M=i; |
---|
191 | hdepth(M); |
---|
192 | hdepth(M,0); |
---|
193 | hdepth(M,1); |
---|
194 | } |
---|