1 /* PLUGIN file to use with plantri.c
2
3 To use this, compile plantri.c using
4 cc -o plantri_ad -O '-DPLUGIN="allowed_deg.c"' plantri.c
5
6 This plug-in deletes those triangulations
7 with vertex degrees which are not allowed.
8 Allowed degrees may be defined by using the -F switch, for example
9
10 plantri_ad -F7F5 14
11
12 makes all triangulations with 14 vertices and degrees 5 or 7.
13
14 authors: Gunnar Brinkmann and Ulrike von Nathusius,
15 contains parts of maxdeg.c by Brendan McKay.
16
17 The upper and lower limits for the number of faces to be used
18 can be given like e.g. plantri_ad -F7_1^3F5F6 14 forcing between
19 one and 3 vertices with valency 7. At the moment this is only
20 implemented as a filter at the end. It could be used for better
21 bounding criteria -- should be implemented once...
22
23 The nonstandard (but common) long long type is required.
24 (Warning: some versions of the Sun "cc" compiler give incorrect
25 results with programs using long long.)
26 */
27
28 #define FAST_FILTER_SIMPLE ((nv<maxnv) || ad_filter())
29 #define FIND_EXTENSIONS_SIMPLE find_extensions_ad
30
31 /************************** Switches ********************************/
32
33 /* The following adds the switch f to those normally present.
34 arg is the address of the command-line argument, and j is the
35 index where 'F' might be. The value of j must be left on the
36 last digit of the switch value. */
37
38 #undef SWITCHES
39 #define SWITCHES "[-F#[_#^#] -uagsh -odG -v]"
40
41 #define PLUGIN_SWITCHES else if (arg[j]=='F') list_of_allowed_degrees(arg, &j);
42
43 #define PLUGIN_INIT if (minconnec >= 0 || polygonsize >= 0 \
44 || minimumdeg >= 0 || pswitch || xswitch || tswitch) \
45 {fprintf(stderr,">E Usage: %s %s n [res/mod] [outfile]\n",cmdname,SWITCHES);\
46 exit(1);}
47
48 /********************************************************************/
49
50 #define INFTY_AD (3*MAXN)
51
52 static int maxdeg = 0; /* the maximum allowed degree */
53 static unsigned long long int LISTE=0; /* a binary representation of the allowed degrees */
54 static unsigned long long int mask[MAXN];
55 static int error_up[MAXN];
56 static int error_down[MAXN];
57 static int error_of_degree[MAXN];
58 static int maxnumber[MAXN],minnumber[MAXN];
59 static int degreelist[MAXN]; /* a list of allowed degrees ended by 0 */
60 static int grenzen=0;
61 static int x_error=0;
62
63 int error_1[MAXN]; /* the error for the special case x=1 */
64 int doextra_x;
65
66 /********************************* Filter **************************/
67
68 static int
ad_filter()69 ad_filter()
70 {
71 int i, run, anzahl[MAXN];
72
73 for(i=0;i<nv ; i++) if (error_of_degree[degree[i]]) return(0);
74
75 if (grenzen)
76 {
77 for (run=0; degreelist[run]; run++) anzahl[degreelist[run]]=0;
78 for(i=0;i<nv ; i++) (anzahl[degree[i]])++;
79 for (run=0; degreelist[run]; run++)
80 if ((anzahl[degreelist[run]]<minnumber[degreelist[run]]) ||
81 (anzahl[degreelist[run]]>maxnumber[degreelist[run]])) return(0);
82 }
83
84 return(1);
85
86 }
87
88
89 /******************* nur_noch_E3 *******************************/
90
91 int
nur_noch_E3()92 nur_noch_E3()
93 {
94 register EDGE *e;
95 unsigned long long int vmask;
96 int a,b,i,ident,deg3;
97
98 /* This function may only be used, if there are exactly two vertices of
99 degree 3. Two vertices of degree three that do not have two common
100 neighbours can only disappear (get a higher degree) by inserting
101 at least two new such vertices with the same property, so only E3
102 can be used from now on:
103 The distance between two such vertices can not decrease with E3,
104 and if E3 is applied to one of them, the new 3-vertex will have
105 (at least) the same distances from the remaining old one.
106
107 The function will return TRUE, if they have no common neighbours
108 so that only E3 can be used further, otherwise FALSE.
109 */
110
111 for (i = deg3 = 0; i < nv ; ++i)
112 if (degree[i] == 3)
113 {
114 if (deg3 == 0) {a = i;deg3=1;}
115 else {b = i;i=nv;}
116 }
117 e = firstedge[a];
118 vmask=mask[e->end]; /* vmask will be the set of neighbours of a */
119 vmask|=mask[e->prev->end];
120 vmask|=mask[e->next->end];
121
122
123 e = firstedge[b]; /* Test if two of the neighbours of b are in vmask*/
124 if (mask[e->end] & vmask)
125 { if ((mask[e->prev->end] & vmask) || (mask[e->next->end] & vmask)) return FALSE; }
126 else
127 { if ((mask[e->prev->end] & vmask) && (mask[e->next->end] & vmask)) return FALSE; }
128
129 return TRUE;
130 }
131
132 /************************** list_of_allowed_degrees ***********************/
133
134 /* The allowed degrees taken from the F-Switch LISTE are described as bits
135 of the binary representation of an integer.
136 The definition of the error_of_degree is explained in the text below.
137 */
138
139 void
list_of_allowed_degrees(char arg[],int * pj)140 list_of_allowed_degrees(char arg[], int *pj )
141 {
142 static int init=1;
143 int i,n,j,x,run,maxdegree, newi, newj, numlargesizes;
144 if (init) {
145 init=0;
146 maxdegree=8*sizeof(unsigned long long)-1;
147 mask[0]=1;
148 for(i=1;(i<=maxdegree) && (i<MAXN);i++) mask[i]=mask[i-1]<<1;
149 for(i=0;i<MAXN;i++)
150 error_up[i]=error_down[i]=error_of_degree[i]=error_1[i]=INFTY_AD;
151 for(i=0;i<MAXN;i++) degreelist[i]=0;
152 }
153 j=*pj;
154 if (!isdigit(arg[j+1]))
155 {fprintf(stderr,"No degrees given! Problem: %s\n",arg+j);exit(0);}
156 else
157 { /* get the switch values of option -F */
158 n = atoi (arg+j+1);
159 if (n>=MAXN) {fprintf(stderr,"Maximum n is %d!\n",MAXN);exit(0);}
160 if (n>maxdegree) {fprintf(stderr,"Maximum degree is %d!\n",maxdegree);exit(0);}
161 LISTE |= mask[n];
162 maxnumber[n]=MAXN; minnumber[n]=0;
163 for (run=0; degreelist[run]; run++); degreelist[run]=n;
164 if(n>maxdeg) maxdeg = n;
165 for (i=3; i<=n; i++){
166 if ( error_up[i] > (n-i)) error_up[i]=n-i;
167 }
168 if(n>=5){
169 for(i=n;i<MAXN;i++){
170 if(error_down[i]> (i-n)){
171 error_down[i]=i-n;
172 }
173 }
174 }
175 if(error_up[3]&&error_up[5]) x=error_up[5] - 1;
176 else x=error_up[5];
177 for(i=3;i<MAXN;i++)
178 { if (2*error_down[i]<error_up[i])
179 error_1[i]=2*error_down[i];
180 else error_1[i]=error_up[i];
181 }
182
183 for(i=3;i<MAXN;i++){
184 if ((1+x)*error_down[i]<error_up[i]){
185 error_of_degree[i]=(1+x)*error_down[i]; }
186 else {
187 error_of_degree[i]=error_up[i];
188 }
189 }
190 j++;
191 while ((arg[j]>='0') && (arg[j]<='9')) j++;
192 j--;
193 /* Two times the same to be independent of the order */
194 if (arg[j+1]=='_') { grenzen=1; j+=2;
195 minnumber[n]=atoi(arg+j);
196 while ((arg[j]>='0') && (arg[j]<='9')) j++; j--; }
197 if (arg[j+1]=='^') { grenzen=1; j+=2;
198 maxnumber[n]=atoi(arg+j);
199 while ((arg[j]>='0') && (arg[j]<='9')) j++; j--; }
200 if (arg[j+1]=='_') { grenzen=1; j+=2;
201 minnumber[n]=atoi(arg+j);
202 while ((arg[j]>='0') && (arg[j]<='9')) j++; j--; }
203 if (arg[j+1]=='^') { grenzen=1; j+=2;
204 maxnumber[n]=atoi(arg+j);
205 while ((arg[j]>='0') && (arg[j]<='9')) j++; j--; }
206 *pj = j;
207 }
208
209 if (error_of_degree[3] && error_of_degree[5]) x_error= error_of_degree[5]-1;
210 else x_error= error_of_degree[5];
211
212 for (i=6,numlargesizes=0;i<=maxdeg;i++) if (error_of_degree[i]==0) numlargesizes++;
213
214 if ((numlargesizes<=2) && (error_of_degree[5]==0))
215 doextra_x=1; else doextra_x=0;
216
217 }
218
219
220
221 /**************************************************************************/
222
223 static void
find_extensions_ad(int numb_total,int numb_pres,EDGE * ext3[],int * numext3,EDGE * ext4[],int * numext4,EDGE * ext5[],int * numext5)224 find_extensions_ad(int numb_total, int numb_pres,
225 EDGE *ext3[], int *numext3,
226 EDGE *ext4[], int *numext4,
227 EDGE *ext5[], int *numext5)
228
229 /* Find all nonequivalent places for extension.
230 These are listed in ext# according to the degree of the future new vertex.
231 The number of cases is returned in numext# (#=3,4,5).
232
233 For explanations of WHY this is correct, see plantri_ad paper.
234
235 */
236 {
237 int d3,d4;
238 register EDGE *e,*e1,*e2,*ex, *start;
239 EDGE **nb0,**nb1,**nbop,**nblim;
240
241 register int i,i1,i2,i3,k,levs,excess;
242 unsigned long long int vmask=0;
243 int a,b,j,x,eintraege,error,ierror,nur_E3,nur_E5,do_E3,do_E4,do_E5,newi,newj;
244 int eckenliste[MAXN];
245 int d[MAXN]; /* number of vertices with degree [i] */
246 int maxd, minimprove, improvebuffer, improvebuffer2;
247 int errors1;
248
249
250 if (nv==4) { *numext3=*numext4=1; *numext5=0;
251 ext3[0]=ext4[0]=firstedge[0];
252 return; }
253
254 *numext3=*numext4=*numext5=0;
255 do_E3=do_E4=do_E5=1;
256 levs = maxnv - nv; /* remaining number of steps */
257 excess = 0;
258
259 for (i = 3; i < nv; ++i) d[i]=0;
260 for (i = 0; i < nv; ++i) (d[degree[i]])++;
261 for (maxd=nv-1;d[maxd]==0;maxd--);
262 for (i=maxd;i>maxdeg;i--) excess+= (d[i]*(i-maxdeg));
263
264
265 d3=d[3]; d4=d[4];
266 nur_E3 = ((d3>2) || ((d3==2) && nur_noch_E3()));
267 if (error_of_degree[3] && nur_E3) { return; }
268
269 i = d3 + d3 + d[4]; /* By (4): */
270
271 if (excess)
272 {
273 if (nur_E3 || (excess > levs)) { return; }
274 if (d3 > 0 && (excess >= levs)) { return; }
275
276 if ((i > 0) && (excess > levs - i + 2)) { return; }
277 }
278
279 if (error_of_degree[3]&&error_of_degree[4]&&(i>levs+1)) { return; }
280
281 if (!nur_E3)
282 {
283 /* In this case it is not possible to exclude any of the three
284 operations. */
285 for (i=3,error=0;i<=maxd;i++) error += (error_of_degree[i]*d[i]);
286
287 if (error_up[3]&&error_up[5])
288 {
289 // if (error>2*levs) return; extensions chosen so that this cannot happen
290 minimprove=error-(2*(levs-1));
291 }
292 else {
293 // if (error>3*levs) return; extensions chosen so that this cannot happen
294 minimprove=error-(3*(levs-1));
295 }
296 do_E3=((minimprove+error_of_degree[3])<=3);
297 do_E4=((minimprove+error_of_degree[4])<=2);
298 do_E5=((minimprove+error_of_degree[5])<=3+x_error);
299
300
301 if (doextra_x)
302 {
303 for (i=3,errors1=0;i<=maxd;i++)
304 { errors1 += (error_1[i]*d[i]); }
305 if (errors1 > (4*levs)) { return; }
306 /* only done if error(5)=0, so m(S)=4 for x=1 */
307 }
308
309
310 } /* end: all operations possible */
311
312 else {
313 /* Only Operation E3 can be used from now on, so that "error_of_degree" may be
314 substituted by "error_up". */
315
316 do_E4=do_E5=0;
317
318 for (i=4,error=0;i<=maxd;i++) /* error_of_degree[3]=0 -- otherwise the routine
319 would already have returned */
320 { error += (error_up[i]*d[i]); }
321 if (error > (3*levs)) { return; }
322 minimprove=error-(3*(levs-1));
323
324 if (error>levs)
325 { /* We try to find a set of independent vertices. E3
326 can decrease the error of only one vertex of the
327 set. So the number of remaining steps cannot be
328 less than the sum of errors of this vertex-set. In
329 order to get an effective condition one starts the
330 computation of neighbours with vertices which have
331 a large error. */
332
333
334
335 /* First rank the vertices with respect to the error of
336 their degrees */
337
338
339 for (i=eintraege=0; i<nv; i++){
340 if (error=error_up[degree[i]])
341 {
342 for (j=eintraege-1;(j>=0) && ((error_up[degree[eckenliste[j]]]) < error); j--)
343 eckenliste[j+1]=eckenliste[j];
344 for ( ;(j>=0) &&
345 (error_up[degree[eckenliste[j]]] == error) && (degree[eckenliste[j]]>degree[i]); j--)
346 eckenliste[j+1]=eckenliste[j];
347
348 eckenliste[j+1]=i;
349 eintraege++;
350 }
351 }
352
353 for(i=ierror=vmask=0;i<eintraege; i++)
354 {
355 a = eckenliste[i]; /* computation of the neighbours */
356 if (!(vmask & mask[a]))
357 {
358 start = e = firstedge[a];
359 vmask|=mask[e->end]; /* vmask: set of forbidden vertices */
360 while((e=e->next) != start){
361 vmask|=mask[e->end];
362 }
363 ierror += error_up[degree[a]];
364 /* error: minimum number of necessary steps */
365 if (ierror>levs) { return; }
366 }
367 }
368 }/* end error>levs -- worth trying to find such an independent set */
369
370 }
371
372
373 nur_E5= (excess==levs);
374 if (nur_E5) { do_E3=do_E4=0;}
375
376 if ((excess>0) && (excess==levs-1)) do_E3=0;
377
378 if (d3 || (d[4] > 2)) do_E5=0;
379
380
381 /* code for trivial group */
382
383 //find_extensions(numb_total, numb_pres,ext3, numext3, ext4, numext4,ext5, numext5); return;
384
385
386 if (numb_total == 1)
387 {
388
389 if (do_E3)
390 { improvebuffer=minimprove+error_of_degree[3];
391 k = 0;
392 for (i = 0; i < nv; ++i)
393 {
394 if (nur_E3) improvebuffer2=(error_up[degree[i]]-error_up[degree[i]+1]);
395 else improvebuffer2=(error_of_degree[degree[i]]-error_of_degree[degree[i]+1]);
396 if (improvebuffer-improvebuffer2<=2)
397 {
398 e = ex = firstedge[i];
399 do
400 {
401 i2=degree[e->end]; i3=degree[e->next->end];
402 if ( (nur_E3 &&
403 (improvebuffer<=improvebuffer2+
404 (error_up[i2]-error_up[i2+1]) +(error_up[i3]-error_up[i3+1])))
405
406 ||
407 (!nur_E3 && (improvebuffer<=improvebuffer2+
408 (error_of_degree[i2]-error_of_degree[i2+1])
409 +(error_of_degree[i3]-error_of_degree[i3+1]))))
410 {
411 e1 = e->invers->prev;
412 if (e1 > e)
413 {
414 e1 = e1->invers->prev;
415 if (e1 > e) ext3[k++] = e;
416 }
417 }
418 e = e->next;
419 } while (e != ex);
420 }
421 }/* end of for() */
422 *numext3 = k;
423 }
424
425 if (do_E4)
426 { improvebuffer=minimprove+error_of_degree[4];
427 k = 0;
428 for (i = 0; i < nv; ++i)
429 {
430 if (degree[i] == 3) continue;
431 e = ex = firstedge[i];
432 do
433 {
434 i2=degree[e->end]; i3=degree[e->next->next->end];
435 if (improvebuffer<=
436 (error_of_degree[i2]-error_of_degree[i2+1])
437 +(error_of_degree[i3]-error_of_degree[i3+1]))
438 {
439 e1 = e->next;
440 if (e1->invers > e1)
441 {
442 e2 = e1->invers->prev;
443 if ((degree[e->end] == 3)
444 + (degree[e2->end] == 3) == d3)
445 ext4[k++] = e;
446 }
447 }
448 e = e->next;
449 } while (e != ex);
450 }
451 *numext4 = k;
452 }
453
454 if (do_E5)
455 { improvebuffer=minimprove+error_of_degree[5];
456 k = 0;
457 for (i = 0; i < nv; ++i)
458 {
459 if (degree[i] < 6) continue;
460 if (nur_E5 && (degree[i]<=maxdeg)) continue;
461 improvebuffer2=(error_of_degree[degree[i]]-error_of_degree[degree[i]-1]);
462 if (improvebuffer-improvebuffer2<=2)
463 {
464 e = ex = firstedge[i];
465 do
466 { i2=degree[e->end]; i3=degree[e->next->next->next->end];
467 if (improvebuffer<=improvebuffer2+
468 (error_of_degree[i2]-error_of_degree[i2+1])
469 +(error_of_degree[i3]-error_of_degree[i3+1]))
470 {
471 e1 = e->next->next->next;
472 if ((degree[e->end] == 4)
473 + (degree[e1->end] == 4) == d4)
474 ext5[k++] = e;
475 }
476 e = e->next;
477 } while (e != ex);
478 }
479 }
480 *numext5 = k;
481 }
482 }
483
484
485
486
487 /* code for nontrivial group */
488
489 else
490 {
491 nb0 = (EDGE**)numbering[0];
492 nbop = (EDGE**)numbering[numb_pres == 0 ? numb_total : numb_pres];
493 nblim = (EDGE**)numbering[numb_total];
494
495 for (i = 0; i < ne; ++i) nb0[i]->index = i;
496
497 if (do_E3)
498 {
499 RESETMARKS;
500 improvebuffer=minimprove+error_of_degree[3];
501 k = 0;
502 for (i = 0; i < ne; ++i)
503 {
504 e = nb0[i];
505 if (ISMARKEDLO(e)) continue;
506 i1=degree[e->start]; i2=degree[e->end]; i3=degree[e->next->end];
507 if (( !nur_E3 && (improvebuffer<=
508 (error_of_degree[i1]-error_of_degree[i1+1])
509 +(error_of_degree[i2]-error_of_degree[i2+1])
510 +(error_of_degree[i3]-error_of_degree[i3+1]))) ||
511 ( nur_E3 && (improvebuffer<=
512 (error_up[i1]-error_up[i1+1])
513 +(error_up[i2]-error_up[i2+1])
514 +(error_up[i3]-error_up[i3+1]))))
515 ext3[k++] = e;
516
517 for (nb1 = nb0 + i; nb1 < nbop; nb1 += MAXE) MARKLO(*nb1);
518
519 for (; nb1 < nblim; nb1 += MAXE) MARKLO((*nb1)->invers);
520
521 e1 = e->invers->prev;
522 for (nb1 = nb0 + e1->index; nb1 < nbop; nb1 += MAXE) MARKLO(*nb1);
523
524 for (; nb1 < nblim; nb1 += MAXE) MARKLO((*nb1)->invers);
525
526 e1 = e1->invers->prev;
527 for (nb1 = nb0 + e1->index; nb1 < nbop; nb1 += MAXE) MARKLO(*nb1);
528
529 for (; nb1 < nblim; nb1 += MAXE) MARKLO((*nb1)->invers);
530 }
531 *numext3 = k;
532 }
533
534 if (do_E4)
535 {
536 RESETMARKS;
537 improvebuffer=minimprove+error_of_degree[4];
538 k = 0;
539 for (i = 0; i < ne; ++i)
540 {
541 e = nb0[i];
542 if (ISMARKEDLO(e)) continue;
543 e1 = e->next->invers->prev;
544 if ((degree[e->end] == 3) + (degree[e1->end] == 3) != d3)
545 continue;
546 i2=degree[e->end]; i3=degree[e->next->next->end];
547 if (improvebuffer<=
548 (error_of_degree[i2]-error_of_degree[i2+1])
549 +(error_of_degree[i3]-error_of_degree[i3+1]))
550 ext4[k++] = e;
551
552 for (nb1 = nb0 + i; nb1 < nbop; nb1 += MAXE) MARKLO(*nb1);
553
554 for (; nb1 < nblim; nb1 += MAXE) MARKLO((*nb1)->prev->prev);
555
556 for (nb1 = nb0 + e1->index; nb1 < nbop; nb1 += MAXE)
557 MARKLO(*nb1);
558
559 for (; nb1 < nblim; nb1 += MAXE) MARKLO((*nb1)->prev->prev);
560 }
561 *numext4 = k;
562 }
563
564 if (do_E5)
565 {
566 RESETMARKS;
567 improvebuffer=minimprove+error_of_degree[5];
568 k = 0;
569 for (i = 0; i < ne; ++i)
570 {
571 e = nb0[i];
572 if (ISMARKEDLO(e) || degree[e->start] < 6) continue;
573
574 if ((degree[e->end] == 4)
575 + (degree[e->next->next->next->end] == 4) != d4)
576 continue;
577 if (nur_E5 && (degree[e->start]<=maxdeg)) continue;
578 i1=degree[e->start]; i2=degree[e->end]; i3=degree[e->next->next->next->end];
579 if (improvebuffer<=
580 (error_of_degree[i1]-error_of_degree[i1-1])
581 +(error_of_degree[i2]-error_of_degree[i2+1])
582 +(error_of_degree[i3]-error_of_degree[i3+1]))
583 ext5[k++] = e;
584
585 for (nb1 = nb0 + i; nb1 < nbop; nb1 += MAXE) MARKLO(*nb1);
586
587 for (; nb1 < nblim; nb1 += MAXE)
588 MARKLO((*nb1)->prev->prev->prev);
589 }
590 *numext5 = k;
591 }
592
593 }
594 }
595
596
597