FORM 4.3
gentopo.cc
1// { ( [
2
3#ifdef HAVE_CONFIG_H
4#ifndef CONFIG_H_INCLUDED
5#define CONFIG_H_INCLUDED
6#include <config.h>
7#endif
8#endif
9
10#include <iostream>
11#include <iomanip>
12#include <sstream>
13#include <cstdio>
14#include <stdlib.h>
15
16#include "gentopo.h"
17
18#define DUMMYUSE(x) (void)(x)
19
20using namespace std;
21
22// The next limitation is imposed by the fact that the latest compilers
23// can give warnings on declarations like "int dtcl1[nNodes]"
24// With a bit of overkil there should be no real problems.
25#define MAXNODES 100
26#define MAXNCLASSES 100
27
28// Generate scalar connected Feynman graphs.
29
30//==============================================================
31
32typedef int Bool;
33const int True = 1;
34const int False = 0;
35
36// compile options
37const int CHECK = True;
38//const int MONITOR = True;
39const int MONITOR = False;
40
41const int OPTPRINT = False;
42
43const int DEBUG0 = False;
44//const int DEBUG1 = False;
45const int DEBUG = False;
46
47// for debugging memory use
48#define DEBUGM False
49
50//==============================================================
51class MGraph;
52class EGraph;
53
54#define Extern
55#define Global
56#define Local
57
58// External functions
59Extern void toForm(EGraph *egraph);
60Extern int countPhiCon(int ex, int lp, int v4);
61
62// Temporal definition: they should be replaced by FROM functions.
63Local BigInt factorial(int n);
64Local BigInt ipow(int n, int p);
65
66// Local functions
67Local int nextPerm(int nelem, int nclass, int *cl, int *r, int *q, int *p, int count);
68
69Local void erEnd(const char *msg);
70Local int *newArray(int size, int val);
71Local void deletArray(int *a);
72Local void copyArray(int size, int *a0, int *a1);
73Local int **newMat(int n0, int n1, int val);
74Local void deleteMat(int **m, int n0, int n1);
75Local void printArray(int n, int *p);
76Local void printMat(int n0, int n1, int **m);
77
78#if DEBUGM
79static int countNMC = 0;
80static int countEG = 0;
81static int countMG = 0;
82#endif
83
84//==============================================================
85// class EGraph
86
87EGraph::EGraph(int nnodes, int nedges, int mxdeg)
88{
89 int j;
90
91 nNodes = nnodes;
92 nEdges = nedges;
93 maxdeg = mxdeg; // maxmum value of degree of nodes
94 nExtern = 0;
95
96 nodes = new ENode[nNodes];
97 edges = new EEdge[nEdges+1];
98 for (j = 0; j < nNodes; j++) {
99 nodes[j].deg = 0;
100 nodes[j].edges = new int[maxdeg];
101 }
102#if DEBUGM
103 printf("+++ new EGraph %d\n", ++countEG);
104 if(countEG > 100) { exit(1); }
105#endif
106}
107
108EGraph::~EGraph()
109{
110 int j;
111
112 for (j = 0; j < nNodes; j++) {
113 delete[] nodes[j].edges;
114 }
115 delete[] nodes;
116 delete[] edges;
117#if DEBUGM
118 printf("+++ delete EGraph %d\n", countEG--);
119#endif
120}
121
122// construct EGraph from adjacency matrix
123void EGraph::init(int pid, long gid, int **adjmat, Bool sopi, BigInt nsm, BigInt esm)
124{
125 int n0, n1, ed, e, eext, eint;
126// Bool ok;
127
128 pId = pid;
129 gId = gid;
130 opi = sopi;
131 nsym = nsm;
132 esym = esm;
133
134 for (n0 = 0; n0 < nNodes; n0++) {
135 nodes[n0].deg = 0;
136 }
137 ed = 1;
138 for (n0 = 0; n0 < nNodes; n0++) {
139 for (e = 0; e < adjmat[n0][n0]/2; e++, ed++) {
140 edges[ed].nodes[0] = n0;
141 edges[ed].nodes[1] = n0;
142 nodes[n0].edges[nodes[n0].deg++] = - ed;
143 nodes[n0].edges[nodes[n0].deg++] = ed;
144 edges[ed].ext = nodes[n0].ext;
145 }
146 for (n1 = n0+1; n1 < nNodes; n1++) {
147 for (e = 0; e < adjmat[n0][n1]; e++, ed++) {
148 edges[ed].nodes[0] = n0;
149 edges[ed].nodes[1] = n1;
150 nodes[n0].edges[nodes[n0].deg++] = - ed;
151 nodes[n1].edges[nodes[n1].deg++] = ed;
152 edges[ed].ext = (nodes[n0].ext || nodes[n1].ext);
153 }
154 }
155 }
156 if (CHECK) {
157 if (ed != nEdges+1) {
158 printf("*** EGraph::init: ed=%d != nEdges=%d\n", ed, nEdges);
159 erEnd("*** EGraph::init: illegal connection\n");
160 }
161 }
162
163 // name of momenta
164 eext = 1;
165 eint = 1;
166 for (ed = 1; ed <= nEdges; ed++) {
167 if (edges[ed].ext) {
168 edges[ed].momn = eext++;
169 edges[ed].momc[0] = 'Q';
170 edges[ed].momc[1] = ((char)0);
171 } else {
172 edges[ed].momn = eint++;
173 edges[ed].momc[0] = 'P';
174 edges[ed].momc[1] = ((char)0);
175 }
176 }
177}
178
179// set external particle to node 'nd'
180void EGraph::setExtern(int nd, Bool val)
181{
182 nodes[nd].ext = val;
183}
184
185// end of calling 'setExtern'
186void EGraph::endSetExtern(void)
187{
188 int n;
189
190 nExtern = 0;
191 for (n = 0; n < nNodes; n++) {
192 if (nodes[n].ext) {
193 nExtern++;
194 }
195 }
196}
197
198// print the EGraph
199void EGraph::print()
200{
201 int nd, lg, ed, nlp;
202
203 nlp = nEdges - nNodes + 1;
204 printf("\n");
205 printf("EGraph: pId=%d gId=%ld nExtern=%d nLoops=%d nNodes=%d nEdges=%d\n",
206 pId, gId, nExtern, nlp, nNodes, nEdges);
207 printf(" sym = (%ld * %ld) maxdeg=%d\n", nsym, esym, maxdeg);
208 printf(" Nodes\n");
209 for (nd = 0; nd < nNodes; nd++) {
210 if (nodes[nd].ext) {
211 printf(" %2d Extern ", nd);
212 } else {
213 printf(" %2d Vertex ", nd);
214 }
215 printf("deg=%d [", nodes[nd].deg);
216 for (lg = 0; lg < nodes[nd].deg; lg++) {
217 printf(" %3d", nodes[nd].edges[lg]);
218 }
219 printf("]\n");
220 }
221 printf(" Edges\n");
222 for (ed = 1; ed <= nEdges; ed++) {
223 if (edges[ed].ext) {
224 printf(" %2d Extern ", ed);
225 } else {
226 printf(" %2d Intern ", ed);
227 }
228 printf("%s%d", (edges[ed].ext)?"Q":"p", edges[ed].momn);
229 printf(" [%3d %3d]\n", edges[ed].nodes[1], edges[ed].nodes[0]);
230 }
231}
232
233//==============================================================
234// class MNode : nodes in MGraph
235
236MNode::MNode(int vid, int vdeg, Bool vext, int vclss)
237{
238 id = vid; // id of the node
239 deg = vdeg; // degree of the node
240 freelg = vdeg; // number of free legs
241 clss = vclss; // class to which the node belogns
242 ext = vext; // external node or not
243}
244
245//===============================================================
246// class of node-classes for MGraph
247//
248class MNodeClass {
249 public:
250 int nNodes; // the number of nodes
251 int nClasses; // the number of classes
252 int *clist; // the number of nodes in each class
253 int *ndcl; // node --> class
254 int **clmat; // matrix used for classification
255 int *flist; // the first node in each class
256 int maxdeg; // maximal value of degree(node)
257 int forallignment;
258
259 MNodeClass(int nnodes, int ncl);
260 ~MNodeClass();
261 void init(int *cl, int mxdeg, int **adjmat);
262 void copy(MNodeClass* mnc);
263 int clCmp(int nd0, int nd1, int cn);
264 void printMat(void);
265
266 void mkFlist(void);
267 void mkNdCl(void);
268 void mkClMat(int **adjmat);
269 void incMat(int nd, int td, int val);
270 Bool chkOrd(int nd, int ndc, MNodeClass *cl, int *dtcl);
271 int cmpArray(int *a0, int *a1, int ma);
272};
273
274MNodeClass::MNodeClass(int nnodes, int ncl)
275{
276 nNodes = nnodes;
277 nClasses = ncl;
278 clist = new int[nClasses];
279 ndcl = new int[nNodes];
280 clmat = newMat(nNodes, nClasses, 0);
281 flist = new int[nClasses+1];
282
283#if DEBUGM
284 printf("+++ new MNodeClass %d\n", ++countNMC);
285 if(countNMC > 100) { exit(1); }
286#endif
287}
288
289MNodeClass::~MNodeClass()
290{
291 delete[] clist;
292 delete[] ndcl;
293 deleteMat(clmat, nNodes, nClasses);
294 delete[] flist;
295
296#if DEBUGM
297 printf("+++ delete MNodeClass %d\n", countNMC--);
298#endif
299}
300
301void MNodeClass::init(int *cl, int mxdeg, int **adjmat)
302{
303 int j;
304
305 for (j = 0; j < nClasses; j++) {
306 clist[j] = cl[j];
307 }
308 maxdeg = mxdeg;
309 mkNdCl();
310 mkClMat(adjmat);
311 mkFlist();
312}
313
314void MNodeClass::copy(MNodeClass* mnc)
315{
316 int j, k;
317
318 for (k = 0; k < nClasses; k++) {
319 clist[k] = mnc->clist[k];
320 }
321 maxdeg = mnc->maxdeg;
322 for (j = 0; j < nNodes; j++) {
323 ndcl[j] = mnc->ndcl[j];
324 for (k = 0; k < nClasses; k++) {
325 clmat[j][k] = mnc->clmat[j][k];
326 }
327 }
328 for (k = 0; k < nClasses+1; k++) {
329 flist[k] = mnc->flist[k];
330 }
331}
332
333// Construct flist
334// The set of nodes in class 'cl' is [flist[cl],...,flist[cl+1]-1]
335void MNodeClass::mkFlist(void)
336{
337 int j, f;
338
339 f = 0;
340 for (j = 0; j < nClasses; j++) {
341 flist[j] = f;
342 f += clist[j];
343 }
344 flist[nClasses] = f;
345}
346
347// Construct ndcl
348// ndcl[nd] = (the class id in which node 'nd' belongs)
349void MNodeClass::mkNdCl(void)
350{
351 int c, k;
352 int nd = 0;
353
354 for (c = 0; c < nClasses; c++) {
355 for (k = 0; k < clist[c]; k++) {
356 ndcl[nd++] = c;
357 }
358 }
359}
360
361// Comparison of two nodes 'nd0' and 'nd1'
362// Ordering is lexicographic (class, connection configuration)
363int MNodeClass::clCmp(int nd0, int nd1, int cn)
364{
365
366 // Wether two nodes are in a same class or not.
367 int cmp = ndcl[nd0] - ndcl[nd1];
368 if (cmp != 0) {
369 return cmp;
370 }
371
372 // Sign '-' signifies the reverse ordering
373 cmp = - cmpArray(clmat[nd0], clmat[nd1], cn);
374 if (cmp != 0) {
375 return cmp;
376 }
377 // for particles ???
378 return cmp;
379}
380
381// Construct a matrix 'clmat[nd][tc]' which is the number
382// of edges connecting 'nd' and all nodes in class 'tc'.
383void MNodeClass::mkClMat(int **adjmat)
384{
385 int nd, td, tc;
386
387 for (nd = 0; nd < nNodes; nd++) {
388 for (td = 0; td < nNodes; td++) {
389 tc = ndcl[td]; // another node
390 clmat[nd][tc] += adjmat[nd][td]; // the number of edges 'nd'--'td'
391 }
392
393 }
394}
395
396// Increase the number of edges by 'val' between 'nd' and 'td'.
397void MNodeClass::incMat(int nd, int td, int val)
398{
399 int tdc = ndcl[td]; // class of 'td'
400 clmat[nd][tdc] += val; // modify matrix 'clmat'.
401}
402
403// Check whether the configuration satisfies the ordering condition or not.
404Bool MNodeClass::chkOrd(int nd, int ndc, MNodeClass *cl, int *dtcl)
405{
406 Bool tcl[MAXNCLASSES];
407// Bool tcl[cl->nClasses];
408 int tn, tc, mxn, cmp, n;
409 DUMMYUSE(nd);
410
411 for (tc = 0; tc < cl->nClasses; tc++) {
412 tcl[tc] = False;
413 }
414 for (tn = 0; tn < nNodes; tn++) {
415 if (dtcl[tn] == 0) {
416 tcl[cl->ndcl[tn]] = False;
417 }
418 }
419 for (tc = 0; tc < cl->nClasses; tc++) {
420 if (tcl[tc] && ndc != tc) {
421 mxn = flist[tc+1];
422 for (n = flist[tc]+1; n < mxn; n++) {
423 cmp = - clmat[n-1][tc] + clmat[n][tc];
424 if (cmp > 0) {
425 return False;
426 }
427 }
428 }
429 }
430 return True;
431}
432
433// Print configuration matrix.
434void MNodeClass::printMat(void)
435{
436 int j1, j2;
437
438 cout << endl;
439
440 // the first line
441 cout << setw(2) << "nd" << ": " << setw(2) << "cl: ";
442 for (j2 = 0; j2 < nClasses; j2++) {
443 cout << setw(2) << j2 << " ";
444 }
445 cout << endl;
446
447 // print raw
448 for (j1 = 0; j1 < nNodes; j1++) {
449 cout << setw(2) << j1 << ": " << setw(2) << ndcl[j1] << ": [";
450 for (j2 = 0; j2 < nClasses; j2++) {
451 cout << " " << setw(2) << clmat[j1][j2];
452 }
453 cout << "] " << endl;
454 }
455}
456
457int MNodeClass::cmpArray(int *a0, int *a1, int ma)
458{
459 for (int j = 0; j < ma; j++) {
460 if (a0[j] < a1[j]) {
461 return -1;
462 } else if (a0[j] > a1[j]) {
463 return 1;
464 }
465 }
466 return 0;
467}
468
469//===============================================================
470// class MGraph : scalar graph expressed by matrix form
471
472MGraph::MGraph(int pid, int ncl, int *cldeg, int *clnum, int *clext, Bool sopi)
473{
474 int nn, ne, j, k;
475
476 // initial conditions
477 nClasses = ncl;
478 clist = new int[ncl];
479
480 mindeg = -1;
481 maxdeg = -1;
482 ne = 0;
483 nn = 0;
484 for (j = 0; j < nClasses; j++) {
485 clist[j] = clnum[j];
486 nn += clnum[j];
487 ne += cldeg[j]*clnum[j];
488 if (mindeg < 0) {
489 mindeg = cldeg[j];
490 } else {
491 mindeg = min(mindeg, cldeg[j]);
492 }
493 if (maxdeg < 0) {
494 maxdeg = cldeg[j];
495 } else {
496 maxdeg = max(maxdeg, cldeg[j]);
497 }
498 }
499 if (ne % 2 != 0) {
500 printf("Sum of degrees are not even\n");
501 for (j = 0; j < nClasses; j++) {
502 printf("class %2d: %2d %2d %2d\n",
503 j, cldeg[j], clnum[j], clext[j]);
504 }
505 erEnd("illegal degrees of nodes");
506 }
507 pId = pid;
508 nNodes = nn;
509 nodes = new MNode*[nNodes];
510
511 selOPI = sopi;
512
513 nEdges = ne / 2;
514 nLoops = nEdges - nNodes + 1;
515
516 egraph = new EGraph(nNodes, nEdges, maxdeg);
517 nn = 0;
518 for (j = 0; j < nClasses; j++) {
519 for (k = 0; k < clist[j]; k++, nn++) {
520 nodes[nn] = new MNode(nn, cldeg[j], clext[j], j);
521 egraph->setExtern(nn, clext[j]);
522 }
523 }
524 egraph->endSetExtern();
525
526 // generated set of graphs
527 ndiag = 0;
528 n1PI = 0;
529
530 // the current graph
531 adjMat = newMat(nNodes, nNodes, 0);
532 nsym = ToBigInt(0);
533 esym = ToBigInt(0);
534 c1PI = 0;
535 wsum = ToFraction(0, 1);
536 wsopi = ToFraction(0, 1);
537
538 // current node classification
539 curcl = new MNodeClass(nNodes, nClasses);
540
541 // measures of efficiency
542 ngen = 0;
543 ngconn = 0;
544
545 if (MONITOR) {
546 nCallRefine = 0;
547 discardOrd = 0;
548 discardRefine = 0;
549 discardDisc = 0;
550 discardIso = 0;
551 }
552
553 // work space for isIsomorphic
554 modmat = newMat(nNodes, nNodes, 0);
555 permp = newArray(nNodes, 0);
556 permq = newArray(nNodes, 0);
557 permr = newArray(nNodes, 0);
558
559 // work space for bisearchM
560 bidef = newArray(nNodes, 0);
561 bilow = newArray(nNodes, 0);
562 bicount = 0;
563 DUMMYUSE(padding);
564
565#if DEBUGM
566 printf("+++ new MGraph %d\n", ++countMG);
567 if(countEG > 100) { exit(1); }
568#endif
569}
570
571MGraph::~MGraph()
572{
573 int j;
574
575 deletArray(bilow);
576 deletArray(bidef);
577 deletArray(permr);
578 deletArray(permq);
579 deletArray(permp);
580
581 deleteMat(modmat, nNodes, nNodes);
582 delete curcl;
583 deleteMat(adjMat, nNodes, nNodes);
584 delete egraph;
585 for (j = 0; j < nNodes; j++) {
586 delete nodes[j];
587 }
588 delete[] nodes;
589 delete[] clist;
590
591#if DEBUGM
592 printf("+++ delete MGraph %d\n", countMG++);
593#endif
594}
595
596
597void MGraph::printAdjMat(MNodeClass *cl)
598{
599 int j1, j2;
600
601 cout << " ";
602 for (j2 = 0; j2 < nNodes; j2++) {
603 cout << " " << setw(2) << j2;
604 }
605 cout << endl;
606 for (j1 = 0; j1 < nNodes; j1++) {
607 cout << setw(2) << j1 << ": [";
608 for (j2 = 0; j2 < nNodes; j2++) {
609 cout << " " << setw(2) << adjMat[j1][j2];
610 }
611 cout << "] " << cl->ndcl[j1] << endl;
612 }
613}
614
615//---------------------------------------------------------------
616// Check graph can be a connected one.
617// If a connected component without free leg is not the whole graph then
618// return False, otherwise return True.
619Bool MGraph::isConnected(void)
620{
621 int j, n, nv;
622
623 for (j = 0; j < nNodes; j++) {
624 nodes[j]->visited = -1;
625 }
626 if (visit(0)) {
627 return True;
628 }
629 nv = 0;
630 for (n = 0; n < nNodes; n++) {
631 if (nodes[n]->visited >= 0) {
632 nv++;
633 }
634 }
635 return (nv == nNodes);
636}
637
638// Visiting connected node used for 'isConnected'
639// If child nodes has free legs, then this function returns True.
640// otherwise it returns False.
641Bool MGraph::visit(int nd)
642{
643 int td;
644
645 // This node has free legs.
646 if (nodes[nd]->freelg > 0) {
647 return True;
648 }
649 nodes[nd]->visited = 0;
650 for (td = 0; td < nNodes; td++) {
651 if ((adjMat[nd][td] > 0) and (nodes[td]->visited < 0)) {
652 if (visit(td)) {
653 return True;
654 }
655 }
656 }
657 // all the child nodes has no free legs.
658 return False;
659}
660
661// Check whether the current graph is the Representative of a isomorphic class.
662// nsym = symmetry factor by the permutation of nodes.
663// esym = symmetry factor by the permutation of edge.
664// If this graph is not a representative, then returns False.
665Bool MGraph::isIsomorphic(MNodeClass *cl)
666{
667 int j1, j2, cmp, count, nself;
668
669 nsym = ToBigInt(0);
670 esym = ToBigInt(1);
671
672 count = 0;
673 while (True) {
674 count = nextPerm(nNodes, cl->nClasses, cl->clist,
675 permr, permq, permp, count);
676
677 if (count < 0) {
678 // calculate permutations of edges
679 esym = ToBigInt(1);
680 for (j1 = 0; j1 < nNodes; j1++) {
681 if (adjMat[j1][j1] > 0) {
682 nself = adjMat[j1][j1]/2;
683 esym *= factorial(nself);
684 esym *= ipow(2, nself);
685 }
686 for (j2 = j1+1; j2 < nNodes; j2++) {
687 if (adjMat[j1][j2] > 0) {
688 esym *= factorial(adjMat[j1][j2]);
689 }
690 }
691 }
692 return True;
693 }
694
695 permMat(nNodes, permp, adjMat, modmat);
696 cmp = compMat(nNodes, adjMat, modmat);
697 if (cmp < 0) {
698 return False;
699 } else if (cmp == 0) {
700 // save permutation ???
701 nsym = nsym + ToBigInt(1);
702 }
703 }
704}
705
706// apply permutation to matrix 'mat0' and obtain 'mat1'
707void MGraph::permMat(int size, int *perm, int **mat0, int **mat1)
708{
709 int j1, j2;
710
711 for (j1 = 0; j1 < size; j1++) {
712 for (j2 = 0; j2 < size; j2++) {
713 mat1[j1][j2] = mat0[perm[j1]][perm[j2]];
714 }
715 }
716}
717
718// comparison of matrix
719int MGraph::compMat(int size, int **mat0, int **mat1)
720{
721 int j1, j2, cmp;
722
723 for (j1 = 0; j1 < size; j1++) {
724 for (j2 = 0; j2 < size; j2++) {
725 cmp = mat0[j1][j2] - mat1[j1][j2];
726 if (cmp != 0) {
727 return cmp;
728 }
729 }
730 }
731 return 0;
732}
733
734
735// Refine the classification
736// cl : the current 'MNodeClass' object
737// cn : the class number
738// Returns (the new class number corresponds to 'cn') if OK, or
739// 'None' if ordering condition is not satisfied
740MNodeClass *MGraph::refineClass(MNodeClass *cl, int cn)
741{
742 MNodeClass *ccl = cl;
743 int ccn = cn;
744 MNodeClass *ncl = NULL;
745 int ucl[MAXNODES];
746 int nucl, nce;
747 int td, cmp;
748
749 nucl = 0;
750 while (ccl->nClasses != nucl) { // repeat refinement.
751 nce = 0;
752 nucl = 0;
753 for (td = 1; td < nNodes; td++) {
754 // 'td' is the next node and the current node is 'td-1'.
755 // Count up the number of the elements in the current class
756 // corresponding to the current node
757 nce++;
758 cmp = ccl->clCmp(td-1, td, ccn);
759
760 // the ordering condition is not satisfied.
761 if (cmp > 0) {
762 if (DEBUG) {
763 cout << "refine: cls = " << cl->clist << endl;
764 cl->printMat();
765 cout << "clmat" << endl;
766 ccl->printMat();
767 cout << "refine: discard: cls = " << ccl->clist
768 << "ucl = " << ucl << endl;
769 }
770 if (ccl != cl) {
771 delete ccl;
772 }
773 return NULL;
774
775 } else if (cmp < 0) {
776 // 'td' is in the next class to the current node.
777 ucl[nucl++] = nce; // close the current class
778
779 // start new class
780 nce = 0;
781 }
782 // nothing to do for the case of 'cmp == 0'.
783
784 }
785
786 // close array 'ucl'.
787 ucl[nucl++] = nce + 1;
788
789 // class is not modified
790 if (nucl == ccl->nClasses) {
791 ncl = ccl;
792 break;
793
794 // inconsistent
795 } else if (nucl < ccl->nClasses) {
796 erEnd("refineClasses : smaller number of classes");
797
798 // preparation of the next repetition.
799 } else {
800 if (cn == ccl->nClasses) {
801 td = nNodes;
802 } else {
803 td = ccl->flist[cn+1]-1;
804 }
805 if (ccl != cl) {
806 delete ccl;
807 }
808 ccl = new MNodeClass(nNodes, nucl);
809 ccl->init(ucl, maxdeg, adjMat);
810 if (td == nNodes) {
811 ccn = ccl->nClasses;
812 } else {
813 ccn = ccl->ndcl[td];
814 }
815 nucl = 0;
816 }
817 }
818
819 if (DEBUG) {
820 cout << "refine: ucl = " << ucl << endl;
821 cout << "refine: ncl = " << ncl->clist << endl;
822 ncl->printMat();
823 }
824
825 return ncl;
826}
827
828// Search biconnected component
829// visit : pd --> nd --> td
830// ne : the number of edges between pd and nd.
831void MGraph::bisearchM(int nd, int pd, int ne)
832{
833 int td;
834
835 bidef[nd] = bicount;
836 bilow[nd] = bicount;
837 bicount++;
838
839 for (td = 0; td < nNodes; td++) {
840 if (nodes[td]->ext) { // ignore external node
841 continue;
842
843 } else if (td == nd) { // pd --> nd --> nd
844 continue;
845
846 } else if (td == pd) { // pd --> nd --> pd
847 if (ne > 1) {
848 bilow[nd] = min(bilow[nd], bidef[pd]);
849 }
850
851 } else if (adjMat[td][nd] < 1) { // td is not adjacent to nd
852 continue;
853
854 } else if (bidef[td] >= 0) { // back edge
855 bilow[nd] = min(bilow[nd], bidef[td]);
856
857 // new node
858 } else {
859
860 bisearchM(td, nd, adjMat[td][nd]);
861
862 // ordinary case
863 if (bilow[td] > bidef[nd]) {
864 nBridges++;
865 }
866 bilow[nd] = min(bilow[nd], bilow[td]);
867 }
868 }
869
870 // nd is the starting point and not an external line
871 // if (pd < 0 && (!nodes[nd]->ext || nodes[nd]->deg > 1)) {
872 // if (pd < 0 && !(nodes[nd]->ext && nodes[nd]->deg ==1)) {
873 if (pd < 0 && nodes[nd]->deg != 1) {
874 // nBridges += nodes[nd]->deg;
875 nBridges++;
876 }
877}
878
879// Count the number of 1PI components.
880// The algorithm is described in
881// A.V. Aho, J.E. Hopcroft and J.D. Ullman
882// 'The Design and Analysis of Computer Algorithms', Chap. 5
883// 1974, Addison-Wesley.
884
885int MGraph::count1PI(void)
886{
887 int j;
888
889 if (nLoops < 0) {
890 return 1;
891 }
892
893 // initialization
894 bicount = 0;
895 nBridges = 0;
896 for (j = 0; j < nNodes; j++) {
897 bidef[j] = -1;
898 bilow[j] = -1;
899 }
900
901 bisearchM(0, -1, 0);
902
903 return nBridges;
904}
905
906//---------------------------------------------------------------
907// Generate graphs
908// the generation process starts from 'connectClass'.
909
910long MGraph::generate(void)
911{
912 MNodeClass *cl;
913 int dscl[MAXNODES];
914 int n;
915
916 for (n = 0; n < nNodes; n++) {
917 dscl[n] = False;
918 }
919
920 // Initial classification of nodes.
921 cl = new MNodeClass(nNodes, nClasses);
922 cl->init(clist, maxdeg, adjMat);
923 connectClass(cl, dscl);
924
925 // Print the result.
926
927 if (MONITOR) {
928 cout << endl;
929 cout << endl;
930 cout << "* Total " << ndiag << " Graphs.";
931 cout << "(" << n1PI << " 1PI)";
932 // cout << " wsum = " << wsum << "(" << wsopi << "1PI)" << endl;
933 cout << endl;
934 }
935
936 if (MONITOR) {
937 cout << "* refine: " << nCallRefine << endl;
938 cout << "* discard for ordering: " << discardOrd << endl;
939 cout << "* discard for refinement: " << discardRefine << endl;
940 cout << "* discard for disconnected: " << discardDisc << endl;
941 cout << "* discard for duplication: " << discardIso << endl;
942 }
943 delete cl;
944
945 return ndiag;
946}
947
948int MGraph::findNextCl(MNodeClass *cl, int *dscl)
949{
950 int mine, cr, c, n, me;
951
952 mine = -1;
953 cr = -1;
954 for (c = 0; c < cl->nClasses; c++) {
955 n = cl->flist[c];
956 if (!dscl[n]) {
957 me = nodes[n]->freelg;
958 if (me > 0) {
959 if (nodes[n]->freelg < nodes[n]->deg) {
960 return c;
961 }
962 if (mine < 0 || mine > me) {
963 mine = me;
964 cr = c;
965 }
966 }
967 }
968 }
969 return cr;
970}
971
972int MGraph::findNextTCl(MNodeClass *cl, int *dtcl)
973{
974 int c, n, mine, cr, me;
975
976 mine = -1;
977 cr = -1;
978 for (c = 0; c < cl->nClasses; c++) {
979 for (n = cl->flist[c]; n < cl->flist[c+1]; n++) {
980 if (!dtcl[n]) {
981 me = nodes[n]->freelg;
982 if (me > 0) {
983 if (nodes[n]->freelg < nodes[n]->deg) {
984 return c;
985 }
986 if (mine < 0 || mine > me) {
987 mine = me;
988 cr = c;
989 }
990 }
991 }
992 }
993 }
994 return cr;
995}
996
997// Connect nodes in a class to others
998void MGraph::connectClass(MNodeClass *cl, int *dscl)
999{
1000 int sc, sn;
1001 MNodeClass *xcl;
1002
1003 if (DEBUG0) {
1004 printf("connectClass:begin:");
1005 printf(" dscl="); printArray(nNodes, dscl);
1006 printf("\n");
1007 }
1008 xcl = refineClass(cl, cl->nClasses);
1009
1010 if (xcl == NULL) {
1011 if (MONITOR) {
1012 discardRefine++;
1013 }
1014 } else {
1015 sc = findNextCl(xcl, dscl);
1016 if (sc < 0) {
1017 newGraph(cl);
1018 } else {
1019 sn = xcl->flist[sc];
1020 connectNode(sc, sn, xcl, dscl);
1021 }
1022 }
1023 if (xcl != cl && xcl != NULL) {
1024 delete xcl;
1025 }
1026 if (DEBUG0) {
1027 printf("connectClass:end\n");
1028 }
1029}
1030
1031//------------------------------------------------------------
1032void MGraph::connectNode(int sc, int ss, MNodeClass *cl, int *dscl)
1033{
1034 int sn;
1035 int dtcl[MAXNODES];
1036
1037 if (DEBUG0) {
1038 printf("connectNode:begin:(%d,%d)", sc, ss);
1039 printf(" dscl="); printArray(nNodes, dscl);
1040 printf("\n");
1041 }
1042 if (ss >= cl->flist[sc+1]) {
1043 connectClass(cl, dscl);
1044 if (DEBUG0) {
1045 printf("connectNode:end1\n");
1046 }
1047 return;
1048 }
1049
1050 copyArray(nNodes, dscl, dtcl);
1051 for (sn = ss; sn < cl->flist[sc+1]; sn++) {
1052 if (!dscl[sn]) {
1053 connectLeg(sc, sn, sc, sn, cl, dscl, dtcl);
1054 if (DEBUG0) {
1055 printf("connectNode:end2\n");
1056 }
1057 return;
1058 }
1059 }
1060 if (DEBUG0) {
1061 printf("connectNode:end3\n");
1062 }
1063}
1064
1065// Add one connection between two legs.
1066// 1. select another node to connect
1067// 2. determine multiplicity of the connection
1068//
1069// Arguments
1070// cn : the current class
1071// nd : the current node to be connected.
1072// nextnd : {nextnd, ...} is the possible target node of the connection.
1073// cl : the current node class
1074
1075void MGraph::connectLeg(int sc, int sn, int tc, int ts, MNodeClass *cl, int *dscl, int* dtcl)
1076{
1077 int tn, maxself, nc2, nc, maxcon, ts1, wc, ncm;
1078 int dtcl1[MAXNODES];
1079
1080 if (DEBUG0) {
1081 printf("connectLeg:begin:(%d,%d,%d,%d)", sc, sn, tc, ts);
1082 printf(" dscl="); printArray(nNodes, dscl);
1083 printf(" dtcl="); printArray(nNodes, dtcl);
1084 printf("\n");
1085 printAdjMat(cl);
1086 }
1087
1088 if (sn >= cl->flist[sc+1]) {
1089 erEnd("*** connectLeg : illegal control");
1090 return;
1091
1092 // There remains no free legs in the node 'sn' : move to next node.
1093 } else if (nodes[sn]->freelg < 1) {
1094
1095 if (!isConnected()) {
1096 if (DEBUG0) {
1097 printf("connectLeg:disconnected\n");
1098 }
1099 if (MONITOR) {
1100 discardDisc++;
1101 }
1102 } else {
1103 if (DEBUG0) {
1104 printf("connectLeg: call conNode\n");
1105 }
1106 dscl[sn] = True;
1107
1108 // next node in the current class.
1109 connectNode(sc, sn+1, cl, dscl);
1110
1111 dscl[sn] = False;
1112 }
1113
1114 // connect a free leg of the current node 'sn'.
1115 } else {
1116 copyArray(nNodes, dtcl, dtcl1);
1117
1118 ts1 = ts;
1119 for (wc = 0; wc < nNodes; wc++) {
1120 if (ts1 >= cl->flist[tc+1]) {
1121 tc = findNextTCl(cl, dtcl1);
1122 if (DEBUG0) {
1123 printf("connectLeg:1:tc=%d\n", tc);
1124 }
1125 if (tc < 0) {
1126 if (DEBUG0) {
1127 printf("connectLeg:end1:(%d,%d,%d,%d)\n", sc, sn, tc, ts);
1128 }
1129 return;
1130 }
1131 if (tc != sc && !cl->chkOrd(sn, sc, cl, dtcl)) {
1132 if (MONITOR) {
1133 discardOrd++;
1134 }
1135 if (DEBUG0) {
1136 printf("connectLeg:end2:(%d,%d,%d,%d)\n", sc, sn, tc, ts);
1137 }
1138 return;
1139 }
1140 }
1141
1142 ts1 = -1;
1143
1144 // repeat for all possible target nodes
1145 // for all tn in tc
1146 if (DEBUG0) {
1147 printf("connectLeg:2:tc=%d, fl:%d-->%d\n", tc, cl->flist[tc], cl->flist[tc+1]);
1148 }
1149 for (tn = cl->flist[tc]; tn < cl->flist[tc+1]; tn++) {
1150 if (DEBUG0) {
1151 printf("connectLeg:2:%d=>try %d\n", sn, tn);
1152 }
1153
1154 if (dtcl1[tn]) {
1155 if (DEBUG0) {
1156 printf("connectLeg:3\n");
1157 }
1158 continue;
1159 }
1160 dtcl1[tn] = True;
1161
1162 if (sc == tc && sn > tn) {
1163 if (DEBUG0) {
1164 printf("connectLeg:4\n");
1165 }
1166 continue;
1167
1168 // self-loop
1169 } else if (sn == tn) {
1170 if (nNodes > 1) {
1171 // there are two or more nodes in the graph :
1172 // avoid disconnected graph
1173 maxself = min((nodes[sn]->freelg)/2, (nodes[sn]->deg-1)/2);
1174 } else {
1175 // there is only one node the graph.
1176 maxself = nodes[sn]->freelg/2;
1177 }
1178
1179 // If we can assume no tadpole, the following line can be used.
1180 if (selOPI and nNodes > 2) {
1181 maxself = min((nodes[sn]->deg-2)/2, maxself);
1182 }
1183
1184 // vary the number of connection.
1185 for (nc2 = maxself; nc2 > 0; nc2--) {
1186 // if nNodes > 1 and ndg == nc2:
1187 // // disconnected
1188 // continue
1189 nc = 2*nc2;
1190 if (DEBUG0) {
1191 printf("connectLeg: call conLeg: (same) %d=>%d(%d)\n", sn, tn,nc);
1192 }
1193 ncm = 5*nc;
1194 adjMat[sn][sn] = nc;
1195 nodes[sn]->freelg -= nc;
1196 cl->incMat(sn, tn, ncm);
1197 ts1 = tn + 1;
1198
1199 // next connection
1200 connectLeg(sc, sn, tc, ts1, cl, dscl, dtcl1);
1201
1202 // restore the configuration
1203 cl->incMat(tn, sn, - ncm);
1204 adjMat[sn][sn] = 0;
1205 nodes[sn]->freelg += nc;
1206 if (DEBUG0) {
1207 printf("connectLeg: ret conLeg: (same) %d=>%d(%d)\n", sn, tn,nc);
1208 }
1209 }
1210
1211 // connections between different nodes.
1212 } else {
1213 // maximum possible connection number
1214 maxcon = min(nodes[sn]->freelg, nodes[tn]->freelg);
1215
1216 // avoid disconnected graphs.
1217 if (nNodes > 2 && nodes[sn]->deg == nodes[tn]->deg) {
1218 maxcon = min(maxcon, nodes[sn]->deg-1);
1219 }
1220
1221 if (CHECK) {
1222 if ((adjMat[sn][tn] != 0) || (adjMat[sn][tn] != adjMat[tn][sn])) {
1223 printf("*** inconsistent connection: sn=%d, tn=%d", sn, tn);
1224 printf(" dscl="); printArray(nNodes, dscl);
1225 printf(" dtcl="); printArray(nNodes, dtcl);
1226 printf("\n");
1227 printAdjMat(cl);
1228 erEnd("*** inconsistent connection ");
1229 }
1230 }
1231
1232 // vary number of connections
1233 for (nc = maxcon; nc > 0; nc--) {
1234 if (DEBUG0) {
1235 printf("connectLeg: call conLeg: (diff) %d=>%d(%d)", sn, tn,nc);
1236 printf(" dtcl="); printArray(nNodes, dtcl);
1237 printf("\n");
1238 }
1239 adjMat[sn][tn] = nc;
1240 adjMat[tn][sn] = nc;
1241 nodes[sn]->freelg -= nc;
1242 nodes[tn]->freelg -= nc;
1243 cl->incMat(sn, tn, nc);
1244 cl->incMat(tn, sn, nc);
1245 ts1 = tn + 1;
1246
1247 // next connection
1248 connectLeg(sc, sn, tc, ts1, cl, dscl, dtcl1);
1249
1250 // restore configuration
1251 cl->incMat(sn, tn, - nc);
1252 cl->incMat(tn, sn, - nc);
1253 adjMat[sn][tn] = 0;
1254 adjMat[tn][sn] = 0;
1255 nodes[sn]->freelg += nc;
1256 nodes[tn]->freelg += nc;
1257 if (DEBUG0) {
1258 printf("connectLeg: ret conLeg: (diff) %d=>%d(%d)\n", sn, tn,nc);
1259 printf(" dtcl="); printArray(nNodes, dtcl);
1260 printAdjMat(cl);
1261 printf("\n");
1262 }
1263 }
1264 }
1265 }
1266 if (ts1 < 0) {
1267 ts1 = cl->flist[tc+1];
1268 }
1269 } // for wc
1270 }
1271 if (DEBUG0) {
1272 printf("connectLeg:end3:(%d,%d,%d,%d)\n", sc, sn, tc, ts);
1273 }
1274}
1275
1276// A new candidate diagram is obtained.
1277// It may be a disconnected graph
1278
1279void MGraph::newGraph(MNodeClass *cl)
1280{
1281 int connected;
1282 MNodeClass *xcl;
1283 Bool sopi;
1284
1285 ngen++;
1286 // printf("newGraph: %d\n", ngen);
1287
1288 // refine class and check ordering condition
1289 xcl = refineClass(cl, cl->nClasses);
1290 if (xcl == NULL) {
1291 // printf("newGraph: fail refine\n");
1292 if (MONITOR) {
1293 discardRefine++;
1294 }
1295
1296 // check whether this is a connected graph or not.
1297 } else {
1298 connected = isConnected();
1299 if (!connected) {
1300 // printf("newGraph: not connected\n");
1301 if (DEBUG0) {
1302 cout << "+++ disconnected graph" << ngen << endl;
1303 xcl->printMat();
1304 }
1305
1306 if (MONITOR) {
1307 discardDisc++;
1308 }
1309
1310 // connected graph : check isomorphism of the graph
1311 } else {
1312 ngconn++;
1313 if (!isIsomorphic(xcl)) {
1314 // printf("newGraph: not isomorphic\n");
1315 if (DEBUG) {
1316 cout << "+++ duplicated graph" << ngen << ngconn << endl;
1317 xcl->printMat();
1318 }
1319
1320 if (MONITOR) {
1321 discardIso++;
1322 }
1323
1324 // We got a new connected and a unique representation of a class.
1325 } else {
1326
1327 c1PI = count1PI();
1328 if (!selOPI || c1PI == 1) {
1329 wsum = wsum + ToFraction(1, nsym*esym);
1330 if (c1PI == 1) {
1331 wsopi = wsopi + ToFraction(1, nsym*esym);
1332 }
1333 curcl->copy(xcl);
1334 ndiag++;
1335 sopi = (c1PI == 1);
1336 if (sopi) {
1337 n1PI++;
1338 }
1339
1340 if (OPTPRINT) {
1341 cout << endl;
1342 cout << "Graph :" << ndiag
1343 << "(" << ngen << ")"
1344 << " 1PI comp. = " << c1PI
1345 << " sym. factor = (" << nsym << "*" << esym << ")"
1346 << endl;
1347 printAdjMat(cl);
1348 // cl->printMat();
1349 if (MONITOR) {
1350 cout << "refine: " << nCallRefine << endl;
1351 cout << "discard for ordering: " << discardOrd << endl;
1352 cout << "discard for refinement: " << discardRefine << endl;
1353 cout << "discard for disconnected: " << discardDisc << endl;
1354 cout << "discard for duplication: " << discardIso << endl;
1355 }
1356 }
1357
1358 // go to next step
1359 egraph->init(pId, ndiag, adjMat, sopi, nsym, esym);
1360 toForm(egraph);
1361 }
1362 }
1363 }
1364 }
1365// printf("in newGraph cl=%p, xcl=%p\n", cl, xcl);
1366 if (xcl != cl && xcl != NULL) {
1367 delete xcl;
1368 }
1369}
1370
1371
1372// go to next step
1373// 1. convert data format which treat the edges as objects.
1374// 2. definition of loop momenta to the edges.
1375// 3. analysis of loop structure in a graph.
1376// 4. assignment of particles to edges and vertices to nodes
1377// 5. recalculate symmetry factor
1378
1379// Permutation
1380Local int nextPerm(int nelem, int nclass, int *cl, int *r, int *q, int *p, int count)
1381{
1382 int j, k, n, e, t;
1383 Bool b;
1384
1385 for (j = 0; j < nelem; j++) {
1386 p[j] = j;
1387 }
1388 if (count < 1) {
1389 for (j = 0; j < nelem; j++) {
1390 q[j] = 0;
1391 r[j] = 0;
1392 }
1393 j = 0;
1394 for (k = 0; k < nclass; k++) {
1395 n = cl[k];
1396 for (e = 0; e < n; e++) {
1397 r[j] = n - e - 1;
1398 j++;
1399 }
1400 }
1401 if (j != nelem) {
1402 erEnd("*** inconsistent # elements");
1403 }
1404 return 1;
1405 }
1406 b = False;
1407 for (j = nelem-1; j >= 0; j--) {
1408 if (q[j] < r[j]) {
1409 for (k = j+1; k < nelem; k++) {
1410 q[k] = 0;
1411 }
1412 q[j]++;
1413 b = True;
1414 break;
1415 }
1416 }
1417 if (!b) {
1418 return (-count);
1419 }
1420
1421 for (j = 0; j < nelem; j++) {
1422 k = j + q[j];
1423 t = p[j];
1424 p[j] = p[k];
1425 p[k] = t;
1426 }
1427 return count + 1;
1428}
1429
1430Local BigInt factorial(int n)
1431/* return (n < 1) ? 1 : n!
1432 */
1433{
1434 int r, j;
1435
1436 r = 1;
1437 for (j = 2; j <= n; j++) {
1438 r *= j;
1439 }
1440 return r;
1441}
1442
1443Local BigInt ipow(int n, int p)
1444{
1445 int r, j;
1446 DUMMYUSE(p);
1447
1448 r = 1;
1449 for (j = 2; j <= n; j++) {
1450 r *= n;
1451 }
1452 return r;
1453}
1454
1455
1456Local void erEnd(const char *msg)
1457{
1458 printf("*** Error : %s\n", msg);
1459 exit(1);
1460}
1461
1462/* memory allocation */
1463Local int *newArray(int size, int val)
1464{
1465 int *a, j;
1466
1467 a = new int[size];
1468 for (j = 0; j < size; j++) {
1469 a[j] = val;
1470 }
1471 return a;
1472}
1473
1474Local void deletArray(int *a)
1475{
1476 delete[] a;
1477}
1478
1479Local void copyArray(int size, int *a0, int *a1)
1480{
1481 int j;
1482 for (j = 0; j < size; j++) {
1483 a1[j] = a0[j];
1484 }
1485}
1486
1487Local int **newMat(int n0, int n1, int val)
1488{
1489 int **m, j;
1490
1491 m = new int*[n0];
1492 for (j = 0; j < n0; j++) {
1493 m[j] = newArray(n1, val);
1494 }
1495 return m;
1496}
1497
1498Local void deleteMat(int **m, int n0, int n1)
1499{
1500 int j;
1501 DUMMYUSE(n1);
1502
1503 for (j = 0; j < n0; j++) {
1504 deletArray(m[j]);
1505 }
1506 delete[] m;
1507}
1508
1509//==============================================================
1510// Functions for testing the program.
1511//
1512//--------------------------------------------------------------
1513// Testing function nextPerm
1514//
1515Local void printArray(int n, int *p)
1516{
1517 int j;
1518
1519 printf("[");
1520 for (j = 0; j < n; j++) {
1521 printf(" %2d", p[j]);
1522 }
1523 printf("]");
1524}
1525
1526Local void printMat(int n0, int n1, int **m)
1527{
1528 int j;
1529
1530 for (j = 0; j < n0; j++) {
1531 printArray(n1, m[j]);
1532 printf("\n");
1533 }
1534}
1535
1536Global void testPerm()
1537{
1538 int nelem, nperm, nclist, n, count;
1539 int *p, *q, *r;
1540 int clist[] = {1, 2, 2, 3};
1541
1542 nclist = sizeof(clist)/sizeof(int);
1543 nelem = 0;
1544 nperm = 1;
1545 for (n = 0; n < nclist; n++) {
1546 nelem += clist[n];
1547 nperm *= factorial(clist[n]);
1548 }
1549
1550 printf("+++ clist = (%d) ", nclist);
1551 printArray(nclist, clist);
1552 printf("\n");
1553 printf("+++ nelem = %d, nperm = %d\n", nelem, nperm);
1554
1555 p = new int[nelem];
1556 q = new int[nelem];
1557 r = new int[nelem];
1558 count = 0;
1559 while (True) {
1560 count = nextPerm(nelem, nclist, clist, r, q, p, count);
1561 if (count < 0) {
1562 count = - count;
1563 break;
1564 }
1565 printf("%4d:", count);
1566 printArray(nelem, p);
1567 printf("\n");
1568
1569 if (count > nperm) {
1570 break;
1571 }
1572 }
1573 if (count != nperm) {
1574 printf("*** %d != %d\n", count, nperm);
1575 }
1576 delete[] p;
1577 delete[] q;
1578 delete[] r;
1579}
1580
1581// } ) ]
1582
#define CHECK(condition)
Definition parallel.c:153