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 
20 using 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 
32 typedef int Bool;
33 const int True = 1;
34 const int False = 0;
35 
36 // compile options
37 const int CHECK = True;
38 //const int MONITOR = True;
39 const int MONITOR = False;
40 
41 const int OPTPRINT = False;
42 
43 const int DEBUG0 = False;
44 //const int DEBUG1 = False;
45 const int DEBUG = False;
46 
47 // for debugging memory use
48 #define DEBUGM False
49 
50 //==============================================================
51 class MGraph;
52 class EGraph;
53 
54 #define Extern
55 #define Global
56 #define Local
57 
58 // External functions
59 Extern void toForm(EGraph *egraph);
60 Extern int countPhiCon(int ex, int lp, int v4);
61 
62 // Temporal definition: they should be replaced by FROM functions.
63 Local BigInt factorial(int n);
64 Local BigInt ipow(int n, int p);
65 
66 // Local functions
67 Local int nextPerm(int nelem, int nclass, int *cl, int *r, int *q, int *p, int count);
68 
69 Local void erEnd(const char *msg);
70 Local int *newArray(int size, int val);
71 Local void deletArray(int *a);
72 Local void copyArray(int size, int *a0, int *a1);
73 Local int **newMat(int n0, int n1, int val);
74 Local void deleteMat(int **m, int n0, int n1);
75 Local void printArray(int n, int *p);
76 Local void printMat(int n0, int n1, int **m);
77 
78 #if DEBUGM
79 static int countNMC = 0;
80 static int countEG = 0;
81 static int countMG = 0;
82 #endif
83 
84 //==============================================================
85 // class EGraph
86 
87 EGraph::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 
108 EGraph::~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
123 void 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'
180 void EGraph::setExtern(int nd, Bool val)
181 {
182  nodes[nd].ext = val;
183 }
184 
185 // end of calling 'setExtern'
186 void 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
199 void 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 
236 MNode::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 //
248 class 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 
274 MNodeClass::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 
289 MNodeClass::~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 
301 void 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 
314 void 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]
335 void 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)
349 void 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)
363 int 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'.
383 void 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'.
397 void 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.
404 Bool 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.
434 void 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 
457 int 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 
472 MGraph::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 
571 MGraph::~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 
597 void 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.
619 Bool 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.
641 Bool 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.
665 Bool 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'
707 void 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
719 int 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
740 MNodeClass *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.
831 void 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 
885 int 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 
910 long 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 
948 int 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 
972 int 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
998 void 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 //------------------------------------------------------------
1032 void 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 
1075 void 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 
1279 void 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
1380 Local 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 
1430 Local 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 
1443 Local 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 
1456 Local void erEnd(const char *msg)
1457 {
1458  printf("*** Error : %s\n", msg);
1459  exit(1);
1460 }
1461 
1462 /* memory allocation */
1463 Local 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 
1474 Local void deletArray(int *a)
1475 {
1476  delete[] a;
1477 }
1478 
1479 Local 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 
1487 Local 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 
1498 Local 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 //
1515 Local 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 
1526 Local 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 
1536 Global 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 
Definition: gentopo.h:61
Definition: gentopo.h:34
Definition: gentopo.h:17
#define CHECK(condition)
Definition: parallel.c:153
Definition: gentopo.h:84
Definition: gentopo.h:24