1 /* genbg.c : version 1.4; B D McKay, Sep 15, 2003. */
2
3 /* TODO: consider colour swaps */
4
5 #define USAGE \
6 "genbg [-c -ugsn -vq -lzF] [-Z#] [-d#|-d#:#] [-D#|-D#:#] n1 n2 \
7 [mine[:maxe]] [res/mod] [file]"
8
9 #define HELPTEXT \
10 " Find all bicoloured graphs of a specified class.\n\
11 \n\
12 n1 : the number of vertices in the first class\n\
13 n2 : the number of vertices in the second class\n\
14 mine:maxe : a range for the number of edges\n\
15 #:0 means '# or more' except in the case 0:0\n\
16 res/mod : only generate subset res out of subsets 0..mod-1\n\
17 file : the name of the output file (defailt stdout)\n\
18 -c : only write connected graphs\n\
19 -z : all the vertices in the second class must have\n\
20 different neighbourhoods\n\
21 -F : the vertices in the second class must have at least two\n\
22 neighbours of degree at least 2\n\
23 -L : there is no vertex in the first class whose removal leaves\n\
24 the vertices in the second class unreachable from each other\n\
25 -Z# : two vertices in the second class may have at most # common nbrs\n\
26 -D# : specify an upper bound for the maximum degree.\n\
27 Example: -D6. You can also give separate maxima for the\n\
28 two parts, for example: -D5:6\n\
29 -d# : specify a lower bound for the minimum degree.\n\
30 Again, you can specify it separately for the two parts: -d1:2\n\
31 -N : use nauty format for output\n\
32 -g : use graph6 format for output (default)\n\
33 -s : use sparse6 format for output\n\
34 -a : use Greechie diagram format for output\n\
35 -u : do not output any graphs, just generate and count them\n\
36 -v : display counts by number of edges to stderr\n\
37 -l : canonically label output graphs (using the 2-part colouring)\n\
38 \n\
39 -q : suppress auxiliary output\n\
40 \n\
41 See program text for much more information.\n"
42
43 /*
44 Output formats.
45
46 If -n is absent, any output graphs are written in graph6 format.
47
48 If -n is present, any output graphs are written in nauty format.
49
50 For a graph of n vertices, the output consists of n+1 long ints
51 (even if a setword is shorter than a long int). The first contains
52 n, and the others contain the adjacency matrix. Long int i of
53 the adjacency matrix (0 <= i <= n-1) is the set of neighbours of
54 vertex i, cast from setword to long int.
55
56 PRUNE feature.
57
58 By defining the C preprocessor variables PRUNE1 and/or PRUNE2 at
59 compile time, you can filter the output of the program efficiently.
60 The value of the variable is a function name with parameter list
61 (graph *g, int *deg, int n1, int n2, int maxn2)
62
63 The function will be called for each intermediate graph generated
64 by the program, including output graphs. The parameters are:
65 g = the graph in nauty format
66 deg = an array giving the degrees of the vertices
67 n1 = the number of vertices in the first colour class
68 (same as the n1 parameter on the command line)
69 n2 = the number of vertices in the second colour class
70 maxn2 = the value of n2 on the command line
71 If n2=maxn2, the graph has the output size.
72
73 If the function returns a non-zero value, neither this graph nor
74 any of its descendants will be written to the output.
75
76 PRUNE1 and PRUNE2 are functionally equivalent, but placed at different
77 points in the program. Essentially, use PRUNE1 for fast tests that
78 eliminate many cases and PRUNE2 for slow tests that eliminate few.
79 If in doubt, try it both ways and choose the fastest.
80 You can use both PRUNE1 and PRUNE2 if you wish, in which case note
81 that the PRUNE1 test has already been passed before PRUNE2 is applied.
82
83 Vertices 0..n1-1 are always present. The program works by successively
84 adding one more vertex to the second colour class. Vertex n1+n2-1 is
85 always the last one that has been added, and (except for n2=1) the
86 subgraph formed by deleting n1+n2-1 has previously been passed by the
87 pruning function.
88
89 If -c is specified, the connectivity test has NOT been performed yet
90 at the time the pruning function is called. However the simplicity
91 test indicated by -z HAS been performed if -z is specified.
92
93 OUTPROC feature.
94
95 By defining the C preprocessor variable OUTPROC at compile time
96 (for Unix the syntax is -DOUTPROC=procname on the cc command),
97 genbg can be made to call a procedure of your manufacture with each
98 output graph instead of writing anything. Your procedure needs to
99 have type void and the argument list (FILE *f, graph *g, int n1, int n2).
100 f is a stream open for writing (in fact, in the current version it
101 is always stdout), g is the graph in nauty format, and n is the number
102 of vertices. Your procedure can be in a separate file so long as it
103 is linked with genbg. The global variables nooutput, nautyformat
104 and canonise (all type boolean) can be used to test for the presence
105 of the flags -u, -n and -l, respectively.
106
107 For backward compatibility, it is possible to instead define OUTPROC1 to
108 be the name of a procedure with argument list (FILE *f, graph *g, int n).
109
110 INSTRUMENT feature.
111
112 If the C preprocessor variable INSTRUMENT is defined at compile time,
113 extra code is inserted to collect statistics during execution, and
114 more information is written to stderr at termination.
115
116 **************************************************************************
117
118 Author: B. D. McKay, Oct 1994. bdm@cs.anu.edu.au
119 Copyright B. McKay (1994-2003). All rights reserved.
120 This software is subject to the conditions and waivers
121 detailed in the file nauty.h.
122 1 May 2003 : fixed PRUNE feature
123 13 Sep 2003 : added Greechie output, all outprocs have n1,n2
124 9 Oct 2003 : changed -l to respect partition
125
126 **************************************************************************/
127
128 #undef MAXN
129 #define MAXN WORDSIZE
130
131 #ifndef MAXN1
132 #define MAXN1 24 /* not more than 30 */
133 #endif
134 #define ONE_WORD_SETS
135 #include "gtools.h" /* which includes nauty.h and stdio.h */
136
137 static void (*outproc)(FILE*,graph*,int,int);
138 #ifdef OUTPROC
139 extern void OUTPROC(FILE*,graph*,int,int);
140 #endif
141
142 static FILE *outfile; /* file for output graphs */
143 static boolean connec; /* presence of -c */
144 static boolean verbose; /* presence of -v */
145 static boolean simple; /* presence of -z */
146 boolean nautyformat; /* presence of -n */
147 boolean nooutput; /* presence of -u */
148 boolean canonise; /* presence of -l */
149 boolean graph6; /* presence of -g */
150 boolean sparse6; /* presence of -s */
151 boolean greout; /* presence of -a */
152 boolean quiet; /* presence of -q */
153 boolean footfree; /* presence of -F */
154 boolean cutfree; /* presence of -L */
155 int class1size; /* same as n1 */
156 int maxcommon; /* -1 or value of -Z */
157 static int maxdeg1,maxdeg2,n1,maxn2,mine,maxe,nprune,mod,res,curres;
158 static int mindeg1,mindeg2;
159 static graph gcan[MAXN];
160 static int xval[MAXN]; /* x-bit version of second class, xval[0..] */
161
162 #if MAXN1 <= 16
163 static int xbit[] = {0x0001,0x0002,0x0004,0x0008,
164 0x0010,0x0020,0x0040,0x0080,
165 0x0100,0x0200,0x0400,0x0800,
166 0x1000,0x2000,0x4000,0x8000};
167
168 #define XNEXTBIT(x) \
169 ((x)&0xFF ? 7-leftbit[(x)&0xFF] : 15-leftbit[((x)>>8)&0xFF])
170 #define XPOPCOUNT(x) (bytecount[((x)>>8)&0xFF] + bytecount[(x)&0xFF])
171 #elif MAXN1 <= 24
172 static int xbit[] = {0x000001,0x000002,0x000004,0x000008,
173 0x000010,0x000020,0x000040,0x000080,
174 0x000100,0x000200,0x000400,0x000800,
175 0x001000,0x002000,0x004000,0x008000,
176 0x010000,0x020000,0x040000,0x080000,
177 0x100000,0x200000,0x400000,0x800000};
178
179 #define XNEXTBIT(x) \
180 ((x)&0xFF ? 7-leftbit[(x)&0xFF] : \
181 (x)&0xFF00 ? 15-leftbit[((x)>>8)&0xFF] : 23-leftbit[((x)>>16)&0xFF])
182 #define XPOPCOUNT(x) (bytecount[((x)>>8)&0xFF] \
183 + bytecount[((x)>>16)&0xFF] + bytecount[(x)&0xFF])
184 #else
185 static int xbit[] = {0x00000001,0x00000002,0x00000004,0x00000008,
186 0x00000010,0x00000020,0x00000040,0x00000080,
187 0x00000100,0x00000200,0x00000400,0x00000800,
188 0x00001000,0x00002000,0x00004000,0x00008000,
189 0x00010000,0x00020000,0x00040000,0x00080000,
190 0x00100000,0x00200000,0x00400000,0x00800000,
191 0x01000000,0x02000000,0x04000000,0x08000000,
192 0x10000000,0x20000000,0x40000000,0x80000000};
193
194 #define XNEXTBIT(x) \
195 ((x)&0xFF ? 7-leftbit[(x)&0xFF] : \
196 (x)&0xFF00 ? 15-leftbit[((x)>>8)&0xFF] : \
197 (x)&0xFF0000 ? 23-leftbit[((x)>>16)&0xFF] : \
198 31-leftbit[((x)>>24)&0xFF])
199 #define XPOPCOUNT(x) (bytecount[((x)>>8)&0xFF] \
200 + bytecount[((x)>>16)&0xFF] + \
201 + bytecount[((x)>>24)&0xFF] + bytecount[(x)&0xFF])
202 #endif
203
204 typedef struct
205 {
206 int ne,dmax; /* values used for xlb,xub calculation */
207 int xlb,xub; /* saved bounds on extension degree */
208 int lo,hi; /* work purposes for orbit calculation */
209 int *xorb; /* min orbit representative */
210 } leveldata;
211
212 typedef struct
213 {
214 long hi,lo;
215 } bigint;
216
217 #define ZEROBIG(big) big.hi = big.lo = 0L
218 #define ISZEROBIG(big) (big.lo == 0 && big.hi == 0)
219 #define SETBIG(big,value) {big.hi = 0L; big.lo = (value);}
220 #define ADDBIG(big,extra) if ((big.lo += (extra)) >= 1000000000L) \
221 { ++big.hi; big.lo -= 1000000000L;}
222 #define PRINTBIG(file,big) if (big.hi == 0) \
223 fprintf(file,"%ld",big.lo); else fprintf(file,"%ld%09ld",big.hi,big.lo)
224 #define BIGTODOUBLE(big) (big.hi * 1000000000.0 + big.lo)
225 #define SUMBIGS(big1,big2) {if ((big1.lo += big2.lo) >= 1000000000L) \
226 {big1.lo -= 1000000000L; big1.hi += big2.hi + 1L;} \
227 else big1.hi += big2.hi;}
228 #define SUBBIGS(big1,big2) {if ((big1.lo -= big2.lo) < 0L) \
229 {big1.lo += 1000000000L; big1.hi -= big2.hi + 1L;} \
230 else big1.hi -= big2.hi;}
231 /* Note: SUBBIGS must not allow the value to go negative.
232 SUMBIGS and SUBBIGS both permit big1 and big2 to be the same bigint. */
233 #define ISEQBIG(big1,big2) (big1.lo == big2.lo && big1.hi == big2.hi)
234 #define ISASBIG(big,value) (big.hi > 0 || big.lo >= (value))
235
236 #define IFLE1BITS(ww) if (!((ww)&((ww)-1)))
237
238 static leveldata data[MAXN]; /* data[n] is data for n -> n+1 */
239 static bigint ecount[1+MAXN*MAXN/4]; /* counts by number of edges */
240 static int xstart[MAXN+1]; /* index into xset[] for each cardinality */
241 static int *xset; /* array of all x-sets in card order */
242 static int *xcard; /* cardinalities of all x-sets */
243 static int *xinv; /* map from x-set to index in xset */
244
245 #ifdef INSTRUMENT
246 static long nodes[MAXN],rigidnodes[MAXN],fertilenodes[MAXN];
247 static long a1calls,a1nauty,a1succs;
248 static long a2calls,a2nauty,a2uniq,a2succs;
249 #endif
250
251 /************************************************************************/
252
253 void
writeny(FILE * f,graph * g,int n1,int n2)254 writeny(FILE *f, graph *g, int n1, int n2)
255 /* write graph g (n1+n2 vertices) to file f in y format */
256 {
257 static char ybit[] = {32,16,8,4,2,1};
258 char s[(MAXN*(MAXN-1)/2 + 5)/6 + 4];
259 register int i,j,k;
260 register char y,*sp;
261 int n;
262
263 n = n1 + n2;
264 sp = s;
265 *(sp++) = 0x40 | n;
266 y = 0x40;
267
268 k = -1;
269 for (j = 1; j < n; ++j)
270 for (i = 0; i < j; ++i)
271 {
272 if (++k == 6)
273 {
274 *(sp++) = y;
275 y = 0x40;
276 k = 0;
277 }
278 if (g[i] & bit[j]) y |= ybit[k];
279 }
280 if (n >= 2) *(sp++) = y;
281 *(sp++) = '\n';
282 *sp = '\0';
283
284 if (fputs(s,f) == EOF || ferror(f))
285 {
286 fprintf(stderr,">E writeny : error on writing file\n");
287 exit(2);
288 }
289 }
290
291 /************************************************************************/
292
293 void
writeg6x(FILE * f,graph * g,int n1,int n2)294 writeg6x(FILE *f, graph *g, int n1, int n2)
295 /* write graph g (n1+n2 vertices) to file f in graph6 format */
296 {
297 writeg6(f,g,1,n1+n2);
298 }
299
300 /************************************************************************/
301
302 void
writes6x(FILE * f,graph * g,int n1,int n2)303 writes6x(FILE *f, graph *g, int n1, int n2)
304 /* write graph g (n1+n2 vertices) to file f in graph6 format */
305 {
306 writes6(f,g,1,n1+n2);
307 }
308
309 /************************************************************************/
310
311 void
writegre(FILE * f,graph * g,int n1,int n2)312 writegre(FILE *f, graph *g, int n1, int n2)
313 /* write graph g (n1+n2 vertices) to file f in Greechie diagram format */
314 {
315 static char atomname[] = "123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ\
316 abcdefghijklmnopqrstuvwxyz!\"#$%&'()*-/:;<=>?@[\\]^_`{|}~";
317 char grestr[MAXN*MAXN+MAXN+5];
318 int i,j,k;
319 setword gi;
320
321 k = 0;
322 for (i = n1; i < n1+n2; ++i)
323 {
324 if (i > n1) grestr[k++] = ',';
325 gi = g[i];
326 while (gi)
327 {
328 TAKEBIT(j,gi);
329 grestr[k++] = atomname[j];
330 }
331 }
332 grestr[k++] = '.';
333 grestr[k++] = '\n';
334 grestr[k] = '\0';
335 if (fputs(grestr,f) == EOF || ferror(f))
336 {
337 fprintf(stderr,">E genbg : error on writing file\n");
338 gt_abort("genbg");
339 }
340 }
341
342 /***********************************************************************/
343
344 static void
nullwrite(FILE * f,graph * g,int n1,int n2)345 nullwrite(FILE *f, graph *g, int n1, int n2)
346 /* don't write graph g (n1+n2 vertices) to file f */
347 {
348 }
349
350 /***********************************************************************/
351
352 #ifdef OUTPROC1
353
354 static void
write12(FILE * f,graph * g,int n1,int n2)355 write12(FILE *f, graph *g, int n1, int n2)
356 /* pass to OUTPROC1 */
357 {
358 OUTPROC1(f,g,1,n1+n2);
359 }
360
361 #endif
362
363 /***********************************************************************/
364
365 void
writenauty(FILE * f,graph * g,int n1,int n2)366 writenauty(FILE *f, graph *g, int n1, int n2)
367 /* write graph g (n1+n2 vertices) to file f in nauty format */
368 {
369 long buffer[MAXN+1];
370 register int i;
371 int n;
372
373 n = n1 + n2;
374 buffer[0] = n;
375 for (i = 0; i < n; ++i) buffer[i+1] = g[i];
376
377 if (fwrite((char*)buffer,sizeof(long),n+1,f) != n+1)
378 {
379 fprintf(stderr,">E writenauty : error on writing file\n");
380 exit(2);
381 }
382 }
383
384 /*********************************************************************/
385
386 static void
fragments(int * x,int nx,int * frag,int * nfrag)387 fragments(int *x, int nx, int *frag, int *nfrag)
388 /* For each v in union(x[0..nx-1]), find the components of the
389 hypergraph x - v and add them to frag if there are more than one. */
390 /* This implementation is shocking. Improve it! */
391 {
392 int allx,i,j,v;
393 int vbit,nw,w[MAXN];
394 boolean done;
395
396 allx = 0;
397 for (i = 0; i < nx; ++i) allx |= x[i];
398
399 *nfrag = 0;
400
401 while (allx)
402 {
403 v = XNEXTBIT(allx);
404 vbit = xbit[v];
405 allx &= ~vbit;
406
407 for (i = 0; i < nx; ++i) w[i] = x[i] & ~vbit;
408 nw = nx;
409
410 done = FALSE;
411 while (!done && nw > 1)
412 {
413 done = TRUE;
414 for (i = 0; i < nw-1; ++i)
415 for (j = i+1; j < nw; )
416 if ((w[i] & w[j]) != 0)
417 {
418 w[i] |= w[j];
419 w[j] = w[nw-1];
420 --nw;
421 done = FALSE;
422 }
423 else
424 ++j;
425 }
426
427 if (nw > 1)
428 for (i = 0; i < nw; ++i)
429 frag[(*nfrag)++] = w[i];
430 }
431 }
432
433 /*********************************************************************/
434
435 static boolean
isconnected(graph * g,int n)436 isconnected(graph *g, int n)
437 /* test if g is connected */
438 {
439 register setword seen,expanded,toexpand;
440 register int i;
441
442 seen = bit[0];
443 expanded = 0;
444
445 while ((toexpand = (seen & ~expanded))) /* not == */
446 {
447 i = FIRSTBIT(toexpand);
448 expanded |= bit[i];
449 seen |= g[i];
450 }
451
452 return POPCOUNT(seen) == n;
453 }
454
455 /**************************************************************************/
456
457 static boolean
distinvar(graph * g,int * invar,int n1,int n2)458 distinvar(graph *g, int *invar, int n1, int n2)
459 /* make distance invariant/
460 exit immediately FALSE if n-1 not maximal else exit TRUE
461 Note: only invar[n1..n1+n2-1] set */
462 {
463 register int w,n;
464 register setword workset,frontier;
465 setword sofar;
466 int inv,d,v;
467
468 n = n1 + n2;
469 for (v = n-1; v >= n1; --v)
470 {
471 inv = 0;
472 sofar = frontier = bit[v];
473 for (d = 1; frontier != 0; ++d)
474 {
475 workset = 0;
476 inv += POPCOUNT(frontier) ^ (0x57 + d);
477 while (frontier)
478 {
479 w = FIRSTBIT(frontier);
480 frontier &= ~bit[w];
481 workset |= g[w];
482 }
483 frontier = workset & ~sofar;
484 sofar |= frontier;
485 }
486 invar[v] = inv;
487 if (v < n-1 && inv > invar[n-1]) return FALSE;
488 }
489 return TRUE;
490 }
491
492 /**************************************************************************/
493
494 static void
makeleveldata(void)495 makeleveldata(void)
496 /* make the level data for each level */
497 {
498 register int i,j,h;
499 int nn,nxsets,tttn;
500 long ncj;
501 leveldata *d;
502 int xw,cw;
503
504 nn = maxdeg2 <= n1 ? maxdeg2 : n1;
505 ncj = nxsets = 1;
506 for (j = 1; j <= nn; ++j)
507 {
508 ncj = (ncj * (n1 - j + 1)) / j;
509 nxsets += ncj;
510 }
511
512 tttn = 1 << n1;
513 xset = (int*) ALLOCS(nxsets,sizeof(int));
514 xcard = (int*) ALLOCS(nxsets,sizeof(int));
515 xinv = (int*) ALLOCS(tttn,sizeof(int));
516 if (xset==NULL || xcard==NULL || xinv==NULL)
517 {
518 fprintf(stderr,">E genbg: malloc failed in makeleveldata()\n");
519 exit(2);
520 }
521
522 j = 0;
523 for (i = 0; i < tttn; ++i)
524 if ((h = XPOPCOUNT(i)) <= maxdeg2)
525 {
526 xset[j] = i;
527 xcard[j] = h;
528 ++j;
529 }
530
531 if (j != nxsets)
532 {
533 fprintf(stderr,">E genbg: j=%d mxsets=%d\n",j,nxsets);
534 exit(2);
535 }
536
537 h = 1;
538 do
539 h = 3 * h + 1;
540 while (h < nxsets);
541
542 do
543 {
544 for (i = h; i < nxsets; ++i)
545 {
546 xw = xset[i];
547 cw = xcard[i];
548 for (j = i; xcard[j-h] > cw ||
549 xcard[j-h] == cw && xset[j-h] > xw; )
550 {
551 xset[j] = xset[j-h];
552 xcard[j] = xcard[j-h];
553 if ((j -= h) < h) break;
554 }
555 xset[j] = xw;
556 xcard[j] = cw;
557 }
558 h /= 3;
559 }
560 while (h > 0);
561
562 for (i = 0; i < nxsets; ++i) xinv[xset[i]] = i;
563
564 xstart[0] = 0;
565 for (i = 1; i < nxsets; ++i)
566 if (xcard[i] > xcard[i-1]) xstart[xcard[i]] = i;
567 xstart[xcard[nxsets-1]+1] = nxsets;
568
569 for (i = 0; i < maxn2; ++i)
570 {
571
572 d = &data[i];
573
574 d->xorb = (int*) ALLOCS(nxsets,sizeof(int));
575
576 if (d->xorb==NULL)
577 {
578 fprintf(stderr,">E genbg: malloc failed in makeleveldata()\n");
579 exit(2);
580 }
581
582 d->ne = d->dmax = d->xlb = d->xub = -1;
583 }
584 }
585
586 /**************************************************************************/
587
588 static UPROC
userautomproc(int count,permutation * p,int * orbits,int numorbits,int stabvertex,int n)589 userautomproc(int count, permutation *p, int *orbits, int numorbits,
590 int stabvertex, int n)
591 /* Automorphism procedure called by nauty
592 Form orbits on powerset of VG
593 Operates on data[n-n1] */
594 {
595 register int i,j1,j2;
596 register int moved,pxi,pi;
597 int w,lo,hi;
598 register int *xorb;
599
600 xorb = data[n-n1].xorb;
601 lo = data[n-n1].lo;
602 hi = data[n-n1].hi;
603
604 if (count == 1) /* first automorphism */
605 for (i = lo; i < hi; ++i) xorb[i] = i;
606
607 moved = 0;
608 for (i = 0; i < n; ++i)
609 if (p[i] != i) moved |= xbit[i];
610
611 for (i = lo; i < hi; ++i)
612 {
613 if ((w = xset[i] & moved) == 0) continue;
614 pxi = xset[i] & ~moved;
615 while (w)
616 {
617 j1 = XNEXTBIT(w);
618 w &= ~xbit[j1];
619 pxi |= xbit[p[j1]];
620 }
621 pi = xinv[pxi];
622
623 j1 = xorb[i];
624 while (xorb[j1] != j1) j1 = xorb[j1];
625 j2 = xorb[pi];
626 while (xorb[j2] != j2) j2 = xorb[j2];
627
628 if (j1 < j2) xorb[j2] = xorb[i] = xorb[pi] = j1;
629 else if (j1 > j2) xorb[j1] = xorb[i] = xorb[pi] = j2;
630 }
631 }
632
633 /*****************************************************************************
634 * *
635 * refinex(g,lab,ptn,level,numcells,count,active,goodret,code,m,n) is a *
636 * custom version of refine() which can exit quickly if required. *
637 * *
638 * Only use at level==0. *
639 * goodret : whether to do an early return for code 1 *
640 * code := -1 for n-1 not max, 0 for maybe, 1 for definite *
641 * *
642 *****************************************************************************/
643
644 static void
refinex(graph * g,int * lab,int * ptn,int level,int * numcells,permutation * count,set * active,boolean goodret,int * code,int m,int n)645 refinex(graph *g, int *lab, int *ptn, int level, int *numcells,
646 permutation *count, set *active, boolean goodret,
647 int *code, int m, int n)
648 {
649 register int i,c1,c2,labc1;
650 register setword x;
651 int split1,split2,cell1,cell2;
652 int cnt,bmin,bmax;
653 set *gptr;
654 setword workset;
655 int workperm[MAXN];
656 int bucket[MAXN+2];
657
658 if (n == 1)
659 {
660 *code = 1;
661 return;
662 }
663
664 *code = 0;
665 split1 = -1;
666 while (*numcells < n && ((split1 = nextelement(active,1,split1)) >= 0
667 || (split1 = nextelement(active,1,-1)) >= 0))
668 {
669 DELELEMENT1(active,split1);
670 for (split2 = split1; ptn[split2] > 0; ++split2)
671 {}
672 if (split1 == split2) /* trivial splitting cell */
673 {
674 gptr = GRAPHROW(g,lab[split1],1);
675 for (cell1 = 0; cell1 < n; cell1 = cell2 + 1)
676 {
677 for (cell2 = cell1; ptn[cell2] > 0; ++cell2) {}
678 if (cell1 == cell2) continue;
679 c1 = cell1;
680 c2 = cell2;
681 while (c1 <= c2)
682 {
683 labc1 = lab[c1];
684 if (ISELEMENT1(gptr,labc1))
685 ++c1;
686 else
687 {
688 lab[c1] = lab[c2];
689 lab[c2] = labc1;
690 --c2;
691 }
692 }
693 if (c2 >= cell1 && c1 <= cell2)
694 {
695 ptn[c2] = 0;
696 ++*numcells;
697 ADDELEMENT1(active,c1);
698 }
699 }
700 }
701
702 else /* nontrivial splitting cell */
703 {
704 workset = 0;
705 for (i = split1; i <= split2; ++i)
706 workset |= bit[lab[i]];
707
708 for (cell1 = 0; cell1 < n; cell1 = cell2 + 1)
709 {
710 for (cell2 = cell1; ptn[cell2] > 0; ++cell2) {}
711 if (cell1 == cell2) continue;
712 i = cell1;
713 if (x = workset & g[lab[i]]) /* not == */
714 cnt = POPCOUNT(x);
715 else
716 cnt = 0;
717 count[i] = bmin = bmax = cnt;
718 bucket[cnt] = 1;
719 while (++i <= cell2)
720 {
721 if (x = workset & g[lab[i]]) /* not == */
722 cnt = POPCOUNT(x);
723 else
724 cnt = 0;
725 while (bmin > cnt) bucket[--bmin] = 0;
726 while (bmax < cnt) bucket[++bmax] = 0;
727 ++bucket[cnt];
728 count[i] = cnt;
729 }
730 if (bmin == bmax) continue;
731 c1 = cell1;
732 for (i = bmin; i <= bmax; ++i)
733 if (bucket[i])
734 {
735 c2 = c1 + bucket[i];
736 bucket[i] = c1;
737 if (c1 != cell1)
738 {
739 ADDELEMENT1(active,c1);
740 ++*numcells;
741 }
742 if (c2 <= cell2) ptn[c2-1] = 0;
743 c1 = c2;
744 }
745 for (i = cell1; i <= cell2; ++i)
746 workperm[bucket[count[i]]++] = lab[i];
747 for (i = cell1; i <= cell2; ++i)
748 lab[i] = workperm[i];
749 }
750 }
751
752 if (ptn[n-2] == 0)
753 {
754 if (lab[n-1] == n-1)
755 {
756 *code = 1;
757 if (goodret) return;
758 }
759 else
760 {
761 *code = -1;
762 return;
763 }
764 }
765 else
766 {
767 i = n - 1;
768 while (1)
769 {
770 if (lab[i] == n-1) break;
771 --i;
772 if (ptn[i] == 0)
773 {
774 *code = -1;
775 return;
776 }
777 }
778 }
779 }
780 }
781
782 /**************************************************************************/
783
784 static void
makecanon(graph * g,graph * gcan,int n1,int n2)785 makecanon(graph *g, graph *gcan, int n1, int n2)
786 /* gcan := canonise(g) */
787 {
788 int lab[MAXN],ptn[MAXN],orbits[MAXN];
789 setword active[1];
790 int i;
791 statsblk stats;
792 static DEFAULTOPTIONS_GRAPH(options);
793 setword workspace[50];
794
795 options.writemarkers = FALSE;
796 options.writeautoms = FALSE;
797 options.getcanon = TRUE;
798 options.defaultptn = FALSE;
799
800 for (i = 0; i < n1+n2; ++i)
801 {
802 lab[i] = i;
803 ptn[i] = 1;
804 }
805 ptn[n1-1] = ptn[n1+n2-1] = 0;
806 EMPTYSET(active,1);
807 ADDELEMENT(active,0);
808 ADDELEMENT(active,n1);
809
810 nauty(g,lab,ptn,active,orbits,&options,&stats,
811 workspace,50,1,n1+n2,gcan);
812 }
813
814 /**************************************************************************/
815
816 static boolean
accept1(graph * g,int n2,int x,graph * gx,int * deg,boolean * rigid)817 accept1(graph *g, int n2, int x, graph *gx, int *deg, boolean *rigid)
818 /* decide if n2 in theta(g+x) -- version for n2+1 < maxn2 */
819 {
820 register int i,n;
821 int lab[MAXN],ptn[MAXN],orbits[MAXN];
822 permutation count[MAXN];
823 graph h[MAXN];
824 int xw;
825 int nx,numcells,code;
826 int i0,i1,degn;
827 set active[MAXM];
828 statsblk stats;
829 static DEFAULTOPTIONS_GRAPH(options);
830 setword workspace[50];
831
832 #ifdef INSTRUMENT
833 ++a1calls;
834 #endif
835
836 n = n1 + n2;
837 nx = n + 1;
838 for (i = 0; i < n; ++i) gx[i] = g[i];
839 gx[n] = 0;
840 deg[n] = degn = XPOPCOUNT(x);
841
842 xw = x;
843 while (xw)
844 {
845 i = XNEXTBIT(xw);
846 xw &= ~xbit[i];
847 gx[i] |= bit[n];
848 gx[n] |= bit[i];
849 ++deg[i];
850 }
851 #ifdef PRUNE1
852 if (PRUNE1(gx,deg,n1,n2+1,maxn2)) return FALSE;
853 #endif
854
855 for (i = 0; i < n1; ++i)
856 {
857 lab[i] = i;
858 ptn[i] = 1;
859 }
860 ptn[n1-1] = 0;
861
862 i0 = n1;
863 i1 = n;
864 for (i = n1; i < nx; ++i)
865 {
866 if (deg[i] == degn) lab[i1--] = i;
867 else lab[i0++] = i;
868 ptn[i] = 1;
869 }
870
871 ptn[n] = 0;
872
873 if (i0 == n1)
874 {
875 numcells = 2;
876 active[0] = bit[0] | bit[n1];
877 }
878 else
879 {
880 numcells = 3;
881 active[0] = bit[0] | bit[n1] | bit[i1+1];
882 ptn[i1] = 0;
883 }
884 refinex(gx,lab,ptn,0,&numcells,count,active,FALSE,&code,1,nx);
885
886 if (code < 0) return FALSE;
887
888 if (numcells == nx)
889 {
890 *rigid = TRUE;
891 #ifdef INSTRUMENT
892 ++a1succs;
893 #endif
894 #ifdef PRUNE2
895 if (PRUNE2(gx,deg,n1,n2+1,maxn2)) return FALSE;
896 #endif
897 return TRUE;
898 }
899
900 options.writemarkers = FALSE;
901 options.writeautoms = FALSE;
902 options.getcanon = TRUE;
903 options.defaultptn = FALSE;
904 options.userautomproc = userautomproc;
905
906 active[0] = 0;
907 #ifdef INSTRUMENT
908 ++a1nauty;
909 #endif
910 nauty(gx,lab,ptn,active,orbits,&options,&stats,workspace,50,1,nx,h);
911
912 if (orbits[lab[n]] == orbits[n])
913 {
914 *rigid = stats.numorbits == nx;
915 #ifdef INSTRUMENT
916 ++a1succs;
917 #endif
918 #ifdef PRUNE2
919 if (PRUNE2(gx,deg,n1,n2+1,maxn2)) return FALSE;
920 #endif
921 return TRUE;
922 }
923 else
924 return FALSE;
925 }
926
927 /**************************************************************************/
928
929 static boolean
accept2(graph * g,int n2,int x,graph * gx,int * deg,boolean nuniq)930 accept2(graph *g, int n2, int x, graph *gx, int *deg, boolean nuniq)
931 /* decide if n in theta(g+x) -- version for n+1 == maxn */
932 {
933 register int i,n;
934 int lab[MAXN],ptn[MAXN],orbits[MAXN];
935 int degx[MAXN],invar[MAXN];
936 setword vmax,gv;
937 int qn,qv;
938 permutation count[MAXN];
939 int xw;
940 int nx,numcells,code;
941 int degn,i0,i1,j,j0,j1;
942 set active[MAXM];
943 statsblk stats;
944 static DEFAULTOPTIONS_GRAPH(options);
945 setword workspace[50];
946
947 #ifdef INSTRUMENT
948 ++a2calls;
949 if (nuniq) ++a2uniq;
950 #endif
951 n = n1 + n2;
952 nx = n + 1;
953 for (i = 0; i < n; ++i)
954 {
955 gx[i] = g[i];
956 degx[i] = deg[i];
957 }
958 gx[n] = 0;
959 degx[n] = degn = XPOPCOUNT(x);
960
961 xw = x;
962 while (xw)
963 {
964 i = XNEXTBIT(xw);
965 xw &= ~xbit[i];
966 gx[i] |= bit[n];
967 gx[n] |= bit[i];
968 ++degx[i];
969 }
970 #ifdef PRUNE1
971 if (PRUNE1(gx,deg,n1,n2+1,maxn2)) return FALSE;
972 #endif
973
974 if (nuniq)
975 {
976 #ifdef INSTRUMENT
977 ++a2succs;
978 #endif
979 #ifdef PRUNE2
980 if (PRUNE2(gx,deg,n1,n2+1,maxn2)) return FALSE;
981 #endif
982 if (canonise) makecanon(gx,gcan,n1,n2+1);
983 return TRUE;
984 }
985
986 for (i = 0; i < n1; ++i)
987 {
988 lab[i] = i;
989 ptn[i] = 1;
990 }
991 ptn[n1-1] = 0;
992
993 i0 = n1;
994 i1 = n;
995 for (i = n1; i < nx; ++i)
996 {
997 if (degx[i] == degn) lab[i1--] = i;
998 else lab[i0++] = i;
999 ptn[i] = 1;
1000 }
1001
1002 ptn[n] = 0;
1003
1004 if (i0 == n1)
1005 {
1006 numcells = 2;
1007 active[0] = bit[0] | bit[n1];
1008
1009 if (!distinvar(gx,invar,n1,n2+1)) return FALSE;
1010 qn = invar[n];
1011 j0 = n1;
1012 j1 = n;
1013 while (j0 <= j1)
1014 {
1015 j = lab[j0];
1016 qv = invar[j];
1017 if (qv < qn)
1018 ++j0;
1019 else
1020 {
1021 lab[j0] = lab[j1];
1022 lab[j1] = j;
1023 --j1;
1024 }
1025 }
1026 if (j0 > n1)
1027 {
1028 if (j0 == n)
1029 {
1030 #ifdef INSTRUMENT
1031 ++a2succs;
1032 #endif
1033 #ifdef PRUNE2
1034 if (PRUNE2(gx,deg,n1,n2+1,maxn2)) return FALSE;
1035 #endif
1036 if (canonise) makecanon(gx,gcan,n1,n2+1);
1037 return TRUE;
1038 }
1039 ptn[j1] = 0;
1040 ++numcells;
1041 active[0] |= bit[j0];
1042 }
1043 }
1044 else
1045 {
1046 numcells = 3;
1047 ptn[i1] = 0;
1048 active[0] = bit[0] | bit[n1] | bit[i1+1];
1049
1050 vmax = 0;
1051 j = MAXN;
1052 for (i = 0; i < n1; ++i)
1053 if (degx[i] < j && degx[i] > 0)
1054 {
1055 j = degx[i];
1056 vmax = bit[i];
1057 }
1058 else if (degx[i] == j)
1059 vmax |= bit[i];
1060
1061 gv = gx[n] & vmax;
1062 qn = POPCOUNT(gv);
1063
1064 j0 = i1+1;
1065 j1 = n;
1066 while (j0 <= j1)
1067 {
1068 j = lab[j0];
1069 gv = gx[j] & vmax;
1070 qv = POPCOUNT(gv);
1071 if (qv > qn)
1072 return FALSE;
1073 else if (qv < qn)
1074 ++j0;
1075 else
1076 {
1077 lab[j0] = lab[j1];
1078 lab[j1] = j;
1079 --j1;
1080 }
1081 }
1082 if (j0 > i1+1)
1083 {
1084 if (j0 == n)
1085 {
1086 #ifdef INSTRUMENT
1087 ++a2succs;
1088 #endif
1089 if (canonise) makecanon(gx,gcan,n1,n2+1);
1090 return TRUE;
1091 }
1092 ptn[j1] = 0;
1093 ++numcells;
1094 active[0] |= bit[j0];
1095 }
1096 }
1097
1098 refinex(gx,lab,ptn,0,&numcells,count,active,TRUE,&code,1,nx);
1099
1100 if (code < 0)
1101 return FALSE;
1102 else if (code > 0 || numcells >= nx-4)
1103 {
1104 #ifdef INSTRUMENT
1105 ++a2succs;
1106 #endif
1107 #ifdef PRUNE2
1108 if (PRUNE2(gx,deg,n1,n2+1,maxn2)) return FALSE;
1109 #endif
1110 if (canonise) makecanon(gx,gcan,n1,n2+1);
1111 return TRUE;
1112 }
1113
1114 options.writemarkers = FALSE;
1115 options.writeautoms = FALSE;
1116 options.getcanon = TRUE;
1117 options.defaultptn = FALSE;
1118
1119 active[0] = 0;
1120 #ifdef INSTRUMENT
1121 ++a2nauty;
1122 #endif
1123 nauty(gx,lab,ptn,active,orbits,&options,&stats,workspace,50,1,nx,gcan);
1124
1125 if (orbits[lab[n]] == orbits[n])
1126 {
1127 #ifdef INSTRUMENT
1128 ++a2succs;
1129 #endif
1130 #ifdef PRUNE2
1131 if (PRUNE2(gx,deg,n1,n2+1,maxn2)) return FALSE;
1132 #endif
1133 if (canonise) makecanon(gx,gcan,n1,n2+1);
1134 return TRUE;
1135 }
1136 else
1137 return FALSE;
1138 }
1139
1140 /**************************************************************************/
1141
1142 static void
xbnds(int n2,int ne,int dmax)1143 xbnds(int n2, int ne, int dmax)
1144 /* find bounds on degree for vertex n2
1145 Store answer in data[*].* */
1146 {
1147 register int xlb,xub,m;
1148
1149 xlb = n2 == 0 ? (connec ? 1 : 0) : dmax;
1150 if (xlb < mindeg2) xlb = mindeg2;
1151 m = mine - ne - (maxn2 - n2 -1)*maxdeg2;
1152 if (m > xlb) xlb = m;
1153
1154 xub = maxdeg2;
1155 m = (maxe - ne) / (maxn2 - n2);
1156 if (m < xub) xub = m;
1157
1158 data[n2].ne = ne;
1159 data[n2].dmax = dmax;
1160 data[n2].xlb = xlb;
1161 data[n2].xub = xub;
1162 }
1163
1164 /**************************************************************************/
1165
1166 static void
genextend(graph * g,int n2,int * deg,int ne,boolean rigid,int xlb,int xub)1167 genextend(graph *g, int n2, int *deg, int ne, boolean rigid, int xlb, int xub)
1168 /* extend from n2 to n2+1 */
1169 {
1170 register int x,y,d;
1171 int *xorb,xc;
1172 int nx,i,j,imin,imax,dmax;
1173 int xlbx,xubx,n;
1174 graph gx[MAXN];
1175 int degx[MAXN];
1176 boolean rigidx;
1177 int dneed,need,nfeet,hideg,deg1,ft[MAXN],nfrag,frag[MAXN];
1178
1179 #ifdef INSTRUMENT
1180 boolean haschild;
1181
1182 haschild = FALSE;
1183 ++nodes[n2];
1184 if (rigid) ++rigidnodes[n2];
1185 #endif
1186
1187 n = n1 + n2;
1188 nx = n2 + 1;
1189 dmax = deg[n-1];
1190
1191 d = 0;
1192 dneed = mindeg1 - maxn2 + n2;
1193 need = 0;
1194 hideg = 0;
1195 deg1 = 0;
1196 for (i = 0; i < n1; ++i)
1197 {
1198 if (deg[i] == maxdeg1) d |= xbit[i];
1199 if (deg[i] <= dneed) need |= xbit[i];
1200 if (deg[i] >= 2) hideg |= xbit[i];
1201 if (deg[i] == 1) deg1 |= xbit[i];
1202 }
1203
1204 if (xlb < XPOPCOUNT(need)) xlb = XPOPCOUNT(need);
1205 if (xlb > xub) return;
1206
1207 imin = xstart[xlb];
1208 imax = xstart[xub+1];
1209 xorb = data[n2].xorb;
1210
1211 if (nx == maxn2)
1212 {
1213 if (footfree)
1214 {
1215 nfeet = 0;
1216 for (j = 0; j < n2; ++j)
1217 {
1218 x = xval[j] & hideg;
1219 IFLE1BITS(x) ft[nfeet++] = xval[j] & deg1;
1220 }
1221 }
1222 if (cutfree) fragments(xval,n2,frag,&nfrag);
1223
1224 for (i = imin; i < imax; ++i)
1225 {
1226 if (!rigid && xorb[i] != i) continue;
1227 x = xset[i];
1228 xc = xcard[i];
1229 if ((x & d) != 0) continue;
1230 if ((need & ~x) != 0) continue;
1231
1232 if (simple)
1233 {
1234 for (j = n2; --j >= 0;)
1235 if (x == xval[j]) break;
1236 if (j >= 0) continue;
1237 }
1238 if (maxcommon >= 0)
1239 {
1240 for (j = n2; --j >= 0;)
1241 {
1242 y = x & xval[j];
1243 if (XPOPCOUNT(y) > maxcommon) break;
1244 }
1245 if (j >= 0) continue;
1246 }
1247 if (footfree)
1248 {
1249 y = x & (hideg | deg1);
1250 IFLE1BITS(y) continue;
1251 for (j = 0; j < nfeet; ++j)
1252 if ((x & ft[j]) == 0) break;
1253 if (j < nfeet) continue;
1254 }
1255 if (cutfree)
1256 {
1257 y = x & (hideg | deg1);
1258 IFLE1BITS(y) continue;
1259 for (j = 0; j < nfrag; ++j)
1260 if ((x & frag[j]) == 0) break;
1261 if (j < nfrag) continue;
1262 }
1263
1264 xval[n2] = x;
1265
1266 if (nx == nprune)
1267 {
1268 if (curres == 0) curres = mod;
1269 if (--curres != 0) continue;
1270 }
1271 if (accept2(g,n2,x,gx,deg,xc > dmax))
1272 if (!connec || isconnected(gx,n+1))
1273 {
1274 ADDBIG(ecount[ne+xc],1);
1275 #ifdef INSTRUMENT
1276 haschild = TRUE;
1277 #endif
1278 (*outproc)(outfile,canonise ? gcan : gx,n1,n2+1);
1279 }
1280 }
1281 }
1282 else
1283 {
1284 for (i = imin; i < imax; ++i)
1285 {
1286 if (!rigid && xorb[i] != i) continue;
1287 x = xset[i];
1288 xc = xcard[i];
1289 if ((x & d) != 0) continue;
1290 if ((need & ~x) != 0) continue;
1291
1292 if (nx == nprune)
1293 {
1294 if (curres == 0) curres = mod;
1295 if (--curres != 0) continue;
1296 }
1297
1298 if (simple)
1299 {
1300 for (j = n2; --j >= 0;)
1301 if (x == xval[j]) break;
1302 if (j >= 0) continue;
1303 }
1304 if (maxcommon >= 0)
1305 {
1306 for (j = n2; --j >= 0;)
1307 {
1308 y = x & xval[j];
1309 if (XPOPCOUNT(y) > maxcommon) break;
1310 }
1311 if (j >= 0) continue;
1312 }
1313 xval[n2] = x;
1314
1315 for (j = 0; j < n; ++j) degx[j] = deg[j];
1316 if (data[nx].ne != ne+xc || data[nx].dmax != xc)
1317 xbnds(nx,ne+xc,xc);
1318 xlbx = data[nx].xlb;
1319 xubx = data[nx].xub;
1320 if (xlbx > xubx) continue;
1321
1322 data[nx].lo = xstart[xlbx];
1323 data[nx].hi = xstart[xubx+1];
1324 if (accept1(g,n2,x,gx,degx,&rigidx))
1325 {
1326 #ifdef INSTRUMENT
1327 haschild = TRUE;
1328 #endif
1329 genextend(gx,nx,degx,ne+xc,rigidx,xlbx,xubx);
1330 }
1331 }
1332 }
1333 #ifdef INSTRUMENT
1334 if (haschild) ++fertilenodes[n2];
1335 #endif
1336 }
1337
1338 /**************************************************************************/
1339 /**************************************************************************/
1340
main(int argc,char * argv[])1341 main(int argc, char *argv[])
1342 {
1343 char *arg;
1344 long ltemp;
1345 boolean badargs,gotD,gote,gotf,gotmr,gotZ,gotd;
1346 long Dval1,Dval2;
1347 long dval1,dval2;
1348 int i,j,imin,imax,argnum,sw;
1349 graph g[MAXN1];
1350 int deg[MAXN1];
1351 bigint nout;
1352 double t1,t2;
1353 char *outfilename;
1354
1355 HELP;
1356 nauty_check(WORDSIZE,1,MAXN,NAUTYVERSIONID);
1357
1358 if (MAXN > WORDSIZE || MAXN1 > 8*sizeof(int)-2)
1359 {
1360 fprintf(stderr,"genbg: incompatible MAXN, MAXN1 or WORDSIZE\n");
1361 fprintf(stderr,"--See notes in program source\n");
1362 exit(1);
1363 }
1364
1365 badargs = FALSE;
1366 connec = FALSE;
1367 verbose = FALSE;
1368 nautyformat = FALSE;
1369 nooutput = FALSE;
1370 canonise = FALSE;
1371 greout = FALSE;
1372 simple = FALSE;
1373 graph6 = FALSE;
1374 sparse6 = FALSE;
1375 quiet = FALSE;
1376 footfree = FALSE;
1377 cutfree = FALSE;
1378
1379 gote = FALSE;
1380 gotf = FALSE;
1381 gotmr = FALSE;
1382 gotD = FALSE;
1383 gotd = FALSE;
1384 gotZ = FALSE;
1385 outfilename = NULL;
1386
1387 maxdeg1 = maxdeg2 = MAXN;
1388 mindeg1 = mindeg2 = 0;
1389
1390 argnum = 0;
1391 for (j = 1; !badargs && j < argc; ++j)
1392 {
1393 arg = argv[j];
1394 if (arg[0] == '-' && arg[1] != '\0')
1395 {
1396 ++arg;
1397 while (*arg != '\0')
1398 {
1399 sw = *arg++;
1400 SWBOOLEAN('n',nautyformat)
1401 else SWBOOLEAN('u',nooutput)
1402 else SWBOOLEAN('q',quiet)
1403 else SWBOOLEAN('v',verbose)
1404 else SWBOOLEAN('z',simple)
1405 else SWBOOLEAN('F',footfree)
1406 else SWBOOLEAN('L',cutfree)
1407 else SWBOOLEAN('l',canonise)
1408 else SWBOOLEAN('c',connec)
1409 else SWBOOLEAN('a',greout)
1410 else SWBOOLEAN('g',graph6)
1411 else SWBOOLEAN('s',sparse6)
1412 else SWINT('Z',gotZ,maxcommon,"genbg -Z")
1413 else SWRANGE('D',":-",gotD,Dval1,Dval2,"genbg -D")
1414 else SWRANGE('d',":-",gotd,dval1,dval2,"genbg -d")
1415 #ifdef PLUGIN_SWITCHES
1416 PLUGIN_SWITCHES
1417 #endif
1418 else badargs = TRUE;
1419 }
1420 }
1421 else if (arg[0] == '-' && arg[1] == '\0')
1422 gotf = TRUE;
1423 else
1424 {
1425 if (argnum == 0)
1426 {
1427 if (sscanf(arg,"%d",&n1) != 1) badargs = TRUE;
1428 ++argnum;
1429 }
1430 else if (argnum == 1)
1431 {
1432 if (sscanf(arg,"%d",&maxn2) != 1) badargs = TRUE;
1433 ++argnum;
1434 }
1435 else if (gotf)
1436 badargs = TRUE;
1437 else
1438 {
1439 if (!gotmr)
1440 {
1441 if (sscanf(arg,"%d/%d",&res,&mod) == 2)
1442 {
1443 gotmr = TRUE;
1444 continue;
1445 }
1446 }
1447 if (!gote)
1448 {
1449 if (sscanf(arg,"%d:%d",&mine,&maxe) == 2
1450 || sscanf(arg,"%d-%d",&mine,&maxe) == 2)
1451 {
1452 gote = TRUE;
1453 if (maxe == 0 && mine > 0) maxe = MAXN*MAXN/4;
1454 continue;
1455 }
1456 else if (sscanf(arg,"%d",&mine) == 1)
1457 {
1458 gote = TRUE;
1459 maxe = mine;
1460 continue;
1461 }
1462 }
1463 if (!gotf)
1464 {
1465 outfilename = arg;
1466 gotf = TRUE;
1467 continue;
1468 }
1469 }
1470 }
1471 }
1472
1473 if (argnum < 2)
1474 badargs = TRUE;
1475 else if (n1 < 1 || maxn2 < 0 || n1 > MAXN1 || n1+maxn2 > MAXN)
1476 {
1477 fprintf(stderr,
1478 ">E genbg: must have n1=1..%d, n1+n2=1..%d\n",MAXN1,MAXN);
1479 badargs = TRUE;
1480 }
1481
1482 if (!gote)
1483 {
1484 mine = 0;
1485 maxe = n1 * maxn2;
1486 }
1487
1488 if (!gotmr)
1489 {
1490 mod = 1;
1491 res = 0;
1492 }
1493 else if (argnum == 5 || argnum > 6)
1494 badargs = TRUE;
1495
1496 if (gotd)
1497 {
1498 mindeg1 = dval1;
1499 mindeg2 = dval2;
1500 }
1501 if (gotD)
1502 {
1503 maxdeg1 = Dval1;
1504 maxdeg2 = Dval2;
1505 }
1506 if (maxdeg1 > maxn2) maxdeg1 = maxn2;
1507 if (maxdeg2 > n1) maxdeg2 = n1;
1508 if (connec && mine < n1+maxn2-1) mine = n1 + maxn2 - 1;
1509 if (connec && mindeg1 == 0) mindeg1 = 1;
1510 if (connec && mindeg2 == 0) mindeg2 = 1;
1511 if (maxe > n1*maxdeg1) maxe = n1*maxdeg1;
1512 if (maxe > maxn2*maxdeg2) maxe = maxn2*maxdeg2;
1513 if (mine < n1*mindeg1) mine = n1*mindeg1;
1514 if (mine < maxn2*mindeg2) mine = maxn2*mindeg2;
1515
1516 if (!badargs && (mine > maxe || maxe < 0 || maxdeg1 < 0 || maxdeg2 < 0))
1517 {
1518 fprintf(stderr,">E genbg: impossible mine,maxe,maxdeg values\n");
1519 badargs = TRUE;
1520 }
1521
1522 if (!gotZ) maxcommon = -1;
1523
1524 if (!badargs && (mine > maxe || maxe < 0 || maxdeg1 < 0 || maxdeg2 < 0))
1525 {
1526 fprintf(stderr,
1527 ">E genbg: impossible mine,maxe,mindeg,maxdeg values\n");
1528 badargs = TRUE;
1529 }
1530
1531 if (badargs)
1532 {
1533 fprintf(stderr,">E Usage: %s\n",USAGE);
1534 GETHELP;
1535 exit(1);
1536 }
1537
1538 if ((nautyformat!=0) + (graph6!=0) + (greout!=0)
1539 + (sparse6!=0) + (nooutput!=0) > 1)
1540 gt_abort(">E genbg: -ungsa are incompatible\n");
1541
1542 #ifdef OUTPROC
1543 outproc = OUTPROC;
1544 #else
1545 #ifdef OUTROC1
1546 outproc = write12;
1547 #endif
1548 if (nautyformat) outproc = writenauty;
1549 else if (nooutput) outproc = nullwrite;
1550 else if (sparse6) outproc = writes6x;
1551 else if (greout) outproc = writegre;
1552 else outproc = writeg6x;
1553 #endif
1554
1555 #ifdef PLUGIN_INIT
1556 PLUGIN_INIT
1557 #endif
1558
1559 #if 0
1560 if (gotd) /* TEMP */
1561 {
1562 gt_abort("genbg -d is not implemented\n");
1563 }
1564 #endif
1565
1566 for (i = 0; i <= maxe; ++i) ZEROBIG(ecount[i]);
1567
1568 if (nooutput)
1569 outfile = stdout;
1570 else if (!gotf || outfilename == NULL)
1571 {
1572 outfilename = "stdout";
1573 outfile = stdout;
1574 }
1575 else if ((outfile = fopen(outfilename,
1576 nautyformat ? "wb" : "w")) == NULL)
1577 {
1578 fprintf(stderr,
1579 ">E genbg: can't open %s for writing\n",outfilename);
1580 gt_abort(NULL);
1581 }
1582
1583 if (!quiet)
1584 {
1585 fprintf(stderr,">A %s n=%d+%d e=%d:%d d=%d:%d D=%d:%d ",
1586 argv[0],n1,maxn2,mine,maxe,
1587 mindeg1,mindeg2,maxdeg1,maxdeg2);
1588 if (simple) fprintf(stderr,"z");
1589 if (footfree) fprintf(stderr,"F");
1590 if (connec) fprintf(stderr,"c");
1591 if (maxcommon >= 0) fprintf(stderr,"Z%d",maxcommon);
1592 if (mod > 1) fprintf(stderr," class=%d/%d",res,mod);
1593 fprintf(stderr,"\n");
1594 }
1595
1596 class1size = n1;
1597
1598 for (i = 0; i < n1; ++i)
1599 {
1600 g[i] = 0;
1601 deg[i] = 0;
1602 }
1603
1604 t1 = CPUTIME;
1605
1606 if (maxn2 == 0)
1607 {
1608 if (res == 0)
1609 {
1610 ADDBIG(ecount[0],1);
1611 (*outproc)(outfile,g,n1,0);
1612 }
1613 }
1614 else
1615 {
1616 makeleveldata();
1617 curres = res;
1618 if (mod <= 1) nprune = 0;
1619 else if (maxn2 >= 6) nprune = maxn2 - 2;
1620 else if (maxn2 >= 3) nprune = maxn2 - 1;
1621 else nprune = maxn2;
1622
1623 xbnds(0,0,0);
1624 imin = xstart[data[0].xlb];
1625 imax = xstart[data[0].xub+1];
1626
1627 for (i = imin; i < imax; ++i)
1628 data[0].xorb[i] = -1;
1629
1630 for (i = data[0].xlb; i <= data[0].xub; ++i)
1631 data[0].xorb[xstart[i]] = xstart[i];
1632
1633 genextend(g,0,deg,0,FALSE,data[0].xlb,data[0].xub);
1634 }
1635 t2 = CPUTIME;
1636
1637 ZEROBIG(nout);
1638 for (i = 0; i <= maxe; ++i) SUMBIGS(nout,ecount[i]);
1639
1640 if (verbose)
1641 for (i = 0; i <= maxe; ++i)
1642 if (!ISZEROBIG(ecount[i]))
1643 {
1644 fprintf(stderr,">C ");
1645 PRINTBIG(stderr,ecount[i]);
1646 fprintf(stderr," graphs with %d edges\n",i);
1647 }
1648
1649 #ifdef INSTRUMENT
1650 fprintf(stderr,"\n>N node counts\n");
1651 for (i = 0; i < maxn2; ++i)
1652 fprintf(stderr," level %2d: %7ld (%ld rigid, %ld fertile)\n",
1653 i,nodes[i],rigidnodes[i],fertilenodes[i]);
1654 fprintf(stderr,">A1 %ld calls to accept1, %ld nauty, %ld succeeded\n",
1655 a1calls,a1nauty,a1succs);
1656 fprintf(stderr,
1657 ">A2 %ld calls to accept2, %ld nuniq, %ld nauty, %ld succeeded\n",
1658 a2calls,a2uniq,a2nauty,a2succs);
1659 fprintf(stderr,"\n");
1660 #endif
1661
1662 if (!quiet)
1663 {
1664 fprintf(stderr,">Z ");
1665 PRINTBIG(stderr,nout);
1666 fprintf(stderr," graphs generated in %3.2f sec\n",t2-t1);
1667 }
1668
1669 exit(0);
1670 }
1671