1 | version="$I@comment d:$"; |
---|
2 | category="Singularities"; |
---|
3 | info=" |
---|
4 | LIBRARY: alexpoly.lib Resolution Graph and Alexander Polynomial |
---|
5 | AUTHOR: Fernando Hernando Carrillo, hernando@agt.uva.es |
---|
6 | Thomas Keilen, keilen@mathematik.uni-kl.de |
---|
7 | |
---|
8 | OVERVIEW: |
---|
9 | A library for computing the resolution graph of a plane curve singularity f, |
---|
10 | the total multiplicities of the total transforms of the branches of f along |
---|
11 | the exceptional divisors of a minimal good resolution of f, the Alexander |
---|
12 | polynomial of f, and the zeta function of its monodromy operator. |
---|
13 | |
---|
14 | PROCEDURES: |
---|
15 | resolutiongraph(f); resolution graph f |
---|
16 | totalmultiplicities(f); resolution graph, total multiplicities and strict multiplicities of f |
---|
17 | alexanderpolynomial(f); Alexander polynomial of f |
---|
18 | semigroup(f); calculates generators for the semigroup of f |
---|
19 | multseq2charexp(v); converts multiplicity sequence to characteristic exponents |
---|
20 | charexp2multseq(v); converts characteristic exponents to multiplicity sequence |
---|
21 | charexp2generators(v); converts characteristic exponents to generators of the semigroup |
---|
22 | charexp2inter(c,e); converts contact matrix and charact. exp. to intersection matrix |
---|
23 | charexp2conductor(v); converts characteristic exponents to conductor |
---|
24 | charexp2poly(v,a); calculates a polynomial f with characteristic exponents v |
---|
25 | tau_es(f); equsingular Tjurina number of f |
---|
26 | |
---|
27 | KEYWORDS: Hamburger-Noether expansion; Puiseux expansion; curve singularities; |
---|
28 | topological invariants; Alexander polynomial; resolution graph; |
---|
29 | total multiplicities; equisingular Tjurina number |
---|
30 | "; |
---|
31 | |
---|
32 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
33 | LIB "hnoether.lib"; |
---|
34 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
35 | |
---|
36 | proc resolutiongraph (def INPUT,list #) |
---|
37 | "USAGE: resolutiongraph(INPUT); INPUT poly or list |
---|
38 | ASSUME: INPUT is either a REDUCED bivariate polynomial defining a plane curve singularity, |
---|
39 | or the output of @code{hnexpansion(f[,\"ess\"])}, or the list @code{hne} in |
---|
40 | the ring created by @code{hnexpansion(f[,\"ess\"])}, or the output of |
---|
41 | @code{develop(f)} resp. of @code{extdevelop(f,n)}, or a list containing |
---|
42 | the contact matrix and a list of integer vectors with the characteristic exponents |
---|
43 | of the branches of a plane curve singularity, or an integer vector containing the |
---|
44 | characteristic exponents of an irreducible plane curve singularity. |
---|
45 | RETURN: intmat, the incidence matrix of the resolution graph of the plane curve |
---|
46 | defined by INPUT, where the entries on the diagonal are the weights of the |
---|
47 | vertices of the graph and a negative entry corresponds to the strict transform |
---|
48 | of a branch of the curve. |
---|
49 | NOTE: In case the Hamburger-Noether expansion of the curve f is needed |
---|
50 | for other purposes as well it is better to calculate this first |
---|
51 | with the aid of @code{hnexpansion} and use it as input instead of |
---|
52 | the polynomial itself. |
---|
53 | @* |
---|
54 | If you are not sure whether the INPUT polynomial is reduced or not, use |
---|
55 | @code{squarefree(INPUT)} as input instead. |
---|
56 | SEE ALSO: develop, hnexpansion, totalmultiplicities, alexanderpolyinomal |
---|
57 | EXAMPLE: example resolutiongraph; shows an example |
---|
58 | " |
---|
59 | { |
---|
60 | return(totalmultiplicities(INPUT,#)[1]); |
---|
61 | } |
---|
62 | example |
---|
63 | { |
---|
64 | "EXAMPLE:"; |
---|
65 | echo=2; |
---|
66 | ring r=0,(x,y),ls; |
---|
67 | poly f1=(y2-x3)^2-4x5y-x7; |
---|
68 | poly f2=y2-x3; |
---|
69 | poly f3=y3-x2; |
---|
70 | resolutiongraph(f1*f2*f3); |
---|
71 | } |
---|
72 | |
---|
73 | proc totalmultiplicities (def INPUT, list #) |
---|
74 | "USAGE: totalmultiplicities(INPUT); INPUT poly or list |
---|
75 | ASSUME: INPUT is either a REDUCED bivariate polynomial defining a plane curve singularity, |
---|
76 | or the output of @code{hnexpansion(f[,\"ess\"])}, or the list @code{hne} in |
---|
77 | the ring created by @code{hnexpansion(f[,\"ess\"])}, or the output of |
---|
78 | @code{develop(f)} resp. of @code{extdevelop(f,n)}, or a list containing |
---|
79 | the contact matrix and a list of integer vectors with the characteristic exponents |
---|
80 | of the branches of a plane curve singularity, or an integer vector containing the |
---|
81 | characteristic exponents of an irreducible plane curve singularity. |
---|
82 | RETURN: list @code{L} of three integer matrices. @code{L[1]} is the incidence matrix of |
---|
83 | the resolution graph of the plane curve defined by INPUT, where the entries on the |
---|
84 | diagonal are the weights of the vertices of the graph and a negative entry corresponds |
---|
85 | to the strict transform of a branch of the curve. @code{L[2]} is an integer matrix, |
---|
86 | which has for each vertex in the graph a row and for each branch of the curve a column. |
---|
87 | The entry in position [i,j] contains the total multiplicity of the j-th branch (i.e. the |
---|
88 | branch with weight -j in @code{L[1]}) along the exceptional divisor corresponding |
---|
89 | to the i-th row in @code{L[1]}. In particular, the i-th row contains |
---|
90 | the total multiplicities of the branches of the plane curve (defined by INPUT) along |
---|
91 | the exceptional divisor which corresponds to the i-th row in the incidence matrix |
---|
92 | @code{L[1]}. @code{L[3]} is an integer matrix which contains the (strict) multiplicities |
---|
93 | of the branches of the curve along the exceptional divisors in the same way as @code{L[2]} |
---|
94 | contains the total multiplicities. |
---|
95 | NOTE: The total multiplicty of a branch along an exceptional divisor is the multiplicity |
---|
96 | with which this exceptional divisor occurs in the total transform of this branch |
---|
97 | under the resolution corresponding to the resolution graph. |
---|
98 | @* |
---|
99 | In case the Hamburger-Noether expansion of the curve f is needed |
---|
100 | for other purposes as well it is better to calculate this first |
---|
101 | with the aid of @code{hnexpansion} and use it as input instead of |
---|
102 | the polynomial itself. |
---|
103 | @* |
---|
104 | If you are not sure whether the INPUT polynomial is reduced or not, use |
---|
105 | @code{squarefree(INPUT)} as input instead. |
---|
106 | SEE ALSO: develop, hnexpansion, alexanderpolynomial, resolutiongraph |
---|
107 | EXAMPLE: example totalmultiplicities; shows an example |
---|
108 | " |
---|
109 | { |
---|
110 | /////////////////////////////////////////////////////////////////////////////////////////////// |
---|
111 | /// The algorithm described in [JP] DeJong, Pfister, Local Analytic Geometry, Vieweg 2000, |
---|
112 | /// is used for the implementation -- the missing case is included by appropriate sorting. |
---|
113 | /////////////////////////////////////////////////////////////////////////////////////////////// |
---|
114 | /// If # is empty, then an additional sorting of the branches will be made, so that |
---|
115 | /// the branches can be split with split_graph in order to get the graphs of the curves |
---|
116 | /// separated on the first exceptional divisor. Otherwise split_graph should not be applied! |
---|
117 | /////////////////////////////////////////////////////////////////////////////////////////////// |
---|
118 | /// If #[1]="tau", then totalmultiplicities is called for the use in tau_es - thus it returns |
---|
119 | /// the prolonged resolutiongraphs of the branches, their multiplicity sequences and their |
---|
120 | /// contact matrix - sorted appropriately. |
---|
121 | /////////////////////////////////////////////////////////////////////////////////////////////// |
---|
122 | int i,j,s,e,ii; |
---|
123 | intmat helpmati,helpmatii; |
---|
124 | intvec helpvec; |
---|
125 | ///////////////////////////////////////////////////////////////////////////////// |
---|
126 | // Decide, which kind of input we have, and define the contact matrix |
---|
127 | ///////////////////////////////////////////////////////////////////////////////// |
---|
128 | // Input: a polynomial defining a plane curve singularity. |
---|
129 | ////////////////////////////////////////////////////////////////////////////// |
---|
130 | if (typeof(INPUT)=="poly") |
---|
131 | { |
---|
132 | /* |
---|
133 | poly f=squarefree(INPUT); |
---|
134 | if ( deg(f)!=deg(INPUT) ) |
---|
135 | { |
---|
136 | dbprint(printlevel-voice+4,"// input polynomial was not reduced"); |
---|
137 | dbprint(printlevel-voice+4,"// we continue with its reduction"); |
---|
138 | } |
---|
139 | list I@N@V=invariants(f); |
---|
140 | */ |
---|
141 | list I@N@V=invariants(INPUT); |
---|
142 | } |
---|
143 | else |
---|
144 | { |
---|
145 | /////////////////////////////////////////////////////////////////////////////////// |
---|
146 | // Input: the vector of characteristic exponents of an irreducible plane curve |
---|
147 | /////////////////////////////////////////////////////////////////////////////////// |
---|
148 | if (typeof(INPUT)=="intvec") |
---|
149 | { |
---|
150 | list charexp; |
---|
151 | charexp[1]=INPUT; |
---|
152 | intmat contact[1][1]=0; |
---|
153 | } |
---|
154 | else |
---|
155 | { |
---|
156 | if (typeof(INPUT)=="list") |
---|
157 | { |
---|
158 | ///////////////////////////////////////////////////////////////////////////////// |
---|
159 | // Input: intersection-matrix and characteristic exponents. |
---|
160 | ////////////////////////////////////////////////////////////////////////////// |
---|
161 | if (typeof(INPUT[1])=="intmat") |
---|
162 | { |
---|
163 | intmat contact=INPUT[1]; |
---|
164 | list charexp=INPUT[2]; |
---|
165 | } |
---|
166 | else |
---|
167 | { |
---|
168 | ///////////////////////////////////////////////////////////////////////////////// |
---|
169 | // Input: output of hnexpansion or hne in the ring created by hnexpansion |
---|
170 | ////////////////////////////////////////////////////////////////////////////// |
---|
171 | if ((typeof(INPUT[1])=="ring") or (typeof(INPUT[1])=="list")) |
---|
172 | { |
---|
173 | list I@N@V=invariants(INPUT); |
---|
174 | } |
---|
175 | else |
---|
176 | { |
---|
177 | //////////////////////////////////////////////////////////////////////////////////// |
---|
178 | // Input: output of reddevelop or extdevelop -- irreducible plane curve singularity |
---|
179 | //////////////////////////////////////////////////////////////////////////////////// |
---|
180 | if (typeof(INPUT[1])=="matrix") |
---|
181 | { |
---|
182 | list charexp=invariants(INPUT)[1]; |
---|
183 | intmat contact[1][1]=0; |
---|
184 | } |
---|
185 | else |
---|
186 | { |
---|
187 | ERROR("The input is invalid!"); |
---|
188 | } |
---|
189 | } |
---|
190 | } |
---|
191 | } |
---|
192 | else |
---|
193 | { |
---|
194 | ERROR("The input is invalid!"); |
---|
195 | } |
---|
196 | } |
---|
197 | } |
---|
198 | /////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
199 | // If the input was a poly or a HN-Expansion, then calculate the contact matrix and char.exponents. |
---|
200 | /////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
201 | if (defined(I@N@V)) |
---|
202 | { |
---|
203 | intmat contact=I@N@V[size(I@N@V)][1]; // contact numbers |
---|
204 | list charexp; // characteristic exponents |
---|
205 | for (i=1;i<=size(I@N@V)-1;i++) |
---|
206 | { |
---|
207 | charexp[i]=I@N@V[i][1]; |
---|
208 | } |
---|
209 | } |
---|
210 | ////////////////////////////////////////////////////////////////////////////// |
---|
211 | // Find the maximal contact numbers of the branches |
---|
212 | ////////////////////////////////////////////////////////////////////////////// |
---|
213 | int r=ncols(contact); // number of branches |
---|
214 | intvec maxcontact; // maximal contactnumber of branch i with other branches |
---|
215 | for (i=1;i<=r;i++) |
---|
216 | { |
---|
217 | maxcontact[i]=max_in_intvec(intvec(contact[i,1..r])); |
---|
218 | } |
---|
219 | /////////////////////////////////////////////////////////////////////////////// |
---|
220 | // Define the graphs of the branches and prolong them if necessary. |
---|
221 | /////////////////////////////////////////////////////////////////////////////// |
---|
222 | intmat gr,gr_red,grp; // a subgraph, a reduced subgraph, a prolonged subgraph |
---|
223 | int omega; // point of highest weight in subgraph |
---|
224 | list graphs; // contains the subgraphs of the C_i |
---|
225 | list totmult; // contains the total multiplicities of the subgraphs |
---|
226 | list multipl; // contains the multiplicities of the subgraphs |
---|
227 | list gr_tm; // takes a single subgraph and tot.mult. |
---|
228 | intvec tm,mt; // total multiplicities and multiplicities |
---|
229 | for (i=1;i<=r;i++) |
---|
230 | { |
---|
231 | gr_tm=irred_resgraph_totmult(charexp[i]); // graph, total multipl. and multipl. of the ith branch |
---|
232 | gr=gr_tm[1]; // the graph of the ith branch |
---|
233 | tm=gr_tm[2]; // the total multiplicities of the ith branch |
---|
234 | mt=gr_tm[3]; // the multiplicities of the ith branch |
---|
235 | omega=nrows(gr)-1; // maximal weight in the graph of the ith branch |
---|
236 | // If the maximal contact of the ith branch with some other branch is larger |
---|
237 | // than the maximal weight omega, then the graph has to be polonged. |
---|
238 | if (omega<maxcontact[i]) |
---|
239 | { |
---|
240 | grp=intmat(intvec(0),maxcontact[i]+1,maxcontact[i]+1); |
---|
241 | // save the graph without the point of the strict transform |
---|
242 | if (omega>=1) // otherwise the graph consists only of the strict transform |
---|
243 | { |
---|
244 | grp[1..omega,1..omega]=gr[1..omega,1..omega]; |
---|
245 | } |
---|
246 | // add the points of multiplicity 1 up to weight maxcontact[i] and the strict transform |
---|
247 | // and connect them |
---|
248 | for (j=omega+1;j<=maxcontact[i]+1;j++) |
---|
249 | { |
---|
250 | // adding the vertex to the graph and adding the total multiplicities |
---|
251 | if (j<=maxcontact[i]) |
---|
252 | { |
---|
253 | grp[j,j]=j; |
---|
254 | if (j>1) |
---|
255 | { |
---|
256 | tm[j]=tm[j-1]+1; |
---|
257 | mt[j]=1; |
---|
258 | } |
---|
259 | else // then there is no previous point in the graph |
---|
260 | { |
---|
261 | tm[j]=1; |
---|
262 | mt[j]=1; |
---|
263 | } |
---|
264 | } |
---|
265 | // connecting the vertex with the previous part of the graph |
---|
266 | if (j>1) // otherwise there is nothing to connect to |
---|
267 | { |
---|
268 | grp[j-1,j]=1; |
---|
269 | grp[j,j-1]=1; |
---|
270 | } |
---|
271 | } |
---|
272 | gr=grp; // replace the subgraph by the prolonged subgraph |
---|
273 | } |
---|
274 | gr[nrows(gr),ncols(gr)]=-i; // give the strict transform in the ith branch weight -i |
---|
275 | graphs[i]=gr; |
---|
276 | totmult[i]=tm; |
---|
277 | multipl[i]=mt; |
---|
278 | } |
---|
279 | /////////////////////////////////////////////////////////////////////////////// |
---|
280 | // Check, if the procedure is called from inside tau_es. |
---|
281 | int check_tau; |
---|
282 | if (size(#)==1) |
---|
283 | { |
---|
284 | if (#[1]=="tau") |
---|
285 | { |
---|
286 | check_tau=1; |
---|
287 | } |
---|
288 | } |
---|
289 | ///////////////////////////////////////////////////////////////////////////////// |
---|
290 | // The procedure tau_es expects the branches to be ordered according to their |
---|
291 | // contact during the resolution process. |
---|
292 | ///////////////////////////////////////////////////////////////////////////////// |
---|
293 | if ((size(#)==0) or (check_tau==1)) |
---|
294 | { |
---|
295 | list SORT; |
---|
296 | for (i=1;i<=ncols(contact)-2;i++) |
---|
297 | { |
---|
298 | SORT=sort_branches(contact,graphs,totmult,multipl,i,ncols(contact)); |
---|
299 | contact=SORT[1]; |
---|
300 | graphs=SORT[2]; |
---|
301 | totmult=SORT[3]; |
---|
302 | multipl=SORT[4]; |
---|
303 | } |
---|
304 | } |
---|
305 | /////////////////////////////////////////////////////////////////////////////////// |
---|
306 | /// If the procedure is called from inside tau_es, the output will be the prolonged |
---|
307 | /// resolutiongraphs of the branches, their multiplicity sequences and their contact matrix. |
---|
308 | if (check_tau==1) |
---|
309 | { |
---|
310 | return(list(graphs,multipl,contact)); |
---|
311 | } |
---|
312 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
313 | // Sort the branches in such a way, that in the blowing up procedure whenever to branches |
---|
314 | // C_i and C_j have contact k (i.e. after k steps they are separated) and one of the following |
---|
315 | // situations occurs: |
---|
316 | // A) some exceptional divisor E_l (l<k) still intersects C_i after k steps of blowing up |
---|
317 | // and no such E_l intersects C_j |
---|
318 | // OR |
---|
319 | // B) some exceptional divisor E_l (l<k-1) still intersects C_i after k steps of blowing up |
---|
320 | // and another one (necessarily E_k-1) intersects C_j, |
---|
321 | // THEN: i<j! |
---|
322 | ///////////////////////////////////////////////////////////////////////////////////////// |
---|
323 | // This means in the algorithm for glueing graphs together described in [JP] p.217-19 |
---|
324 | // the Case 3 never occurs -- and Case 1 only occurs, if it is unavoidable, that is when |
---|
325 | // C_1 there meets an E_l with l<k. |
---|
326 | // (Otherwise we say that (C_i,C_j) has "bad contact".) |
---|
327 | //////////////////////////////////////////////////////////////////////////////////////// |
---|
328 | // We can detect this by looking at the graphs of C_i and C_j. If E_l stays together with |
---|
329 | // C_i and not with C_j in the k-th step of blowing up, then in the graph of C_i |
---|
330 | // k and l are NOT connected, while in the graph of C_j they are! |
---|
331 | ////////////////////////////////////////////////////////////////////////////////////////// |
---|
332 | // Note also, that for a fixed pair (i,j) this situation can only occur for ONE l (in Case A: l<k |
---|
333 | // and possibly l=k-1; in Case B: l<k-1). |
---|
334 | ////////////////////////////////////////////////////////////////////////////////////////// |
---|
335 | // Our algorithm now works as follows. Start with i=1 and j=2. |
---|
336 | // While (i<r=number of branches) do as follows: Check C_i with C_j. |
---|
337 | // If a bad contact (C_i,C_j) occurs, then MOVE C_j to the position of C_i (and move all C_t |
---|
338 | // with t>=i, t!=j one position further). Else do nothing particular. |
---|
339 | // Then, if j<r set j=j+1. Else set j=i+2 and i=i+1. End of While. |
---|
340 | // To see that this works we note that if (C_i,C_j) has bad contact for some j>i, then |
---|
341 | // (C_j,C_t) does NOT have bad contact for t=i+1,...,j-1. For this we consider the 3 cases: |
---|
342 | // 1) contact(C_t,C_j)=contact(C_i,C_j)=k: E_l stays with C_j in k-th step and C_t and C_j |
---|
343 | // separate there, so (C_t,C_j) has bad contact |
---|
344 | // and hence (C_j,C_t) does NOT have bad contact |
---|
345 | // 2) contact(C_t,C_j)>contact(C_i,C_j)=k: CANNOT HAPPEN - C_i had no bad contact with C_t, |
---|
346 | // hence E_l did not stay with C_t in the k-th step; |
---|
347 | // however, E_l stays with C_j there, so C_t and C_j |
---|
348 | // cannot have a contact higher than k |
---|
349 | // 3) contact(C_t,C_j)<contact(C_i,C_j)=k: when C_t and C_j separate, C_j and C_i are still |
---|
350 | // together; and since (C_i,C_t) did not have bad |
---|
351 | // contact also (C_j,C_t) cannot have bad contact. |
---|
352 | // This ensures that when we do the moving of the graphs, we do NOT create any new bad contacts. |
---|
353 | ////////////////////////////////////////////////////////////////////////////// |
---|
354 | i=1; |
---|
355 | j=2; |
---|
356 | intvec connections_i,connections_j; |
---|
357 | int n_i,n_j; |
---|
358 | int bad_contact; |
---|
359 | while (i<r) |
---|
360 | { |
---|
361 | // Test graphs[i] with graphs[j] for bad contact. |
---|
362 | connections_i=find_lower_connection_points(graphs[i],contact[i,j]); |
---|
363 | if (connections_i[1]==0) // no connection point of lower weight there |
---|
364 | { |
---|
365 | n_i=0; |
---|
366 | } |
---|
367 | else |
---|
368 | { |
---|
369 | n_i=size(connections_i); |
---|
370 | } |
---|
371 | connections_j=find_lower_connection_points(graphs[j],contact[i,j]); |
---|
372 | if (connections_j[1]==0) // no connection point of lower weight there |
---|
373 | { |
---|
374 | n_j=0; |
---|
375 | } |
---|
376 | else |
---|
377 | { |
---|
378 | n_j=size(connections_j); |
---|
379 | } |
---|
380 | bad_contact=0; |
---|
381 | ii=1; |
---|
382 | // The points of weight k=contact(C_i,C_j) in C_i and C_j are connectd to n_i respectively |
---|
383 | // n_j points of weight less than k, where n_i and n_j might take values among 0, 1 and 2. |
---|
384 | // n_j=2: Then E_k is intersected by E_k-1 and E_l (l<k-1) in the k-th step, however |
---|
385 | // C_j does not stay with anyone of those, so (C_i,C_j) does not have bad contact. |
---|
386 | // n_i=0: Then C_i stays with E_k-1 and there is no E_l with l<k-1 in the k-th step, |
---|
387 | // hence (C_i,C_j) does not have bad contact. |
---|
388 | // n_i=1, n_j=0: Analogously, then C_j stays with E_k-1 and there is no E_l (l<k-1). This means |
---|
389 | // (C_i,C_j) has bad contact. |
---|
390 | // n_i=1, n_j=1: EITHER there is no E_l (l<k-1) and C_i and C_j both separate from E_k-1 |
---|
391 | // in the k-th step of blowing up (hence (C_i,C_j) does not have bad contact). |
---|
392 | // OR there is an E_l (l<k-1) with which one stays and the other one stays |
---|
393 | // with E_k-1 (bad contact, if E_l stays with C_j; good contact otherwise). |
---|
394 | // n_i=2, E_k is intersected by E_k-1 and E_l in the k-th step and C_i does not |
---|
395 | // stay together with any of those. |
---|
396 | // n_j=0: This CANNOT occur, since then C_j would have to stay together with E_l |
---|
397 | // and E_k-1, which is impossible. |
---|
398 | // n_j=1: Then C_j stays together with either E_l or with E_k-1, however both cases |
---|
399 | // imply that (C_i,C_j) has bad contact |
---|
400 | // ALTERNATIVE CONSIDERATION FOR DISTINGUISHING THE CASES: |
---|
401 | // n_i<n_j : C_i stays with "more" E's than C_j does (only one is possible of course!), |
---|
402 | // hence (C_i,C_j) does not have bad contact. |
---|
403 | // n_i>n_j : C_i stays with "less" E's than C_j does (i.e. with none) hence (C_i,C_j) has bad contact. |
---|
404 | // n_i=n_j=0 : Not possible, since then both would stay with E_k-1. |
---|
405 | // n_i=n_j=2 : Both separate from E_k-1 as well as from some E_l (l<k) - hence NO bad contact. |
---|
406 | // n_i=n_j=1 : One stays with E_k-1 and one with E_l (l<k). Bad contact, if C_j stays with E_l, |
---|
407 | // i.e. if in the graph of C_i the point k is connected to l, and no bad contact otherwise. |
---|
408 | if ((n_i>n_j) or ((n_i==1) and (n_j==1) and (connections_i[1]<contact[i,j]-1))) |
---|
409 | { |
---|
410 | bad_contact=1; |
---|
411 | // Move the graph (etc.) of C_j into the position i. |
---|
412 | graphs=insert(graphs,graphs[j],i-1); |
---|
413 | graphs=delete(graphs,j+1); |
---|
414 | totmult=insert(totmult,totmult[j],i-1); |
---|
415 | totmult=delete(totmult,j+1); |
---|
416 | multipl=insert(multipl,multipl[j],i-1); |
---|
417 | multipl=delete(multipl,j+1); |
---|
418 | contact=move_row_col(contact,i,j); |
---|
419 | } |
---|
420 | if (j<r) // There are still some C_j against which C_i has to be tested. |
---|
421 | { |
---|
422 | j++; |
---|
423 | } |
---|
424 | else // Now C_i has been tested against all C_j, and we may continue with C_i+1. |
---|
425 | { |
---|
426 | j=i+2; |
---|
427 | i=i+1; |
---|
428 | } |
---|
429 | } |
---|
430 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
431 | // Glue the graphs together. |
---|
432 | ////////////////////////////////////////////////////////////////////////////// |
---|
433 | /////////////////////////////////////////////////////////////////////////////////// |
---|
434 | intmat rgraph=graphs[1]; // keeps the resolution graph |
---|
435 | intmat rtm=intmat(totmult[1],nrows(rgraph),1); // keeps the tot.mult. at the vertices of the graph as rows |
---|
436 | intmat rmt=intmat(multipl[1],nrows(rgraph),1); // keeps the mult. at the vertices of the graph as rows |
---|
437 | intvec stricttransforms=0,ncols(rgraph); // keeps the position of the ith strict |
---|
438 | // transform in rgraph at position i+1 !!! |
---|
439 | intvec k,kp,p,q,o; // highest contact numbers, num. of branch with hgt contact, separation points |
---|
440 | int maxc,maxcp; |
---|
441 | for (i=2;i<=r;i++) |
---|
442 | { |
---|
443 | ////////////////////////////////////////////////////////////////////////////////// |
---|
444 | // Find j<i minimal s.t. contact[i,j] is maximal. |
---|
445 | maxcp=i; |
---|
446 | for (j=1;j<i;j++) |
---|
447 | { |
---|
448 | if (contact[i,j]>contact[i,maxcp]){maxcp=j;} |
---|
449 | } |
---|
450 | kp[i]=maxcp; // the j such that C_i has its maximal contact to C_j |
---|
451 | k[i]=contact[i,maxcp]; // the maximal contact of C_i with C_1,...,C_i-1 |
---|
452 | /////////////////////////////////////////////////////////////////////////////////// |
---|
453 | // Find in the graph of C_1,...,C_i-1 the points p of wgt k and q of wgt k-1 |
---|
454 | // connected to C_maxcp. |
---|
455 | // Since non of C_j for j<maxcp has contact k with C_i, the point p lies in |
---|
456 | // the remaining part of the graph of C_maxcp. |
---|
457 | s=rgraph[stricttransforms[maxcp]+1,stricttransforms[maxcp]+1]; |
---|
458 | p[i]=stricttransforms[maxcp]+1+k[i]-s; // pt. to which reduced subgraph of C_i is glued |
---|
459 | // If s<k[i], then q also lies in this part, otherwise it lies in the remaining part |
---|
460 | // of the graph of the C_j to which C_maxcp is connected, i.e. j=kp[maxcp], since |
---|
461 | // the contact of C_i and of C_maxcp to this C_j is strictly less than k. |
---|
462 | // If s=k[i]=1, then p=1 and there is no q! We may thus set q=0. |
---|
463 | if ((s<k[i]) or ((s==k[i]) and (s==1))) // i.e. q is on the same subgraph as p, or q does not exist |
---|
464 | { |
---|
465 | q[i]=p[i]-1; |
---|
466 | } |
---|
467 | else // i.e. q is on the subgraph to which the subgraph of p has been glued |
---|
468 | { |
---|
469 | s=rgraph[stricttransforms[kp[maxcp]]+1,stricttransforms[kp[maxcp]]+1]; |
---|
470 | q[i]=stricttransforms[kp[maxcp]]+k[i]-s; |
---|
471 | } |
---|
472 | ////////////////////////////////////////////////////////////////////////////////////// |
---|
473 | // Reduce the graph of C_i and add it to the graph of C_1,...,C_i-1. |
---|
474 | gr=graphs[i]; |
---|
475 | s=nrows(rgraph); |
---|
476 | // Delete in the graph of C_i everything of weight <=k. |
---|
477 | gr_red=intmat(intvec(gr[k[i]+1..nrows(gr),k[i]+1..ncols(gr)]),nrows(gr)-k[i],ncols(gr)-k[i]); |
---|
478 | // Add this part to the graph of C_1,...,C_i-1. |
---|
479 | rgraph=addmat(rgraph,gr_red); |
---|
480 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
481 | // Connect the two parts of the graph. |
---|
482 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
483 | // Connect the points connected to the point of wgt k in the graph of C_i to p[i]. |
---|
484 | for (j=k[i]+1;j<=ncols(gr);j++) |
---|
485 | { |
---|
486 | if(gr[k[i],j]==1) |
---|
487 | { |
---|
488 | rgraph[s+j-k[i],p[i]]=1; |
---|
489 | rgraph[p[i],s+j-k[i]]=1; |
---|
490 | } |
---|
491 | } |
---|
492 | // If pt. of wgt k is not connected to pt of wgt k-1 in graph of C_i, then points |
---|
493 | // connected to the one of wgt k-1 have to be connected to q[i]. |
---|
494 | if (k[i]>1) |
---|
495 | { |
---|
496 | if (gr[k[i],k[i]-1]!=1) |
---|
497 | { |
---|
498 | for (j=k[i]+1;j<=ncols(gr);j++) |
---|
499 | { |
---|
500 | if(gr[k[i]-1,j]==1) |
---|
501 | { |
---|
502 | rgraph[s+j-k[i],q[i]]=1; |
---|
503 | rgraph[q[i],s+j-k[i]]=1; |
---|
504 | } |
---|
505 | } |
---|
506 | // Cut the connection from p[i] to q[i]. |
---|
507 | rgraph[p[i],q[i]]=0; |
---|
508 | rgraph[q[i],p[i]]=0; |
---|
509 | } |
---|
510 | } |
---|
511 | stricttransforms[i+1]=ncols(rgraph); |
---|
512 | //////////////////////////////////////////////////////////////////////////////// |
---|
513 | // Adjust the total multiplicities |
---|
514 | //////////////////////////////////////////////////////////////////////////////// |
---|
515 | // Add the total multiplicities for the added subgraph to rtm |
---|
516 | tm=totmult[i]; |
---|
517 | mt=multipl[i]; |
---|
518 | if (k[i]<size(tm)) // if the reduced subgraph of C_i has more than one point |
---|
519 | { |
---|
520 | rtm=addmat(rtm,intmat(intvec(tm[k[i]+1..size(tm)]),nrows(gr_red),1)); |
---|
521 | rmt=addmat(rmt,intmat(intvec(mt[k[i]+1..size(mt)]),nrows(gr_red),1)); |
---|
522 | } |
---|
523 | else // if the reduced subgraph of C_i consists only of the strict transform |
---|
524 | { |
---|
525 | rtm=addmat(rtm,0); |
---|
526 | rmt=addmat(rmt,0); |
---|
527 | } |
---|
528 | // Adjust the total multiplicities at the places where the subgraph has been glued. |
---|
529 | e=k[i]; // the highest weight of a point that has not yet been assigned its tot. mult. |
---|
530 | while(e>=1) |
---|
531 | { |
---|
532 | s=stricttransforms[maxcp]+1; // Line in the graph of the start. pt. of the subgraph of C_maxcp. |
---|
533 | for (j=rgraph[s,s];j<=e;j++) // Adjust the multiplicities. |
---|
534 | { |
---|
535 | rtm[s+j-rgraph[s,s],i]=tm[j]; |
---|
536 | rmt[s+j-rgraph[s,s],i]=mt[j]; |
---|
537 | } |
---|
538 | maxcp=kp[maxcp]; // To which subgraph has the subgraph of C_maxcp been glued? |
---|
539 | e=rgraph[s,s]-1; // What is the highest weight of a pt. that has not yet been assigned its tot.mult.? |
---|
540 | } |
---|
541 | e=nrows(rtm); // Number of rows in the matrix of totalmultiplicities. |
---|
542 | // The total multiplicities of the C_s for s=1,...,i-1 along the exceptional divisors |
---|
543 | // which are introduced after the strict transform of C_s has separated (i.e. the entries |
---|
544 | // in rows stricttransform[i]+1,...,stricttransform[i+1]-1 in the s-th column of the matrix |
---|
545 | // of total multiplicities still have to be calculated. |
---|
546 | for (s=1;s<i;s++) |
---|
547 | { |
---|
548 | rtm[1..e,s]=adjust_tot_mult(intvec(rtm[1..e,i]),intvec(rtm[1..e,s]),intvec(rmt[1..e,i]),intvec(rmt[1..e,s]),p,q,stricttransforms,i); |
---|
549 | } |
---|
550 | // The total multiplicities of the C_i along the exceptional divisors |
---|
551 | // which are introduced for the sake of C_s, s=1,...,i-1, after the strict transform |
---|
552 | // of C_i has separated (i.e. the entries in rows stricttransform[s]+1,...,stricttransform[s+1]-1 |
---|
553 | // in the i-th column of the matrix of total multiplicities still have to be calculated. |
---|
554 | for (s=1;s<i;s++) |
---|
555 | { |
---|
556 | rtm[1..e,i]=adjust_tot_mult(intvec(rtm[1..e,s]),intvec(rtm[1..e,i]),intvec(rmt[1..e,s]),intvec(rmt[1..e,i]),p,q,stricttransforms,s); |
---|
557 | } |
---|
558 | } |
---|
559 | list result=rgraph,rtm,rmt; |
---|
560 | return(result); |
---|
561 | } |
---|
562 | example |
---|
563 | { |
---|
564 | "EXAMPLE:"; |
---|
565 | echo=2; |
---|
566 | ring r=0,(x,y),ls; |
---|
567 | poly f1=(y2-x3)^2-4x5y-x7; |
---|
568 | poly f2=y2-x3; |
---|
569 | poly f3=y3-x2; |
---|
570 | totalmultiplicities(f1*f2*f3); |
---|
571 | } |
---|
572 | |
---|
573 | |
---|
574 | |
---|
575 | proc alexanderpolynomial (def INPUT) |
---|
576 | "USAGE: alexanderpolynomial(INPUT); INPUT poly or list |
---|
577 | ASSUME: INPUT is either a REDUCED bivariate polynomial defining a plane curve singularity, |
---|
578 | or the output of @code{hnexpansion(f[,\"ess\"])}, or the list @code{hne} in |
---|
579 | the ring created by @code{hnexpansion(f[,\"ess\"])}, or the output of |
---|
580 | @code{develop(f)} resp. of @code{extdevelop(f,n)}, or a list containing |
---|
581 | the contact matrix and a list of integer vectors with the characteristic exponents |
---|
582 | of the branches of a plane curve singularity, or an integer vector containing the |
---|
583 | characteristic exponents of an irreducible plane curve singularity. |
---|
584 | CREATE: a ring with variables t, t(1), ..., t(r) (where r is the number of branches of |
---|
585 | the plane curve singularity f defined by INPUT) and ordering ls over the |
---|
586 | ground field of the basering. @* |
---|
587 | Moreover, the ring contains the Alexander polynomial of f in variables t(1), ..., t(r) |
---|
588 | (@code{alexpoly}), the zeta function of the monodromy operator of f in the variable t |
---|
589 | (@code{zeta_monodromy}), and a list containing the factors of the Alexander |
---|
590 | polynomial with multiplicities (@code{alexfactors}). |
---|
591 | RETURN: a list, say @code{ALEX}, where @code{ALEX[1]} is the created ring |
---|
592 | NOTE: to use the ring type: @code{def ALEXring=ALEX[i]; setring ALEXring;}. |
---|
593 | @* |
---|
594 | Alternatively you may use the procedure sethnering and type: sethnering(ALEX,\"ALEXring\"); |
---|
595 | @* |
---|
596 | To access the Alexander polynomial resp. the zeta function resp. the |
---|
597 | factors of the Alexander polynomial type: @code{alexpoly} resp. @code{zeta_monodromy} |
---|
598 | resp. @code{alexfactors}.@* |
---|
599 | In case the Hamburger-Noether expansion of the curve f is needed |
---|
600 | for other purposes as well it is better to calculate this first |
---|
601 | with the aid of @code{hnexpansion} and use it as input instead of |
---|
602 | the polynomial itself. |
---|
603 | @* |
---|
604 | If you are not sure whether the INPUT polynomial is reduced or not, use |
---|
605 | @code{squarefree(INPUT)} as input instead. |
---|
606 | SEE ALSO: resolutiongraph, totalmultiplicities |
---|
607 | EXAMPLE: example alexanderpolynomial; shows an example |
---|
608 | " |
---|
609 | { |
---|
610 | // Get the resolution graph and the total multiplicities. |
---|
611 | list gr_tm=totalmultiplicities(INPUT); |
---|
612 | intmat gr=gr_tm[1]; |
---|
613 | intmat tm=gr_tm[2]; |
---|
614 | int r=ncols(tm); |
---|
615 | int e=ncols(gr); |
---|
616 | // Define the Ring for the Alexander Polynomial and the Zeta Function of the Monodromy. |
---|
617 | execute("ring ALEXring=("+charstr(basering)+"),(t,t(1..r)),dp;"); |
---|
618 | poly hilfspoly=1; |
---|
619 | poly alexnumerator=1; // numerator of the Alexander polynomial |
---|
620 | poly alexdenominator=1; // denominator of the Alexander polynomial |
---|
621 | list alexfactors; // the factors of the Alexanderpolynomial with multiplicities |
---|
622 | list alexfactor; |
---|
623 | int chi=2; |
---|
624 | int i,j,k; |
---|
625 | int s=1; |
---|
626 | // Consider every vertex of the resolution graph. |
---|
627 | for (i=1;i<=e;i++) |
---|
628 | { |
---|
629 | if (gr[i,i]>0) // If it belongs to an exceptional curve. |
---|
630 | { |
---|
631 | for (j=1;j<=e;j++) // Calculate the Euler charateristik of the smooth locus of the exc. curve. |
---|
632 | { |
---|
633 | if ((gr[i,j]==1) and (i!=j)) |
---|
634 | { |
---|
635 | chi=chi-1; |
---|
636 | } |
---|
637 | } |
---|
638 | if (chi!=0) // If the Euler characteristik is not zero, then it gives a factor in the AP. |
---|
639 | { |
---|
640 | for (k=1;k<=r;k++) |
---|
641 | { |
---|
642 | hilfspoly=hilfspoly*t(k)^tm[i,k]; |
---|
643 | } |
---|
644 | hilfspoly=1-hilfspoly; |
---|
645 | if (chi<0) // ... either in the numerator ... |
---|
646 | { |
---|
647 | alexnumerator=alexnumerator * hilfspoly^-chi; |
---|
648 | } |
---|
649 | else // ... or in the denominator. |
---|
650 | { |
---|
651 | alexdenominator=alexdenominator * hilfspoly^chi; |
---|
652 | } |
---|
653 | alexfactor=hilfspoly,-chi; |
---|
654 | alexfactors[s]=alexfactor; |
---|
655 | s++; |
---|
656 | } |
---|
657 | chi=2; |
---|
658 | hilfspoly=1; |
---|
659 | } |
---|
660 | } |
---|
661 | // Calculate the Alexander Polynomial. |
---|
662 | if (ncols(tm)==1) // If we have only one branch, then the numerator has to be multiplied by 1-t. |
---|
663 | { |
---|
664 | alexnumerator=alexnumerator*(1-t(1)); |
---|
665 | alexfactor=1-t(1),1; |
---|
666 | alexfactors[size(alexfactors)+1]=alexfactor; |
---|
667 | } |
---|
668 | poly alexpoly=alexnumerator / alexdenominator; |
---|
669 | // Calculate the Zeta Function of the Monodromy Operator. |
---|
670 | poly zeta_monodromy=alexpoly; |
---|
671 | for (i=1;i<=r;i++) |
---|
672 | { |
---|
673 | zeta_monodromy=subst(zeta_monodromy,t(i),t); |
---|
674 | } |
---|
675 | export zeta_monodromy; |
---|
676 | export alexnumerator; |
---|
677 | export alexdenominator; |
---|
678 | export alexfactors; |
---|
679 | export alexpoly; |
---|
680 | list ALEX=ALEXring; |
---|
681 | return(ALEX); |
---|
682 | } |
---|
683 | example |
---|
684 | { |
---|
685 | "EXAMPLE:"; |
---|
686 | echo=2; |
---|
687 | ring r=0,(x,y),ls; |
---|
688 | poly f1=(y2-x3)^2-4x5y-x7; |
---|
689 | poly f2=y2-x3; |
---|
690 | poly f3=y3-x2; |
---|
691 | list ALEX=alexanderpolynomial(f1*f2*f3); |
---|
692 | def ALEXring=ALEX[1]; |
---|
693 | setring ALEXring; |
---|
694 | alexfactors; |
---|
695 | alexpoly; |
---|
696 | zeta_monodromy; |
---|
697 | } |
---|
698 | |
---|
699 | proc semigroup (def INPUT) |
---|
700 | "USAGE: semigroup(INPUT); INPUT poly or list |
---|
701 | ASSUME: INPUT is either a REDUCED bivariate polynomial defining a plane curve singularity, |
---|
702 | or the output of @code{hnexpansion(f[,\"ess\"])}, or the list @code{hne} in |
---|
703 | the ring created by @code{hnexpansion(f[,\"ess\"])}, or the output of |
---|
704 | @code{develop(f)} resp. of @code{extdevelop(f,n)}, or a list containing |
---|
705 | the contact matrix and a list of integer vectors with the characteristic exponents |
---|
706 | of the branches of a plane curve singularity, or an integer vector containing |
---|
707 | the characteristic exponents of an irreducible plane curve singularity. |
---|
708 | RETURN: a list with three entries. The first and the second are lists @code{v_1,...,v_s} |
---|
709 | and @code{w_1,...,w_r} respectively of integer vectors such that the semigroup |
---|
710 | of the plane curve defined by the INPUT is generated by the vectors |
---|
711 | @code{v_1,...,v_s,w_1+k*e_1,...,w_r+k*e_r}, where e_i denotes the i-th standard |
---|
712 | basis vector and k runs through all non-negative integers. The thrid entry is the conductor |
---|
713 | of the plane curve singularity. Note that r is the number of branches of the plane curve |
---|
714 | singularity and integer vectors thus have size r. |
---|
715 | NOTE: If the output is zero this means that the curve has one branch and is regular. |
---|
716 | In the reducible case the set of generators may not be minimal. |
---|
717 | @* |
---|
718 | If you are not sure whether the INPUT polynomial is reduced or not, use |
---|
719 | @code{squarefree(INPUT)} as input instead. |
---|
720 | SEE ALSO: resolutiongraph, totalmultiplicities |
---|
721 | EXAMPLE: example semigroup; shows an example |
---|
722 | " |
---|
723 | { |
---|
724 | //////////////////////////////////////////////////////////////////////////////////////// |
---|
725 | // Uses the algorithm in [CDG99] A. Campillo, F. Delgado, S.M. Gusein-Zade, On the |
---|
726 | // generators of the semigroup of a plane curve singularity, |
---|
727 | // J. London Math. Soc. (2) 60 (1999), 420-430. |
---|
728 | //////////////////////////////////////////////////////////////////////////////////////// |
---|
729 | intvec conductor; // conductor of the singularity |
---|
730 | list charexp; // characteristic exponents of the singularity |
---|
731 | int i,j; |
---|
732 | ///////////////////////////////////////////////////////////////////////////////// |
---|
733 | // Decide, which kind of input we have, and define the contact matrix |
---|
734 | ///////////////////////////////////////////////////////////////////////////////// |
---|
735 | // Input: a polynomial defining a plane curve singularity. |
---|
736 | ////////////////////////////////////////////////////////////////////////////// |
---|
737 | if (typeof(INPUT)=="poly") |
---|
738 | { |
---|
739 | /* |
---|
740 | poly FFF=squarefree(INPUT); |
---|
741 | if ( deg(FFF)!=deg(INPUT) ) |
---|
742 | { |
---|
743 | dbprint(printlevel-voice+3,"// input polynomial was not reduced"); |
---|
744 | dbprint(printlevel-voice+3,"// we continue with its reduction"); |
---|
745 | } |
---|
746 | list I@N@V=invariants(FFF); |
---|
747 | */ |
---|
748 | list I@N@V=invariants(INPUT); |
---|
749 | } |
---|
750 | else |
---|
751 | { |
---|
752 | /////////////////////////////////////////////////////////////////////////////////// |
---|
753 | // Input: the vector of characteristic exponents of an irreducible plane curve |
---|
754 | /////////////////////////////////////////////////////////////////////////////////// |
---|
755 | if (typeof(INPUT)=="intvec") |
---|
756 | { |
---|
757 | // Calculate the semigroup and the conductor directly. |
---|
758 | conductor[1]=charexp2conductor(INPUT); |
---|
759 | intvec gener=charexp2generators(INPUT); |
---|
760 | list genera; |
---|
761 | for (i=1;i<=size(gener);i++) |
---|
762 | { |
---|
763 | genera[i]=gener[i]; |
---|
764 | } |
---|
765 | return(list(genera,list(),conductor)); |
---|
766 | } |
---|
767 | else |
---|
768 | { |
---|
769 | if (typeof(INPUT)=="list") |
---|
770 | { |
---|
771 | ///////////////////////////////////////////////////////////////////////////////// |
---|
772 | // Input: intersection-matrix and characteristic exponents. |
---|
773 | ////////////////////////////////////////////////////////////////////////////// |
---|
774 | if (typeof(INPUT[1])=="intmat") |
---|
775 | { |
---|
776 | intmat contact=INPUT[1]; |
---|
777 | charexp=INPUT[2]; |
---|
778 | int n=ncols(contact); // to know how many branches we have |
---|
779 | if (n==1) // Only one branch! |
---|
780 | { |
---|
781 | return(semigroup(charexp[1])); |
---|
782 | } |
---|
783 | intmat intersecmult=charexp2inter(contact,charexp); |
---|
784 | for(i=1;i<=ncols(contact);i++) |
---|
785 | { |
---|
786 | conductor[i]=charexp2conductor(charexp[i]);//list with the characteristic exponents |
---|
787 | } |
---|
788 | } |
---|
789 | else |
---|
790 | { |
---|
791 | ///////////////////////////////////////////////////////////////////////////////// |
---|
792 | // Input: output of hnexpansion or hne in the ring created by hnexpansion |
---|
793 | ////////////////////////////////////////////////////////////////////////////// |
---|
794 | if ((typeof(INPUT[1])=="ring") or (typeof(INPUT[1])=="list")) |
---|
795 | { |
---|
796 | list I@N@V=invariants(INPUT); |
---|
797 | } |
---|
798 | else |
---|
799 | { |
---|
800 | //////////////////////////////////////////////////////////////////////////////////// |
---|
801 | // Input: output of develop or extdevelop -- irreducible plane curve singularity |
---|
802 | //////////////////////////////////////////////////////////////////////////////////// |
---|
803 | if (typeof(INPUT[1])=="matrix") |
---|
804 | { |
---|
805 | return(semigroup(invariants(INPUT)[1])); |
---|
806 | } |
---|
807 | else |
---|
808 | { |
---|
809 | ERROR("The input is invalid!"); |
---|
810 | } |
---|
811 | } |
---|
812 | } |
---|
813 | } |
---|
814 | else |
---|
815 | { |
---|
816 | ERROR("The input is invalid!"); |
---|
817 | } |
---|
818 | } |
---|
819 | } |
---|
820 | /////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
821 | // If the input was a poly or an HN-Expansion, then calculate the contact matrix and char.exponents. |
---|
822 | /////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
823 | if (defined(I@N@V)) |
---|
824 | { |
---|
825 | int n =size(I@N@V)-1;// number of branches |
---|
826 | // If the singularity is irreducible, then we calculate the semigroup directly. |
---|
827 | if (n==1) |
---|
828 | { |
---|
829 | return(semigroup(I@N@V[1][1])); |
---|
830 | } |
---|
831 | // If the singularity is not irreducible, then we go on. |
---|
832 | intmat contact=I@N@V[size(I@N@V)][1]; // contact matrix |
---|
833 | intmat intersecmult=I@N@V[size(I@N@V)][2]; // intersection multiplicities |
---|
834 | for(i=1;i<=n;i++) |
---|
835 | { |
---|
836 | conductor[i]=I@N@V[i][5]; |
---|
837 | charexp[i]=I@N@V[i][1]; |
---|
838 | } |
---|
839 | } |
---|
840 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
841 | // If we have come so far, the curve is reducible! |
---|
842 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
843 | // Compute the conductor of the curve. |
---|
844 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
845 | for(i=1;i<=size(conductor);i++) |
---|
846 | { |
---|
847 | for(j=1;j<=size(conductor);j++) |
---|
848 | { |
---|
849 | conductor[i]=conductor[i]+intersecmult[i,j]; |
---|
850 | } |
---|
851 | } |
---|
852 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
853 | /// The total multiplicity vectors of the exceptional divisors generate the semigroup, |
---|
854 | /// however, this system of generators is not minimal. Theorem 2 in [CDG99] leads |
---|
855 | /// to the following algorithm for minimizing the set of generators. |
---|
856 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
857 | list resgr_totmult=totalmultiplicities(list(contact,charexp)); // resolution graph and tot. mult. |
---|
858 | intmat resgr=resgr_totmult[1]; // resolution graph |
---|
859 | intmat totmult=resgr_totmult[2]; // total multiplicities |
---|
860 | list conpts; // ith entry = points to which i is connected and their weights |
---|
861 | intvec series; // ith entry = row in totmult which corresponds to w_i |
---|
862 | intvec deadarc; // Star point and the point connected to the star point of a dead arc. |
---|
863 | list deadarcs; // Saves deadarc for all dead arcs. |
---|
864 | int stop,ctp,ctpstop; |
---|
865 | list ordinary_generators, series_generators; // the v_i and the w_i from the description |
---|
866 | // Find for each branch of the curve all the points in the graph to which it is connected |
---|
867 | // and save the numbers of the corresponding rows in the graph as well as the weight of the point. |
---|
868 | for (i=1;i<=ncols(resgr);i++) |
---|
869 | { |
---|
870 | conpts[i]=find_connection_points(resgr,i); |
---|
871 | } |
---|
872 | ////////////////////////////////////////////////////////////////////////////////// |
---|
873 | // Check for each point in the graph, whether it contributes to the semigroup. |
---|
874 | // If it does not contribute to the semigroup, then set the weight of the point |
---|
875 | // to zero, i.e. resgr[i,i]=0. Note that the vertex 1 always contributes! |
---|
876 | ////////////////////////////////////////////////////////////////////////////////// |
---|
877 | // Find first all the points corresponding to the branches and the end points of dead arcs. |
---|
878 | // The points to which the branch k is connected contributes to w_k. The end points of |
---|
879 | // dead arcs contribute to the v_i, while the remaining points of the dead arcs are to be removed. |
---|
880 | j=1; // Counter for dead arcs - the star points of dead arcs have to be saved for future use. |
---|
881 | for (i=2;i<=ncols(resgr);i++) |
---|
882 | { |
---|
883 | // The end points of dead arcs and the end points corresponding to the branches are |
---|
884 | // recognized by the fact that they are only connected to one other point. |
---|
885 | if (size(conpts[i][1])==1) // row i in the graph corresponds to an end point |
---|
886 | { |
---|
887 | // For an end point corresponding to a branch k the last point E_alpha(k), to which it is |
---|
888 | // connected, contributes to the generator w_k. |
---|
889 | if (resgr[i,i]<0) |
---|
890 | { |
---|
891 | series[-resgr[i,i]]=conpts[i][1][1]; // Find E_alpha(k), where k=-resgr[i,i] |
---|
892 | ctp=conpts[i][1][1]; |
---|
893 | resgr[ctp,ctp]=0; // So E_alph(k) does not contribute to the v_i. |
---|
894 | } |
---|
895 | // For an end point of a dead arc the end point contributes, but all the other points |
---|
896 | // of the dead arc - including the star point = separation pt of the dead arc - do not contribute. |
---|
897 | if (resgr[i,i]>0) |
---|
898 | { |
---|
899 | stop=0; // Stop checks whether the star point of the dead arc has been reached. |
---|
900 | ctp=conpts[i][1][1]; // The point to which the end point is connected. |
---|
901 | // Set the weights of all the other points in the dead arc to zero. |
---|
902 | while ((stop==0) and (ctp!=1)) // If the star point or the vertex 1 has been reached, stop. |
---|
903 | { |
---|
904 | deadarc[2]=i; // i might be the point connected to the star point of the dead arc. |
---|
905 | resgr[ctp,ctp]=0; |
---|
906 | if (size(conpts[ctp][1])==2) // If ctp was an ordinary vertex, we go on. |
---|
907 | { |
---|
908 | deadarc[2]=ctp; // ctp might be the point connectd to the star point of the dead arc. |
---|
909 | ctp=conpts[ctp][1][2]; // This is the second point (of higher weight) to which ctp is connected. |
---|
910 | } |
---|
911 | else // If ctp is the star point, we stop. |
---|
912 | { |
---|
913 | deadarc[1]=ctp; // Star point of this dead arc. |
---|
914 | deadarcs[j]=deadarc; // Save the j-th dead arc. |
---|
915 | j++; |
---|
916 | stop=1; |
---|
917 | } |
---|
918 | } |
---|
919 | } |
---|
920 | } |
---|
921 | } |
---|
922 | ////////////////////////////////////////////////////////////////////////////////// |
---|
923 | // Points (!=1) which are on the geodesics of every branch don't contribute. |
---|
924 | // (The geodesic of a branch is the shortest connection of the strict transform to 1 in the graph.) |
---|
925 | stop=0; // Checks whether the points on all geodesics have been found. |
---|
926 | ctp=1; // Point (=row in resgr) to be considered. |
---|
927 | int prev_ctp; // Previous point in the graph, to which ctp was connected. |
---|
928 | int dead_arc_ctp; // If ctp is the star pt of a dead arc, this is the connection pt of ctp in the d.a. |
---|
929 | int dastarptcheck; // Checks whether a point is a star point of a dead arc. |
---|
930 | if (size(conpts[1][1])>=2) // The graphs separate already at point 1. |
---|
931 | { |
---|
932 | stop=1; |
---|
933 | } |
---|
934 | else // The graphs do not separate at point 1. |
---|
935 | { |
---|
936 | prev_ctp=1; |
---|
937 | ctp=conpts[1][1][1]; // Next point to be considered. |
---|
938 | } |
---|
939 | // Pass on along the graph until we reach the first point where some branches separate. |
---|
940 | while (stop==0) |
---|
941 | { |
---|
942 | if (size(conpts[ctp][1])==2) // Point ctp is connected to 2 other points, hence is a normal vertex. |
---|
943 | { |
---|
944 | resgr[ctp,ctp]=0; // Point ctp is a normal vertex. |
---|
945 | prev_ctp=ctp; // Save the position of ctp for future use. |
---|
946 | ctp=conpts[ctp][1][2]; // Next point to which ctp is connected. |
---|
947 | } |
---|
948 | if (size(conpts[ctp][1])>3) // More than three points are connected to ctp. |
---|
949 | { |
---|
950 | resgr[ctp,ctp]=0; |
---|
951 | stop=1; // The graphs separate at point ctp. |
---|
952 | } |
---|
953 | if (size(conpts[ctp][1])==3) // At ctp a dead arc might depart or some branch(es)! |
---|
954 | { // If a dead arc departs, then the branches stay together. |
---|
955 | resgr[ctp,ctp]=0; |
---|
956 | // Check if a dead arc departs at point ctp (i.e. if ctp is the star point of a dead arc), |
---|
957 | // then the branches do not separate at ctp. |
---|
958 | dastarptcheck=0; |
---|
959 | i=1; |
---|
960 | while ((i<=size(deadarcs)) and (dastarptcheck==0)) |
---|
961 | { |
---|
962 | if (ctp==deadarcs[i][1]) // ctp is the star point of a dead arc. |
---|
963 | { |
---|
964 | dastarptcheck=1; |
---|
965 | dead_arc_ctp=deadarcs[i][2]; // The point in the dead arc to which ctp is connected. |
---|
966 | } |
---|
967 | i++; |
---|
968 | } |
---|
969 | if (dastarptcheck==0) // ctp is not the star point of a dead arc, hence the graphs separate at ctp. |
---|
970 | { |
---|
971 | stop=1; |
---|
972 | } |
---|
973 | else |
---|
974 | { |
---|
975 | // Set ctp to the point connected to ctp which is not in the dead arc and is not prev_ctp. |
---|
976 | i=1; |
---|
977 | ctpstop=0; |
---|
978 | while ((i<=3) and (ctpstop==0)) |
---|
979 | { |
---|
980 | if ((conpts[ctp][1][i]!=prev_ctp) and (conpts[ctp][1][i]!=dead_arc_ctp)) |
---|
981 | { |
---|
982 | prev_ctp=ctp; |
---|
983 | ctp=conpts[ctp][1][i]; |
---|
984 | ctpstop=1; |
---|
985 | } |
---|
986 | i++; |
---|
987 | } |
---|
988 | } |
---|
989 | } |
---|
990 | } |
---|
991 | ///////////////////////////////////////////////////////////////////////////////////////////// |
---|
992 | // Collect the generators v_j by checking which points in the graph still have |
---|
993 | // a positive weight! These points contribute their total multiplicity vector as |
---|
994 | // generator v_j. |
---|
995 | j=1; |
---|
996 | for (i=1;i<=ncols(resgr);i++) |
---|
997 | { |
---|
998 | if (resgr[i,i]>0) |
---|
999 | { |
---|
1000 | ordinary_generators[j]=intvec(totmult[i,1..ncols(totmult)]); |
---|
1001 | j++; |
---|
1002 | } |
---|
1003 | } |
---|
1004 | // The "exceptional" generators w_i, for which we have to include w_i+ke_i (for all k) |
---|
1005 | // are the total multiplicity vectors of the points saved in series. |
---|
1006 | for (i=1;i<=ncols(totmult);i++) |
---|
1007 | { |
---|
1008 | series_generators[i]=intvec(totmult[series[i],1..ncols(totmult)]); |
---|
1009 | } |
---|
1010 | return(list(ordinary_generators,series_generators,conductor)); |
---|
1011 | } |
---|
1012 | example |
---|
1013 | { |
---|
1014 | "EXAMPLE:"; |
---|
1015 | echo=2; |
---|
1016 | ring r=0,(x,y),ls; |
---|
1017 | // Irreducible Case |
---|
1018 | semigroup((x2-y3)^2-4x5y-x7); |
---|
1019 | // In the irreducible case, invariants() also calculates a minimal set of |
---|
1020 | // generators of the semigroup. |
---|
1021 | invariants((x2-y3)^2-4x5y-x7)[1][2]; |
---|
1022 | // Reducible Case |
---|
1023 | poly f=(y2-x3)*(y2+x3)*(y4-2x3y2-4x5y+x6-x7); |
---|
1024 | semigroup(f); |
---|
1025 | } |
---|
1026 | |
---|
1027 | |
---|
1028 | |
---|
1029 | proc charexp2generators (intvec charexp) |
---|
1030 | "USAGE: charexp2generators(v), v intvec |
---|
1031 | ASSUME: v contains the characteristic exponents of an irreducible plane |
---|
1032 | curve singularity |
---|
1033 | RETURN: intvec, the minimal set of generators of the semigroup of the plane curve singularity |
---|
1034 | SEE ALSO: invariants, resolutiongraph, totalmultiplicities, alexanderpolynomial |
---|
1035 | KEYWORDS: generators; semigroup; characteristic exponents; topological invariants |
---|
1036 | EXAMPLE: example charexp2generators; shows an example" |
---|
1037 | { |
---|
1038 | int end=size(charexp); |
---|
1039 | // If the singularity is smooth! |
---|
1040 | if (end==1) |
---|
1041 | { |
---|
1042 | return(1); |
---|
1043 | } |
---|
1044 | int i,j; |
---|
1045 | intvec gener; |
---|
1046 | intvec GGT; |
---|
1047 | for (i=1;i<=end;i++) |
---|
1048 | { |
---|
1049 | // Calculate the sequence of gcd's of the characteristic exponents. |
---|
1050 | if (i==1) |
---|
1051 | { |
---|
1052 | GGT[1]=charexp[1]; |
---|
1053 | } |
---|
1054 | else |
---|
1055 | { |
---|
1056 | GGT[i]=gcd(GGT[i-1],charexp[i]); |
---|
1057 | } |
---|
1058 | // Calculate the generators. |
---|
1059 | gener[i]=charexp[i]; |
---|
1060 | for (j=2;j<=i-1;j++) |
---|
1061 | { |
---|
1062 | gener[i]=gener[i]+((GGT[j-1]-GGT[j])/GGT[i-1])*charexp[j]; |
---|
1063 | } |
---|
1064 | } |
---|
1065 | return(gener); |
---|
1066 | } |
---|
1067 | example |
---|
1068 | { |
---|
1069 | "EXAMPLE:"; |
---|
1070 | echo=2; |
---|
1071 | charexp2generators(intvec(28,64,66,77)); |
---|
1072 | } |
---|
1073 | |
---|
1074 | |
---|
1075 | proc charexp2multseq (intvec charexp) |
---|
1076 | "USAGE: charexp2multseq(v), v intvec |
---|
1077 | ASSUME: v contains the characteristic exponents of an irreducible plane |
---|
1078 | curve singularity |
---|
1079 | RETURN: intvec, the multiplicity sequence of the plane curve singularity |
---|
1080 | NOTE: If the curve singularity is smooth, then the multiplicity sequence is empty. |
---|
1081 | This is expressed by returning zero. |
---|
1082 | SEE ALSO: invariants, resolutiongraph, totalmultiplicities, alexanderpolynomial |
---|
1083 | KEYWORDS: characteristic exponents; multiplicity sequence; topological invariants |
---|
1084 | EXAMPLE: example charexp2multseq; shows an example" |
---|
1085 | { |
---|
1086 | int end=size(charexp); |
---|
1087 | // If the singularity is smooth! |
---|
1088 | if (end==1) |
---|
1089 | { |
---|
1090 | return(1); // ERROR: Should be 0, but for the time being, Singular returns 1. |
---|
1091 | } |
---|
1092 | intvec multseq=euclidseq(charexp[2],charexp[1]); |
---|
1093 | for (int i=3;i<=end;i++) |
---|
1094 | { |
---|
1095 | multseq=multseq,euclidseq(charexp[i]-charexp[i-1],multseq[size(multseq)]); |
---|
1096 | } |
---|
1097 | return(multseq); |
---|
1098 | } |
---|
1099 | example |
---|
1100 | { |
---|
1101 | "EXAMPLE:"; |
---|
1102 | echo=2; |
---|
1103 | charexp2multseq(intvec(28,64,66,77)); |
---|
1104 | } |
---|
1105 | |
---|
1106 | proc multseq2charexp(def v) // Procedure written by Fernando. |
---|
1107 | "USAGE: multseq2charexp(v), v intvec |
---|
1108 | ASSUME: The input is an intvec, which contains the mutiplicity sequence |
---|
1109 | of an irreducible plane curve singularity . |
---|
1110 | RETURN: An intvec, which contains the sequence of characteristic |
---|
1111 | exponents of the irreducible plane curve singularity defined by v. |
---|
1112 | EXAMPLE: example multseq2charexp; shows an example. |
---|
1113 | " |
---|
1114 | { |
---|
1115 | /////// Preamble which reduces the input of an intvec to ///////////// |
---|
1116 | /////// the originally assumed input of a list of intvecs ///////////// |
---|
1117 | /////// and tests the input. ///////////// |
---|
1118 | if (typeof(v)=="intvec") |
---|
1119 | { |
---|
1120 | list RESULT=multseq2charexp(list(v)); |
---|
1121 | return(RESULT[1]); |
---|
1122 | } |
---|
1123 | if (typeof(v)!="list") |
---|
1124 | { |
---|
1125 | ERROR("Invalid Input!"); |
---|
1126 | } |
---|
1127 | if (typeof(v)=="list") |
---|
1128 | { |
---|
1129 | int TESTV; |
---|
1130 | for (int III=1;III<=size(v);III++) |
---|
1131 | { |
---|
1132 | if (typeof(v[III])!="intvec") |
---|
1133 | { |
---|
1134 | TESTV=1; |
---|
1135 | } |
---|
1136 | } |
---|
1137 | if (TESTV==1) |
---|
1138 | { |
---|
1139 | ERROR("Invalid Input!"); |
---|
1140 | } |
---|
1141 | } |
---|
1142 | /////////////////////////////////////////////////////////// |
---|
1143 | list L=v; |
---|
1144 | int n =size(L); |
---|
1145 | // /////////////////////////////////////////////////////// |
---|
1146 | // we look the size of each vector |
---|
1147 | intvec mm; |
---|
1148 | for(int j=1;j<=n;j++) |
---|
1149 | { |
---|
1150 | mm[j]=size(L[j]); |
---|
1151 | } |
---|
1152 | // /////////////////////////////////////////////////////// |
---|
1153 | // we define some variables |
---|
1154 | list LL; |
---|
1155 | int temp, temp1,temp2; |
---|
1156 | int ind,r,l,boolean; |
---|
1157 | int old,new; |
---|
1158 | int contador; |
---|
1159 | list EUCLI,EUCLI1; |
---|
1160 | intvec exponent,exponentes1; |
---|
1161 | int new1,old1; |
---|
1162 | int contador1; |
---|
1163 | int a,b,f; |
---|
1164 | //with the for we round each branch. |
---|
1165 | for(int k=1;k<=n;k++) |
---|
1166 | { |
---|
1167 | l=1; |
---|
1168 | old=L[k][l]; |
---|
1169 | //if the vertor has more than one element |
---|
1170 | if(mm[k]<>1) |
---|
1171 | { |
---|
1172 | // /////////////////////////////////////////////////////////////////////////////// |
---|
1173 | // the first step is special because is easy to know the two first characteristic exponents |
---|
1174 | new=L[k][l+1]; |
---|
1175 | contador=1; |
---|
1176 | while(old==new)//we check how many consecutives elements are equal, starting in the first. |
---|
1177 | { |
---|
1178 | contador=contador+1; |
---|
1179 | old=new; |
---|
1180 | new=L[k][contador+1]; |
---|
1181 | } |
---|
1182 | exponent=L[k][l],contador*L[k][l]+L[k][l+contador];// those are the first two characteristic exponents. |
---|
1183 | LL[k]=exponent;// we insert them in the final matrix |
---|
1184 | EUCLI=euclides(LL[k][2],LL[k][1]);// compute the euclides algorithm for the two first characteristic exponents. |
---|
1185 | temp=size(EUCLI[1]); |
---|
1186 | // measure how many elements of the multiplicity sequence belong to the first euclidean algorithm. |
---|
1187 | for(ind=1;ind<=temp;ind=ind+1) |
---|
1188 | { |
---|
1189 | l=l+EUCLI[1][ind]; |
---|
1190 | } |
---|
1191 | l=l-1;//this is the number we are looking for. |
---|
1192 | /////////////////////////////////////////////////////////////// |
---|
1193 | r=1; |
---|
1194 | //repeat the same process until the end of the multiplicity sequence. |
---|
1195 | if(mm[k]-1>l) |
---|
1196 | { |
---|
1197 | while( l<mm[k]-1) |
---|
1198 | { |
---|
1199 | r=r+1; |
---|
1200 | old1=L[k][l]; |
---|
1201 | new1=L[k][l+1]; |
---|
1202 | contador1=0; |
---|
1203 | boolean=1; |
---|
1204 | if(old1==new1) |
---|
1205 | { |
---|
1206 | while(old1==new1 and boolean==1) |
---|
1207 | { |
---|
1208 | contador1=contador1+1; |
---|
1209 | old1=new1; |
---|
1210 | new1=L[k][l+contador1+1]; |
---|
1211 | if(size(L[k])<=l+contador1+1) |
---|
1212 | { |
---|
1213 | boolean =0; |
---|
1214 | } |
---|
1215 | } |
---|
1216 | } |
---|
1217 | temp1=size(LL[k]); |
---|
1218 | exponentes1=LL[k],LL[k][temp1]+(contador1*L[k][l])+L[k][contador1+l+1]; |
---|
1219 | LL[k]=exponentes1; |
---|
1220 | EUCLI1=euclides(LL[k][temp1+1]-LL[k][temp1],L[k][l]); |
---|
1221 | temp2=size(EUCLI1[1]); |
---|
1222 | for(ind=1;ind<=temp2;ind=ind+1) |
---|
1223 | { |
---|
1224 | l=l+EUCLI1[1][ind]; |
---|
1225 | } |
---|
1226 | } |
---|
1227 | } |
---|
1228 | } |
---|
1229 | // if the vector has only one element then the charexp is only 1. |
---|
1230 | else |
---|
1231 | { |
---|
1232 | LL[k]=1; |
---|
1233 | } |
---|
1234 | } |
---|
1235 | return(LL); |
---|
1236 | } |
---|
1237 | example |
---|
1238 | { |
---|
1239 | "EXAMPLE:";echo=2; |
---|
1240 | intvec v=2,1,1; |
---|
1241 | multseq2charexp(v); |
---|
1242 | intvec v1=4,2,2,1,1; |
---|
1243 | multseq2charexp(v1); |
---|
1244 | } |
---|
1245 | |
---|
1246 | proc charexp2inter (intmat contact, list charexp) |
---|
1247 | "USAGE: charexp2inter(contact,charexp), contact matrix, charexp list |
---|
1248 | ASSUME: charexp contains the integer vectors of characteristic exponents |
---|
1249 | of the branches of a plane curve singularity, and contact is their |
---|
1250 | contact matrix |
---|
1251 | RETURN: intmat, the matrix intersection multiplicities of the branches |
---|
1252 | SEE ALSO: invariants, resolutiongraph, totalmultiplicities, semigroup |
---|
1253 | KEYWORDS: contact matrix; characteristic exponents; intersection multiplicity; topological invariants |
---|
1254 | EXAMPLE: example charexp2inter; shows an example" |
---|
1255 | { |
---|
1256 | int n=ncols(contact); |
---|
1257 | int i,j,k; |
---|
1258 | list multpl; |
---|
1259 | int max=0; |
---|
1260 | intvec helpvect; |
---|
1261 | intmat inters[n][n]; |
---|
1262 | // Calculate the multiplicity sequences of the branches. |
---|
1263 | for (i=1;i<=n;i++) |
---|
1264 | { |
---|
1265 | multpl[i]=charexp2multseq(charexp[i]); |
---|
1266 | // Find the maximal length of a multiplicity sequence. |
---|
1267 | if (max<size(multpl[i])) |
---|
1268 | { |
---|
1269 | max=size(multpl[i]); |
---|
1270 | } |
---|
1271 | } |
---|
1272 | // If the contact of certain branches is higher than the maximal length of the |
---|
1273 | // multiplicity sequence, then max should be the maximal contact! |
---|
1274 | helpvect=max,intvec(contact); |
---|
1275 | max=max_in_intvec(helpvect)[1]; |
---|
1276 | // Prolong them by 1s, in order to take care of higher contact. |
---|
1277 | for (i=1;i<=n;i++) |
---|
1278 | { |
---|
1279 | helpvect=multpl[i]; |
---|
1280 | for (j=size(multpl[i]+1);j<=max;j++) |
---|
1281 | { |
---|
1282 | helpvect[j]=1; |
---|
1283 | } |
---|
1284 | multpl[i]=helpvect; |
---|
1285 | } |
---|
1286 | // Calculate the intersection numbers of the branches: for two branches f_i and f_j |
---|
1287 | // this is the sum over mult_k(f_i)*mult_k(f_j), where k runs over all infinitely |
---|
1288 | // near points which f_i and f_j share, i.e. from 1 to contact[i,j]. |
---|
1289 | for (i=1;i<=n;i++) |
---|
1290 | { |
---|
1291 | for (j=i+1;j<=n;j++) |
---|
1292 | { |
---|
1293 | for (k=1;k<=contact[i,j];k++) |
---|
1294 | { |
---|
1295 | inters[i,j]=inters[i,j]+multpl[i][k]*multpl[j][k]; |
---|
1296 | } |
---|
1297 | inters[j,i]=inters[i,j]; |
---|
1298 | } |
---|
1299 | } |
---|
1300 | return(inters); |
---|
1301 | } |
---|
1302 | example |
---|
1303 | { |
---|
1304 | "EXAMPLE:";echo=2; |
---|
1305 | ring r=0,(x,y),ds; |
---|
1306 | list INV=invariants((x2-y3)*(x3-y2)*((x2-y3)^2-4x5y-x7)); |
---|
1307 | intmat contact=INV[4][1]; |
---|
1308 | list charexp=INV[1][1],INV[2][1],INV[3][1]; |
---|
1309 | // The intersection matrix is INV[4][2]. |
---|
1310 | print(INV[4][2]); |
---|
1311 | // And it is calulated as ... |
---|
1312 | print(charexp2inter(contact,charexp)); |
---|
1313 | } |
---|
1314 | |
---|
1315 | |
---|
1316 | proc charexp2conductor(intvec B) // Procedure written by Fernando |
---|
1317 | "USAGE: charexp2conductor(v), v intvec |
---|
1318 | ASSUME: v contains the characteristic exponents of an irreducible plane |
---|
1319 | curve singularity |
---|
1320 | RETURN: int, the conductor of the plane curve singularity |
---|
1321 | NOTE: If the curve singularity is smooth, the conductor is zero. |
---|
1322 | SEE ALSO: invariants, resolutiongraph, totalmultiplicities, semigroup |
---|
1323 | KEYWORDS: conductor; characteristic exponents; multiplicity sequence; topological invariants |
---|
1324 | EXAMPLE: example charexp2conductor; shows an example" |
---|
1325 | { |
---|
1326 | intvec E; |
---|
1327 | int i,conductor; |
---|
1328 | E[1]=B[1]; |
---|
1329 | for(i=2;i<=size(B);i++) |
---|
1330 | { |
---|
1331 | E[i]=gcd(E[i-1],B[i]); |
---|
1332 | } |
---|
1333 | conductor=1-E[1]; |
---|
1334 | for(i=2;i<=size(B);i++) |
---|
1335 | { |
---|
1336 | conductor=conductor+(B[i]*(E[i-1]-E[i])); |
---|
1337 | } |
---|
1338 | return(conductor); |
---|
1339 | } |
---|
1340 | example |
---|
1341 | { |
---|
1342 | "EXAMPLE:"; |
---|
1343 | echo=2; |
---|
1344 | charexp2conductor(intvec(2,3)); // A1-Singularity |
---|
1345 | charexp2conductor(intvec(28,64,66,77)); |
---|
1346 | } |
---|
1347 | |
---|
1348 | |
---|
1349 | proc charexp2poly(intvec v, vector a) // Procedure written by Fernando. |
---|
1350 | "USAGE: charexp2poly(v,a); intvec v, vector a. |
---|
1351 | ASSUME: v an intvec containing the characterictic exponents of an irreducible plane curve singularity. |
---|
1352 | a a vector containing the coefficients of a parametrization given by x(t)=x^v[1], |
---|
1353 | y(t)=a(1)t^v[2]+...+a[n-1]t^v[n], i.e. the entries of a are of type number. |
---|
1354 | RETURN: A polynomial f in the first two variables of the basering, such that f defines an |
---|
1355 | irreducible plane curve singularity with characteristic exponents v. |
---|
1356 | NOTE: The entries in a should be of type number and the vector v should be the sequence of |
---|
1357 | characteristic exponents of an irreducible plane curve singularity in order to |
---|
1358 | get a sensible result, |
---|
1359 | SEE ALSO: charexp2multseq, multseq2charexp. |
---|
1360 | EXAMPLE: example charexp2poly; shows an example |
---|
1361 | " |
---|
1362 | { |
---|
1363 | int n=size(v); |
---|
1364 | vector expo; |
---|
1365 | int i,j,s; |
---|
1366 | for(i=1;i<=v[1];i++) |
---|
1367 | { |
---|
1368 | expo[i]=0;//initialize to 0. |
---|
1369 | } |
---|
1370 | for(i=2;i<=n;i++) |
---|
1371 | { |
---|
1372 | s=v[i] mod v[1];//calculate the position. |
---|
1373 | expo=expo-a[i-1]*var(1)^((v[i]-s)/v[1])*gen(s+1);//save in expo -var(1) to the power the corresponding |
---|
1374 | //but only in the right positions |
---|
1375 | } |
---|
1376 | matrix M[v[1]][v[1]]; |
---|
1377 | //construct the matrix that generates the polynomial |
---|
1378 | for(i=1;i<=v[1];i++) |
---|
1379 | { |
---|
1380 | M[i,i]=var(2)+expo[1];//The diagonal |
---|
1381 | for(j=1;j<=v[1]-i;j++) |
---|
1382 | { |
---|
1383 | M[j,j+i]=expo[i+1];//over diagonal |
---|
1384 | } |
---|
1385 | for(j=1;j<=v[1]-i;j++) |
---|
1386 | { |
---|
1387 | M[j+i,j]=var(1)*expo[1+v[1]-i];//under diagonal |
---|
1388 | } |
---|
1389 | } |
---|
1390 | //the poynomial is the determinant of the matrix |
---|
1391 | poly irredpoly=det(M); |
---|
1392 | return(irredpoly) |
---|
1393 | } |
---|
1394 | example |
---|
1395 | { |
---|
1396 | "EXAMPLE:";echo=2; |
---|
1397 | ring r=0,(x,y),dp; |
---|
1398 | intvec v=8,12,14,17; |
---|
1399 | vector a=[1,1,1]; |
---|
1400 | poly f=charexp2poly(v,a); |
---|
1401 | f; |
---|
1402 | invariants(f)[1][1]; // The characteristic exponents of f. |
---|
1403 | } |
---|
1404 | |
---|
1405 | proc tau_es (def INPUT, list #) |
---|
1406 | "USAGE: tau_es(INPUT); INPUT poly or list |
---|
1407 | ASSUME: INPUT is either a REDUCED bivariate polynomial defining a plane curve singularity, |
---|
1408 | or the output of @code{hnexpansion(f[,\"ess\"])}, or the list @code{hne} in |
---|
1409 | the ring created by @code{hnexpansion(f[,\"ess\"])}, or the output of |
---|
1410 | @code{develop(f)} resp. of @code{extdevelop(f,n)}, or a list containing |
---|
1411 | the contact matrix and a list of integer vectors with the characteristic exponents |
---|
1412 | of the branches of a plane curve singularity, or an integer vector containing the |
---|
1413 | characteristic exponents of an irreducible plane curve singularity. |
---|
1414 | RETURN: int, the equisingular Tjurina number of f, i. e. the codimension of the mu-constant |
---|
1415 | stratum in the semiuniversal deformation of f, where mu is the Milnor number of f. |
---|
1416 | NOTE: The equisingular Tjurina number is calculated with the aid of a Hamburger-Noether |
---|
1417 | expansion, which is the hard part of the calculation. |
---|
1418 | In case the Hamburger-Noether expansion of the curve f is needed |
---|
1419 | for other purposes as well it is better to calculate this first |
---|
1420 | with the aid of @code{hnexpansion} and use it as input instead of |
---|
1421 | the polynomial itself. |
---|
1422 | @* |
---|
1423 | If you are not sure whether the INPUT polynomial is reduced or not, use |
---|
1424 | @code{squarefree(INPUT)} as input instead. |
---|
1425 | SEE ALSO: develop, hnexpansion, totalmultiplicities, equising.lib |
---|
1426 | EXAMPLE: example tau_es; shows an example |
---|
1427 | " |
---|
1428 | { |
---|
1429 | // If the input is a weighted homogeneous polnomial, then use a direct algorithm to |
---|
1430 | // calculate the equisingular Tjurina number, by caluclating a K-basis of the |
---|
1431 | // Tjurina algebra and omitting those elements with weighted degree at least the |
---|
1432 | // weighted degree of the polynomial. -- If an additional input # is given (which is |
---|
1433 | // not just the number 1 !!!), the procedure always uses the recursive algorithm. |
---|
1434 | if ((typeof(INPUT)=="poly") and (size(#)==0)) |
---|
1435 | { |
---|
1436 | if (qhweight(INPUT)[1]!=0) |
---|
1437 | { |
---|
1438 | /* |
---|
1439 | poly f=squarefree(INPUT); |
---|
1440 | if ( deg(f)!=deg(INPUT) ) |
---|
1441 | { |
---|
1442 | dbprint(printlevel-voice+3,"// input polynomial was not reduced"); |
---|
1443 | dbprint(printlevel-voice+3,"// we continue with its reduction"); |
---|
1444 | } |
---|
1445 | return(tau_es_qh(f)); |
---|
1446 | */ |
---|
1447 | return(tau_es_qh(INPUT)); |
---|
1448 | } |
---|
1449 | } |
---|
1450 | // Else apply the recursive algorithm from Eugenii Shustin, On Manifolds of Singular |
---|
1451 | // Algebraic Curves, Selecta Math. Sov. Vol 10, No. 1 (1991), p. 31. |
---|
1452 | int i,j,k; // Laufvariablen |
---|
1453 | int tau,tau_i; // Variable for es-Tjurina no., multiplicitiy, and es-Tjurina at blow ups |
---|
1454 | intmat contact; // contact matrix |
---|
1455 | list graphs, multipl; // resolution graphs and multiplicity sequences |
---|
1456 | list graphs_ed, multipl_ed; // res.graphs and mult.seq. of one curve on 1st ex. div. |
---|
1457 | intmat contact_ed; // contact matrix of one curve on 1st exceptional divisor |
---|
1458 | int d_i; // correction term |
---|
1459 | intvec connections; |
---|
1460 | intmat graph; |
---|
1461 | intvec multseq; |
---|
1462 | int s,e; |
---|
1463 | int multi; // multiplicity of the curve defined by the input. |
---|
1464 | ///////////////////////////////////////////////////////// |
---|
1465 | // Check what type the input is, and act accordingly. |
---|
1466 | if (typeof(INPUT)=="list") |
---|
1467 | { |
---|
1468 | if (size(INPUT)==3) |
---|
1469 | { |
---|
1470 | // If the INPUT is the output of totalmultiplicities(INPUT,"tau") |
---|
1471 | if ((typeof(INPUT[1])=="list") and (typeof(INPUT[2])=="list") and (typeof(INPUT[3])=="intmat")) |
---|
1472 | { |
---|
1473 | graphs=INPUT[1]; |
---|
1474 | multipl=INPUT[2]; |
---|
1475 | contact=INPUT[3]; |
---|
1476 | } |
---|
1477 | // Otherwise call the procedure with the output of totalmultiplicities(INPUT,"tau"). |
---|
1478 | else |
---|
1479 | { |
---|
1480 | return(tau_es(totalmultiplicities(INPUT,"tau"))); |
---|
1481 | } |
---|
1482 | } |
---|
1483 | else |
---|
1484 | { |
---|
1485 | return(tau_es(totalmultiplicities(INPUT,"tau"))); |
---|
1486 | } |
---|
1487 | } |
---|
1488 | // Otherwise call the procedure with the output of totalmultiplicities(INPUT,"tau"). |
---|
1489 | else |
---|
1490 | { |
---|
1491 | return(tau_es(totalmultiplicities(INPUT,"tau"))); |
---|
1492 | } |
---|
1493 | /// End of checking the input |
---|
1494 | /////////////////////////////////////////////////////// |
---|
1495 | // If the singularity is smooth, the equisingular Tjurina number is zero. |
---|
1496 | if ((ncols(contact)==1) and (multipl[1][1]<=1)) |
---|
1497 | { |
---|
1498 | return(0); |
---|
1499 | } |
---|
1500 | // Otherwise calculate the multiplicity of the singularity. |
---|
1501 | for (i=1;i<=size(multipl);i++) |
---|
1502 | { |
---|
1503 | multi=multi+multipl[i][1]; |
---|
1504 | } |
---|
1505 | // Find the branches which stay together after blowing up once, and group them together. |
---|
1506 | k=0; |
---|
1507 | intvec curves_on_first_ex_div=k; // If this is 0=i_0,i_1,...i_s=ncols(contact), then the branches |
---|
1508 | while (k<ncols(contact)) // (i_j)+1,...,i_(j+1) stay together after the first blowing up, |
---|
1509 | { // and s is the number of infinitely near points of first order. |
---|
1510 | k=find_last_non_one(intvec(contact[k+1,1..ncols(contact)]),k+2); |
---|
1511 | curves_on_first_ex_div=curves_on_first_ex_div,k; |
---|
1512 | } |
---|
1513 | // And then calculate the equisingular Tjurina numbers in the infinitely near points of |
---|
1514 | // first order plus some correction term (= the number of blow ups needed to separate the |
---|
1515 | // strict transform from the exceptional divisor in the corresponding infinitely near point) |
---|
1516 | // and add them. |
---|
1517 | for (i=1;i<=size(curves_on_first_ex_div)-1;i++) |
---|
1518 | { |
---|
1519 | s=curves_on_first_ex_div[i]+1; |
---|
1520 | e=curves_on_first_ex_div[i+1]; |
---|
1521 | graphs_ed=list(graphs[s..e]); |
---|
1522 | multipl_ed=list(multipl[s..e]); |
---|
1523 | contact_ed=intmat_minus_one(intmat(intvec(contact[s..e,s..e]),e-s+1,e-s+1)); |
---|
1524 | connections=0; |
---|
1525 | // Adjust the graphs and multiplicity sequences by cutting the first level. |
---|
1526 | // And find in each graph the point to which the point 1 is connected = 1 + number of |
---|
1527 | // blow ups needed to separate the strict transform from the first exceptional divisor. |
---|
1528 | for (j=1;j<=size(graphs_ed);j++) |
---|
1529 | { |
---|
1530 | // Adjust the graphs and find the connection points. |
---|
1531 | graph=graphs_ed[j]; |
---|
1532 | connections[j]=find_connection_point(intvec(graph[1,1..ncols(graph)]),1); |
---|
1533 | graph=delete_row(delete_col(graph,1),1); |
---|
1534 | graphs_ed[j]=graph; |
---|
1535 | // Adjust the multiplicity sequences. |
---|
1536 | multseq=multipl_ed[j]; |
---|
1537 | if (size(multseq)==1) // If multseq has only one entry, then their is only |
---|
1538 | { // one branch left and it is smooth! So set multipl_ed[j]=1. |
---|
1539 | multipl_ed[j]=intvec(1); |
---|
1540 | } |
---|
1541 | else // Otherwise just cut the first entry. |
---|
1542 | { |
---|
1543 | multipl_ed[j]=intvec(multseq[2..size(multseq)]); |
---|
1544 | } |
---|
1545 | } |
---|
1546 | // Calculate the equisingular Tjurina number of the strict transform at the i-th |
---|
1547 | // intersection point with the first exceptional divisor. |
---|
1548 | tau_i=tau_es(list(graphs_ed,multipl_ed,contact_ed)); |
---|
1549 | // Calculate the number no. of blow ups needed to separate the strict transform of |
---|
1550 | // the curve defined by the input from the first exceptional divisor at their i-th |
---|
1551 | // intersection point. |
---|
1552 | d_i=max_in_intvec(connections)-1; |
---|
1553 | // If d_i is negative, then there was only one branch left and it was smooth, so d_i should be 1. |
---|
1554 | if (d_i<0) |
---|
1555 | { |
---|
1556 | d_i=1; |
---|
1557 | } |
---|
1558 | // If the curve defined by graphs_ed was smooth, then the summand has to be reduced by 1. |
---|
1559 | if (tau_i==0) |
---|
1560 | { |
---|
1561 | tau_i=-1; |
---|
1562 | } |
---|
1563 | tau=tau+tau_i+d_i; |
---|
1564 | } |
---|
1565 | // The equisingular Tjurina number is then calculated by adding the following term. |
---|
1566 | tau=tau+((multi*(multi+1))/2)-2; |
---|
1567 | return(tau); |
---|
1568 | } |
---|
1569 | example |
---|
1570 | { |
---|
1571 | "EXAMPLE:"; |
---|
1572 | echo=2; |
---|
1573 | ring r=0,(x,y),ls; |
---|
1574 | poly f1=y2-x3; |
---|
1575 | poly f2=(y2-x3)^2-4x5y-x7; |
---|
1576 | poly f3=y3-x2; |
---|
1577 | tau_es(f1); |
---|
1578 | tau_es(f2); |
---|
1579 | tau_es(f1*f2*f3); |
---|
1580 | } |
---|
1581 | |
---|
1582 | |
---|
1583 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
1584 | // Static procedures. |
---|
1585 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
1586 | |
---|
1587 | static proc euclidseq(int a,int b) |
---|
1588 | "USAGE: euclidseq(a,b), a,b integers |
---|
1589 | RETURN: intvec, the divisors in the euclidean alogrithm with multiplicities |
---|
1590 | KEYWORDS: Euclidean algorithm; multiplicity sequence |
---|
1591 | NOTE : This procedure is for internal use only; it is called by charexp2multseq. |
---|
1592 | " |
---|
1593 | { |
---|
1594 | int i; |
---|
1595 | // multseq saves in each step of the Euclidean algorithm q-times the divisor b |
---|
1596 | intvec multseq; |
---|
1597 | int q=a div b; |
---|
1598 | int r=a mod b; |
---|
1599 | for (i=1;i<=q;i++) |
---|
1600 | { |
---|
1601 | multseq[i]=b; |
---|
1602 | } |
---|
1603 | int s=q; // size of multseq |
---|
1604 | a=b; |
---|
1605 | b=r; |
---|
1606 | while(r<>0) |
---|
1607 | { |
---|
1608 | q=a div b; |
---|
1609 | r=a mod b; |
---|
1610 | for (i=1;i<=q;i++) |
---|
1611 | { |
---|
1612 | multseq[s+i]=b; |
---|
1613 | } |
---|
1614 | s=s+q; // size of multseq |
---|
1615 | a=b; |
---|
1616 | b=r; |
---|
1617 | } |
---|
1618 | return(multseq); |
---|
1619 | } |
---|
1620 | |
---|
1621 | static proc puiseuxchainpart (int piij, int muij, intvec multpl, int initial_tm, int end_tm, int j) |
---|
1622 | "USAGE: puiseuxchainpart(piij,muiij,multpl,initial_tm,end_tm,j); |
---|
1623 | RETURN: list L, L[1] incidence matrix of a part of the Puiseux chain, L[2] the total |
---|
1624 | multiplicities of this part of the Puiseux chain. |
---|
1625 | NOTE: This procedure is only for internal use; it is called by puiseuxchain. |
---|
1626 | " |
---|
1627 | { |
---|
1628 | int delta=1; |
---|
1629 | if (j==1){delta=0;} // Delta measures whether j is 1 or not. |
---|
1630 | intvec totalmultiplicity; // Keeps the total multiplicities. |
---|
1631 | intmat pcp[muij][muij]; // Keeps the incidence matrix of the Puiseuxchainpart S_i,j. |
---|
1632 | // Calculate the total multiplicities and the Puiseuxchainpart S_i,j. |
---|
1633 | totalmultiplicity[1]=initial_tm+end_tm+multpl[1]; |
---|
1634 | pcp[1,1]=piij+1; |
---|
1635 | for (int k=2;k<=muij;k++) |
---|
1636 | { |
---|
1637 | pcp[k,k]=piij+k; |
---|
1638 | pcp[k-1,k]=1; |
---|
1639 | pcp[k,k-1]=1; |
---|
1640 | totalmultiplicity[k]=totalmultiplicity[k-1]+delta*initial_tm+multpl[k]; |
---|
1641 | } |
---|
1642 | list result=pcp,totalmultiplicity; |
---|
1643 | return(result); |
---|
1644 | } |
---|
1645 | |
---|
1646 | static proc puiseuxchain (int initial, intvec divseq, intvec multpl, int initial_tm) |
---|
1647 | "USAGE: puiseuxchain(initial,divseq,multpl,initial_tm); int initial, initial_tm, intvec divseq, multpl |
---|
1648 | RETURN: list L, L[1] incidence matrix of a Puiseux chain, L[2] the weight of the point to which the |
---|
1649 | previous Puiseux chain has to be connected, L[3] the sequence of total multiplicities of |
---|
1650 | the points in this Puiseux chain. |
---|
1651 | NOTE: This procedure is only for internal use; it is called by irred_resgraph_totmult. |
---|
1652 | " |
---|
1653 | { |
---|
1654 | int j,k,l,connectpoint; |
---|
1655 | intvec multpli; |
---|
1656 | int pc_tm=initial_tm; // Keeps the total multipl. of the endpoint of P_i-1. |
---|
1657 | int end_tm=0; |
---|
1658 | int start=1; |
---|
1659 | int omega=size(divseq); |
---|
1660 | // Keep the endpoints of the puiseuxchainparts (minus initial) s_i,j in divseq_sum. |
---|
1661 | intvec divseq_sum=divseq[1]; |
---|
1662 | for (j=2;j<=omega;j++) |
---|
1663 | { |
---|
1664 | divseq_sum[j]=divseq_sum[j-1]+divseq[j]; |
---|
1665 | } |
---|
1666 | // Define the connecting point of the Puiseuxchain P_i-1 with P_i. |
---|
1667 | if (divseq[1]==0) |
---|
1668 | { |
---|
1669 | // If divseq[1]=mu_i,1=0, then the Puiseuxchainpart S_i,1 is empty. |
---|
1670 | // We may start building the Puiseuxchain with part 2, i.e. set start=2. |
---|
1671 | start=2; |
---|
1672 | if (omega>=3) |
---|
1673 | { |
---|
1674 | connectpoint=initial+divseq_sum[2]+1; // startpoint of s_i+1,3 |
---|
1675 | } |
---|
1676 | else |
---|
1677 | { |
---|
1678 | connectpoint=initial+divseq_sum[2]; // endpoint of s_i+1,2 |
---|
1679 | } |
---|
1680 | } |
---|
1681 | else |
---|
1682 | { |
---|
1683 | connectpoint=initial+1; |
---|
1684 | } |
---|
1685 | // Build the Puiseuxchainparts s_i,j and put them as blocks into pc=P_i. |
---|
1686 | multpli=multpl[initial+1..initial+divseq_sum[start]]; |
---|
1687 | list pcp=puiseuxchainpart(initial,divseq[start],multpli,initial_tm,end_tm,start); |
---|
1688 | intmat pc=pcp[1]; |
---|
1689 | intvec tm=pcp[2]; |
---|
1690 | for (j=start+1;j<=omega;j++) |
---|
1691 | { |
---|
1692 | // Keep the final total multipl. of the puiseuxchainpart S_i,j-2 resp. P_i-1, if S_i,1 empty. |
---|
1693 | if (j>2){end_tm=initial_tm;} |
---|
1694 | // Calculate the endpoint of S_i,j-1. |
---|
1695 | initial=initial+divseq[j-1]; |
---|
1696 | // Keep the total multiplicity of the endpoint of S_i,j-1 |
---|
1697 | initial_tm=pcp[2][size(pcp[2])]; |
---|
1698 | // Build the new puiseuxchainpart S_i,j |
---|
1699 | multpli=multpl[initial+1..initial+divseq[j]]; |
---|
1700 | pcp=puiseuxchainpart(initial,divseq[j],multpli,initial_tm,end_tm,j); |
---|
1701 | pc=addmat(pc,pcp[1]); |
---|
1702 | tm=tm,pcp[2]; |
---|
1703 | } |
---|
1704 | // Connect the Puiseuxchainparts s_i,j. |
---|
1705 | for (j=start;j<=omega-2;j++) |
---|
1706 | { |
---|
1707 | k=divseq_sum[j]; // endpoint of s_i,j |
---|
1708 | l=divseq_sum[j+1]+1; // startpoint of s_i,j+2 |
---|
1709 | pc[k,l]=1; // connecting these |
---|
1710 | pc[l,k]=1; |
---|
1711 | } |
---|
1712 | if (omega>=start+1) // If there are at least two non-empty s_i,j. |
---|
1713 | { |
---|
1714 | k=divseq_sum[omega-1]; // endpoint of s_i,omega-1 |
---|
1715 | l=divseq_sum[omega]; // endpoint of s_i,omega |
---|
1716 | pc[k,l]=1; // connecting these |
---|
1717 | pc[l,k]=1; |
---|
1718 | } |
---|
1719 | list ergebnis=pc,connectpoint,tm; |
---|
1720 | return(ergebnis); |
---|
1721 | } |
---|
1722 | |
---|
1723 | static proc irred_resgraph_totmult (intvec charexp) |
---|
1724 | "USAGE: irred_resgraph_totmult(charexp); charexp intvec |
---|
1725 | ASSUME: charexp is an integer vector containg the characteristic exponents of |
---|
1726 | an irreducible plane curve singularity |
---|
1727 | RETURN: list L, L[1] is the incidence matrix of the resolution graph of the plane curve |
---|
1728 | singularity defined by INPUT, and L[2] is its sequence of total multiplicities |
---|
1729 | NOTE: This procedure is only for internal use; it is called by resgraph. |
---|
1730 | " |
---|
1731 | { |
---|
1732 | int k,l; |
---|
1733 | intvec multpl=charexp2multseq(charexp); // multiplicity sequence of the singularity |
---|
1734 | // Do first the case where the singularity is actually smooth. |
---|
1735 | if (size(charexp)==1) |
---|
1736 | { |
---|
1737 | intmat resgraph[1][1]=0; |
---|
1738 | intvec tm=1; // there is no exceptional divisor in the resolution - ERROR: should be 0, |
---|
1739 | // but for the time being, Singular returns 1 as multiplicity of smooth curves |
---|
1740 | list result=resgraph,tm,multpl; |
---|
1741 | return(result); |
---|
1742 | } |
---|
1743 | // End of the smooth case |
---|
1744 | int initial_tm=0; // total multipl. of the endpoint of Puiseux chain P_i-1 |
---|
1745 | int g=size(charexp); |
---|
1746 | list es=divsequence(charexp[2],charexp[1]); // keeps the lenghts of the Puiseuxchainparts s_i,j |
---|
1747 | intvec divseq=es[1]; |
---|
1748 | int r=es[2]; |
---|
1749 | int initial=0; |
---|
1750 | // Build the Puiseuxchains P_i and put them as blocks into a matrix. |
---|
1751 | list pc=puiseuxchain(initial,divseq,multpl,initial_tm); |
---|
1752 | intmat resgraph=pc[1]; |
---|
1753 | intvec endpoints=resgraph[nrows(resgraph),ncols(resgraph)]; |
---|
1754 | intvec connectpoints=pc[2]; |
---|
1755 | intvec tm=pc[3]; |
---|
1756 | for (int i=3;i<=g;i++) |
---|
1757 | { |
---|
1758 | initial_tm=tm[size(tm)]; |
---|
1759 | es=divsequence(charexp[i]-charexp[i-1],r); |
---|
1760 | divseq=es[1]; |
---|
1761 | r=es[2]; |
---|
1762 | initial=endpoints[size(endpoints)]; |
---|
1763 | pc=puiseuxchain(initial,divseq,multpl,initial_tm); |
---|
1764 | resgraph=addmat(resgraph,pc[1]); |
---|
1765 | endpoints=endpoints,resgraph[nrows(resgraph),ncols(resgraph)]; |
---|
1766 | connectpoints=connectpoints,pc[2]; |
---|
1767 | tm=tm,pc[3]; |
---|
1768 | } |
---|
1769 | // Adding the * for the strict transform to the Graph. |
---|
1770 | resgraph=addmat(resgraph,0); |
---|
1771 | // The connecting point of the * with the graph. |
---|
1772 | connectpoints=connectpoints,nrows(resgraph); |
---|
1773 | // Connect the P_i with each other and with *. |
---|
1774 | for (i=2;i<=g;i++) |
---|
1775 | { |
---|
1776 | k=endpoints[i-1]; // endpoint of P_i-1 |
---|
1777 | l=connectpoints[i]; // conncting point of P_i resp. of * |
---|
1778 | resgraph[k,l]=1; // connecting these |
---|
1779 | resgraph[l,k]=1; |
---|
1780 | } |
---|
1781 | list result=resgraph,tm,multpl; //HIER GEAENDERT!!!! |
---|
1782 | return(result); |
---|
1783 | } |
---|
1784 | |
---|
1785 | |
---|
1786 | static proc max_in_intvec (intvec v, list #) |
---|
1787 | "USAGE: max_in_intvec(v); v intvec |
---|
1788 | RETURN: int m, maximum of the integers in v |
---|
1789 | USAGE: max_in_intvec(v,1); v intvec |
---|
1790 | RETURN: intvec m, m[1] maximum of the integers in v, m[2] position of the |
---|
1791 | last occurence of the maximum in v |
---|
1792 | NOTE: This procedure is only for internal use; this procedure is called by |
---|
1793 | totalmultiplicities and semigroup. |
---|
1794 | " |
---|
1795 | { |
---|
1796 | int max=v[1]; |
---|
1797 | int maxpos=1; |
---|
1798 | for (int i=2;i<=size(v);i++) |
---|
1799 | { |
---|
1800 | if (v[i]>max) |
---|
1801 | { |
---|
1802 | max=v[i]; |
---|
1803 | maxpos=i; |
---|
1804 | } |
---|
1805 | } |
---|
1806 | if (size(#)==0) |
---|
1807 | { |
---|
1808 | return(max); |
---|
1809 | } |
---|
1810 | else |
---|
1811 | { |
---|
1812 | return(intvec(max,maxpos)); |
---|
1813 | } |
---|
1814 | } |
---|
1815 | |
---|
1816 | static proc addmat (intmat A,intmat B) |
---|
1817 | "USAGE: max_in_intvec(A,B); A, B integer matrices |
---|
1818 | RETURN: intmat C, block matrix with left-upper block A, right-lower block B |
---|
1819 | NOTE: This procedure is only for internal use; this procedure is called several times. |
---|
1820 | " |
---|
1821 | { |
---|
1822 | int nc=ncols(A); |
---|
1823 | int nr=nrows(A); |
---|
1824 | int mc=ncols(B); |
---|
1825 | int mr=nrows(B); |
---|
1826 | int i,j; |
---|
1827 | intmat AB[nr+mr][nc+mc]; |
---|
1828 | for (i=1;i<=nr;i++) |
---|
1829 | { |
---|
1830 | for (j=1;j<=nc;j++) |
---|
1831 | { |
---|
1832 | AB[i,j]=A[i,j]; |
---|
1833 | } |
---|
1834 | } |
---|
1835 | for (i=1;i<=mr;i++) |
---|
1836 | { |
---|
1837 | for (j=1;j<=mc;j++) |
---|
1838 | { |
---|
1839 | AB[i+nr,j+nc]=B[i,j]; |
---|
1840 | } |
---|
1841 | } |
---|
1842 | return(AB); |
---|
1843 | } |
---|
1844 | |
---|
1845 | static proc divsequence(int a,int b) |
---|
1846 | "USAGE: divsequence(a,b); a,b integers |
---|
1847 | RETURN: list l, l[1] the multiplicities of the divisors in the Euclidean algorithm |
---|
1848 | and l[2] the last non-zero remainder in the Euclidean algorithm |
---|
1849 | NOTE: This procedure is only for internal use; it is called in irred_res_graph. |
---|
1850 | " |
---|
1851 | { |
---|
1852 | int q=a div b; |
---|
1853 | int r=a mod b; |
---|
1854 | intvec divseq=q; |
---|
1855 | while(r<>0) |
---|
1856 | { |
---|
1857 | a=b; |
---|
1858 | b=r; |
---|
1859 | q=a div b; |
---|
1860 | r=a mod b; |
---|
1861 | divseq = divseq,q; |
---|
1862 | } |
---|
1863 | list result=divseq,b; |
---|
1864 | return(result); |
---|
1865 | } |
---|
1866 | |
---|
1867 | |
---|
1868 | |
---|
1869 | static proc adjust_tot_mult (intvec rtm_fix, rtm_var, rmt_fix, rmt_var, p, q, stricttransforms, int k) |
---|
1870 | "USAGE: adjust_tot_mult(v1,v2,v3,v4,p,q,st,k); v1,...,st intvecs and k an integer |
---|
1871 | RETURN: intvec rtm_var, adjusted total multiplicities |
---|
1872 | NOTE: This procedure is only for internal use; it is called in totalmultiplicities. |
---|
1873 | " |
---|
1874 | { |
---|
1875 | int j,l,store; // Help variables. |
---|
1876 | // Recalculate the entries in rtm_var from stricttransforms[k]+1,...,stricttransforms[k+1]-1. |
---|
1877 | for (j=stricttransforms[k]+1;j<stricttransforms[k+1];j++) |
---|
1878 | { |
---|
1879 | if (rtm_var[j]==0) // If the entry is non-zero, we know that it is already correct. |
---|
1880 | { |
---|
1881 | // Check if the vertex in the fixed part is connected to one or to two vertices of lower weight. |
---|
1882 | if (j==stricttransforms[k]+1) // The vertex of weight 1 less is p[k], to which the subgraph is glued. |
---|
1883 | { |
---|
1884 | store=rtm_fix[j]-rmt_fix[j]-rtm_fix[p[k]]; |
---|
1885 | } |
---|
1886 | else // The vertex of weight 1 less belongs to the subgraph. |
---|
1887 | { |
---|
1888 | store=rtm_fix[j]-rmt_fix[j]-rtm_fix[j-1]; |
---|
1889 | } |
---|
1890 | // It is connected to two vertices V (which has weight one less) and W. |
---|
1891 | if (store>0) |
---|
1892 | { |
---|
1893 | if (j==stricttransforms[k]+1) // V is p[k] (to which the subgraph was glued) and W is q[k], the |
---|
1894 | { // vertex of weight one less, to which p[k] is connected. |
---|
1895 | rtm_var[j]=rtm_var[p[k]]+rtm_var[q[k]]; // In this case the subgraph separates p[k] and q[k]! |
---|
1896 | } |
---|
1897 | if (j==stricttransforms[k]+2) // V belongs to the subgraph (it is the vertex considerd in the |
---|
1898 | { // previous case!), and W is p[k] or q[k]. In this case the subgraph separates p[k] and q[k]. |
---|
1899 | if (store==rtm_fix[p[k]]) // If W=p[k], ... |
---|
1900 | { |
---|
1901 | rtm_var[j]=rtm_var[j-1]+rtm_var[p[k]]; |
---|
1902 | } |
---|
1903 | else // else W=q[k] ... . |
---|
1904 | { |
---|
1905 | rtm_var[j]=rtm_var[j-1]+rtm_var[q[k]]; // separates p[k] and q[k]. |
---|
1906 | } |
---|
1907 | } |
---|
1908 | if (j>stricttransforms[k]+2) // V belongs to the subgraph and W either does as well or is p[k]. |
---|
1909 | { |
---|
1910 | l=j-2; |
---|
1911 | while (store<rtm_fix[l]) // Find the second vertex W which is connected to which it is |
---|
1912 | { // connected. It has total multipl. = store! |
---|
1913 | if (l>stricttransforms[k]+1) // If l-1 belongs to the graph, then reduce l. |
---|
1914 | { |
---|
1915 | l=l-1; |
---|
1916 | } |
---|
1917 | else // If l-1 is stricttransform[k], hence isn't in the graph, then reducing l gives p[k]. |
---|
1918 | { // If l=p[k] and still store<rtm_fix[l], then j must be connected to q[k]! |
---|
1919 | if (l==stricttransforms[k]+1) |
---|
1920 | { |
---|
1921 | l=p[k]; |
---|
1922 | } |
---|
1923 | else |
---|
1924 | { |
---|
1925 | l=q[k]; |
---|
1926 | } |
---|
1927 | } |
---|
1928 | } |
---|
1929 | rtm_var[j]=rtm_var[j-1]+rtm_var[l]; |
---|
1930 | } |
---|
1931 | } |
---|
1932 | // It is only connected to one vertex V, which then must be the one of weight one less. |
---|
1933 | else |
---|
1934 | { |
---|
1935 | if (j==stricttransforms[k]+1) // V is p[k], the vertex, to which the subgraph was glued. |
---|
1936 | { |
---|
1937 | rtm_var[j]=rtm_var[p[k]]; |
---|
1938 | } |
---|
1939 | else |
---|
1940 | { |
---|
1941 | rtm_var[j]=rtm_var[j-1]; // V is belongs already to the subgraph. |
---|
1942 | } |
---|
1943 | } |
---|
1944 | } |
---|
1945 | } |
---|
1946 | return(rtm_var); |
---|
1947 | } |
---|
1948 | |
---|
1949 | |
---|
1950 | static proc find_connection_point (intvec v, int k) |
---|
1951 | "USAGE: find_connection_point(v,c); where v is an intvec, and k is an integer |
---|
1952 | RETURN: The largest integer i>k, such that v[i]=1, or 0 if no such i exists. |
---|
1953 | NOTE: This procedure is only for internal use; it is called in totalmultiplicities. |
---|
1954 | " |
---|
1955 | { |
---|
1956 | for (int i=size(v)-1;i>=k+1;i--) |
---|
1957 | { |
---|
1958 | if (v[i]==1) |
---|
1959 | { |
---|
1960 | return(i); |
---|
1961 | } |
---|
1962 | } |
---|
1963 | return(0); |
---|
1964 | } |
---|
1965 | |
---|
1966 | static proc find_connection_points (intmat resgr, int k) |
---|
1967 | "USAGE: find_connection_points(resgr,k); where resgr is an intmat, and k is an integer |
---|
1968 | RETURN: list of two intvecs ctps and ctpswts, where ctps contains all integers i!=k, such |
---|
1969 | that resgr[k,i]=1, and ctpswts contains for each such i the weight resgr[i,i]. |
---|
1970 | NOTE: This procedure is only for internal use; it is called in totalmultiplicities. |
---|
1971 | " |
---|
1972 | { |
---|
1973 | intvec ctps; |
---|
1974 | intvec ctpswts; |
---|
1975 | int j=1; |
---|
1976 | for (int i=1;i<=ncols(resgr);i++) |
---|
1977 | { |
---|
1978 | if ((resgr[k,i]==1) and (i!=k)) |
---|
1979 | { |
---|
1980 | ctps[j]=i; |
---|
1981 | ctpswts[j]=resgr[i,i]; |
---|
1982 | j++; |
---|
1983 | } |
---|
1984 | } |
---|
1985 | return(list(ctps,ctpswts)); |
---|
1986 | } |
---|
1987 | |
---|
1988 | static proc find_lower_connection_points (intmat resgr, int k) |
---|
1989 | "USAGE: find_lower_connection_points(resgr,k); where resgr is an intmat, and k is an integer |
---|
1990 | ASSUME: resgr is the resolutiongraph of an IRREDUCIBLE curve singularity and k<ncols(resgr). |
---|
1991 | RETURN: intvec, which contains the weights of points of lower weight than k, to which |
---|
1992 | the point of weight k in resgr is connected, and 0 if no such point exists. |
---|
1993 | NOTE: This procedure is only for internal use; it is called in totalmultiplicities. |
---|
1994 | " |
---|
1995 | { |
---|
1996 | intvec ctps=find_connection_points(resgr,k)[2]; |
---|
1997 | intvec lower_ctps; |
---|
1998 | int i=1; |
---|
1999 | while ((ctps[i]<k) and (ctps[i]>0)) |
---|
2000 | { |
---|
2001 | lower_ctps[i]=ctps[i]; |
---|
2002 | i++; |
---|
2003 | } |
---|
2004 | return(lower_ctps); |
---|
2005 | } |
---|
2006 | |
---|
2007 | |
---|
2008 | static proc euclides(int a,int b) // Procedure of Fernando. |
---|
2009 | "USAGE: euclides(m,n);where m,n are integers. |
---|
2010 | RETURN: The divisor, the quotients and the remainders of the euclidean algorithm performing for m and n. |
---|
2011 | NOTE: This procedure is only for internal use; it is called in multseq2charexp. |
---|
2012 | " |
---|
2013 | { |
---|
2014 | int c=a div b;//we compute in c the integer division between a and b. |
---|
2015 | int r=a mod b;//in r the remainer of the division between a and b |
---|
2016 | intvec dividendo=c;// we define the intvec of the dividens and we initialize it to c |
---|
2017 | intvec remain=r;// we define the intvec of the remainders and we initialize it to r |
---|
2018 | a=b;//change a to b |
---|
2019 | b=r;// and b to r |
---|
2020 | |
---|
2021 | while(r<>0)// while the current remainder is diferent to 0 we make the same af before |
---|
2022 | { |
---|
2023 | c=a div b; |
---|
2024 | r=a mod b; |
---|
2025 | dividendo =dividendo,c; |
---|
2026 | if(r<>0) |
---|
2027 | { |
---|
2028 | remain=remain,r; |
---|
2029 | } |
---|
2030 | a=b; |
---|
2031 | b=r; |
---|
2032 | } |
---|
2033 | list L=dividendo,remain;//we put in a list all the dividens and all the remainders |
---|
2034 | return(L);// and return the list |
---|
2035 | } |
---|
2036 | |
---|
2037 | |
---|
2038 | |
---|
2039 | static proc tau_es_qh (poly f) |
---|
2040 | "USAGE: tau_es_qh(f), poly f |
---|
2041 | RETURN: int, equisingular Tjurina number |
---|
2042 | NOTE: This procedure is only for internal use; it is called in Tau_es. |
---|
2043 | " |
---|
2044 | { |
---|
2045 | intvec qh=qhweight(f); |
---|
2046 | if (qh[1]==0) |
---|
2047 | { |
---|
2048 | ERROR("Input is not quasi-homogenous."); |
---|
2049 | } |
---|
2050 | else |
---|
2051 | { |
---|
2052 | int d_f = deg(f,qh); |
---|
2053 | list Tj=Tjurina(f,1); |
---|
2054 | ideal tj=Tj[2]; |
---|
2055 | int Taues=size(tj); |
---|
2056 | for (int i=1;i<=size(tj);i++) |
---|
2057 | { |
---|
2058 | if (deg(tj[i],qh)>=d_f) |
---|
2059 | { |
---|
2060 | Taues--; |
---|
2061 | } |
---|
2062 | } |
---|
2063 | } |
---|
2064 | return(Taues); |
---|
2065 | } |
---|
2066 | |
---|
2067 | |
---|
2068 | static proc move_row (intmat M, int i,j) |
---|
2069 | "USAGE: move_row(M,i,j), intmat M, int i,j |
---|
2070 | RETURN: intmat, the matrix M with j-th row now i-th row and the remaining rows moved accordingly. |
---|
2071 | NOTE: This procedure is only for internal use; it is called in move_row_col. |
---|
2072 | " |
---|
2073 | { |
---|
2074 | if ((i>nrows(M)) or (j>nrows(M))) |
---|
2075 | { |
---|
2076 | ERROR("The matrix has not enough rows."); |
---|
2077 | } |
---|
2078 | if (i==j) |
---|
2079 | { |
---|
2080 | return(M); |
---|
2081 | } |
---|
2082 | if (i>1) |
---|
2083 | { |
---|
2084 | intmat N[nrows(M)+1][ncols(M)]=intvec(M[1..i-1,1..ncols(M)]),intvec(M[j,1..ncols(M)]),intvec(M[i..nrows(M),1..ncols(M)]); |
---|
2085 | } |
---|
2086 | if (i==1) |
---|
2087 | { |
---|
2088 | intmat N[nrows(M)+1][ncols(M)]=intvec(M[j,1..ncols(M)]),intvec(M[i..nrows(M),1..ncols(M)]); |
---|
2089 | } |
---|
2090 | if (i<j) |
---|
2091 | { |
---|
2092 | N=delete_row(N,j+1); |
---|
2093 | } |
---|
2094 | if (i>j) |
---|
2095 | { |
---|
2096 | N=delete_row(N,j); |
---|
2097 | } |
---|
2098 | return(N); |
---|
2099 | } |
---|
2100 | |
---|
2101 | static proc move_col (intmat M, int i,j) |
---|
2102 | "USAGE: move_col(M,i,j), intmat M, int i,j |
---|
2103 | RETURN: intmat, the matrix M with j-th column now i-th column and the remaining columns moved accordingly. |
---|
2104 | NOTE: This procedure is only for internal use; it is called in move_row_col. |
---|
2105 | " |
---|
2106 | { |
---|
2107 | return(transpose(move_row(transpose(M),i,j))); |
---|
2108 | } |
---|
2109 | |
---|
2110 | static proc move_row_col (intmat M,int i,j) |
---|
2111 | "USAGE: move_row(M,i,j), intmat M, int i,j |
---|
2112 | RETURN: intmat, the matrix M with j-th row/column now i-th row/column and the remaining |
---|
2113 | rows and columns moved accordingly. |
---|
2114 | NOTE: This procedure is only for internal use; it is called in totalmultiplicities. |
---|
2115 | " |
---|
2116 | { |
---|
2117 | return(move_col(move_row(M,i,j),i,j)); |
---|
2118 | } |
---|
2119 | |
---|
2120 | |
---|
2121 | static proc delete_row (intmat M,int i) |
---|
2122 | "USAGE: delete_row(M,i); M intmat, i int |
---|
2123 | RETURN: intmat, which is derived from M by removing the ith row |
---|
2124 | NOTE: This procedure is only for internal use; it is called in move_row and tau_es. |
---|
2125 | " |
---|
2126 | { |
---|
2127 | if ((i!=1) and (i!=nrows(M))) |
---|
2128 | { |
---|
2129 | return(intmat(intvec(M[1..i-1,1..ncols(M)],M[i+1..nrows(M),1..ncols(M)]),nrows(M)-1,ncols(M))); |
---|
2130 | } |
---|
2131 | if (i==1) |
---|
2132 | { |
---|
2133 | return(intmat(intvec(M[2..nrows(M),1..ncols(M)]),nrows(M)-1,ncols(M))); |
---|
2134 | } |
---|
2135 | else |
---|
2136 | { |
---|
2137 | return(intmat(intvec(M[1..nrows(M)-1,1..ncols(M)]),nrows(M)-1,ncols(M))); |
---|
2138 | } |
---|
2139 | } |
---|
2140 | |
---|
2141 | static proc delete_col (intmat M, int i) |
---|
2142 | "USAGE: delete_col(M,i); M intmat, i int |
---|
2143 | RETURN: intmat, which is derived from M by removing the ith column |
---|
2144 | NOTE: This procedure is only for internal use; it is called in tau_es. |
---|
2145 | " |
---|
2146 | { |
---|
2147 | return(transpose(delete_row(transpose(M),i))); |
---|
2148 | } |
---|
2149 | |
---|
2150 | |
---|
2151 | |
---|
2152 | static proc sort_branches (intmat contact, list graphs, list totmult, list multpl, int k,l) |
---|
2153 | "USAGE: sort_branches(K,L,M,N,k,l); K intmat, L,M,N lists, k,l integers |
---|
2154 | ASSUME: K = contact matrix of the branches of a curve, L = their resolutiongraphs, |
---|
2155 | M = their totalmultiplicities, N = their multiplicities |
---|
2156 | RETURN: list LL, LL[1] transformed K, LL[2..4] transformed L,M,N. |
---|
2157 | The procedure sorts the branches from k+1 to l according to descending contact with |
---|
2158 | with the branch k. |
---|
2159 | NOTE: This procedure is only for internal use; it is called in totalmultiplicities. |
---|
2160 | " |
---|
2161 | { |
---|
2162 | intvec max; |
---|
2163 | for (int i=k+1;i<=l;i++) |
---|
2164 | { |
---|
2165 | // Find the last graph max between i and l which has maximal contact with k. |
---|
2166 | max=max_in_intvec(intvec(contact[k,i..l]),1); |
---|
2167 | max[2]=max[2]+i-1; |
---|
2168 | if (i!=max[2]) // If max is not i, then move max to position i! |
---|
2169 | { |
---|
2170 | graphs=insert(graphs,graphs[max[2]],i-1); |
---|
2171 | graphs=delete(graphs,max[2]+1); |
---|
2172 | totmult=insert(totmult,totmult[max[2]],i-1); |
---|
2173 | totmult=delete(totmult,max[2]+1); |
---|
2174 | multpl=insert(multpl,multpl[max[2]],i-1); |
---|
2175 | multpl=delete(multpl,max[2]+1); |
---|
2176 | contact=move_row_col(contact,i,max[2]); |
---|
2177 | } |
---|
2178 | } |
---|
2179 | return(list(contact,graphs,totmult,multpl)); |
---|
2180 | } |
---|
2181 | |
---|
2182 | |
---|
2183 | static proc find_last_non_one (intvec v,int k) |
---|
2184 | "USAGE: find_last_non_one (v,k); intvec v, int k |
---|
2185 | RETURN: int i, the last index i>=k such that v[i]!=1, or k-1 if no such i exists. |
---|
2186 | NOTE: This procedure is only for internal use; it is called in tau_es. |
---|
2187 | " |
---|
2188 | { |
---|
2189 | int i=size(v); |
---|
2190 | while (i>=k) |
---|
2191 | { |
---|
2192 | if (v[i]!=1) |
---|
2193 | { |
---|
2194 | return(i); |
---|
2195 | } |
---|
2196 | else |
---|
2197 | { |
---|
2198 | i--; |
---|
2199 | } |
---|
2200 | } |
---|
2201 | return(i); |
---|
2202 | } |
---|
2203 | |
---|
2204 | static proc intmat_minus_one (intmat M) |
---|
2205 | "USAGE: intmat_minus_one(M); intmat M |
---|
2206 | RETURN: intmat, all non-zero entries of M deminuished by 1. |
---|
2207 | NOTE: This procedure is only for internal use; it is called in tau_es. |
---|
2208 | " |
---|
2209 | { |
---|
2210 | int i,j; |
---|
2211 | for (i=1;i<=nrows(M);i++) |
---|
2212 | { |
---|
2213 | for (j=1;j<=ncols(M);j++) |
---|
2214 | { |
---|
2215 | if (M[i,j]!=0) |
---|
2216 | { |
---|
2217 | M[i,j]=M[i,j]-1; |
---|
2218 | } |
---|
2219 | } |
---|
2220 | } |
---|
2221 | return(M); |
---|
2222 | } |
---|
2223 | |
---|