1 /* cubhamg.c : pick those inputs that are nonhamiltonian and
2 have max degree <= 3.
3
4 Usage:
5 cubhamg [-#] [-v|-V] [-n#-#|-y#-#|-i|-I|-o|-x|-e|-E] [-b|-t] [infile [outfile]]
6
7 infile is the name of the input file in graph6/sparse6 format
8 outfile is the name of the output file in the same format
9
10 stdin and stdout are the defaults for infile and outfile
11
12 The output file will have a header >>graph6<< or >>sparse6<<
13 if and only if the input file does.
14
15 Optional switches:
16
17 -# A parameter useful for tuning (default 100)
18 -v Report nonhamiltonian graphs and noncubic graphs
19 -V .. in addition give a cycle for the hamiltonian ones
20 -n#-# If the two numbers are v and i, then the i-th edge
21 out of vertex v is required to be not in the cycle.
22 It must be that i=1..3 and v=0..n-1.
23 -y#-# If the two numbers are v and i, then the i-th edge
24 out of vertex v is required to be in the cycle.
25 It must be that i=1..3 and v=0..n-1.
26 You can use any number of -n/-y switches to force
27 edges. Out of range first arguments are ignored.
28 If -y and -n give same edge, -y wins.
29 -i Test + property: for each edge e, there is a hamiltonian
30 cycle using e.
31 -I Test ++ property: for each pair of edges e,e', there is
32 a hamiltonian cycle which uses both e and e'.
33 -o Test - property: for each edge e, there is a hamiltonian
34 cycle avoiding e.
35 -x Test +- property: for each pair of edges e,e', there is
36 a hamiltonian cycle which uses e but avoids e'.
37 -e Test 3/4 property: for each edge e, at least 3 of the 4
38 paths of length 3 passing through e lie on hamiltonian cycles.
39 -E Test 3/4+ property: for each edge e failing the 3/4 property,
40 all three ways of joining e to the rest of the graph are
41 hamiltonian avoiding e.
42 -T# Specify a timeout, being a limit on how many search tree
43 nodes are made. If the timeout occurs, the graph is
44 written to the output as if it is nonhamiltonian.
45 -R# Specify the number of repeat attempts for each stage.
46 -F Analyze covering paths from 2 or 4 vertices of degree 2.
47
48 -b Require biconnectivity
49 -t Require triconnectivity (note: quadratic algorithm)
50
51 -y, -n, -#, -R and -T are ignored for -i, -I, -x, -o, -e, -E, -F
52
53 B. D. McKay, Nov 1995 + Aug 1996 + Feb 2002 + Jul 2008 + Nov 2015
54
55 **************************************************************************/
56
57 #ifndef MAXN
58 #define MAXN 30002 /* 2 more than largest graph size! */
59 #endif
60
61 #if MAXN==0
62 #error MAXN cannot be zero for cubhamg
63 #endif
64
65 #include "gtools.h"
66 #include "naurng.h"
67
68 /**************************************************************************/
69
70 #define RANPERM TRUE
71
72 /* cubham.h */
73
74 #define MAXNE ((3 * MAXN) / 2)
75 #define FALSE 0
76 #define TRUE 1
77 #define YES 1
78 #define DUNNO 0
79 #define NO (-1)
80 #define HABORT 5
81
82 #define BADLIM 100 /* limit on number of bad things to report per graph */
83 #define MAXA 100 /* max number of -y or -n switches */
84
85 typedef int cubgraph[MAXN][4];
86 typedef int vertvec[MAXN];
87 typedef int edgevec[MAXNE+1];
88 static long nodecount,maxnodes,totalnodes;
89 static long timeout;
90 static long repeats;
91 static int verbose;
92
93 #define NO_LIMIT 0x7FFFFFFFL
94
95 static long standard[]
96 = {30,40,50,60,100,200,300,400,500,1000,2000,3000,5000,
97 10000,20000,30000,100000,300000,1000000};
98 #define NUMMAXNODES (sizeof(standard)/sizeof(long))
99 static nauty_counter numtries[NUMMAXNODES+1];
100
101 /* cubham.c */
102
103 /***************************************************************************
104
105 ***** Instructions for using cubham. *****
106
107 cubham can find hamiltonian cycles in graphs of maximum
108 degree at most 3, with specifed edges either required or forbidden.
109 To use it:
110
111 (1) Check that MAXN (defined in cubham.h) is at least equal to the
112 number of vertices.
113
114 (2) #include cubham.h.
115
116 (3) The cubic graph must be stored in an object of type cubgraph. Say
117 you call it g. The vertices are numbered 0,1,..., and the neighbours of
118 vertex i are g[i][0], g[i][1] and g[i][2]. Entries of the form g[i][3]
119 are unused. If the degree is less than 3, -1 is used as padding.
120
121 (4) Call cubinit(g,eno,v1,v2,nv,&ne), where
122 g = the cubic graph (input)
123 eno = another object of type cubgraph (output)
124 v1,v2 = objects of type edgevec (output)
125 nv = the number of vertices (input)
126 ne = the number of edges (output)
127
128 This numbers the edges 0,1,...ne-1, and defines
129 eno[i][j] = the number of the edge {i, g[i][j]}
130 {v1[j], v2[j]} = the vertices of the j-th edge.
131
132 (5) Call cubham(g,eno,initclass,v1,v2,cycle,outclass,nv,ne), where
133 g,eno,v1,v2,nv = as above (input)
134 initclass = initial edge classification (input)
135 outclass = final edge classifiction (output)
136 cycle = the hamiltonian cycle found, if any (output)
137
138 The value returned by cubham is either YES (hamiltonian cycle found)
139 or NO (there isn't any).
140
141 The initial edge classification is specified by you in the edgevec
142 initclass[]. For each edge j, (0 <= j < 3*nv/2), set
143 initclass[j] = NO, YES or DUNNO, if edge j must not be in the cycle,
144 must be in the cycle, or don't care, respectively. Passing NULL as
145 the initclass parameter is equivalent to setting each edge to DUNNO.
146 All initial classifications are legal, even those clearly impossible.
147
148 The final edge classification is set by cubham in the edgevec
149 outclass[], if a hamiltonian cycle is found. Each entry should be
150 either NO or YES. No final classification is provided if NULL is passed
151 as the outclass parameter.
152
153 The hamiltonian cycle itself, if any, is returned as cycle[0], cycle[1],
154 ..., cycle[nv-1], cycle[0]. If the cycle is not needed, you can pass
155 NULL for this parameter.
156
157 Step (4) need not be repeated if the same graph is processed again with
158 a different initial classification.
159
160 ****************************************************************************/
161
162 #define POP(x) (onstack[x = *(--stackptr)] = 0)
163 #define PUSH(x) if(onstack[x]!=stacklev){onstack[x]=stacklev;*(stackptr++)=x;}
164 #define RESETSTACK {stacklev++; stackptr = stack;}
165
166 typedef struct
167 {
168 edgevec class;
169 vertvec din,dout,farend;
170 } nodedata;
171
172 static nodedata hcnodat;
173 static cubgraph eno;
174 static vertvec onstack,stack; /* stack contains vertex numbers */
175 static int *stackptr,stacklev; /* stackptr points above top */
176 static int classstack[4*MAXNE]; /* stack of classifications */
177 /* x >= 0 : edge number
178 (x < 0 above y : farend[-x-1] = y */
179 static int *classstackptr; /* points above top of classstack */
180 static int classout(cubgraph,nodedata*,int,int,int);
181 static int classin(cubgraph,cubgraph,nodedata*,int,int,int,int*,int);
182
183 #define MAXES 0
184
185 #if MAXES
186 static int maxlevel,maxclassstack;
187 #endif
188
189 static void
dummy(void)190 dummy(void)
191 {
192 }
193
194 static void
check_it(int index,cubgraph g,cubgraph eno,edgevec v1,edgevec v2,int * din,int * dout,int * class,int * farend,int nin,int nv,int stable)195 check_it(int index, cubgraph g, cubgraph eno, edgevec v1, edgevec v2,
196 int *din, int *dout, int *class, int *farend, int nin, int nv,
197 int stable)
198 /* Check some things */
199 {
200 int xdin[MAXN],xdout[MAXN],xnin,v,i,j,k,l,*gv,xin,xout,has1;
201
202 for (i = 0; i < nv; ++i) xdin[i] = xdout[i] = 0;
203 xnin = has1 = 0;
204
205 for (v = 0; v < nv; ++v)
206 {
207 gv = g[v];
208 xin = xout = 0;
209 for (i = 0; i < 3; ++i)
210 if (gv[i] >= 0)
211 {
212 j = eno[v][i];
213 if (class[j] == NO) ++xout;
214 if (class[j] == YES) ++xin;
215 }
216 if (xout != dout[v] || xin != din[v])
217 {
218 fprintf(stderr,">E%d degrees of %d: din,dout=%d,%d really %d,%d\n",
219 index,v,din[v],dout[v],xin,xout);
220 dummy();
221 }
222 if (xin == 1) ++has1;
223 xnin += xin;
224 }
225
226 xnin /= 2;
227 if (xnin != nin)
228 {
229 fprintf(stderr,">E%d nin=%d actually %d\n",index,nin,xnin);
230 dummy();
231 }
232 if (nin != 0 && !has1)
233 {
234 fprintf(stderr,">E%d nin=%d has no in=1\n",index,nin);
235 dummy();
236 }
237
238 for (i = 0; i < nv; ++i)
239 if (din[i] == 0)
240 {
241 if (farend[i] != i)
242 {
243 fprintf(stderr,">E%d farend[isolate %d]=%d\n",
244 index,i,farend[i]);
245 dummy();
246 }
247 }
248 else if (din[i] == 1)
249 {
250 k = -1;
251 j = i;
252 do
253 {
254 for (l = 0; l < 3; ++l)
255 if (g[j][l] >= 0 && g[j][l] != k && class[eno[j][l]] == YES)
256 break;
257 k = j;
258 if (l < 3) j = g[j][l];
259 } while (l < 3);
260 if (farend[i] != j)
261 {
262 fprintf(stderr,">E%d farend[%d]=%d really %d\n",
263 index,i,farend[i],j);
264 dummy();
265 }
266 }
267
268 if (stable)
269 for (i = 0; i < nv; ++i)
270 if ((dout[i] == 1 && din[i] != 2) || (din[i] == 2 && dout[i] != 1)
271 || dout[i] > 1 || din[i] > 2)
272 {
273 fprintf(stderr,">E%d din[%d]=%d dout[%d]=%d\n",
274 index,i,din[i],i,dout[i]);
275 dummy();
276 }
277 }
278
279 static void
cubinit(cubgraph g,cubgraph eno,edgevec v1,edgevec v2,int nv,int ne)280 cubinit(cubgraph g, cubgraph eno, edgevec v1, edgevec v2, int nv, int ne)
281 /* initialise edge numbers, etc. */
282 {
283 int *gpx,*gpy,*enop,x,y,i,j,n,en;
284
285 n = nv;
286 en = 0;
287 for (x = 0; x < n; ++x)
288 {
289 gpx = g[x];
290 enop = eno[x];
291 for (i = 0; i < 3; ++i)
292 if ((y = gpx[i]) < 0)
293 enop[i] = ne;
294 else if (y > x)
295 {
296 v1[en] = x;
297 v2[en] = y;
298 enop[i] = en++;
299 }
300 else
301 {
302 gpy = g[y];
303 for (j = 0; gpy[j] != x; j++)
304 {}
305 enop[i] = eno[y][j];
306 }
307 }
308
309 if (en != ne)
310 fprintf(stderr,"%% cubinit got en=%d when ne=%d\n",en,ne);
311 }
312
313 static int
propagate(cubgraph g,cubgraph eno,nodedata * ndptr,int * nin,int nv)314 propagate(cubgraph g, cubgraph eno, nodedata *ndptr, int *nin, int nv)
315 /* propagate classifications: */
316 /* ans = YES, NO or DUNNO */
317 {
318 int v,w,i,status;
319 nodedata *np;
320 int *gp,*enop,*class,*din,*dout;
321
322 status = DUNNO;
323 np = ndptr;
324 class = np->class;
325 din = np->din;
326 dout = np->dout;
327
328 while (status == DUNNO && stackptr > stack)
329 {
330 POP(v);
331 gp = g[v];
332 enop = eno[v];
333 if (dout[v] == 0)
334 {
335 if (din[v] == 2)
336 {
337 if (class[enop[0]] == DUNNO) i = 0;
338 else if (class[enop[1]] == DUNNO) i = 1;
339 else i = 2;
340 w = gp[i];
341 status = classout(g,np,v,w,enop[i]);
342 PUSH(w);
343 }
344 else if (din[v] == 3)
345 status = NO;
346 }
347 else if (dout[v] == 1)
348 {
349 for (i = 0; i < 3; ++i)
350 if (class[enop[i]] == DUNNO)
351 {
352 w = gp[i];
353 if ((status = classin(g,eno,np,v,w,enop[i],nin,nv))
354 != DUNNO)
355 break;
356 else
357 PUSH(w);
358 }
359 }
360 else
361 status = NO;
362 }
363
364 if (status != NO && *nin == nv) return YES;
365 else return status;
366 }
367
368 static int
classout(cubgraph g,nodedata * nodat,int v,int w,int en)369 classout(cubgraph g, nodedata *nodat, int v, int w, int en)
370 /* classify edge en = vw out */
371 {
372 nodedata *np;
373
374 np = nodat;
375 ++np->dout[v];
376 ++np->dout[w];
377 np->class[en] = NO;
378 *classstackptr++ = en;
379 #if MAXES
380 if (classstackptr-classstack > maxclassstack)
381 maxclassstack = classstackptr-classstack;
382 #endif
383
384 return DUNNO;
385 }
386
387 static int
classin(cubgraph g,cubgraph eno,nodedata * nodat,int v,int w,int en,int * nin,int nv)388 classin(cubgraph g, cubgraph eno, nodedata *nodat,
389 int v, int w, int en, int *nin, int nv)
390 /* classify edge en = vw in */
391 {
392 nodedata *np;
393 int *farend,*gp,fv,fw,i;
394
395 np = nodat;
396 ++np->din[v];
397 ++np->din[w];
398 np->class[en] = YES;
399 *classstackptr++ = en;
400 #if MAXES
401 if (classstackptr-classstack > maxclassstack)
402 maxclassstack = classstackptr-classstack;
403 #endif
404
405 ++*nin;
406 if (*nin == nv)
407 {
408 return DUNNO;
409 }
410
411 farend = np->farend;
412 fv = farend[v];
413 fw = farend[w];
414 *classstackptr++ = farend[fv];
415 *classstackptr++ = -fv-1;
416 *classstackptr++ = farend[fw];
417 *classstackptr++ = -fw-1;
418 #if MAXES
419 if (classstackptr-classstack > maxclassstack)
420 maxclassstack = classstackptr-classstack;
421 #endif
422 farend[fv] = fw;
423 farend[fw] = fv;
424
425 gp = g[fv];
426 if (gp[0] == fw) i = 0;
427 else if (gp[1] == fw) i = 1;
428 else if (gp[2] == fw) i = 2;
429 else return DUNNO;
430
431 i = eno[fv][i];
432 if (np->class[i] == DUNNO)
433 {
434 PUSH(fv);
435 PUSH(fw);
436 if (*nin == nv - 1)
437 return classin(g,eno,np,fv,fw,i,nin,nv);
438 else
439 return classout(g,np,fv,fw,i);
440 }
441
442 return DUNNO;
443 }
444
445 static int
hamnode(cubgraph g,cubgraph eno,edgevec v1,edgevec v2,nodedata * nodat,int level,int nin,int nv)446 hamnode(cubgraph g, cubgraph eno, edgevec v1, edgevec v2,
447 nodedata *nodat, int level, int nin, int nv)
448 /* main node for recursion */
449 {
450 int i,p,q,status;
451 int v,w,en,*gv,*enov;
452 int *csptr;
453
454 #if MAXES
455 if (level > maxlevel) maxlevel = level;
456 #endif
457
458 if (++nodecount > maxnodes && maxnodes != NO_LIMIT) return HABORT;
459 status = propagate(g,eno,nodat,&nin,nv);
460
461 if (status != DUNNO) return status;
462
463 for (v = nv; --v >= 0;)
464 if (nodat->din[v] == 1) break;
465
466 if (v < 0) v = 0;
467
468 gv = g[v];
469 enov = eno[v];
470
471 for (i = 0; i < 3; ++i)
472 {
473 en = enov[i];
474 if (nodat->class[en] == DUNNO)
475 {
476 w = gv[i];
477 csptr = classstackptr;
478 status = classout(g,nodat,v,w,en);
479 if (status == YES) break;
480 if (status == NO)
481 {
482 while (classstackptr > csptr)
483 {
484 p = *--classstackptr;
485 if (p >= 0)
486 {
487 if (nodat->class[p] == YES)
488 {
489 --nodat->din[v1[p]];
490 --nodat->din[v2[p]];
491 }
492 else
493 {
494 --nodat->dout[v1[p]];
495 --nodat->dout[v2[p]];
496 }
497 nodat->class[p] = DUNNO;
498 }
499 else
500 {
501 q = *--classstackptr;
502 nodat->farend[-p-1] = q;
503 }
504 }
505 continue;
506 }
507 RESETSTACK;
508 PUSH(v);
509 PUSH(w);
510 status = hamnode(g,eno,v1,v2,nodat,level+1,nin,nv);
511 if (status == YES) break;
512 while (classstackptr > csptr)
513 {
514 p = *--classstackptr;
515 if (p >= 0)
516 {
517 if (nodat->class[p] == YES)
518 {
519 --nodat->din[v1[p]];
520 --nodat->din[v2[p]];
521 }
522 else
523 {
524 --nodat->dout[v1[p]];
525 --nodat->dout[v2[p]];
526 }
527 nodat->class[p] = DUNNO;
528 }
529 else
530 {
531 q = *--classstackptr;
532 nodat->farend[-p-1] = q;
533 }
534 }
535 if (status == HABORT) return HABORT;
536 }
537 }
538
539 if (status == DUNNO)
540 fprintf(stderr,"hamnode returning DUNNO, this can't happen\n");
541 return status;
542 }
543
544 static int
cubham(cubgraph g,cubgraph eno,edgevec initclass,edgevec v1,edgevec v2,vertvec cycle,edgevec outclass,int nv,int ne)545 cubham(cubgraph g, cubgraph eno, edgevec initclass, edgevec v1, edgevec v2,
546 vertvec cycle, edgevec outclass, int nv, int ne)
547 /* external interface */
548 {
549 int i,j,status,nin,v,w;
550
551 for (i = ne; --i >= 0;)
552 hcnodat.class[i] = DUNNO;
553 if (3*nv > 2*ne) hcnodat.class[ne] = NO;
554
555 for (i = nv; --i >= 0;)
556 {
557 hcnodat.din[i] = hcnodat.dout[i] = 0;
558 hcnodat.farend[i] = i;
559 onstack[i] = 0;
560 }
561 nin = 0;
562 stacklev = 0;
563 RESETSTACK;
564
565 for (i = nv; --i >= 0;)
566 {
567 if (g[i][1] < 0) return NO;
568 if (g[i][2] < 0)
569 {
570 hcnodat.dout[i] = 1;
571 PUSH(i);
572 }
573 }
574
575 status = DUNNO;
576 classstackptr = classstack;
577 if (initclass)
578 for (i = 0; i < ne; ++i)
579 if (initclass[i] != DUNNO)
580 {
581 v = v1[i];
582 w = v2[i];
583 if (initclass[i] == NO)
584 {
585 if (hcnodat.class[i] == YES)
586 status = NO;
587 else if (hcnodat.class[i] == DUNNO)
588 {
589 if (hcnodat.dout[v] == 0)
590 {
591 status = classout(g,&hcnodat,v,w,i);
592 PUSH(v);
593 PUSH(w);
594 }
595 else
596 status = NO;
597 }
598 }
599 else if (initclass[i] == YES)
600 {
601 if (hcnodat.class[i] == NO)
602 status = NO;
603 else if (hcnodat.class[i] == DUNNO)
604 {
605 if (hcnodat.din[v] < 2)
606 {
607 status = classin(g,eno,&hcnodat,v,w,i,&nin,nv);
608 PUSH(v);
609 PUSH(w);
610 }
611 else
612 status = NO;
613 }
614 }
615
616 if (status != DUNNO) break;
617 }
618
619 if (status == DUNNO)
620 status = hamnode(g,eno,v1,v2,&hcnodat,0,nin,nv);
621
622 if (status == YES && cycle)
623 {
624 w = -1;
625 v = 0;
626 cycle[0] = 0;
627 for (i = 1; i < nv; ++i)
628 {
629 for (j = 0; g[v][j] == w || hcnodat.class[eno[v][j]] != YES; ++j)
630 {}
631 w = v;
632 v = g[v][j];
633 cycle[i] = v;
634 }
635 }
636 if (status == YES && outclass)
637 for (i = 0; i < ne; ++i)
638 outclass[i] = hcnodat.class[i];
639
640 return status;
641 }
642
643 /********************************************************************/
644
645 static int
isham(cubgraph cub,int n,int ne,int weight,int * vv,int * vi,int nvv,int * yy,int * yi,int nyy,int * cyc)646 isham(cubgraph cub,
647 int n, int ne, int weight,
648 int *vv, int *vi, int nvv,
649 int *yy, int *yi, int nyy, int *cyc)
650 /* test if hamiltonian; optionally return a cycle
651 Forbid the vi[i]-th nbr of vv[i], for i=0..nvv-1
652 Force the yi[i]-th nbr of yy[i], for i=0..nyy-1
653 WARNING: vi[i]/yi[i] is numbered starting at 1 */
654 {
655 int i,j,k;
656 int nmax,ch;
657 cubgraph cubcopy;
658 edgevec v1,v2,initclass,outclass;
659 int perm[MAXN],pinv[MAXN];
660 double tmp;
661
662 #if !RANPERM
663
664 maxnodes = NO_LIMIT;
665 nodecount = 0;
666 cubinit(cub,eno,v1,v2,n,ne);
667
668 for (i = 0; i < ne; ++i)
669 initclass[i] = DUNNO;
670
671 for (i = 0; i < nvv; ++i)
672 if (vv[i] < n) initclass[eno[vv[i]][vi[i]-1]] = NO;
673 for (i = 0; i < nyy; ++i)
674 if (yy[i] < n) initclass[eno[yy[i]][yi[i]-1]] = YES;
675
676 ch = cubham(cub,eno,initclass,v1,v2,cyc,outclass,n,ne);
677
678 totalnodes += nodecount;
679 ++numtries[0];
680
681 #else
682 ch = HABORT;
683 maxnodes = -1;
684 for (nmax = 0; ch == HABORT && maxnodes != timeout; ++nmax)
685 {
686 if (nmax/repeats < NUMMAXNODES)
687 {
688 tmp = (double)standard[nmax/repeats] * (double)weight
689 * (double)n/ 10000.0;
690 if (tmp >= (double)NO_LIMIT) maxnodes = NO_LIMIT;
691 else maxnodes = tmp;
692 if (timeout > 0 && timeout < maxnodes) maxnodes = timeout;
693 }
694 else if (timeout > 0)
695 maxnodes = timeout;
696 else
697 maxnodes = NO_LIMIT;
698
699 if (nmax != 0)
700 {
701 for (i = n; --i > 0;)
702 {
703 k = KRAN(i+1);
704 j = perm[i];
705 perm[i] = perm[k];
706 perm[k] = j;
707 }
708 }
709 else
710 {
711 for (i = 0; i < n; ++i)
712 perm[i] = i;
713 }
714
715 for (i = 0; i < n; ++i)
716 {
717 j = perm[i];
718 cubcopy[j][0] = cub[i][0] < 0 ? -1 : perm[cub[i][0]];
719 cubcopy[j][1] = cub[i][1] < 0 ? -1 : perm[cub[i][1]];
720 cubcopy[j][2] = cub[i][2] < 0 ? -1 : perm[cub[i][2]];
721 }
722
723 cubinit(cubcopy,eno,v1,v2,n,ne);
724 nodecount = 0;
725
726 for (i = 0; i < ne; ++i)
727 initclass[i] = DUNNO;
728
729 for (i = 0; i < nvv; ++i)
730 if (vv[i] < n) initclass[eno[perm[vv[i]]][vi[i]-1]] = NO;
731 for (i = 0; i < nyy; ++i)
732 if (yy[i] < n) initclass[eno[perm[yy[i]]][yi[i]-1]] = YES;
733
734 ch = cubham(cubcopy,eno,initclass,v1,v2,cyc,outclass,n,ne);
735 totalnodes += nodecount;
736 ++numtries[nmax/repeats];
737 }
738
739 if (cyc != NULL && ch == YES)
740 {
741 for (i = 0; i < n; ++i)
742 pinv[perm[i]] = i;
743 for (i = 0; i < n; ++i)
744 cyc[i] = pinv[cyc[i]];
745 }
746 #endif
747
748 return ch;
749 }
750
751 /**************************************************************************/
752
753 static int
optadd(cubgraph cub,int v1,int v2)754 optadd(cubgraph cub, int v1, int v2)
755 /* v1 and v2 must have degree 2 and be distinct.
756 Add edge v1-v2 if not already present; return index of edge in cub[v1]. */
757 {
758 if (cub[v1][0] == v2) return 0;
759 if (cub[v1][1] == v2) return 1;
760 cub[v1][2] = v2;
761 cub[v2][2] = v1;
762 return 2;
763 }
764
765 /**************************************************************************/
766
767 static void
dofragment(nauty_counter id,cubgraph cub,int n,int ne,int weight)768 dofragment(nauty_counter id, cubgraph cub, int n, int ne, int weight)
769 /* Test for coverage by one or two paths between vertices of degree 2 */
770 {
771 int i,i1,i2,i3,i4;
772 int v1,v2,v3,v4,j1,j3;
773 int deg2[MAXN],ndeg2;
774 int yy[3],yi[3],newne;
775 int cyc[MAXN];
776 int status;
777
778 ndeg2 = 0;
779 for (i = 0; i < n; ++i)
780 {
781 if (cub[i][0] < 0 || cub[i][1] < 0)
782 gt_abort(">E -F forbids degree 0,1\n");
783
784 if (cub[i][2] < 0) deg2[ndeg2++] = i;
785 }
786
787 printf("Input " COUNTER_FMT ":",id);
788 for (i = 0; i < ndeg2; ++i) printf(" %d",deg2[i]);
789 printf("\n");
790
791 printf(" Pairs: ");
792 for (i1 = 0; i1 < ndeg2; ++i1)
793 for (i2 = i1+1; i2 < ndeg2; ++i2)
794 {
795 v1 = deg2[i1]; v2 = deg2[i2];
796 j1 = optadd(cub,v1,v2);
797 yy[0] = v1; yi[0] = j1+1;
798 newne = ne + (j1==2);
799 status = isham(cub,n,newne,weight,NULL,NULL,0,yy,yi,1,cyc);
800 if (status == HABORT)
801 printf(" T%d-%d",v1,v2);
802 if (status == NO)
803 printf(" N%d-%d",v1,v2);
804 else
805 {
806 printf(" Y%d-%d",v1,v2);
807 if (verbose > 1)
808 {
809 printf("[");
810 for (i = 0; i < n; ++i) printf(" %d",cyc[i]);
811 printf("]\n");
812 }
813 }
814 cub[v1][2] = cub[v2][2] = -1;
815 }
816 printf("\n");
817
818 printf(" Quartets: ");
819 for (i1 = 0; i1 < ndeg2; ++i1)
820 for (i2 = i1+1; i2 < ndeg2; ++i2)
821 for (i3 = i1+1; i3 < ndeg2; ++i3)
822 for (i4 = i3+1; i4 < ndeg2; ++i4)
823 {
824 if (i3 == i2 || i4 == i2) continue;
825
826 v1 = deg2[i1]; v2 = deg2[i2];
827 j1 = optadd(cub,v1,v2);
828 v3 = deg2[i3]; v4 = deg2[i4];
829 j3 = optadd(cub,v3,v4);
830 yy[0] = v1; yi[0] = j1+1;
831 yy[1] = v3; yi[1] = j3+1;
832 newne = ne + (j1==2) + (j3==2);
833 status = isham(cub,n,newne,weight,NULL,NULL,0,yy,yi,2,cyc);
834 if (status == HABORT)
835 printf(" T%d-%d,%d-%d",v1,v2,v3,v4);
836 if (status == NO)
837 printf(" N%d-%d,%d-%d",v1,v2,v3,v4);
838 else
839 {
840 printf(" Y%d-%d,%d-%d",v1,v2,v3,v4);
841 if (verbose > 1)
842 {
843 printf("[");
844 for (i = 0; i < n; ++i) printf(" %d",cyc[i]);
845 printf("]\n");
846 }
847 }
848 cub[v1][2] = cub[v2][2] = -1;
849 cub[v3][2] = cub[v4][2] = -1;
850 }
851 printf("\n");
852 }
853
854 /********************************************************************/
855
856 static int
hasinout(cubgraph cub,int n,int ne,int * x0,int * x1,int * y0,int * y1,int limit)857 hasinout(cubgraph cub,
858 int n, int ne, int *x0, int *x1, int *y0, int *y1, int limit)
859 /* test if cub has in-out (+-) property */
860 {
861 edgevec v1,v2,initclass,outclass;
862 set *d0,*di,*dii;
863 int i,ii,j,jj;
864 int me,nbad;
865 DYNALLSTAT(graph,done,done_sz);
866
867 me = (ne + WORDSIZE - 1) / WORDSIZE;
868
869 DYNALLOC2(graph,done,done_sz,ne,me,"hasinout");
870
871 d0 = (set*)done;
872 EMPTYSET(d0,me);
873 for (j = 0; j < ne; ++j)
874 ADDELEMENT(d0,j);
875
876 for (i = 1, di = d0 + me; i < ne; ++i, di += me)
877 {
878 for (j = 0; j < me; ++j)
879 di[j] = d0[j];
880 DELELEMENT(di,i);
881 }
882 DELELEMENT(d0,0);
883
884 cubinit(cub,eno,v1,v2,n,ne);
885 for (i = 0; i < ne; ++i)
886 initclass[i] = DUNNO;
887
888 maxnodes = NO_LIMIT;
889 nbad = 0;
890
891 for (i = 0, di = (set*)done; i < ne; ++i, di += me)
892 for (j = -1; (j = nextelement(di,me,j)) >= 0;)
893 {
894 initclass[i] = NO;
895 initclass[j] = YES;
896 ++numtries[0];
897 if (cubham(cub,eno,initclass,v1,v2,NULL,
898 outclass,n,ne) == NO)
899 {
900 x0[nbad] = v1[i]; x1[nbad] = v2[i];
901 y0[nbad] = v1[j]; y1[nbad] = v2[j];
902 ++nbad;
903 if (nbad >= limit) return nbad;
904 }
905 else
906 {
907 for (ii = i, dii = di; ii < ne; ++ii, dii += me)
908 if (outclass[ii] == NO)
909 for (jj = 0; jj < ne; ++jj)
910 if (outclass[jj] == YES)
911 DELELEMENT(dii,jj);
912 }
913
914 initclass[i] = DUNNO;
915 initclass[j] = DUNNO;
916 }
917
918 return nbad;
919 }
920
921 /**************************************************************************/
922
923 static int
hasinin(cubgraph cub,int n,int ne,int * x0,int * x1,int * y0,int * y1,int limit)924 hasinin(cubgraph cub,
925 int n, int ne, int *x0, int *x1, int *y0, int *y1, int limit)
926 /* test if cub has in-in (++) property */
927 {
928 edgevec v1,v2,initclass,outclass;
929 set *d0,*di,*dii;
930 int i,ii,j,jj;
931 int me,nbad;
932 DYNALLSTAT(graph,done,done_sz);
933
934 me = (ne + WORDSIZE - 1) / WORDSIZE;
935
936 DYNALLOC2(graph,done,done_sz,ne,me,"hasinin");
937
938 d0 = (set*)done;
939 EMPTYSET(d0,me);
940 for (j = 0; j < ne; ++j)
941 ADDELEMENT(d0,j);
942
943 for (i = 1, di = d0 + me; i < ne; ++i, di += me)
944 {
945 for (j = 0; j < me; ++j)
946 di[j] = d0[j];
947 DELELEMENT(di,i);
948 }
949 DELELEMENT(d0,0);
950
951 cubinit(cub,eno,v1,v2,n,ne);
952 for (i = 0; i < ne; ++i)
953 initclass[i] = DUNNO;
954
955 maxnodes = NO_LIMIT;
956 nbad = 0;
957
958 for (i = 0, di = (set*)done; i < ne; ++i, di += me)
959 for (j = i; (j = nextelement(di,me,j)) >= 0;)
960 {
961 initclass[i] = YES;
962 initclass[j] = YES;
963 ++numtries[0];
964 if (cubham(cub,eno,initclass,v1,v2,NULL,outclass,n,ne) == NO)
965 {
966 x0[nbad] = v1[i]; x1[nbad] = v2[i];
967 y0[nbad] = v1[j]; y1[nbad] = v2[j];
968 ++nbad;
969 if (nbad >= limit) return nbad;
970 }
971 else
972 {
973 for (ii = i, dii = di; ii < ne; ++ii, dii += me)
974 if (outclass[ii] == YES)
975 for (jj = ii; jj < ne; ++jj)
976 if (outclass[jj] == YES)
977 DELELEMENT(dii,jj);
978 }
979
980 initclass[i] = DUNNO;
981 initclass[j] = DUNNO;
982 }
983 return nbad;
984 }
985
986 /**************************************************************************/
987
988 static int
hasin(cubgraph cub,int n,int ne,int * x0,int * x1,int limit)989 hasin(cubgraph cub, int n, int ne, int *x0, int *x1, int limit)
990 /* test if cub has "in" property */
991 {
992 edgevec v1,v2,initclass,outclass;
993 boolean done[MAXNE];
994 int i,ii;
995 int nbad;
996
997 cubinit(cub,eno,v1,v2,n,ne);
998 for (i = 0; i < ne; ++i)
999 {
1000 initclass[i] = DUNNO;
1001 done[i] = FALSE;
1002 }
1003
1004 maxnodes = NO_LIMIT;
1005 nbad = 0;
1006
1007 for (i = 0; i < ne; ++i)
1008 if (!done[i])
1009 {
1010 initclass[i] = YES;
1011 ++numtries[0];
1012 if (cubham(cub,eno,initclass,v1,v2,NULL,outclass,n,ne) == NO)
1013 {
1014 x0[nbad] = v1[i]; x1[nbad] = v2[i];
1015 ++nbad;
1016 if (nbad >= limit) return nbad;
1017 }
1018 else
1019 {
1020 for (ii = i; ii < ne; ++ii)
1021 if (outclass[ii] == YES) done[ii] = TRUE;
1022 }
1023
1024 initclass[i] = DUNNO;
1025 }
1026 return nbad;
1027 }
1028
1029 /**************************************************************************/
1030
1031 static boolean
eplus(cubgraph acub,int n,int ne,int x,int y,int * pwhy)1032 eplus(cubgraph acub, int n, int ne, int x, int y, int *pwhy)
1033 {
1034 cubgraph cub;
1035 edgevec v1,v2,initclass;
1036 int i,a,b,c,d,xy,why;
1037
1038 if (3*n != 2*ne)
1039 {
1040 fprintf(stderr,
1041 "cubhamg: eplus() not implemented for noncubic graphs\n");
1042 exit(1);
1043 }
1044
1045 for (i = 0; i < n; ++i)
1046 {
1047 cub[i][0] = acub[i][0];
1048 cub[i][1] = acub[i][1];
1049 cub[i][2] = acub[i][2];
1050 }
1051
1052 if (cub[x][0] == y)
1053 {
1054 xy = 0;
1055 a = cub[x][1];
1056 b = cub[x][2];
1057 }
1058 else if (cub[x][1] == y)
1059 {
1060 xy = 1;
1061 a = cub[x][0];
1062 b = cub[x][2];
1063 }
1064 else
1065 {
1066 xy = 2;
1067 a = cub[x][0];
1068 b = cub[x][1];
1069 }
1070
1071 if (cub[y][0] == x)
1072 {
1073 c = cub[y][1];
1074 d = cub[y][2];
1075 }
1076 else if (cub[y][1] == x)
1077 {
1078 c = cub[y][0];
1079 d = cub[y][2];
1080 }
1081 else
1082 {
1083 c = cub[y][0];
1084 d = cub[y][1];
1085 }
1086
1087 for (i = 0; i < ne; ++i)
1088 initclass[i] = DUNNO;
1089
1090 why = 0;
1091 cubinit(cub,eno,v1,v2,n,ne);
1092 initclass[eno[x][xy]] = NO;
1093 if (cubham(cub,eno,initclass,v1,v2,NULL,NULL,n,ne) == YES)
1094 why |= 1;
1095 initclass[eno[x][xy]] = DUNNO;
1096
1097 #define CHANGE(z,p,q) if (cub[z][0]==p) cub[z][0]=q;\
1098 else if (cub[z][1]==p) cub[z][1]=q; else cub[z][2]=q;
1099
1100 CHANGE(x,b,c);
1101 CHANGE(y,c,b);
1102 CHANGE(b,x,y);
1103 CHANGE(c,y,x);
1104
1105 cubinit(cub,eno,v1,v2,n,ne);
1106 initclass[eno[x][xy]] = NO;
1107 if (cubham(cub,eno,initclass,v1,v2,NULL,NULL,n,ne) == YES)
1108 why |= 2;
1109 initclass[eno[x][xy]] = DUNNO;
1110
1111 CHANGE(x,c,d);
1112 CHANGE(y,d,c);
1113 CHANGE(c,x,y);
1114 CHANGE(d,y,x);
1115
1116 cubinit(cub,eno,v1,v2,n,ne);
1117 initclass[eno[x][xy]] = NO;
1118 if (cubham(cub,eno,initclass,v1,v2,NULL,NULL,n,ne) == YES)
1119 why |= 4;
1120
1121 *pwhy = why;
1122 return why == 7;
1123 }
1124
1125 /**************************************************************************/
1126
1127 static int
hase34(cubgraph cub,int n,int ne,int * x0,int * x1,int * why,boolean plus,int limit)1128 hase34(cubgraph cub,
1129 int n, int ne, int *x0, int *x1, int *why, boolean plus, int limit)
1130 /* test if cub has "e34" property */
1131 {
1132 edgevec v1,v2,initclass,outclass;
1133 int ea[4*MAXNE],eb[4*MAXNE],ec[4*MAXNE];
1134 boolean done[4*MAXNE];
1135 int count[MAXNE];
1136 int i,ii;
1137 int vi,nbad,pluswhy;
1138 static int pop[] = {0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4};
1139
1140 if (3*n != 2*ne)
1141 {
1142 fprintf(stderr,
1143 "cubhamg: hase34() not implemented for noncubic graphs\n");
1144 exit(1);
1145 }
1146 cubinit(cub,eno,v1,v2,n,ne);
1147 for (i = 0; i < ne; ++i)
1148 {
1149 initclass[i] = DUNNO;
1150 count[i] = 0;
1151 }
1152
1153 for (i = ii = 0; i < ne; ++i)
1154 {
1155 eb[ii] = eb[ii+1] = eb[ii+2] = eb[ii+3] = i;
1156 done[ii] = done[ii+1] = done[ii+2] = done[ii+3] = FALSE;
1157
1158 vi = v1[i];
1159 if (eno[vi][0] != i) ea[ii++] = eno[vi][0];
1160 if (eno[vi][1] != i) ea[ii++] = eno[vi][1];
1161 if (eno[vi][2] != i) ea[ii++] = eno[vi][2];
1162 ea[ii] = ea[ii-2];
1163 ea[ii+1] = ea[ii-1];
1164
1165 vi = v2[i];
1166 if (eno[vi][0] != i) ec[ii++] = eno[vi][0];
1167 if (eno[vi][1] != i) ec[ii++] = eno[vi][1];
1168 if (eno[vi][2] != i) ec[ii++] = eno[vi][2];
1169 ec[ii-4] = ec[ii-3] = ec[ii-2];
1170 ec[ii-2] = ec[ii-1];
1171 }
1172
1173 maxnodes = NO_LIMIT;
1174 nbad = 0;
1175
1176 for (i = 0; i < 4*ne; ++i)
1177 if (!done[i])
1178 {
1179 initclass[ea[i]] = YES;
1180 initclass[eb[i]] = YES;
1181 initclass[ec[i]] = YES;
1182 ++numtries[0];
1183 if (cubham(cub,eno,initclass,v1,v2,NULL,
1184 outclass,n,ne) == YES)
1185 for (ii = i; ii < 4*ne; ++ii)
1186 if (!done[ii] && outclass[ea[ii]] == YES &&
1187 outclass[eb[ii]] == YES && outclass[ec[ii]] == YES)
1188 {
1189 done[ii] = TRUE;
1190 count[eb[ii]] |= 1 << (ii & 3);
1191 }
1192 ++numtries[0];
1193 initclass[ea[i]] = DUNNO;
1194 initclass[eb[i]] = DUNNO;
1195 initclass[ec[i]] = DUNNO;
1196 }
1197
1198 pluswhy = 0;
1199 for (i = 0; i < ne; ++i)
1200 if (pop[count[i]] < 3
1201 && (!plus || !eplus(cub,n,ne,v1[i],v2[i],&pluswhy)))
1202 {
1203 x0[nbad] = v1[i]; x1[nbad] = v2[i];
1204 why[nbad] = (pluswhy << 4) | count[i];
1205 ++nbad;
1206 if (nbad >= limit) return nbad;
1207 }
1208
1209 return nbad;
1210 }
1211
1212 /**************************************************************************/
1213
1214 #if 0
1215 static int
1216 oldhasout(cubgraph cub, int n, int ne, int *x0, int *x1, int limit)
1217 /* test if cub has "out" property */
1218 {
1219 edgevec v1,v2,initclass,outclass;
1220 boolean done[MAXNE];
1221 int i,ii;
1222 int nbad;
1223
1224 cubinit(cub,eno,v1,v2,n,ne);
1225 for (i = 0; i < ne; ++i)
1226 {
1227 initclass[i] = DUNNO;
1228 done[i] = FALSE;
1229 }
1230
1231 maxnodes = NO_LIMIT;
1232 nbad = 0;
1233
1234 for (i = 0; i < ne; ++i)
1235 if (!done[i])
1236 {
1237 initclass[i] = NO;
1238 ++numtries[0];
1239 if (cubham(cub,eno,initclass,v1,v2,NULL,outclass,n,ne) == NO)
1240 {
1241 x0[nbad] = v1[i]; x1[nbad] = v2[i];
1242 ++nbad;
1243 if (nbad >= limit) return nbad;
1244 }
1245 else
1246 {
1247 for (ii = i; ii < ne; ++ii)
1248 if (outclass[ii] == NO) done[ii] = TRUE;
1249 }
1250
1251 initclass[i] = DUNNO;
1252 }
1253 return nbad;
1254 }
1255 #endif
1256
1257 /**************************************************************************/
1258
1259 static int
hasout(cubgraph cub,int n,int ne,int weight,int * x0,int * x1,int limit)1260 hasout(cubgraph cub, int n, int ne, int weight, int *x0, int *x1, int limit)
1261 /* test if cub has "out" property
1262 Returns are -2 for timeout, -1 for nonhamiltonian, otherwise
1263 number of edges not present in any cycle */
1264 {
1265 cubgraph done;
1266 int cyc[MAXN];
1267 int vv[2],vi[2];
1268 int nbad,ch;
1269 int i,j,ii,jj,x,y,z;
1270
1271 nbad = 0;
1272
1273 for (i = 0; i < n; ++i) done[i][0] = done[i][1] = done[i][2] = 0;
1274
1275 ch = isham(cub,n,ne,weight,vv,vi,0,NULL,NULL,0,cyc);
1276 if (ch == HABORT) return -2;
1277 if (ch == NO) return -1;
1278
1279 for (i = 0; i < n; ++i)
1280 {
1281 x = cyc[i];
1282 y = cyc[i==n-1?0:i+1];
1283 z = cyc[i==0?n-1:i-1];
1284 for (j = 0; j < 3; ++j)
1285 if (cub[x][j] >= 0 &&
1286 cub[x][j] != y && cub[x][j] != z) break;
1287 if (j < 3) done[x][j] = 1;
1288 }
1289
1290 for (ii = 0; ii < n; ++ii)
1291 for (jj = 0; jj < 3; ++jj)
1292 {
1293 if (cub[ii][jj] > ii && !done[ii][jj])
1294 {
1295 vv[0] = ii;
1296 vi[0] = jj+1;
1297 ch = isham(cub,n,ne,weight,vv,vi,1,NULL,NULL,0,cyc);
1298 if (ch == HABORT) return -2;
1299 if (ch == NO)
1300 {
1301 x0[nbad] = ii;
1302 x1[nbad] = cub[ii][jj];
1303 ++nbad;
1304 }
1305 else
1306 {
1307 for (i = 0; i < n; ++i)
1308 {
1309 x = cyc[i];
1310 y = cyc[i==n-1?0:i+1];
1311 z = cyc[i==0?n-1:i-1];
1312 for (j = 0; j < 3; ++j)
1313 if (cub[x][j] >= 0 &&
1314 cub[x][j] != y && cub[x][j] != z) break;
1315 if (j < 3) done[x][j] = 1;
1316 }
1317
1318 }
1319 }
1320 }
1321
1322 return nbad;
1323 }
1324
1325 /**************************************************************************/
1326
1327 static boolean
biconnected_cub(cubgraph cub,int n)1328 biconnected_cub(cubgraph cub, int n)
1329 /* test whether cub is biconnected */
1330 {
1331 int i,sp,v,w,x;
1332 set visited[MAXM];
1333 int numvis,num[MAXN],lp[MAXN],stack[MAXN];
1334 int m,*gv;
1335
1336 if (n <= 2) return FALSE;
1337
1338 m = (n + WORDSIZE - 1) / WORDSIZE;
1339 EMPTYSET(visited,m);
1340 ADDELEMENT(visited,0);
1341
1342 stack[0] = 0;
1343 num[0] = 0;
1344 lp[0] = 0;
1345 numvis = 1;
1346 sp = 0;
1347 v = 0;
1348 gv = (int*)cub[v];
1349
1350 for (;;)
1351 {
1352 for (i = 0; i < 3; ++i)
1353 if (gv[i] >= 0 && !ISELEMENT(visited,gv[i])) break;
1354
1355 if (i < 3)
1356 {
1357 w = v;
1358 v = gv[i]; /* visit next child */
1359 stack[++sp] = v;
1360 gv = (int*)cub[v];
1361 ADDELEMENT(visited,v);
1362 lp[v] = num[v] = numvis++;
1363
1364 for (i = 0; i < 3; ++i)
1365 {
1366 x = gv[i];
1367 if (x >= 0 && x != w && ISELEMENT(visited,x))
1368 if (num[x] < lp[v]) lp[v] = num[x];
1369 }
1370 }
1371 else
1372 {
1373 w = v; /* back up to parent */
1374 if (sp <= 1) return numvis == n;
1375 v = stack[--sp];
1376 gv = (int*)cub[v];
1377 if (lp[w] >= num[v]) return FALSE;
1378 if (lp[w] < lp[v]) lp[v] = lp[w];
1379 }
1380 }
1381 }
1382
1383 /**************************************************************************/
1384
1385 static boolean
biconnected_cub_v(cubgraph cub,int vv,int n)1386 biconnected_cub_v(cubgraph cub, int vv, int n)
1387 /* test whether cub-vv is biconnected */
1388 {
1389 int i,sp,v,w,x,start;
1390 set visited[MAXM];
1391 int numvis,num[MAXN],lp[MAXN],stack[MAXN];
1392 int m,*gv;
1393
1394 if (n <= 3) return FALSE;
1395 start = (vv == 0 ? 1 : 0);
1396
1397 m = (n + WORDSIZE - 1) / WORDSIZE;
1398 EMPTYSET(visited,m);
1399 ADDELEMENT(visited,start);
1400 ADDELEMENT(visited,vv);
1401
1402 stack[0] = start;
1403 num[start] = 0;
1404 lp[start] = 0;
1405 numvis = 1;
1406 sp = 0;
1407 v = start;
1408 gv = (int*)cub[v];
1409
1410 for (;;)
1411 {
1412 for (i = 0; i < 3; ++i)
1413 if (gv[i] >= 0 && !ISELEMENT(visited,gv[i])) break;
1414
1415 if (i < 3)
1416 {
1417 w = v;
1418 v = gv[i]; /* visit next child */
1419 stack[++sp] = v;
1420 gv = (int*)cub[v];
1421 ADDELEMENT(visited,v);
1422 lp[v] = num[v] = numvis++;
1423
1424 for (i = 0; i < 3; ++i)
1425 {
1426 x = gv[i];
1427 if (x >= 0 && x != w && x != vv && ISELEMENT(visited,x))
1428 if (num[x] < lp[v]) lp[v] = num[x];
1429 }
1430 }
1431 else
1432 {
1433 w = v; /* back up to parent */
1434 if (sp <= 1) return numvis == n-1;
1435 v = stack[--sp];
1436 gv = (int*)cub[v];
1437 if (lp[w] >= num[v]) return FALSE;
1438 if (lp[w] < lp[v]) lp[v] = lp[w];
1439 }
1440 }
1441 }
1442
1443 /**************************************************************************/
1444
1445 static boolean
biconnected_v(graph * g,int vv,int m,int n)1446 biconnected_v(graph *g, int vv, int m, int n)
1447 /* test whether g-vv is biconnected */
1448 /* version for arbitrary sizes */
1449 {
1450 int i,sp,v,w;
1451 setword ww;
1452 set sw[MAXM],visited[MAXM];
1453 int numvis,num[MAXN],lp[MAXN],stack[MAXN];
1454 int start;
1455 set *gv;
1456
1457 if (n <= 3) return FALSE;
1458
1459 start = vv == 0 ? 1 : 0;
1460 EMPTYSET(visited,m);
1461 ADDELEMENT(visited,start);
1462 ADDELEMENT(visited,vv);
1463
1464 stack[0] = start;
1465 num[start] = 0;
1466 lp[start] = 0;
1467 numvis = 1;
1468 sp = 0;
1469 v = start;
1470 gv = (set*)g + m*v;
1471
1472 for (;;)
1473 {
1474 for (i = 0; i < m; ++i)
1475 if ((ww = gv[i] & ~visited[i])) break; /* = */
1476
1477 if (i < m)
1478 {
1479 w = v;
1480 v = TIMESWORDSIZE(i) + FIRSTBIT(ww); /* visit next child */
1481 stack[++sp] = v;
1482 gv = (set*)g + m*v;
1483 ADDELEMENT(visited,v);
1484 lp[v] = num[v] = numvis++;
1485 for (i = 0; i < m; ++i)
1486 sw[i] = gv[i] & visited[i];
1487 DELELEMENT(sw,w);
1488 DELELEMENT(sw,vv);
1489
1490 for (i = 0; i < m; ++i)
1491 {
1492 ww = sw[i];
1493 while (ww)
1494 {
1495 TAKEBIT(w,ww);
1496 w += TIMESWORDSIZE(i);
1497 if (num[w] < lp[v]) lp[v] = num[w];
1498 }
1499 }
1500 }
1501 else
1502 {
1503 w = v; /* back up to parent */
1504 if (sp <= 1) return numvis == n-1;
1505 v = stack[--sp];
1506 gv = (set*)g + m*v;
1507 if (lp[w] >= num[v]) return FALSE;
1508 if (lp[w] < lp[v]) lp[v] = lp[w];
1509 }
1510 }
1511 }
1512
1513 /**************************************************************************/
1514
1515 static boolean
triconnected_cub(cubgraph cub,int n)1516 triconnected_cub(cubgraph cub, int n)
1517 /* test whether triconnected; awfully inefficient */
1518 {
1519 int vv;
1520
1521 if (n <= 3) return FALSE;
1522
1523 for (vv = n-1; --vv >= 0;)
1524 if (!biconnected_cub_v(cub,vv,n)) return FALSE;
1525
1526 return TRUE;
1527 }
1528
1529 /**************************************************************************/
1530
1531 static boolean
gtocub(graph * g,int m,int n,cubgraph cub,int * pne)1532 gtocub(graph *g, int m, int n, cubgraph cub, int *pne)
1533 /* Convert nauty-format graph into cubgraph. Returns FALSE if there
1534 * are any vertices of degree 4 or more.
1535 */
1536 {
1537 int i,j,nde;
1538 set *gi;
1539
1540 nde = 0;
1541 for (i = 0, gi = (set*)g; i < n; ++i, gi += m)
1542 {
1543 cub[i][0] = cub[i][1] = cub[i][2] = -1;
1544 j = nextelement(gi,m,-1);
1545 if (j < 0) continue;
1546 cub[i][0] = j;
1547 ++nde;
1548 j = nextelement(gi,m,j);
1549 if (j < 0) continue;
1550 cub[i][1] = j;
1551 ++nde;
1552 j = nextelement(gi,m,j);
1553 if (j < 0) continue;
1554 cub[i][2] = j;
1555 ++nde;
1556 if (nextelement(gi,m,j) >= 0) return FALSE;
1557 }
1558
1559 *pne = nde / 2;
1560 return TRUE;
1561 }
1562
1563 /**************************************************************************/
1564
1565 static boolean
sgtocub(sparsegraph * sg,cubgraph cub,int * pne)1566 sgtocub(sparsegraph *sg, cubgraph cub, int *pne)
1567 /* Convert sparse-format graph into cubgraph. Returns FALSE if there
1568 * are any vertices of degree 4 or more.
1569 */
1570 {
1571 int *d,*e;
1572 size_t *v;
1573 int n,i,j,vi;
1574 unsigned long int nde;
1575
1576 nde = 0;
1577 n = sg->nv;
1578 SG_VDE(sg,v,d,e);
1579 for (i = 0; i < n; ++i)
1580 {
1581 if (d[i] >= 4) return FALSE;
1582 vi = v[i];
1583
1584 cub[i][0] = cub[i][1] = cub[i][2] = -1;
1585 for (j = 0; j < d[i]; ++j)
1586 cub[i][j] = e[vi+j];
1587 nde += d[i];
1588 }
1589
1590 *pne = nde / 2;
1591 return TRUE;
1592 }
1593
1594 /**************************************************************************/
1595
1596 int
main(int argc,char * argv[])1597 main(int argc, char *argv[])
1598 {
1599 char *infilename,*outfilename;
1600 FILE *infile,*outfile,*msgfile;
1601 boolean badargs,biconn,triconn,fragment;
1602 boolean in,out,inin,inout,e34,e34plus,testing;
1603 nauty_counter numread,noncub,nonham,nonconn,numto;
1604 int ch,weight,m,n,ne,i,namesgot,nl;
1605 int nbad,limit,x0[BADLIM],x1[BADLIM],y0[BADLIM],y1[BADLIM];
1606 sparsegraph sg;
1607 int vv[MAXA],vi[MAXA],nvv,cyc[MAXN];
1608 int yy[MAXA],yi[MAXA],nyy;
1609 double t0,t1;
1610 cubgraph cub;
1611 char *arg;
1612 int codetype;
1613
1614 infilename = outfilename = NULL;
1615 badargs = FALSE;
1616 e34plus = e34 = in = out = inin = inout = FALSE;
1617 fragment = biconn = triconn = testing = FALSE;
1618 verbose = 0;
1619 weight = 100;
1620 nvv = nyy = 0;
1621 timeout = 0;
1622 repeats = 1;
1623
1624 /* parse argument list */
1625
1626 namesgot = 0;
1627 for (i = 1; i < argc && !badargs; ++i)
1628 {
1629 arg = argv[i];
1630 if (arg[0] == '-' && arg[1] != '\0')
1631 {
1632 if (arg[1] == 'v')
1633 verbose = (verbose == 0 ? 1 : verbose);
1634 else if (arg[1] == 'V')
1635 verbose = 2;
1636 else if (arg[1] == 'i')
1637 in = TRUE;
1638 else if (arg[1] == 'I')
1639 inin = TRUE;
1640 else if (arg[1] == 'o')
1641 out = TRUE;
1642 else if (arg[1] == 'x')
1643 inout = TRUE;
1644 else if (arg[1] == 'e')
1645 e34 = TRUE;
1646 else if (arg[1] == 'E')
1647 e34 = e34plus = TRUE;
1648 else if (arg[1] == 'b')
1649 biconn = TRUE;
1650 else if (arg[1] == 't')
1651 triconn = TRUE;
1652 else if (arg[1] == 'Q')
1653 testing = TRUE;
1654 else if (arg[1] == 'F')
1655 fragment = TRUE;
1656 else if (arg[1] == 'n')
1657 {
1658 if (nvv == MAXA)
1659 {
1660 fprintf(stderr,">E cubhamg : too many -n switches\n");
1661 exit(1);
1662 }
1663 if (sscanf(arg+2,"%d-%d",&vv[nvv],&vi[nvv]) != 2)
1664 badargs = TRUE;
1665 else if (vi[nvv] < 1 || vi[nvv] > 3)
1666 {
1667 fprintf(stderr,
1668 ">E cubhamg : second arg of -n must be 1..3\n");
1669 exit(1);
1670 }
1671 ++nvv;
1672 }
1673 else if (arg[1] == 'y')
1674 {
1675 if (nyy == MAXA)
1676 {
1677 fprintf(stderr,">E cubhamg : too many -y switches\n");
1678 exit(1);
1679 }
1680 if (sscanf(arg+2,"%d-%d",&yy[nyy],&yi[nyy]) != 2)
1681 badargs = TRUE;
1682 else if (yi[nyy] < 1 || yi[nyy] > 3)
1683 {
1684 fprintf(stderr,
1685 ">E cubhamg : second arg of -y must be 1..3\n");
1686 exit(1);
1687 }
1688 ++nyy;
1689 }
1690 else if (arg[1] >= '0' && arg[1] <= '9')
1691 sscanf(arg+1,"%d",&weight);
1692 else if (arg[1] == 'T')
1693 {
1694 if (sscanf(arg+2,"%ld",&timeout) != 1 || timeout < 0)
1695 {
1696 fprintf(stderr,">E cubhamg : bad -T value\n");
1697 exit(1);
1698 }
1699 else if (timeout > NO_LIMIT-1)
1700 timeout = NO_LIMIT - 1;
1701 }
1702 else if (arg[1] == 'R')
1703 {
1704 if (sscanf(arg+2,"%ld",&repeats) != 1 || repeats < 1)
1705 {
1706 fprintf(stderr,">E cubhamg : bad -R value\n");
1707 exit(1);
1708 }
1709 }
1710 else
1711 badargs = TRUE;
1712 }
1713 else
1714 {
1715 if (namesgot == 0)
1716 {
1717 namesgot = 1;
1718 infilename = arg;
1719 }
1720 else if (namesgot == 1)
1721 {
1722 namesgot = 2;
1723 outfilename = arg;
1724 }
1725 else
1726 badargs = TRUE;
1727 }
1728 }
1729
1730 if (badargs)
1731 {
1732 fprintf(stderr,
1733 ">E Usage: cubhamg [-#] [-v | -V] [-n#-#] [-y#-#] [infile [outfile]]\n");
1734 exit(1);
1735 }
1736
1737 /* open input file */
1738
1739 if (infilename && infilename[0] == '-') infilename = NULL;
1740 infile = opengraphfile(infilename,&codetype,FALSE,1);
1741 if (!infile) exit(1);
1742 if (!infilename) infilename = "stdin";
1743
1744 if (infilename == NULL) infilename = "stdin";
1745
1746 NODIGRAPHSYET(codetype);
1747
1748 /* open output file */
1749
1750 if (outfilename == NULL || outfilename[0] == '-')
1751 {
1752 outfile = stdout;
1753 outfilename = "stdout";
1754 msgfile = stderr;
1755 }
1756 else
1757 {
1758 msgfile = stdout;
1759 if ((outfile = fopen(outfilename,"w")) == NULL)
1760 {
1761 fprintf(stderr,
1762 ">E cubhamg: can't open %s for writing\n",outfilename);
1763 ABORT(">E cubhamg");
1764 }
1765 }
1766
1767 if (triconn) biconn = FALSE;
1768 if (e34) in = out = inin = inout = FALSE;
1769 if (inout) in = out = inin = FALSE;
1770 if (inin) in = out = FALSE;
1771 if (out) in = FALSE;
1772 if (testing) out = TRUE;
1773
1774 if (codetype&HAS_HEADER)
1775 {
1776 if (codetype&SPARSE6) writeline(outfile,SPARSE6_HEADER);
1777 else writeline(outfile,GRAPH6_HEADER);
1778 }
1779
1780 //if (codetype&SPARSE6) outproc = writes6;
1781 //else outproc = writeg6;
1782
1783 numread = noncub = nonham = nonconn = 0;
1784 numto = 0;
1785 INITRANBYTIME;
1786 t0 = CPUTIME;
1787
1788 limit = verbose ? BADLIM : 1;
1789
1790 SG_INIT(sg);
1791 while (read_sg(infile,&sg))
1792 {
1793 ++numread;
1794 n = sg.nv;
1795
1796 if (n >= MAXN-1)
1797 {
1798 fprintf(stderr,
1799 ">E cubhamg: input " COUNTER_FMT " too big\n",numread);
1800 exit(1);
1801 }
1802 #if MAXM==1
1803 if (3*n >= 2*WORDSIZE)
1804 {
1805 fprintf(stderr,
1806 ">E cubhamg: must compile with MAXM>1 for ne>WORDSIZE\n");
1807 exit(1);
1808 }
1809 #endif
1810
1811 if (!sgtocub(&sg,cub,&ne))
1812 {
1813 if (verbose)
1814 fprintf(msgfile,"Input " COUNTER_FMT
1815 " has maxdeg>3.\n",numread);
1816 ++noncub;
1817 }
1818 else if (biconn && !biconnected_cub(cub,n))
1819 ++nonconn;
1820 else if (triconn && !triconnected_cub(cub,n))
1821 ++nonconn;
1822 else if (e34)
1823 {
1824 if ((nbad = hase34(cub,n,ne,x0,x1,y0,e34plus,limit)) > 0)
1825 {
1826 if (verbose)
1827 {
1828 fprintf(msgfile,"Input " COUNTER_FMT
1829 " fails property e34%s:",
1830 numread,e34plus ? "+" : "");
1831 if (e34plus)
1832 for (i = 0; i < nbad; ++i)
1833 fprintf(msgfile," %d-%d[%02x]",
1834 x0[i],x1[i],y0[i]);
1835 else
1836 for (i = 0; i < nbad; ++i)
1837 fprintf(msgfile," %d-%d",x0[i],x1[i]);
1838 fprintf(msgfile,"\n");
1839 }
1840 ++nonham;
1841 writelast(outfile);
1842 }
1843 }
1844 else if (inout)
1845 {
1846 if ((nbad = hasinout(cub,n,ne,x0,x1,y0,y1,limit)) > 0)
1847 {
1848 if (verbose)
1849 {
1850 fprintf(msgfile,"Input " COUNTER_FMT
1851 " fails property -+:",numread);
1852 for (i = 0; i < nbad; ++i)
1853 fprintf(msgfile," %d-%d/%d-%d",
1854 x0[i],x1[i],y0[i],y1[i]);
1855 fprintf(msgfile,"\n");
1856 }
1857 ++nonham;
1858 writelast(outfile);
1859 }
1860 }
1861 else if (inin)
1862 {
1863 if ((nbad = hasinin(cub,n,ne,x0,x1,y0,y1,limit)) > 0)
1864 {
1865 if (verbose)
1866 {
1867 fprintf(msgfile,"Input " COUNTER_FMT
1868 " fails property ++:",numread);
1869 for (i = 0; i < nbad; ++i)
1870 fprintf(msgfile," %d-%d/%d-%d",
1871 x0[i],x1[i],y0[i],y1[i]);
1872 fprintf(msgfile,"\n");
1873 }
1874 ++nonham;
1875 writelast(outfile);
1876 }
1877 }
1878 else if (out)
1879 {
1880 if ((nbad = hasout(cub,n,ne,weight,x0,x1,limit)) != 0)
1881 {
1882 if (verbose)
1883 {
1884 if (nbad == -2)
1885 fprintf(msgfile,"Input " COUNTER_FMT
1886 " timed out\n",numread);
1887 else if (nbad == -1)
1888 fprintf(msgfile,"Input " COUNTER_FMT
1889 " is nonhamiltonian\n",numread);
1890 else
1891 {
1892 fprintf(msgfile,"Input " COUNTER_FMT
1893 " fails property -:",numread);
1894 for (i = 0; i < nbad; ++i)
1895 fprintf(msgfile," %d-%d",x0[i],x1[i]);
1896 fprintf(msgfile,"\n");
1897 }
1898 }
1899 if (nbad != -1)
1900 {
1901 ++nonham;
1902 writelast(outfile);
1903 }
1904 }
1905 }
1906 else if (in)
1907 {
1908 if ((nbad = hasin(cub,n,ne,x0,x1,limit)) > 0)
1909 {
1910 if (verbose)
1911 {
1912 fprintf(msgfile,"Input " COUNTER_FMT
1913 " fails property +:",numread);
1914 for (i = 0; i < nbad; ++i)
1915 fprintf(msgfile," %d-%d",x0[i],x1[i]);
1916 fprintf(msgfile,"\n");
1917 }
1918 ++nonham;
1919 writelast(outfile);
1920 }
1921 }
1922 else if (fragment)
1923 {
1924 dofragment(numread,cub,n,ne,weight);
1925 }
1926 else if ((ch = isham(cub,n,ne,weight,vv,vi,nvv,yy,yi,nyy,cyc)) == NO)
1927 {
1928 if (verbose)
1929 fprintf(msgfile,"Input " COUNTER_FMT
1930 " is not hamiltonian.\n",numread);
1931 ++nonham;
1932 writelast(outfile);
1933 }
1934 else if (ch == HABORT)
1935 {
1936 if (verbose)
1937 fprintf(msgfile,"Input " COUNTER_FMT " timed out.\n",numread);
1938 ++numto;
1939 writelast(outfile);
1940 }
1941 else if (verbose >= 2)
1942 {
1943 fprintf(msgfile,"Cycle in input " COUNTER_FMT ":\n",numread);
1944 if (n <= 100) nl = 26;
1945 else nl = 19;
1946 for (i = 0; i < n; ++i)
1947 {
1948 if (i > 0 && i % nl == 0) fprintf(msgfile,"\n ");
1949 fprintf(msgfile," %d",cyc[i]);
1950 }
1951 fprintf(msgfile,"\n");
1952 }
1953 }
1954 t1 = CPUTIME;
1955
1956 fprintf(msgfile,">C " COUNTER_FMT
1957 " graphs read from %s\n",numread,infilename);
1958 if (noncub > 0)
1959 fprintf(msgfile,">C " COUNTER_FMT " graphs with maxdeg > 3\n",noncub);
1960 if (nonconn > 0)
1961 fprintf(msgfile,">C " COUNTER_FMT " graphs not %sconnected\n",
1962 nonconn,biconn ? "bi" : "tri");
1963 if (numto > 0)
1964 fprintf(msgfile,">C " COUNTER_FMT " graphs timed out\n",numto);
1965 fprintf(msgfile,
1966 ">Z " COUNTER_FMT " %s graphs written to %s; %3.2f sec\n",
1967 nonham+numto, e34 ? (e34plus ? "non-e34+" : "non-e34") :
1968 inout ? "non-inout" :
1969 out ? "non-out" :
1970 in ? "non-in" : "nonhamiltonian",
1971 outfilename,t1-t0);
1972 if (verbose)
1973 {
1974 fprintf(msgfile,"Tries:");
1975 for (i = 0; i <= NUMMAXNODES && numtries[i] > 0; ++i)
1976 fprintf(msgfile," " COUNTER_FMT,numtries[i]);
1977 fprintf(msgfile,"\n");
1978 }
1979 #if MAXES
1980 fprintf(msgfile,"Maximum level = %d; Maximum classstack = %d\n",
1981 maxlevel,maxclassstack);
1982 #endif
1983
1984 return 0;
1985 }
1986