87#define CLWIGHTD(x) (5*(x))
88#define CLWIGHTO(x) (3*(x)-2)
91#define MASK(n) ((1ul) << (n))
93#define BOOLSTR(x) ((x) ? "True" : "False")
99 {
"Step",
"Generate particle assigned graphs", GRCC_AGraph, 0},
100 {
"Outgrf",
"Output to file (out.grf)", False, 0},
101 {
"Outgrp",
"Output to file (out.grp)", False, 0},
102 {
"OPI",
"Only 1PI graphs", True , 0},
103 {
"NoSelfLoop",
"Exclude graphs with loops consist of 1 edge", True , 0},
104 {
"NoTadpole",
"Exclude graphs with tadpoles (2 edge connected)", True , 0},
105 {
"No1PtBlock",
"Exclude graphs with tadpole blocks (2 node connected)", False, 0},
106 {
"No2PtL1PI",
"Exclude graphs with 2-point subgraphs", False, 0},
107 {
"NoExtSelf",
"Exclude graphs with 2-pt subgraphs at ext. particles", False, 0},
108 {
"NoAdj2PtV",
"Exclude graphs with an edge connecting 2-pt vertices", False, 0},
109 {
"Block",
"Exclude graphs with more than one block", False, 0},
110 {
"NoMultiEdge",
"Exclude graphs with multi-edges", False, 0},
111 {
"SymmInitial",
"Symmetrize initial particles", False, 0},
112 {
"SymmFinal",
"Symmetrize final particles", False, 0},
114static OptDef optDef1[] = {
115 {
"Step",
"Generate particle assigned graphs", GRCC_AGraph, 0},
116 {
"Outgrf",
"Ouput to file (out.grf)", False, 0},
117 {
"Outgrp",
"Ouput to file (out.grp)", False, 0},
118 {
"OPI",
"Only 1PI graphs", False, 0},
119 {
"NoSelfLoop",
"Exclude graphs with loops consist of 1 edge", False, 0},
120 {
"NoTadpole",
"Exclude graphs with tadpoles (2 edge connected)", False, 0},
121 {
"No1PtBlock",
"Exclude graphs with tadpole blocks (2 node connected)", False, 0},
122 {
"No2PtL1PI",
"Exclude graphs with 2-point subgraphs", False, 0},
123 {
"NoExtSelf",
"Exclude graphs with 2-pt subgraphs at ext. particles", False, 0},
124 {
"NoAdj2PtV",
"Exclude graphs with an edge connecting 2-pt vertices", False, 0},
125 {
"Block",
"Exclude graphs with more than one block", False, 0},
126 {
"NoMultiEdge",
"Exclude graphs with multi-edges", False, 0},
127 {
"SymmInitial",
"Symmetrize initial particles", False, 0},
128 {
"SymmFinal",
"Symmetrize final particles", False, 0},
131static OptDef *optDef = &(optDef0[0]);
132static int nOptDef =
sizeof(optDef0)/
sizeof(
OptDef);
134static OptDef *optDef = &(optDef1[0]);
135static int nOptDef =
sizeof(optDef1)/
sizeof(
OptDef);
140 {
"onepi",
"onepr",
"one-particle illreducible"},
141 {
"onshell",
"offshell",
"without self-energy part at external particles"},
142 {
"nosigma",
"sigma",
"no 2-point functions"},
143 {
"nosnail",
"snail",
"without snail"},
144 {
"notadpole",
"tadpole",
"without tadpole"},
145 {
"simple",
"notsimple",
"without self-loop nor multi-edge"},
146 {
"bipart",
"nonbipart",
"only bipartite graph"},
147 {
"cycli",
"cyclr",
"???"},
148 {
"floop",
"",
"without fermion loops of odd length"},
149#ifdef GRCC_QGRAF_OPT_TOPOL
150 {
"topol",
"?",
"topology"},
154static int nOptQGDef =
sizeof(optQGDef)/
sizeof(
OptQGDef);
156static int prlevel = 2;
157static ErExit *erExit = NULL;
158static void *erExitArg = NULL;
164static void grcc_fprintf(FILE* out,
const char* fmt, ...);
165static void erEnd(
const char *msg);
166static Bool nextPart(
int nelem,
int nclist,
int *clist,
int *nl,
int *r);
167static void prilist(
int n,
const int *a,
const char *msg);
168static int *newArray(
int size,
int val);
169static int *deleteArray(
int *a);
170static int **newMat(
int n0,
int n1,
int val);
171static int **deleteMat(
int **m,
int n0);
172static void prIntArray(
int n,
int *p,
const char *msg);
173static void prIntArrayErr(
int n,
int *p,
const char *msg);
174static int *intdup(
int n,
int *v);
175static int *delintdup(
int *v);
176static void bsort(
int n,
int *ord,
int *a);
177static void sorti(
int n,
int *a);
178static int sortb(
int n,
int *a);
179static int toSList(
int n,
int *a);
180static int intSetAdd(
int n,
int *a,
int v,
const int size);
181static int intSListAdd(
int n,
int *a,
int v,
const int size);
182static int cmpArray(
int na,
int *a,
int nb,
int *b);
183static Bool leqArray(
int n,
int *a,
int *b);
184static void prMomStr(
int mom,
const char *ms,
int mn);
185static void listDiff(
int na,
int *a,
int nb,
int *b,
int *np,
int *p,
int *nq,
int *q,
int *nr,
int *r);
186static Bool isSubSList(
int na,
int *a,
int nb,
int *b);
187static int subtrSet(
int na,
int *a,
int nb,
int *b,
int *c,
int size);
188static int nextPerm(
int nelem,
int nclass,
int *cl,
int *r,
int *q,
int *p,
int count);
189static BigInt ipow(
int n,
int p);
190static BigInt factorial(
int n);
194static Bool isIn(
int n,
int *a,
int v);
200Options::Options(
void)
202 if (nOptDef != GRCC_OPT_Size) {
203 grcc_fprintf(GRCC_Stderr,
"*** Options: inconsistent default values\n");
204 grcc_fprintf(GRCC_Stderr,
"nOptDef=%d, GRCC_OPT_Size=%d\n",
205 nOptDef, GRCC_OPT_Size);
224 if (nOptQGDef != GRCC_QGRAF_OPT_Size) {
225 grcc_fprintf(GRCC_Stderr,
"*** Options: inconsistent default values\n");
226 grcc_fprintf(GRCC_Stderr,
"nOptQGDef=%d, GRCC_QGRAF_OPT_Size=%d\n",
227 nOptQGDef, GRCC_QGRAF_OPT_Size);
231 for (
int j = 0; j < nOptQGDef; j++) {
232 qgref[nqgopt].name = optQGDef[j].name;
233 qgref[nqgopt].index = j;
234 qgref[nqgopt].sign = +1;
236 if (strlen(optQGDef[j].cname) > 0) {
237 qgref[nqgopt].name = optQGDef[j].cname;
238 qgref[nqgopt].index = j;
239 qgref[nqgopt].sign = -1;
245 for (
int j = 0; j < GRCC_QGRAF_OPT_Size; j++) {
261Options::~Options(
void)
273void Options::setDefaultValues(
void)
278 for (j = 0; j < GRCC_OPT_Size; j++) {
279 values[j] = optDef[j].defaultv;
282 values[GRCC_OPT_Step] = GRCC_AGraph;
289void Options::setOldDefaultValues(
void)
294 for (j = 0; j < GRCC_OPT_Size; j++) {
295 values[j] = optDef0[j].defaultv;
298 values[GRCC_OPT_Step] = GRCC_AGraph;
305void Options::setOutMG(OutEGB *omg,
void *pt)
312void Options::setEndMG(OutEGB *emg,
void *pt)
319void Options::setOutAG(OutEG *oag,
void *pt)
326void Options::setErExit(ErExit *ere,
void *erearg)
333void Options::printLevel(
int l)
339void Options::setValue(
int ind,
int val)
341 if (0 <= ind && ind < GRCC_OPT_Size) {
345 grcc_fprintf(GRCC_Stderr,
"*** Options::setValue : invalid index=%d ",
347 grcc_fprintf(GRCC_Stderr,
"(val=%d)\n", val);
349 erEnd(
"Options::setValue : invalid index");
354int Options::getValue(
int ind)
356 if (0 <= ind && ind < GRCC_OPT_Size) {
360 grcc_fprintf(GRCC_Stderr,
"*** Options::getValue : invalid index=%d\n", ind);
362 erEnd(
"Options::getValue : invalid index");
368void Options::setQGrafOpt(
int *qg)
373 values[GRCC_OPT_1PI] = False;
374 values[GRCC_OPT_NoSelfLoop] = False;
375 values[GRCC_OPT_NoTadpole] = False;
376 values[GRCC_OPT_No1PtBlock] = False;
377 values[GRCC_OPT_No2PtL1PI] = False;
378 values[GRCC_OPT_NoExtSelf] = False;
379 values[GRCC_OPT_NoAdj2PtV] = False;
380 values[GRCC_OPT_Block] = False;
381 values[GRCC_OPT_SymmInitial] = False;
382 values[GRCC_OPT_SymmFinal] = False;
385 for (
int j = 0; j < GRCC_QGRAF_OPT_Size; j++) {
390 if (qgopt[GRCC_QGRAF_OPT_ONEPI] > 0) {
391 values[GRCC_OPT_1PI] = True;
392 }
else if (qgopt[GRCC_QGRAF_OPT_ONEPI] < 0) {
393 values[GRCC_OPT_1PI] = False;
397 if (qgopt[GRCC_QGRAF_OPT_ONSHELL] > 0) {
398 values[GRCC_OPT_NoExtSelf] = 1;
399 }
else if (qgopt[GRCC_QGRAF_OPT_ONSHELL] < 0) {
400 values[GRCC_OPT_NoExtSelf] = 0;
406 if (qgopt[GRCC_QGRAF_OPT_NOSNAIL] > 0) {
407 values[GRCC_OPT_NoTadpole] = True;
408 values[GRCC_OPT_No1PtBlock] = True;
409 }
else if (qgopt[GRCC_QGRAF_OPT_NOSNAIL] < 0) {
410 values[GRCC_OPT_NoTadpole] = False;
411 values[GRCC_OPT_No1PtBlock] = False;
415 if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] > 0) {
416 values[GRCC_OPT_NoTadpole] = True;
417 }
else if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] < 0) {
418 values[GRCC_OPT_NoTadpole] = False;
422 if (qgopt[GRCC_QGRAF_OPT_SIMPLE] > 0) {
423 values[GRCC_OPT_NoSelfLoop] = True;
424 values[GRCC_OPT_NoMultiEdge] = True;
425 }
else if (qgopt[GRCC_QGRAF_OPT_SIMPLE] < 0) {
426 values[GRCC_OPT_NoSelfLoop] = False;
427 values[GRCC_OPT_NoMultiEdge] = False;
433#ifdef GRCC_QGRAF_OPT_TOPOL
435 if (qgopt[GRCC_QGRAF_OPT_TOPOL] > 0) {
436 values[GRCC_OPT_SymmInitial] = True;
437 values[GRCC_OPT_SymmFinal] = True;
443void Options::print(
void)
447 grcc_fprintf(GRCC_Stdout,
"Options\n");
448 grcc_fprintf(GRCC_Stdout,
"+++ GRCC_OPT_Size=%d, print level=%d: ",
449 GRCC_OPT_Size, prlevel);
450 grcc_fprintf(GRCC_Stdout,
"symbol = value (default)\n");
451 for (j=0; j < GRCC_OPT_Size; j++) {
452 grcc_fprintf(GRCC_Stdout,
" %4d GRCC_OPT_%-15s = %2d (%2d)\n",
453 j, optDef[j].name, values[j], optDef[j].defaultv);
455 grcc_fprintf(GRCC_Stdout,
" outgrf=%s, outgrp=%s\n", out->outgrf, out->outgrp);
456 grcc_fprintf(GRCC_Stdout,
" GRCC_QGRAF_OPT_Size=%d:\n", GRCC_QGRAF_OPT_Size);
457 for (j=0; j < GRCC_QGRAF_OPT_Size; j++) {
458 grcc_fprintf(GRCC_Stdout,
" %4d %-10s = %2d\n", j, optQGDef[j].name, qgopt[j]);
463const OptDef *Options::getDef(
void)
469const OptDef *Options::getOldDef(
void)
475const OptQGDef *Options::getQGDef(
void)
481void Options::setOutputF(Bool outf,
const char *fname)
483 values[GRCC_OPT_Outgrf] = outf;
484 out->setOutgrf(fname);
488void Options::setOutputP(Bool outp,
const char *fname)
490 values[GRCC_OPT_Outgrp] = outp;
491 out->setOutgrp(fname);
495void Options::printModel(
void)
501 grcc_fprintf(GRCC_Stdout,
"*** model is not defined\n");
507void Options::outModel(
void)
509 if (out != NULL && model != NULL) {
516void Options::begin(
Model *mdl)
520 if (values[GRCC_OPT_Outgrf]) {
521 out->outBeginF(model, values[GRCC_OPT_Outgrf]);
526 if (values[GRCC_OPT_Outgrp]) {
527 out->outBeginP(model, values[GRCC_OPT_Outgrp]);
536void Options::end(
void)
539 grcc_fprintf(GRCC_Stdout,
"Optimization: ");
541 grcc_fprintf(GRCC_Stdout,
"SIMPSEL=1 ");
543 grcc_fprintf(GRCC_Stdout,
"SIMPSEL=0 ");
546 grcc_fprintf(GRCC_Stdout,
"MINMAXLEG=1 ");
548 grcc_fprintf(GRCC_Stdout,
"MINMAXLEG=0 ");
551 grcc_fprintf(GRCC_Stdout,
"OPTEXTONLY=1 ");
553 grcc_fprintf(GRCC_Stdout,
"OPTEXTONLY=0 ");
555 grcc_fprintf(GRCC_Stdout,
"\n");
558 if (out != NULL && values[GRCC_OPT_Outgrf]) {
561 if (out != NULL && values[GRCC_OPT_Outgrp]) {
568void Options::beginProc(
Process *prc)
571 if (out != NULL && values[GRCC_OPT_Outgrf]) {
572 out->outProcBeginF(prc);
574 if (out != NULL && values[GRCC_OPT_Outgrp]) {
575 out->outProcBeginP(prc);
584void Options::endProc(
void)
589 if (out != NULL && values[GRCC_OPT_Outgrf]) {
592 if (out != NULL && values[GRCC_OPT_Outgrp]) {
598 grcc_fprintf(GRCC_Stdout,
"\n");
599 grcc_fprintf(GRCC_Stdout,
"+++ Proc %d: ext=%d, loop=%d, ",
600 proc->id, proc->nExtern, proc->loop);
602 grcc_fprintf(GRCC_Stdout,
"order=");
603 prIntArray(model->ncouple, proc->clist,
": ");
604 model->prParticleArray(proc->ninitl, proc->initlPart,
"-->");
605 model->prParticleArray(proc->nfinal, proc->finalPart,
"");
607 grcc_fprintf(GRCC_Stdout,
" (%8.2f sec)\n", proc->sec);
610 grcc_fprintf(GRCC_Stdout,
" Proc %d: Total M-Graphs=%ld, M-Graphs=",
611 proc->id, proc->nMGraphs);
612 proc->wMGraphs.print(
" (Conn)\n");
614 grcc_fprintf(GRCC_Stdout,
" Proc %d: Total M-Graphs=%ld, M-Graphs=",
615 proc->id, proc->nMOPI);
616 proc->wMOPI.print(
" (1PI)\n");
618 grcc_fprintf(GRCC_Stdout,
" Proc %d: Total A-Graphs=%ld, A-Graphs=",
619 proc->id, proc->nAGraphs);
620 proc->wAGraphs.print(
" (Conn)\n");
622 grcc_fprintf(GRCC_Stdout,
" Proc %d: Total A-Graphs=%ld, A-Graphs=",
623 proc->id, proc->nAOPI);
624 proc->wAOPI.print(
" (1PI)\n");
626 grcc_fprintf(GRCC_Stdout,
"# { %d,{", proc->ninitl);
627 for (k = 0; k < proc->ninitl; k++) {
629 grcc_fprintf(GRCC_Stdout,
", ");
631 if (proc->model != NULL) {
632 grcc_fprintf(GRCC_Stdout,
"\"%s\"", proc->model->particleName(proc->initlPart[k]));
634 grcc_fprintf(GRCC_Stdout,
"%d", proc->initlPart[k]);
637 grcc_fprintf(GRCC_Stdout,
"}, %d,{", proc->nfinal);
638 for (k = 0; k < proc->nfinal; k++) {
640 grcc_fprintf(GRCC_Stdout,
", ");
642 if (proc->model != NULL) {
643 grcc_fprintf(GRCC_Stdout,
"\"%s\"", proc->model->particleName(proc->finalPart[k]));
645 grcc_fprintf(GRCC_Stdout,
"%d", proc->initlPart[k]);
648 grcc_fprintf(GRCC_Stdout,
"}, ");
650 grcc_fprintf(GRCC_Stdout,
"{");
651 for (k = 0; k < model->ncouple; k++) {
653 grcc_fprintf(GRCC_Stdout,
", ");
655 grcc_fprintf(GRCC_Stdout,
"%d", proc->clist[k]);
657 grcc_fprintf(GRCC_Stdout,
"},");
659 grcc_fprintf(GRCC_Stdout,
"},%6ldL,%6ldL,%3ldL, -1.0, %4.2f},\n",
660 proc->nAOPI, proc->wAOPI.num, proc->wAOPI.den, proc->sec);
663 if (out != NULL && values[GRCC_OPT_Outgrf]) {
666 if (out != NULL && values[GRCC_OPT_Outgrp]) {
673void Options::beginSubProc(
SProcess *sprc)
678 if (out != NULL && values[GRCC_OPT_Outgrf]) {
679 out->outSProcBeginF(sprc);
681 if (out != NULL && values[GRCC_OPT_Outgrp]) {
682 out->outSProcBeginP(sprc);
694void Options::endSubProc(
void)
699 if (out != NULL && values[GRCC_OPT_Outgrf]) {
702 if (out != NULL && values[GRCC_OPT_Outgrp]) {
707 mgraph = sproc->mgraph;
709 sproc->resultMGraph(mgraph->cDiag, mgraph->wscon,
710 mgraph->c1PI, mgraph->wsopi);
714 grcc_fprintf(GRCC_Stdout,
"\n");
715 grcc_fprintf(GRCC_Stdout,
"+++ Subproc %d: ext=%d, loop=%d, nodes=%d, edges=%d\n",
716 sproc->id, sproc->nExtern, sproc->loop,
717 sproc->nNodes, sproc->nEdges);
719 grcc_fprintf(GRCC_Stdout,
" Subproc %d: Total M-Graphs=%ld, M-Wsum=",
720 sproc->id, sproc->nMGraphs);
721 sproc->wMGraphs.print(
" (Conn)\n");
723 grcc_fprintf(GRCC_Stdout,
" Subproc %d: Total M-Graphs=%ld, M-Wsum=",
724 sproc->id, sproc->nMOPI);
725 sproc->wMOPI.print(
" (1PI)\n");
727 grcc_fprintf(GRCC_Stdout,
" Subproc %d: Total A-Graphs=%ld, A-Wsum=",
728 sproc->id, sproc->nAGraphs);
729 sproc->wAGraphs.print(
" (Conn)\n");
731 grcc_fprintf(GRCC_Stdout,
" Subproc %d: Total A-Graphs=%ld, A-Wsum=",
732 sproc->id, sproc->nAOPI);
733 sproc->wAOPI.print(
" (1PI)\n");
736 grcc_fprintf(GRCC_Stdout,
"\n");
737 grcc_fprintf(GRCC_Stdout,
"* Total %ld MGraphs; %ld 1PI", mgraph->cDiag, mgraph->c1PI);
738 grcc_fprintf(GRCC_Stdout,
" wscon = ");
739 mgraph->wscon.print(
" ( ");
740 mgraph->wsopi.print(
" 1PI)\n");
743 grcc_fprintf(GRCC_Stdout,
"* refine: %ld\n", mgraph->nCallRefine);
744 grcc_fprintf(GRCC_Stdout,
"* discarded for refinement: %ld\n", mgraph->discardRefine);
745 grcc_fprintf(GRCC_Stdout,
"* discarded for disconnected: %ld\n", mgraph->discardDisc);
746 grcc_fprintf(GRCC_Stdout,
"* discarded for duplication: %ld\n", mgraph->discardIso);
753 if (out != NULL && values[GRCC_OPT_Outgrf]) {
756 if (out != NULL && values[GRCC_OPT_Outgrp]) {
762void Options::newMGraph(
MGraph *mgr)
768 mgr->egraph->mId = proc->mgrcount;
769 }
else if (sproc != NULL) {
771 mgr->egraph->mId = sproc->mgrcount;
775 && values[GRCC_OPT_Step] == GRCC_MGraph) {
776 if (values[GRCC_OPT_Outgrf]) {
777 out->outEGraphF(mgraph->egraph);
779 if (values[GRCC_OPT_Outgrp]) {
780 out->outEGraphP(mgraph->egraph);
784 sproc->endMGraph(mgr);
787 (*endmg)(mgr->egraph, argemg);
792void Options::newAGraph(
EGraph *egraph)
798 egraph->sId = proc->agrcount;
799 }
else if (sproc != NULL) {
801 egraph->sId = sproc->agrcount;
806 if ((egraph->nsym)*(egraph->esym) != 0) {
807 sf =
Fraction(egraph->extperm, egraph->nsym*egraph->esym);
811 if (egraph->mgraph->opi) {
812 sproc->resultAGraph(1, sf, 1, sf);
814 sproc->resultAGraph(1, sf, 0, zr);
818 if (out != NULL && values[GRCC_OPT_Outgrf]) {
819 out->outEGraphF(egraph);
821 if (out != NULL && values[GRCC_OPT_Outgrp]) {
822 out->outEGraphP(egraph);
825 (*outag)(egraph, argag);
828 sproc->endAGraph(egraph);
850 if (outgrfp != NULL) {
852 grcc_fprintf(GRCC_Stderr,
"*** file has not been closed : \"%s\"\n",
857 erEnd(
"file has not been closed");
859 if (outgrf != NULL) {
863 if (outgrpp != NULL) {
865 grcc_fprintf(GRCC_Stderr,
"*** file has not been closed : \"%s\"\n",
870 erEnd(
"file has not been closed");
872 if (outgrp != NULL) {
879void Output::setOutgrf(
const char *fname)
881 if (outgrf != NULL && strcmp(outgrf, fname) == 0) {
885 if (outgrf != NULL) {
889 if (fname == NULL || strlen(fname) < 1) {
892 outgrf = strdup(fname);
897void Output::setOutgrp(
const char *fname)
899 if (outgrp != NULL && strcmp(outgrp, fname) == 0) {
903 if (outgrp != NULL) {
906 if (fname == NULL || strlen(fname) < 1) {
909 outgrp = strdup(fname);
914Bool Output::outBeginF(
Model *mdl, Bool pr)
918 char sdate[MAXSTRLEN];
922 if (outgrf == NULL || strlen(outgrf) < 1) {
925 }
else if (outgrfp != NULL) {
927 grcc_fprintf(GRCC_Stderr,
"*** outBegin: \"%s\" is already opened\n",
930 erEnd(
"outBegin: file is already opened\n");
936 if ((outgrfp = fopen(outgrf,
"w")) == NULL) {
938 grcc_fprintf(GRCC_Stderr,
"*** outBegin: cannot open %s\n", outgrf);
944 strftime(sdate, MAXSTRLEN,
"%Y/%m/%d %H:%M:%S", tm);
946 fprintf(outgrfp,
"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
947 fprintf(outgrfp,
"%% Generated by grcc at \"%s\"\n", sdate);
948 fprintf(outgrfp,
"Version={2,2,0,0};\n");
950 fprintf(outgrfp,
"Model=\"./%s.mdl\";\n", model->name);
952 fprintf(outgrfp,
"Model=\"./Phi34.mdl\";\n");
960Bool Output::outBeginP(
Model *mdl, Bool pr)
964 char sdate[MAXSTRLEN];
968 if (outgrp == NULL || strlen(outgrp) < 1) {
971 }
else if (outgrpp != NULL) {
973 grcc_fprintf(GRCC_Stderr,
"*** outBegin: \"%s\" is already opened\n",
976 erEnd(
"outBegin: file is already opened\n");
982 if ((outgrpp = fopen(outgrp,
"w")) == NULL) {
984 grcc_fprintf(GRCC_Stderr,
"*** outBegin: cannot open %s\n", outgrp);
990 strftime(sdate, MAXSTRLEN,
"%Y/%m/%d %H:%M:%S", tm);
992 fprintf(outgrpp,
"################################\n");
993 fprintf(outgrpp,
"# Generated by 'grcc' at \"%s\"\n", sdate);
999void Output::outEndF(
void)
1001 if (outgrfp == NULL) {
1004 fprintf(outgrfp,
"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
1005 fprintf(outgrfp,
"End=1;\n");
1006 fprintf(outgrfp,
"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
1015void Output::outEndP(
void)
1017 if (outgrpp == NULL) {
1020 fprintf(outgrpp,
"################################\n");
1021 fprintf(outgrpp,
"# End\n");
1022 fprintf(outgrpp,
"################################\n");
1031void Output::outProcBeginF(
Process *prc)
1043 if (outgrfp == NULL) {
1051 fprintf(outgrfp,
"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
1052 fprintf(outgrfp,
"Process=%d;\n", proc->id);
1053 fprintf(outgrfp,
"External=%d;\n", proc->ninitl+proc->nfinal);
1054 for (k = 0, ex = 0; k < proc->ninitl; k++, ex++) {
1055 fprintf(outgrfp,
"%4d= initial %s;\n", ex,
1056 proc->model->particleName(proc->initlPart[k]));
1058 for (k = 0; k < proc->nfinal; k++, ex++) {
1059 fprintf(outgrfp,
"%4d= final %s;\n", ex,
1060 proc->model->particleName(proc->finalPart[k]));
1062 fprintf(outgrfp,
"Eend;\n");
1063 for (k = 0; k < proc->model->ncouple; k++) {
1064 fprintf(outgrfp,
"%s=%d; ", proc->model->cnlist[k], proc->clist[k]);
1066 fprintf(outgrfp,
"Loop=%d;\n", proc->loop);
1067 fprintf(outgrfp,
"OPI=%s;\n", (opt->values[GRCC_OPT_1PI] > 0?
"Yes" :
"No"));
1068 fprintf(outgrfp,
"Assign=%s;\n", (opt->values[GRCC_OPT_Step]==GRCC_AGraph ?
"Yes" :
"No"));
1069 fprintf(outgrfp,
"%% Options : dummy values\n");
1070 fprintf(outgrfp,
"Expand=Yes; ExpMin=0; ExpMax=-1;\n");
1071 fprintf(outgrfp,
"Block=No; Tadpole=Yes; Extself=Yes;\n");
1072 fprintf(outgrfp,
"Selfe=Yes; Countert=No; AnyCT=No;Undefp=No;\n");
1076void Output::outProcBeginP(
Process *prc)
1086 if (outgrpp == NULL) {
1094 fprintf(outgrpp,
"# Process=%d;\n", proc->id);
1098void Output::outProcBegin0(
int next,
int couple,
int loop)
1104 if (outgrfp == NULL) {
1112 fprintf(outgrfp,
"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
1113 fprintf(outgrfp,
"Process=%d;\n", procId);
1114 fprintf(outgrfp,
"External=%d;\n", next);
1115 int ninit = (next + 1) / 2;
1116 for (k = 0, ex = 0; k < ninit; k++, ex++) {
1117 fprintf(outgrfp,
"%4d= initial undef;\n", ex);
1119 for (k = 0; k < next - ninit; k++, ex++) {
1120 fprintf(outgrfp,
"%4d= final undef;\n", ex);
1122 fprintf(outgrfp,
"Eend;\n");
1123 fprintf(outgrfp,
"GRCC_PHI=%d; ", couple);
1125 fprintf(outgrfp,
"Loop=%d;\n", loop);
1126 fprintf(outgrfp,
"OPI=%s;\n", (opt->values[GRCC_OPT_1PI] > 0?
"Yes" :
"No"));
1127 fprintf(outgrfp,
"Assign=No;\n");
1128 fprintf(outgrfp,
"%% Options : dummy values\n");
1129 fprintf(outgrfp,
"Expand=Yes; ExpMin=0; ExpMax=-1;\n");
1130 fprintf(outgrfp,
"Block=No; Tadpole=Yes; Extself=Yes;\n");
1131 fprintf(outgrfp,
"Selfe=Yes; Countert=No; AnyCT=No;Undefp=No;\n");
1135void Output::outSProcBeginF(
SProcess *sprc)
1146 pnc = sproc->pnclass;
1148 if (outgrfp == NULL) {
1155 fprintf(outgrfp,
"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
1156 fprintf(outgrfp,
"Process=%d;\n", sproc->id);
1157 fprintf(outgrfp,
"External=%d;\n", sproc->nExtern);
1160 for (j = 0; j < pnc->nclass; j++) {
1161 if (pnc->type[j] == GRCC_AT_Initial) {
1162 for (k = 0; k < pnc->count[j]; k++) {
1163 fprintf(outgrfp,
"%4d= initial %s;\n", ex,
1164 model->particleName(pnc->particle[j]));
1167 }
else if (pnc->type[j] == GRCC_AT_Final) {
1168 for (k = 0; k < pnc->count[j]; k++) {
1169 fprintf(outgrfp,
"%4d= final %s;\n", ex,
1170 model->particleName(-pnc->particle[j]));
1176 fprintf(outgrfp,
"Eend;\n");
1177 for (k = 0; k < sproc->model->ncouple; k++) {
1178 fprintf(outgrfp,
"%s=%d; ", sproc->model->cnlist[k], sproc->clist[k]);
1180 fprintf(outgrfp,
"Loop=%d;\n", sproc->loop);
1181 fprintf(outgrfp,
"OPI=%s;\n", (opt->values[GRCC_OPT_1PI] > 0?
"Yes" :
"No"));
1182 fprintf(outgrfp,
"Assign=%s;\n", (opt->values[GRCC_OPT_Step]==GRCC_AGraph ?
"Yes" :
"No"));
1183 fprintf(outgrfp,
"%% Options : dummy values\n");
1184 fprintf(outgrfp,
"Expand=Yes; ExpMin=0; ExpMax=-1;\n");
1185 fprintf(outgrfp,
"Block=No; Tadpole=Yes; Extself=Yes;\n");
1186 fprintf(outgrfp,
"Selfe=Yes; Countert=No; AnyCT=No;Undefp=No;\n");
1190void Output::outSProcBeginP(
SProcess *sprc)
1199 if (outgrpp == NULL) {
1206 fprintf(outgrpp,
"# SProcess=%d;\n", sproc->id);
1210void Output::outProcEndF(
void)
1212 if (outgrfp == NULL) {
1216 fprintf(outgrfp,
"%%------------------------------\n");
1217 fprintf(outgrfp,
"Pend=%d;\n", proc->id);
1218 }
else if (sproc != NULL) {
1219 fprintf(outgrfp,
"%%------------------------------\n");
1220 fprintf(outgrfp,
"Pend=%d;\n", sproc->id);
1222 fprintf(outgrfp,
"%%------------------------------\n");
1223 fprintf(outgrfp,
"Pend=1;\n");
1229void Output::outProcEndP(
void)
1231 if (outgrpp == NULL) {
1238void Output::outEGraphF(
EGraph *egraph)
1240 int nd, lg, ptcl, ed, intr, j, k, loop;
1241 Model *mdl = egraph->model;
1245 if (outgrfp == NULL) {
1248 fprintf(outgrfp,
"%%------------------------------\n");
1249 if (egraph->assigned) {
1250 fprintf(outgrfp,
"Graph=%ld;\n", egraph->sId);
1251 fprintf(outgrfp,
"%% AGraph=%ld;\n", egraph->aId);
1253 fprintf(outgrfp,
"Graph=%ld;\n", egraph->mId);
1255 fprintf(outgrfp,
"Gtype=%ld;\n", egraph->mId);
1256 if (egraph->assigned) {
1257 fprintf(outgrfp,
"Sfactor=%ld;\n",
1258 egraph->nsym * egraph->esym * egraph->fsign);
1260 fprintf(outgrfp,
"Sfactor=%ld;\n", egraph->nsym * egraph->esym);
1262 fprintf(outgrfp,
"%% ExtPerm=%ld; NSym=%ld; ESym=%ld; NSym1=%ld;"
1264 egraph->extperm, egraph->nsym, egraph->esym, egraph->nsym1,
1267 fprintf(outgrfp,
"Vertex=%d;\n", egraph->nNodes - egraph->nExtern);
1269 for (nd = 0; nd < egraph->nNodes; nd++) {
1270 fprintf(outgrfp,
"%4d", nd);
1271 if (mdl != NULL && egraph->assigned && !egraph->isExternal(nd)) {
1272 intr = egraph->nodes[nd]->intrct;
1273 loop = mdl->interacts[intr]->loop;
1278 fprintf(outgrfp,
"[loop=%d", loop);
1283 if (mdl->ncouple > 1) {
1285 fprintf(outgrfp,
";order={");
1287 fprintf(outgrfp,
"[order={");
1290 for (j = 0; j < mdl->interacts[intr]->nclist; j++) {
1292 fprintf(outgrfp,
",");
1294 fprintf(outgrfp,
"%d", mdl->interacts[intr]->clist[j]);
1296 fprintf(outgrfp,
"}");
1299 fprintf(outgrfp,
"]");
1301 }
else if (!egraph->isExternal(nd)) {
1302 loop = egraph->nodes[nd]->extloop;
1305 fprintf(outgrfp,
"[loop=%d]", loop);
1308 fprintf(outgrfp,
"={");
1311 for (lg = 0; lg < egraph->nodes[nd]->deg; lg++) {
1313 fprintf(outgrfp,
", ");
1315 ed = V2Iedge(egraph->nodes[nd]->edges[lg]);
1317 fprintf(outgrfp,
"%4d", abs(egraph->nodes[nd]->edges[lg]));
1318 ptcl = egraph->edges[ed]->ptcl;
1319 if (mdl != NULL && ptcl != 0 && egraph->assigned) {
1320 if (egraph->nodes[nd]->edges[lg] < 0) {
1321 ptcl = mdl->normalParticle(-ptcl);
1323 ptcl = mdl->normalParticle(ptcl);
1325 fprintf(outgrfp,
"[%s]", mdl->particleName(ptcl));
1327 fprintf(outgrfp,
"[undef]");
1330 fprintf(outgrfp,
"};\n");
1333 if (egraph->nFlines >= 0) {
1334 fprintf(outgrfp,
"%% FLines=%d; FSign=%d; sId=%ld;\n",
1335 egraph->nFlines, egraph->fsign, egraph->sId);
1337 for (j = 0; j < egraph->nFlines; j++) {
1338 fl = egraph->flines[j];
1339 fprintf(outgrfp,
"%% %4d", j);
1340 if (fl->ftype == FL_Open) {
1341 fprintf(outgrfp,
"[FOpen]=[");
1342 }
else if (fl->ftype == FL_Closed) {
1343 fprintf(outgrfp,
"[FLoop]=[");
1345 fprintf(outgrfp,
" ?%d", fl->ftype);
1347 for (k = 0; k < fl->nlist; k++) {
1349 fprintf(outgrfp,
", ");
1351 fprintf(outgrfp,
"%d", fl->elist[k]);
1353 fprintf(outgrfp,
"];\n");
1355 fprintf(outgrfp,
"Vend;\n");
1356 fprintf(outgrfp,
"Gend;\n");
1361void Output::outEGraphP(
EGraph *eg)
1363 if (outgrpp == NULL) {
1368 grcc_fprintf(GRCC_Stdout,
"outEGraphP:egraph:egraph[%ld]\n", eg->mId-1);
1369 eg->printPy(outgrpp, eg->mId);
1372 grcc_fprintf(GRCC_Stdout,
"outEGraphP:mgraph:egraph[%ld]\n", eg->mId-1);
1373 eg->mgraph->printPy(outgrpp, eg->mId);
1378void Output::outEGraphP(
EGraph *eg)
1380 int j, nd, lg, ed, k;
1382 static char undef[] =
"Undef";
1386 if (outgrpp == NULL) {
1391 fprintf(outgrpp,
"egraph[%ld] = {\n", eg->sId-1);
1393 fprintf(outgrpp,
"egraph[%ld] = {\n", eg->mId-1);
1398 fprintf(outgrpp,
" \"Process\": {\n");
1399 proc->outProcP(outgrpp);
1400 fprintf(outgrpp,
" },\n");
1404 fprintf(outgrpp,
" \"GId\": [%ld, %ld, %ld],\n",
1405 eg->mId, eg->aId, eg->sId);
1406 fprintf(outgrpp,
" \"GParam\": {\"Assigned\":%d, "
1407 "\"NNodes\":%d, \"NEdges\":%d, \"NFLines\":%d},\n",
1408 eg->assigned, eg->nNodes, eg->nEdges, eg->nFlines);
1409 fprintf(outgrpp,
" \"Sym\": {\"FSign\":%d, \"NSym\":%ld, "
1410 "\"ESym\":%ld, \"NSym1\":%ld, \"Multp\":%ld},\n",
1411 eg->fsign, eg->nsym, eg->esym, eg->nsym1, eg->multp);
1414 fprintf(outgrpp,
" \"Nodes\": [\n");
1415 for (nd = 0; nd < eg->nNodes; nd++) {
1416 node = eg->nodes[nd];
1419 }
else if(!eg->assigned) {
1420 if (eg->isExternal(nd)) {
1421 ed = Abs(eg->nodes[nd]->edges[0])-1;
1422 s = model->particleName(eg->edges[ed]->ptcl);
1426 }
else if (eg->isExternal(nd)) {
1427 s = model->particleName(
node->intrct);
1429 s = model->interacts[
node->intrct]->name;
1431 fprintf(outgrpp,
" [%d, \"%s\", %d, [",
1432 nd, s,
node->extloop);
1433 for (lg = 0; lg < eg->nodes[nd]->deg; lg++) {
1435 fprintf(outgrpp,
", ");
1437 fprintf(outgrpp,
"%d", eg->nodes[nd]->edges[lg]);
1439 fprintf(outgrpp,
"]],\n");
1441 fprintf(outgrpp,
" ],\n");
1444 fprintf(outgrpp,
" \"Edges\": [\n");
1445 for (ed = 0; ed < eg->nEdges; ed++) {
1446 if (model != NULL) {
1447 s = model->particleName(eg->edges[ed]->ptcl);
1448 fprintf(outgrpp,
" [%d, \"%s\", [%d, %d]],\n",
1449 ed+1, s, eg->edges[ed]->nodes[0],
1450 eg->edges[ed]->nodes[1]);
1452 fprintf(outgrpp,
" [%d, \"Undef\", [%d, %d]],\n",
1453 ed+1, eg->edges[ed]->nodes[0],
1454 eg->edges[ed]->nodes[1]);
1458 fprintf(outgrpp,
" ],\n");
1460 if (eg->nFlines >= 0) {
1462 fprintf(outgrpp,
" \"FLines\": [\n");
1463 for (j = 0; j < eg->nFlines; j++) {
1465 fprintf(outgrpp,
" [%d, ", j);
1466 if (fl->ftype == FL_Open) {
1467 fprintf(outgrpp,
"\"FOpen\", ");
1468 }
else if (fl->ftype == FL_Closed) {
1469 fprintf(outgrpp,
"\"FLoop\", ");
1471 fprintf(outgrpp,
"\"Undef%d\", ", fl->ftype);
1473 fprintf(outgrpp,
"[");
1474 for (k = 0; k < fl->nlist; k++) {
1476 fprintf(outgrpp,
", ");
1478 fprintf(outgrpp,
"%d", fl->elist[k]);
1480 fprintf(outgrpp,
"]],\n");
1482 fprintf(outgrpp,
" ]\n");
1486 fprintf(outgrpp,
"};\n");
1491void Output::outModelF(
void)
1500 if (model == NULL) {
1502 grcc_fprintf(GRCC_Stderr,
"*** Output::outModel : model is not defined\n");
1506 ln = strlen(model->name)+5;
1507 mdlfn =
new char[ln];
1508 snprintf(mdlfn, ln,
"%s.mdl", model->name);
1510 if ((mdlfp = fopen(mdlfn,
"w")) == NULL) {
1512 grcc_fprintf(GRCC_Stderr,
"*** Output::outModel : cannot open \"%s\"\n",
1518 for (l = 0; l < 60; l++) { putc(
'%', mdlfp); } putc(
'\n', mdlfp);
1519 fprintf(mdlfp,
"%% File \"%s\"\n", mdlfn);
1520 fprintf(mdlfp,
"%% Generated by graph generater\n");
1521 fprintf(mdlfp,
" Version={2,2,0};\n");
1524 fprintf(mdlfp,
" Order={");
1525 for (j = 0; j < model->ncouple; j++) {
1527 fprintf(mdlfp,
", ");
1529 fprintf(mdlfp,
"%s", model->cnlist[j]);
1531 fprintf(mdlfp,
"};\n");
1533 for (l = 0; l < 60; l++) { putc(
'%', mdlfp); } putc(
'\n', mdlfp);
1534 fprintf(mdlfp,
"%% particles\n");
1535 for (l = 0; l < 60; l++) { putc(
'%', mdlfp); } putc(
'\n', mdlfp);
1536 for (j = 1; j < model->nParticles; j++) {
1537 pt = model->particles[j];
1538 fprintf(mdlfp,
" Particle=%s; Antiparticle=%s;\n",
1539 pt->name, pt->aname);
1540 ptn = pt->typeGName();
1541 fprintf(mdlfp,
" PType=%s; Charge=0; Color=1; Mass=0; Width=0;\n", ptn);
1542 fprintf(mdlfp,
" PCode=%d; Massless; \n", pt->pcode);
1543 fprintf(mdlfp,
" Pend;\n");
1544 fprintf(mdlfp,
"%%\n");
1546 for (l = 0; l < 60; l++) { putc(
'%', mdlfp); } putc(
'\n', mdlfp);
1547 fprintf(mdlfp,
"%% interactions\n");
1548 for (l = 0; l < 60; l++) { putc(
'%', mdlfp); } putc(
'\n', mdlfp);
1549 for (j = 0; j < model->nInteracts; j++) {
1550 vt = model->interacts[j];
1551 fprintf(mdlfp,
" Vertex={");
1552 for (k = 0; k < vt->nplist; k++) {
1554 fprintf(mdlfp,
", ");
1556 fprintf(mdlfp,
"%s", model->particleName(vt->plist[k]));
1558 fprintf(mdlfp,
"}; ");
1559 for (k = 0; k < model->ncouple; k++) {
1560 fprintf(mdlfp,
"%s=%d; ", model->cnlist[k], vt->clist[k]);
1562 fprintf(mdlfp,
"FName=%s; Vend;\n", vt->name);
1564 for (l = 0; l < 60; l++) { putc(
'%', mdlfp); } putc(
'\n', mdlfp);
1565 fprintf(mdlfp,
" Mend;\n");
1566 for (l = 0; l < 60; l++) { putc(
'%', mdlfp); } putc(
'\n', mdlfp);
1572void Output::outModelP(
void)
1581 if (model == NULL) {
1583 grcc_fprintf(GRCC_Stderr,
"*** Output::outModel : ");
1584 grcc_fprintf(GRCC_Stderr,
"model is not defined\n");
1588 ln = strlen(model->name)+5;
1589 mdlfn =
new char[ln];
1590 snprintf(mdlfn, ln,
"%s.mdp", model->name);
1592 if ((mdlfp = fopen(mdlfn,
"w")) == NULL) {
1594 grcc_fprintf(GRCC_Stderr,
"*** Output::outModel : cannot open \"%s\"\n",
1600 for (l = 0; l < 60; l++) { putc(
'#', mdlfp); } putc(
'\n', mdlfp);
1601 fprintf(mdlfp,
"# File \"%s\"\n", mdlfn);
1602 fprintf(mdlfp,
"# Generated by graph generater grcc\n");
1605 fprintf(mdlfp,
"Model={\n");
1606 fprintf(mdlfp,
" \"Model\": [");
1607 fprintf(mdlfp,
"\"%s\", %d, [",
1608 model->name, model->ncouple);
1609 for (j = 0; j < model->ncouple; j++) {
1611 fprintf(mdlfp,
", ");
1613 fprintf(mdlfp,
"\"%s\"", model->cnlist[j]);
1615 fprintf(mdlfp,
"], ");
1616 if (model->defpart == GRCC_DEFBYCODE) {
1617 fprintf(mdlfp,
"\"ByCode\"");
1619 fprintf(mdlfp,
"\"ByName\"");
1621 fprintf(mdlfp,
"],\n");
1623 fprintf(mdlfp,
" \"Particles\": [\n");
1624 for (j = 1; j < model->nParticles; j++) {
1625 pt = model->particles[j];
1626 ptn = pt->typeName();
1627 fprintf(mdlfp,
" [\"%s\", %d, \"%s\", %d, \"%s\", %d],\n",
1628 pt->name, pt->pcode, pt->aname, pt->acode,
1631 fprintf(mdlfp,
" ],\n");
1632 fprintf(mdlfp,
" \"Interactions\": [\n");
1633 for (j = 0; j < model->nInteracts; j++) {
1634 vt = model->interacts[j];
1635 fprintf(mdlfp,
" [\"%s\", %d, %d, [", vt->name, vt->id, vt->nplist);
1636 for (k = 0; k < vt->nplist; k++) {
1638 fprintf(mdlfp,
", ");
1640 if (model->defpart == GRCC_DEFBYCODE) {
1641 fprintf(mdlfp,
"%d", vt->plist[k]);
1643 fprintf(mdlfp,
"\"%s\"", model->particleName(vt->plist[k]));
1646 fprintf(mdlfp,
"], [");
1647 for (k = 0; k < model->ncouple; k++) {
1649 fprintf(mdlfp,
",");
1651 fprintf(mdlfp,
"%d", vt->clist[k]);
1653 fprintf(mdlfp,
"]],\n");
1655 fprintf(mdlfp,
" ],\n");
1656 fprintf(mdlfp,
"};\n");
1666static const char *ptypenames[] = GRCC_PT_Table;
1667static const char *pgtypenames[] = GRCC_PT_GTable;
1672Particle::Particle(
Model *modl,
int pid,
PInput *pinp)
1674 static char buff[100];
1675 static int nbuff=100;
1680 if (pinp->name == NULL || strlen(pinp->name) < 1) {
1681 if (pinp->pcode != 0) {
1682 snprintf(buff, nbuff,
"pa%d", pinp->pcode);
1683 name = strdup(buff);
1685 snprintf(buff, nbuff,
"pi%d", pid);
1686 name = strdup(buff);
1689 name = strdup(pinp->name);
1691 if (pinp->aname == NULL) {
1692 if (pinp->acode != 0) {
1693 snprintf(buff, nbuff,
"ap%d", Abs(pinp->acode));
1694 aname = strdup(buff);
1696 snprintf(buff, nbuff,
"ai%d", pid);
1697 aname = strdup(buff);
1700 aname = strdup(pinp->aname);
1702 pcode = pinp->pcode;
1703 acode = pinp->acode;
1704 if (Abs(mdl->defpart) == GRCC_DEFBYCODE) {
1705 neutral = (pcode == acode);
1706 }
else if (Abs(mdl->defpart) == GRCC_DEFBYNAME) {
1707 neutral = ((strcmp(name, aname)==0) ? True : False);
1711 if (Abs(mdl->defpart) == GRCC_DEFBYCODE) {
1712 if (pinp->ptypec < GRCC_PT_Undef || pinp->ptypec > GRCC_PT_Size) {
1714 grcc_fprintf(GRCC_Stderr,
"*** particle type is not defined: %d\n",
1717 erEnd(
"particle type is not defined (illegal code)");
1719 ptype = pinp->ptypec;
1720 }
else if (Abs(mdl->defpart) == GRCC_DEFBYNAME) {
1721 if (pinp->ptypen == NULL) {
1722 erEnd(
"particle type is not defined (NULL)");
1725 for (j = 0; j < GRCC_PT_Size; j++) {
1726 if (strcmp(pinp->ptypen, ptypenames[j]) == 0) {
1733 grcc_fprintf(GRCC_Stderr,
"*** particle type \"%s\" is not defined\n",
1736 erEnd(
"particle type is not defined (name)");
1739 if (ptype == GRCC_PT_Undef && ptype != 0) {
1740 erEnd(
"illegal ptype");
1743 extonly = pinp->extonly;
1749Particle::~Particle(
void)
1756char *Particle::particleName(
int p)
1766int Particle::particleCode(
int p)
1776int Particle::isNeutral(
void)
1782const char *Particle::typeName(
void)
1784 return ptypenames[ptype];
1788const char *Particle::typeGName(
void)
1790 return pgtypenames[ptype];
1794char *Particle::aparticle(
void)
1796 static char prtcl[] =
"Particle";
1806void Particle::prParticle(
void)
1809 grcc_fprintf(GRCC_Stdout,
"pid=%d, name=\"%s\", real_field, ",
id, name);
1811 grcc_fprintf(GRCC_Stdout,
"pid=%d, name=\"%s\", aname=\"%s\", ",
id, name, aname);
1813 grcc_fprintf(GRCC_Stdout,
"ptype=%s, pcode=%d, acode=%d, extonly=%d, cdeg=(%d,%d)\n",
1814 ptypenames[ptype], pcode, acode, extonly, cmindeg, cmaxdeg);
1820Interaction::Interaction(
Model *modl,
int iid,
const char *nam,
int icd,
int *cpl,
int nlgs,
int *plst,
int csm,
int lp)
1822 static char buff[100];
1823 static int nbuff=100;
1824 static Bool prmsg = True;
1826 int j, ndir, sdir, jdir, nmaj, jmaj, ngho, sgho, jgho, ptcl;
1836 if (nam == NULL || strlen(nam) < 1) {
1837 snprintf(buff, nbuff,
"vt%d", icd);
1838 name = strdup(buff);
1842 nclist = mdl->ncouple;
1843 clist =
new int[nclist];
1845 for (j = 0; j < mdl->ncouple; j++) {
1848 plist =
new int[nlegs];
1850 for (j = 0; j < nlegs; j++) {
1866 for (j = 0; j < nplist; j++) {
1868 p = mdl->particles[plist[j]];
1871 p = mdl->particles[-plist[j]];
1875 if (p->ptype == GRCC_PT_Dirac) {
1882 if (Abs(sdir) == 1) {
1884 }
else if ((Abs(sdir) == 0 && j != jdir + 1) || Abs(sdir) > 1) {
1885 grcc_fprintf(GRCC_Stderr,
"*** Interaction: Dirac and anti-Dirac should "
1886 "arranged in pairs in an interaction.\n");
1889 }
else if (p->ptype == GRCC_PT_Ghost) {
1896 if (Abs(sgho) == 1) {
1898 }
else if ((Abs(sgho) == 0 && j != jgho + 1) || Abs(sgho) > 1) {
1899 grcc_fprintf(GRCC_Stderr,
"*** Interaction: Ghost and anti-Ghost should "
1900 "arranged in pairs in an interaction.\n");
1903 }
else if (p->ptype == GRCC_PT_Majorana) {
1905 if (nmaj % 2 == 1) {
1907 }
else if (j != jmaj + 1) {
1908 grcc_fprintf(GRCC_Stderr,
"*** Interaction: two Majoranas should "
1909 "arranged in pairs in an interaction.\n");
1915 grcc_fprintf(GRCC_Stderr,
"*** Interaction: Dirac number is not conserved\n");
1919 grcc_fprintf(GRCC_Stderr,
"*** Interaction: Ghost number is not conserved\n");
1922 if (nmaj % 2 != 0) {
1923 grcc_fprintf(GRCC_Stderr,
"*** Interaction: odd number of Majorana particles\n");
1926 if (ndir + nmaj + ngho > 2 && prmsg) {
1927 grcc_fprintf(GRCC_Stderr,
"+++ Interaction: more than 2 "
1928 "Dirac/Majorana/Ghost "
1929 "particles are interacting.\n");
1930 grcc_fprintf(GRCC_Stderr,
" Interaction has %d Dirac, %d Ghost "
1931 "and %d Majorana particles.\n",
1933 grcc_fprintf(GRCC_Stderr,
" Sign factors related to fermions may "
1934 "be inconsistent with your calculation method.\n");
1939 erEnd(
"illegal Dirac/Majorana/Ghost interaction");
1944Interaction::~Interaction(
void)
1946 if (slist != NULL) {
1949 if (plist != NULL) {
1952 if (clist != NULL) {
1959void Interaction::prInteraction(
void)
1963 grcc_fprintf(GRCC_Stdout,
"vid=%d, icode=%d, name=\"%s\", loop=%d, csum=%d, cpl=",
1964 id, icode, name, loop, csum);
1965 prIntArray(mdl->ncouple, clist,
", legs=[");
1966 for (j = 0; j < nplist; j++) {
1968 grcc_fprintf(GRCC_Stdout,
", ");
1970 grcc_fprintf(GRCC_Stdout,
"%s", mdl->particleName(plist[j]));
1972 grcc_fprintf(GRCC_Stdout,
"];\n");
1979Model::Model(
MInput *minp)
1981 static PInput pundef = {
"undef",
"undef", 0, 0,
"undef", 0, 0};
1984 if (minp->name == NULL || strlen(minp->name) < 1) {
1985 erEnd(
"model should have a name");
1987 name = strdup(minp->name);
1988 ncouple = minp->ncouple;
1989 if (ncouple >= GRCC_MAXNCPLG) {
1990 erEnd(
"too many coupling constants in the model (GRCC_MAXNCPLG)");
1993 cnlist =
new char*[ncouple];
1994 for (j = 0; j < ncouple; j++) {
1995 cnlist[j] = strdup(minp->cnamlist[j]);
1999 particles =
new Particle*[GRCC_MAXMPARTICLES];
2000 for (j = 0; j < GRCC_MAXMPARTICLES; j++) {
2001 particles[j] = NULL;
2007 for (j = 0; j < GRCC_MAXMINTERACT; j++) {
2008 interacts[j] = NULL;
2012 if (particles == NULL || interacts == NULL) {
2013 erEnd(
"memory allocation failed");
2017 if (minp->defpart != GRCC_DEFBYNAME
2018 && minp->defpart != GRCC_DEFBYCODE) {
2019 erEnd(
"illegal value of defpart");
2022 if (Abs(minp->defpart) != GRCC_DEFBYNAME
2023 && Abs(minp->defpart) != GRCC_DEFBYCODE) {
2024 erEnd(
"illegal value of defpart");
2028 defpart = minp->defpart;
2036 addParticle(&pundef);
2045 for (j = ncplgcp-1; j >= 0; j--) {
2046 cplgvl[j] = delintdup(cplgvl[j]);
2049 cplgnvl = delintdup(cplgnvl);
2050 cplglg = delintdup(cplglg);
2051 cplgcp = delintdup(cplgcp);
2054 for (j = GRCC_MAXMINTERACT-1; j >= 0; j--) {
2055 if (interacts[j] != NULL) {
2056 interacts[j]->slist = delintdup(interacts[j]->slist);
2057 delete interacts[j];
2058 interacts[j] = NULL;
2064 for (j = GRCC_MAXMPARTICLES-1; j >= 0; j--) {
2065 if (particles[j] != NULL) {
2066 delete particles[j];
2067 particles[j] = NULL;
2073 for (j=ncouple-1; j>=0; j--) {
2082void Model::prModel(
void)
2084 static char hdr[] =
"#=================================================\n";
2085 static char hd1[] =
"#-------------------------------------------------\n";
2088 grcc_fprintf(GRCC_Stdout,
"%s", hdr);
2089 grcc_fprintf(GRCC_Stdout,
"Model=\"%s\", ", name);
2090 grcc_fprintf(GRCC_Stdout,
"coupling=[");
2091 for (j = 0; j < ncouple; j++) {
2093 grcc_fprintf(GRCC_Stdout,
", ");
2095 grcc_fprintf(GRCC_Stdout,
"\"%s\"", cnlist[j]);
2097 grcc_fprintf(GRCC_Stdout,
"];\n");
2098 grcc_fprintf(GRCC_Stdout,
"%s", hd1);
2099 grcc_fprintf(GRCC_Stdout,
"Particles=%d;\n", nParticles);
2100 for (j = 1; j < nParticles; j++) {
2101 particles[j]->prParticle();
2103 grcc_fprintf(GRCC_Stdout,
"EndParticle;\n");
2104 grcc_fprintf(GRCC_Stdout,
"allPart (%d) = [", nallPart);
2105 for (j = 0; j < nallPart; j++) {
2107 grcc_fprintf(GRCC_Stdout,
", ");
2109 grcc_fprintf(GRCC_Stdout,
"%d", allPart[j]);
2111 grcc_fprintf(GRCC_Stdout,
"]\n");
2113 grcc_fprintf(GRCC_Stdout,
"%s", hd1);
2114 grcc_fprintf(GRCC_Stdout,
"Interactions=%d;\n", nInteracts);
2115 for (j = 0; j < nInteracts; j++) {
2116 interacts[j]->prInteraction();
2118 grcc_fprintf(GRCC_Stdout,
"EndInteraction;\n");
2119 grcc_fprintf(GRCC_Stdout,
"%s", hd1);
2120 grcc_fprintf(GRCC_Stdout,
"InteractionTable;\n");
2121 grcc_fprintf(GRCC_Stdout,
"# %d class in (total coupling, degree)\n", ncplgcp);
2122 for (j = 0; j < ncplgcp; j++) {
2123 grcc_fprintf(GRCC_Stdout,
" class=%d: cp=%d, lg=%d, vl=", j, cplgcp[j], cplglg[j]);
2124 prIntArray(cplgnvl[j], cplgvl[j],
"\n");
2126 grcc_fprintf(GRCC_Stdout,
"EndInteractionTable;\n");
2127 grcc_fprintf(GRCC_Stdout,
"%s", hdr);
2131void Model::addParticle(
PInput *pinp)
2133 int nid, aid, pid, pcd, acd;
2135 if (nParticles >= GRCC_MAXMPARTICLES) {
2136 erEnd(
"particle table overflow (GRCC_MAXMPARTICLES)");
2140 if (Abs(defpart) == GRCC_DEFBYCODE) {
2141 pcd = findParticleCode(pinp->pcode);
2142 acd = findParticleCode(pinp->acode);
2143 if (pcd > 0 || acd > 0) {
2145 grcc_fprintf(GRCC_Stderr,
"*** particle code [%d, %d] ",
2146 pinp->pcode, pinp->acode);
2147 grcc_fprintf(GRCC_Stderr,
"is already used.\n");
2149 erEnd(
"particle code is already defined");
2151 }
else if (Abs(defpart) == GRCC_DEFBYNAME) {
2152 if (pinp->name == NULL || pinp->aname == NULL ||
2153 strlen(pinp->name) < 1 || strlen(pinp->aname) < 1) {
2154 erEnd(
"no name of particle or anti-particle.");
2156 nid = findParticleName(pinp->name);
2157 aid = findParticleName(pinp->aname);
2158 if (nid > 0 || aid > 0) {
2160 grcc_fprintf(GRCC_Stderr,
"*** particle name [%s, %s] ",
2161 pinp->name, pinp->aname);
2162 grcc_fprintf(GRCC_Stderr,
"is already used.\n");
2164 erEnd(
"particle name is already defined.");
2169 particles[pid] =
new Particle(
this, pid, pinp);
2173void Model::addParticleEnd(
void)
2181 for (p = nParticles-1; p >= 1; p--) {
2182 if (!particles[p]->extonly) {
2183 ap = antiParticle(p);
2185 allPart[nallPart++] = ap;
2189 for (p = 1; p < nParticles; p++) {
2190 if (!particles[p]->extonly) {
2191 allPart[nallPart++] = p;
2194 sorti(nallPart, allPart);
2197 for (p = nParticles-1; p >= 1; p--) {
2198 ap = antiParticle(p);
2200 allPart[nallPart++] = ap;
2203 for (p = 1; p < nParticles; p++) {
2204 allPart[nallPart++] = p;
2206 sorti(nallPart, allPart);
2211void Model::addInteraction(
IInput *iinp)
2213 int vid, nlegs, j, k, c;
2215 int plist[GRCC_MAXLEGS];
2218 erEnd(
"call 'addParticleEnd' before 'addInteraction'");
2220 if (nInteracts >= GRCC_MAXMINTERACT) {
2221 erEnd(
"too many interactions in the model (GRCC_MAXMINTERACT)");
2223 nlegs = iinp->nplistn;
2224 if (nlegs >= GRCC_MAXLEGS) {
2225 erEnd(
"too many legs in an interaction (GRCC_MAXLEGS)");
2229 if (defpart == GRCC_DEFBYCODE) {
2230 if (iinp->icode < 0) {
2232 grcc_fprintf(GRCC_Stderr,
"*** vertex code %d should be positive\n",
2235 erEnd(
"vertex code should be positive");
2237 vid = findInteractionCode(iinp->icode);
2240 grcc_fprintf(GRCC_Stderr,
"*** vertex code %d is already used\n",
2243 erEnd(
"vertex code is already used");
2246 }
else if (defpart == GRCC_DEFBYNAME) {
2247 if (iinp->name == NULL || strlen(iinp->name) < 1) {
2248 erEnd(
"no name of vertex\n");
2250 vid = findInteractionName(iinp->name);
2253 grcc_fprintf(GRCC_Stderr,
"*** vertex name %s is already used\n",
2255 erEnd(
"vertex name is already used");
2261 erEnd(
"illegal value of defpart");
2269 for (j = 0; j < nlegs; j++) {
2270 if (Abs(defpart) == GRCC_DEFBYCODE) {
2271 plist[j] = findParticleCode(iinp->plistc[j]);
2272 if (plist[j] == 0) {
2274 grcc_fprintf(GRCC_Stderr,
"*** particle code %d ",
2276 grcc_fprintf(GRCC_Stderr,
"is not defined\n");
2278 erEnd(
"particle is not defined (code)");
2280 }
else if (Abs(defpart) == GRCC_DEFBYNAME) {
2281 plist[j] = findParticleName(iinp->plistn[j]);
2282 if (plist[j] == 0) {
2284 grcc_fprintf(GRCC_Stderr,
"*** particle name %s ",
2286 grcc_fprintf(GRCC_Stderr,
"is not defined\n");
2288 erEnd(
"particle is not defined (name)");
2294 for (j = 0; j < ncouple; j++) {
2295 c = iinp->cvallist[j];
2298 grcc_fprintf(GRCC_Stderr,
"*** illegal value of c-constants \n");
2299 for (k = 0; k < ncouple; k++) {
2300 grcc_fprintf(GRCC_Stderr,
" %s = %d\n",
2301 cnlist[k], iinp->cvallist[k]);
2304 erEnd(
"illegal value of c-constants");
2310 grcc_fprintf(GRCC_Stderr,
"*** illegal total coupling-constants %d\n",
2312 for (k = 0; k < ncouple; k++) {
2313 grcc_fprintf(GRCC_Stderr,
" %s = %d\n",
2314 cnlist[k], iinp->cvallist[k]);
2317 erEnd(
"illegal value of c-constants");
2321 lp2 = csum - nlegs + 2;
2322 if (lp2 % 2 != 0 || lp2 < 0) {
2324 grcc_fprintf(GRCC_Stderr,
"*** illegal coupling const : ");
2325 grcc_fprintf(GRCC_Stderr,
"nlegs - 2 + 2*loop: 2*loop = %d\n", lp2);
2327 erEnd(
"illegal value of c-constants");
2331 maxnlegs = Max(maxnlegs, nlegs);
2332 maxloop = Max(maxloop , lp);
2333 maxcpl = Max(maxcpl , csum);
2335 interacts[vid] =
new Interaction(
this, vid, iinp->name, iinp->icode,
2336 iinp->cvallist, nlegs, plist, csum, lp);
2340void Model::addInteractionEnd(
void)
2342 int lst[GRCC_MAXMINTERACT];
2343 int tcp[GRCC_MAXMINTERACT], tlg[GRCC_MAXMINTERACT];
2344 int tnvl[GRCC_MAXMINTERACT], *tvl[GRCC_MAXMINTERACT];
2346 int cp, lg, j, k, p;
2350 for (cp = 0; cp <= maxcpl; cp++) {
2352 for (lg = 0; lg <= maxnlegs; lg++) {
2354 for (j = 0; j < nInteracts; j++) {
2356 if (vt->csum == cp && vt->nlegs == lg) {
2364 tvl[ncplgcp] = intdup(k, lst);
2369 cplgcp = intdup(ncplgcp, tcp);
2370 cplglg = intdup(ncplgcp, tlg);
2371 cplgnvl = intdup(ncplgcp, tnvl);
2372 cplgvl =
new int*[ncplgcp];
2373 for (j = 0; j < ncplgcp; j++) {
2374 if (cplgnvl[j] > 0) {
2381 for (j = 0; j < nInteracts; j++) {
2383 if (vt->slist != NULL) {
2384 erEnd(
"addInteractionEnd: vt->nslist != NULL");
2386 vt->slist = intdup(vt->nplist, vt->plist);
2387 vt->nslist = toSList(vt->nplist, vt->slist);
2390 for (p = 0; p < nParticles; p++) {
2391 particles[p]->cmindeg = 0;
2392 particles[p]->cmaxdeg = 0;
2394 for (j = 0; j < nInteracts; j++) {
2396 for (k = 0; k < vt->nplist; k++) {
2397 p = abs(vt->plist[k]);
2398 if (particles[p]->cmindeg == 0) {
2399 particles[p]->cmindeg = vt->nlegs;
2400 particles[p]->cmaxdeg = vt->nlegs;
2402 particles[p]->cmindeg = Min(particles[p]->cmindeg, vt->nlegs);
2403 particles[p]->cmaxdeg = Max(particles[p]->cmaxdeg, vt->nlegs);
2410int Model::findParticleName(
const char *name)
2414 for (j = 0; j < nParticles; j++) {
2415 if (strcmp(name, particles[j]->name) == 0) {
2419 for (j = 0; j < nParticles; j++) {
2420 if (strcmp(name, particles[j]->aname) == 0) {
2428int Model::findParticleCode(
int pcd)
2432 for (j = 0; j < nParticles; j++) {
2433 if (pcd == particles[j]->pcode) {
2435 }
else if (pcd == particles[j]->acode) {
2443char *Model::particleName(
int p)
2447 if (Abs(p) >= nParticles) {
2449 grcc_fprintf(GRCC_Stderr,
"\n*** Model::particleName: ");
2450 grcc_fprintf(GRCC_Stderr,
"illegal particle id=%d\n", p);
2452 erEnd(
"Model::particleName: illegal particle");
2459 return particles[q]->particleName(p);
2463int Model::particleCode(
int p)
2467 if (Abs(p) >= nParticles) {
2468 grcc_fprintf(GRCC_Stderr,
"\n*** Model::particleCode: illegal particle id=%d\n",
2470 erEnd(
"Model::particleCode: illegal particle");
2477 return particles[q]->particleCode(p);
2481void Model::prParticleArray(
int n,
int *a,
const char *msg)
2485 grcc_fprintf(GRCC_Stdout,
"[");
2486 for (j = 0; j < n; j++) {
2488 grcc_fprintf(GRCC_Stdout,
", ");
2490 grcc_fprintf(GRCC_Stdout,
"%s", particleName(a[j]));
2492 grcc_fprintf(GRCC_Stdout,
"]%s", msg);
2496int Model::findMClass(
const int cpl,
const int dgr)
2500 for (j = 0; j < ncplgcp; j++) {
2501 if (cplgcp[j] == cpl && cplglg[j] == dgr) {
2509int Model::normalParticle(
int pt)
2514 if (ptc >= nParticles) {
2518 if (p->isNeutral()) {
2526int Model::antiParticle(
int pt)
2532 if (p->isNeutral()) {
2540int *Model::allParticles(
int *len)
2547int Model::findInteractionName(
const char *name)
2551 for (j = 0; j < nInteracts; j++) {
2552 if (strcmp(name, interacts[j]->name) == 0) {
2560int Model::findInteractionCode(
int icd)
2564 for (j = 0; j < nInteracts; j++) {
2565 if (icd == interacts[j]->icode) {
2573void Model::printMInput(
MInput *min)
2575 grcc_fprintf(GRCC_Stdout,
" \"Model\": [\"%s\", %d, [",
2576 min->name, min->ncouple);
2577 for (
int j = 0; j < min->ncouple; j++) {
2579 grcc_fprintf(GRCC_Stdout,
", ");
2581 grcc_fprintf(GRCC_Stdout,
"\"%s\"", min->cnamlist[j]);
2583 grcc_fprintf(GRCC_Stdout,
"], ");
2584 if (min->defpart == GRCC_DEFBYCODE) {
2585 grcc_fprintf(GRCC_Stdout,
"\"ByCode\"");
2587 grcc_fprintf(GRCC_Stdout,
"\"ByName\"");
2589 grcc_fprintf(GRCC_Stdout,
"],\n");
2593void Model::printPInput(
PInput *pin)
2595 grcc_fprintf(GRCC_Stdout,
" [\"%s\"(%d), \"%s\"(%d), \"%s\"(%d), %d],\n",
2596 pin->name, pin->pcode, pin->aname, pin->acode,
2597 pin->ptypen, pin->ptypec, pin->extonly);
2601void Model::printIInput(
IInput *iin)
2605 grcc_fprintf(GRCC_Stdout,
" [\"%s\"(%d), %d, [", iin->name, iin->icode, iin->nplistn);
2606 for (j = 0; j < iin->nplistn; j++) {
2608 grcc_fprintf(GRCC_Stdout,
", ");
2610 grcc_fprintf(GRCC_Stdout,
"\"%s\"(%d)", iin->plistn[j], iin->plistc[j]);
2612 grcc_fprintf(GRCC_Stdout,
"], [");
2613 for (j = 0; j < GRCC_MAXNCPLG; j++) {
2615 grcc_fprintf(GRCC_Stdout,
", ");
2617 grcc_fprintf(GRCC_Stdout,
"%d", iin->cvallist[j]);
2619 grcc_fprintf(GRCC_Stdout,
"],\n");
2628PNodeClass::PNodeClass(
SProcess *spc,
int nnods,
int nclss,
NCInput *cls)
2641 deg =
new int[nclass];
2642 type =
new int[nclass];
2643 particle =
new int[nclass];
2644 count =
new int[nclass];
2645 couple =
new int[nclass];
2646 cmindeg =
new int[nclass];
2647 cmaxdeg =
new int[nclass];
2648 for (j = 0; j < nclass; j++) {
2649 deg[j] = cls[j].cldeg;
2650 count[j] = cls[j].clnum;
2651 type[j] = cls[j].cltyp;
2652 cmindeg[j] = cls[j].cmind;
2653 cmaxdeg[j] = cls[j].cmaxd;
2654 particle[j] = cls[j].ptcl;
2655 couple[j] = cls[j].cple;
2656 lp2 = (couple[j]-deg[j]+2);
2657 if (!isATExternal(type[j]) && (lp2 % 2 != 0 || lp2 < 0)) {
2659 grcc_fprintf(GRCC_Stderr,
"*** PNodeClass: illegal loop: "
2660 ": 2*loop=cpl[%d](%d)-deg[%d](%d)+2 =2*loop = %d\n",
2661 j, couple[j], j, deg[j], lp2);
2662 for (k = 0; k < nclss; k++) {
2663 grcc_fprintf(GRCC_Stderr,
"k=%d, deg=%d, typ=%d, ptcl=%d, ",
2664 k, deg[k], type[k], particle[k]);
2665 grcc_fprintf(GRCC_Stderr,
"cpl=%d, cnt=%d\n",
2666 couple[k], count[k]);
2673 erEnd(
"PNodeClass: illegal loop");
2675 cl2nd =
new int[nclass+1];
2676 nd2cl =
new int[nnodes];
2677 cl2mcl =
new int[nnodes];
2680 for (j = 0; j < nclass; j++) {
2682 for (k = 0; k < count[j]; k++, nn++) {
2685 if (isATExternal(type[j])) {
2688 cl2mcl[j] = spc->model->findMClass(couple[j], deg[j]);
2689 if (cl2mcl[j] < 0) {
2691 grcc_fprintf(GRCC_Stderr,
"*** PNodeClass : no vertex : ");
2692 grcc_fprintf(GRCC_Stderr,
"coupling=%d, degree=%d\n",
2695 erEnd(
"PNodeClass : no vertex");
2699 cl2nd[nclass] = nnodes;
2702PNodeClass::PNodeClass(
SProcess *spc,
int nnods,
int nclss,
int *dgs,
int *typ,
int *ptcl,
int *cpl,
int *cnt,
int *cmind,
int *cmaxd)
2719 deg =
new int[nclass];
2720 type =
new int[nclass];
2721 particle =
new int[nclass];
2722 count =
new int[nclass];
2723 couple =
new int[nclass];
2724 cmindeg =
new int[nclass];
2725 cmaxdeg =
new int[nclass];
2726 for (j = 0; j < nclass; j++) {
2729 particle[j] = ptcl[j];
2732 cmindeg[j] = cmind[j];
2733 cmaxdeg[j] = cmaxd[j];
2734 lp2 = (cpl[j]-dgs[j]+2);
2735 if (!isATExternal(type[j]) && (lp2 % 2 != 0 || lp2 < 0)) {
2737 grcc_fprintf(GRCC_Stderr,
"*** PNodeClass: illegal loop: "
2738 ": 2*loop=cpl[%d](%d)-deg[%d](%d)+2 =2*loop = %d\n",
2739 j, cpl[j], j, dgs[j], lp2);
2740 for (k = 0; k < nclss; k++) {
2741 grcc_fprintf(GRCC_Stderr,
"k=%d, dgs=%d, typ=%d, ptcl=%d, ",
2742 k, dgs[k], typ[k], ptcl[k]);
2743 grcc_fprintf(GRCC_Stderr,
"cpl=%d, cnt=%d\n", cpl[k], cnt[k]);
2750 erEnd(
"PNodeClass: illegal loop");
2752 cl2nd =
new int[nclass+1];
2753 nd2cl =
new int[nnodes];
2754 cl2mcl =
new int[nnodes];
2757 for (j = 0; j < nclass; j++) {
2759 for (k = 0; k < count[j]; k++, nn++) {
2762 if (isATExternal(type[j])) {
2765 cl2mcl[j] = spc->model->findMClass(couple[j], deg[j]);
2766 if (cl2mcl[j] < 0) {
2768 grcc_fprintf(GRCC_Stderr,
"*** PNodeClass : no vertex : ");
2769 grcc_fprintf(GRCC_Stderr,
"coupling=%d, degree=%d\n",
2772 erEnd(
"PNodeClass : no vertex");
2776 cl2nd[nclass] = nnodes;
2780PNodeClass::~PNodeClass(
void)
2795void PNodeClass::prPNodeClass(
void)
2799 grcc_fprintf(GRCC_Stdout,
"+++ PNodeClass: nclass=%d, nnodes=%d\n", nclass, nnodes);
2800 for (j = 0; j < nclass; j++) {
2803 grcc_fprintf(GRCC_Stdout,
"\n");
2807void PNodeClass::prElem(
int j)
2811 grcc_fprintf(GRCC_Stdout,
"%3d: %-7s(%2d), ", j, GRCC_AT_NdStr(type[j]), type[j]);
2812 grcc_fprintf(GRCC_Stdout,
"deg=%d, count=%d, nodes[%d--%d], couple=%d, cmindeg=%d, cmaxdeg=%d",
2813 deg[j], count[j], cl2nd[j], cl2nd[j+1], couple[j], cmindeg[j], cmaxdeg[j]);
2815 if (isATExternal(tp)) {
2816 if (sproc->model != NULL) {
2817 grcc_fprintf(GRCC_Stdout,
", ptcl=%s ", sproc->model->particleName(particle[j]));
2819 grcc_fprintf(GRCC_Stdout,
", ptcl=%d ", particle[j]);
2822 grcc_fprintf(GRCC_Stdout,
"\n");
2829SProcess::SProcess(
Model *mdl,
Process *prc,
Options *opts,
int sid,
int *clst,
int ncls,
int *cdeg,
int *ctyp,
int *ptcl,
int *cpl,
int *cnum,
int *cmind,
int *cmaxd)
2845 int j, cp, lp2, ndeg, nvrt, tcpl0;
2849 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: no class %d\n", ncls);
2850 erEnd(
"SProcess::SProcess: no Class");
2874 wMGraphs.setValue(0,1);
2875 wMOPI.setValue(0,1);
2879 wAGraphs.setValue(0,1);
2880 wAOPI.setValue(0,1);
2891 for (j = 0; j < model->ncouple; j++) {
2895 for (j = 0; j < nclass; j++) {
2898 lp2 = cpl[j] - cdeg[j] + 2;
2899 if (lp2 % 2 != 0 || lp2 < 0) {
2901 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: illegal loop:"
2902 "class %d, cp=%d, deg=%d, lp2=%d\n",
2903 j, cp, cdeg[j], lp2);
2908 tCouple += cp*cnum[j];
2909 ndeg += cdeg[j]*cnum[j];
2911 if (!isATExternal(ctyp[j])) {
2913 }
else if (isATExternal(ctyp[j])) {
2914 if (model->normalParticle(ptcl[j]) == 0) {
2916 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: ");
2917 grcc_fprintf(GRCC_Stderr,
"illegal particle code: ");
2918 grcc_fprintf(GRCC_Stderr,
"code: class %d, ptcl=%d\n", j, ptcl[j]);
2924 extperm *= factorial(cnum[j]);
2928 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: illegal type: ");
2929 grcc_fprintf(GRCC_Stderr,
"class %d, type=%d\n", j, ctyp[j]);
2934 if (tcpl0 != tCouple) {
2936 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: ");
2937 grcc_fprintf(GRCC_Stderr,
"illegal coupling constants:");
2938 grcc_fprintf(GRCC_Stderr,
" %d != %d\n", tcpl0, tCouple);
2942 if (ndeg == nExtern) {
2944 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: ");
2945 grcc_fprintf(GRCC_Stderr,
"no vertices: ");
2946 grcc_fprintf(GRCC_Stderr,
"nExtern=%d, ndeg=%d\n", nExtern, ndeg);
2950 if (ndeg % 2 == 0) {
2954 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: ");
2955 grcc_fprintf(GRCC_Stderr,
"illegal total deg = %d (not even)\n", ndeg);
2959 nNodes = nvrt + nExtern;
2960 if (nNodes >= GRCC_MAXNODES) {
2962 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: ");
2963 grcc_fprintf(GRCC_Stderr,
"too many nodes = %d\n", nNodes);
2964 grcc_fprintf(GRCC_Stderr,
" nExtern=%d, nvert=%d (GRCC_MAXNODES)\n",
2971 erEnd(
"SProcess::SProcess: illegal input");
2974 opt->beginProc(NULL);
2977 lp2 = tCouple - nExtern + 2;
2978 if (lp2 % 2 != 0 || lp2 < 0) {
2980 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: illegal loop: "
2981 "tCouple=%d, nExtern=%d, 2*loop=%d\n",
2982 tCouple, nExtern, lp2);
2989 if (opt->values[GRCC_OPT_Step] == GRCC_AGraph) {
2990 astack =
new AStack(GRCC_MAXNSTACK, GRCC_MAXESTACK);
2994 pnclass =
new PNodeClass(
this, nNodes, nclass, cdeg, ctyp, ptcl, cpl, cnum, cmind, cmaxd);
3008 int j, cp, lp2, ndeg, nvrt, tcpl0;
3012 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: no class %d\n", ncls);
3013 erEnd(
"SProcess::SProcess: no Class");
3037 wMGraphs.setValue(0,1);
3038 wMOPI.setValue(0,1);
3042 wAGraphs.setValue(0,1);
3043 wAOPI.setValue(0,1);
3054 for (j = 0; j < model->ncouple; j++) {
3058 for (j = 0; j < nclass; j++) {
3061 lp2 = cls[j].cple - cls[j].cldeg + 2;
3062 if (lp2 % 2 != 0 || lp2 < 0) {
3064 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: illegal loop: "
3065 "class %d, cp=%d, deg=%d, lp2=%d\n",
3066 j, cp, cls[j].cldeg, lp2);
3071 tCouple += cp*cls[j].clnum;
3072 ndeg += cls[j].cldeg*cls[j].clnum;
3074 if (!isATExternal(cls[j].cltyp)) {
3075 nvrt += cls[j].clnum;
3076 }
else if (isATExternal(cls[j].cltyp)) {
3077 if (model->normalParticle(cls[j].ptcl) == 0) {
3079 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: ");
3080 grcc_fprintf(GRCC_Stderr,
"illegal particle code: ");
3081 grcc_fprintf(GRCC_Stderr,
"code: class %d, ptcl=%d\n", j, cls[j].ptcl);
3085 nExtern += cls[j].clnum;
3087 extperm *= factorial(cls[j].clnum);
3091 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: illegal type: ");
3092 grcc_fprintf(GRCC_Stderr,
"class %d, type=%d\n", j, cls[j].cltyp);
3097 if (tcpl0 != tCouple) {
3099 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: ");
3100 grcc_fprintf(GRCC_Stderr,
"illegal coupling constants:");
3101 grcc_fprintf(GRCC_Stderr,
" %d != %d\n", tcpl0, tCouple);
3105 if (ndeg == nExtern) {
3107 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: ");
3108 grcc_fprintf(GRCC_Stderr,
"no vertices: ");
3109 grcc_fprintf(GRCC_Stderr,
"nExtern=%d, ndeg=%d\n", nExtern, ndeg);
3113 if (ndeg % 2 == 0) {
3117 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: ");
3118 grcc_fprintf(GRCC_Stderr,
"illegal total deg = %d (not even)\n", ndeg);
3122 nNodes = nvrt + nExtern;
3123 if (nNodes >= GRCC_MAXNODES) {
3125 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: ");
3126 grcc_fprintf(GRCC_Stderr,
"too many nodes = %d\n", nNodes);
3127 grcc_fprintf(GRCC_Stderr,
" nExtern=%d, nvert=%d (GRCC_MAXNODES)\n",
3133 erEnd(
"SProcess::SProcess: illegal input");
3136 opt->beginProc(NULL);
3139 lp2 = tCouple - nExtern + 2;
3140 if (lp2 % 2 != 0 || lp2 < 0) {
3142 grcc_fprintf(GRCC_Stderr,
"*** SProcess::SProcess: illegal loop: "
3143 "tCouple=%d, nExtern=%d, 2*loop=%d\n",
3144 tCouple, nExtern, lp2);
3151 if (opt->values[GRCC_OPT_Step] == GRCC_AGraph) {
3152 astack =
new AStack(GRCC_MAXNSTACK, GRCC_MAXESTACK);
3156 pnclass =
new PNodeClass(
this, nNodes, nclass, cls);
3160SProcess::~SProcess(
void)
3169void SProcess::prSProcess(
void)
3171 grcc_fprintf(GRCC_Stdout,
"\n");
3172 grcc_fprintf(GRCC_Stdout,
"+++ Subprocess %d, class=%d\n",
id, nclass);
3173 if (pnclass != NULL) {
3174 pnclass->prPNodeClass();
3176 grcc_fprintf(GRCC_Stdout,
" pnclass = NULL\n");
3181BigInt SProcess::generate(
void)
3184 int cldeg[GRCC_MAXNODES], clnum[GRCC_MAXNODES], cltyp[GRCC_MAXNODES];
3185 int cmind[GRCC_MAXNODES], cmaxd[GRCC_MAXNODES];
3189 ncl = toMNodeClass(cltyp, cldeg, clnum, cmind, cmaxd);
3191 opt->beginSubProc(
this);
3193 mgraph =
new MGraph(
id, ncl, cldeg, clnum, cltyp, cmind, cmaxd, opt);
3195 ng = mgraph->generate();
3206int SProcess::toMNodeClass(
int *cltyp,
int *cldeg,
int *clnum,
int *cmind,
int *cmaxd)
3210 for (j = 0; j < nclass; j++) {
3211 if (isATExternal(pnclass->type[j])) {
3214 cltyp[j] = (pnclass->couple[j]-pnclass->deg[j]+2)/2;
3216 cldeg[j] = pnclass->deg[j];
3217 clnum[j] = pnclass->count[j];
3218 cmind[j] = pnclass->cmindeg[j];
3219 cmaxd[j] = pnclass->cmaxdeg[j];
3225void SProcess::assign(
MGraph *mgr)
3231 agraph =
new Assign(
this, mgr, pnc);
3234 agraph->checkAG(
"SProcess::assign");
3246 int typ[GRCC_MAXNODES], dgs[GRCC_MAXNODES];
3247 int cnt[GRCC_MAXNODES], ptcl[GRCC_MAXNODES];
3248 int cpl[GRCC_MAXNODES];
3249 int n2m[GRCC_MAXNODES];
3250 int cmind[GRCC_MAXNODES], cmaxd[GRCC_MAXNODES];
3251 int j, k, l, nd, mc, mcl, mclss;
3253 if (mgr->nNodes != nNodes) {
3255 grcc_fprintf(GRCC_Stderr,
"*** SProcess::match: different nNodes=%d != %d\n",
3256 nNodes, mgr->nNodes);
3258 erEnd(
"SProcess::match: different nNodes");
3261 mcl = toMNodeClass(typ, dgs, cnt, cmind, cmaxd);
3264 for (j = 0; j < mcl; j++) {
3265 for (k = 0; k < cnt[j]; k++, nd++) {
3271 mclss = mgr->curcl->nClasses;
3272 for (j = 0; j < mclss; j++) {
3273 for (k = 0; k < mgr->curcl->clist[j]; k++, nd++) {
3274 mc = mgr->nodes[nd]->clss;
3275 if (mc != n2m[nd]) {
3277 grcc_fprintf(GRCC_Stderr,
"*** SProcess::match: ");
3278 grcc_fprintf(GRCC_Stderr,
"inconsistent class\n");
3279 for (l = 0; l < nNodes; l++) {
3280 grcc_fprintf(GRCC_Stderr,
" %d: %d, %d\n",
3281 j, mgr->nodes[l]->clss, n2m[l]);
3284 erEnd(
"SProcess::match: inconsistent class");
3288 for (j = 0; j < mcl; j++) {
3289 dgs[j] = pnclass->deg[j];
3290 typ[j] = pnclass->type[j];
3291 ptcl[j] = pnclass->particle[j];
3292 cnt[j] = pnclass->count[j];
3293 cpl[j] = pnclass->couple[j];
3294 cmind[j] = pnclass->cmindeg[j];
3295 cmaxd[j] = pnclass->cmaxdeg[j];
3297 pnc =
new PNodeClass(
this, nNodes, mcl, dgs, typ, ptcl, cpl, cnt, cmind, cmaxd);
3303void SProcess::resultMGraph(BigInt nmgraphs,
Fraction mwsum, BigInt nmopi,
Fraction mwopi)
3305 nMGraphs += nmgraphs;
3307 wMGraphs.add(mwsum);
3312void SProcess::resultAGraph(BigInt nmgraphs,
Fraction mwsum, BigInt nmopi,
Fraction mwopi)
3314 nAGraphs += nmgraphs;
3316 wAGraphs.add(mwsum);
3321void SProcess::endMGraph(
MGraph *mgr)
3323 egraph = mgr->egraph;
3325 if (opt->values[GRCC_OPT_Step] != GRCC_MGraph) {
3331void SProcess::endAGraph(
EGraph *egr)
3340Process::Process(
int pid,
Model *modl,
Options *optn,
int nini,
int *intlPrt,
int nfin,
int *finlPrt,
int *coupling)
3358 initlPart = intdup(ninitl, intlPrt);
3359 finalPart = intdup(nfinal, finlPrt);
3363 for (j = 0; j < GRCC_MAXNCPLG; j++) {
3364 if (j < model->ncouple) {
3365 clist[j] = coupling[j];
3366 ctotal += coupling[j];
3372 erEnd(
"Process: coupling constant = 0");
3388 wMGraphs.setValue(0,1);
3389 wMOPI.setValue(0,1);
3393 wAGraphs.setValue(0,1);
3394 wAOPI.setValue(0,1);
3399 nExtern = ninitl + nfinal;
3400 maxnlegs = Max(nExtern, model->maxnlegs);
3403 lp2 = ctotal - nExtern + 2;
3406 grcc_fprintf(GRCC_Stderr,
"*** cannot generate : 2*loop is odd : "
3407 "2*loop = %d, ctotal=%d, nExtern=%d\n",
3408 lp2, ctotal, nExtern);
3411 }
else if (lp2 < 0) {
3413 grcc_fprintf(GRCC_Stderr,
"*** cannot make a connected graph : "
3414 "2*loop=%d, ctotal=%d, nExtern=%d\n",
3415 lp2, ctotal, nExtern);
3424 grcc_fprintf(GRCC_Stderr,
"*** Process: illegal input: lp2 = %d\n", lp2);
3426 erEnd(
"Process: illegal input");
3429 if (opt->out->outgrp != NULL) {
3430 snprintf(buff, MAXSTR,
"out%d.prp", pid);
3439Process::~Process(
void)
3441 initlPart = delintdup(initlPart);
3442 finalPart = delintdup(finalPart);
3443 for (
int j = 0; j < nSubproc; j++) {
3450void Process::prProcess(
void)
3452 grcc_fprintf(GRCC_Stdout,
"+++ process options : OPI = %d, (Step) = %d, coupling=%d: ",
3453 opt->values[GRCC_OPT_1PI], opt->values[GRCC_OPT_Step], ctotal);
3454 prIntArray(model->ncouple, clist,
"\n");
3458void Process::prProcessP(
const char *fname)
3462 if ((fp = fopen(fname,
"w")) == NULL) {
3463 grcc_fprintf(GRCC_Stderr,
"*** cannot open \"%s\"\n", fname);
3466 fprintf(fp,
"Process = {\n");
3468 fprintf(fp,
"};\n");
3473void Process::outProcP(FILE *fp)
3477 if (model == NULL) {
3480 fprintf(fp,
" \"Model\": \"%s\",\n", model->name);
3481 fprintf(fp,
" \"Initial\": [");
3482 for (j = 0; j < ninitl; j++) {
3486 fprintf(fp,
"\"%s\"", model->particleName(initlPart[j]));
3488 fprintf(fp,
"],\n");
3489 fprintf(fp,
" \"Final\": [");
3490 for (j = 0; j < nfinal; j++) {
3494 fprintf(fp,
"\"%s\"", model->particleName(finalPart[j]));
3496 fprintf(fp,
"],\n");
3497 fprintf(fp,
" \"Couple\": [");
3498 for (j = 0; j < model->ncouple; j++) {
3502 fprintf(fp,
"%d", clist[j]);
3504 fprintf(fp,
"],\n");
3505 fprintf(fp,
" \"Options\": {\n");
3506 for (j = 0; j < GRCC_OPT_Size; j++) {
3507 if (j == GRCC_OPT_Outgrf || j == GRCC_OPT_Outgrp) {
3510 fprintf(fp,
" \"%s\":%d,\n", optDef[j].name, opt->values[j]);
3512 fprintf(fp,
" \"%s\":", optDef[GRCC_OPT_Outgrf].name);
3513 if (opt->out->outgrf != NULL) {
3514 fprintf(fp,
"\"%s\",\n", opt->out->outgrf);
3516 fprintf(fp,
"0,\n");
3518 fprintf(fp,
" \"%s\":", optDef[GRCC_OPT_Outgrp].name);
3519 if (opt->out->outgrp != NULL) {
3520 fprintf(fp,
"\"%s\",\n", opt->out->outgrp);
3522 fprintf(fp,
"0,\n");
3524 fprintf(fp,
" }\n");
3528void Process::mkSProcess(
void)
3536 NCInput cls[GRCC_MAXMINTERACT];
3537 int r, j, k, lin2, nvtx, nleg, nc, n0;
3538 int nl[GRCC_MAXMINTERACT];
3539 double proct0 = 0, proct1 = 0;
3541 if (model->ncplgcp < 1) {
3542 erEnd(
"function 'addInteractionEnd' has not been called");
3545 opt->beginProc(
this);
3562 while (nextPart(ctotal, model->ncplgcp, model->cplgcp, nl, &r)) {
3566 for (j = 0; j < model->ncplgcp; j++) {
3568 nleg += nl[j]*model->cplglg[j];
3572 lin2 = nExtern + nleg;
3573 if (lin2 % 2 != 0) {
3576 loop = lin2/2 - nExtern - nvtx + 1;
3582 if (opt->values[GRCC_OPT_SymmInitial]) {
3584 for (j = 0; j < ninitl; j++) {
3585 for (k = n0; k < nc; k++) {
3586 if (cls[k].ptcl == initlPart[j]) {
3592 cls[nc].ptcl = initlPart[j];
3597 cls[nc].cltyp = GRCC_AT_Initial;
3598 cls[nc].ptcl = initlPart[j];
3600 = model->particles[Abs(cls[nc].ptcl)]->cmindeg;
3602 = model->particles[Abs(cls[nc].ptcl)]->cmaxdeg;
3603 if (opt->values[GRCC_OPT_NoExtSelf] > 0 && nExtern > 2) {
3604 cls[nc].cmind = Max(cls[nc].cmind, 3);
3610 for (j = 0; j < ninitl; j++) {
3611 cls[nc].ptcl = initlPart[j];
3615 cls[nc].cltyp = GRCC_AT_Initial;
3616 cls[nc].cmind = model->particles[Abs(cls[nc].ptcl)]->cmindeg;
3617 cls[nc].cmaxd = model->particles[Abs(cls[nc].ptcl)]->cmaxdeg;
3618 if (opt->values[GRCC_OPT_NoExtSelf] > 0 && nExtern > 2) {
3619 cls[nc].cmind = Max(cls[nc].cmind, 3);
3624 if (opt->values[GRCC_OPT_SymmFinal]) {
3626 for (j = 0; j < nfinal; j++) {
3627 for (k = n0; k < nc; k++) {
3628 if (cls[k].ptcl == model->antiParticle(finalPart[j])) {
3634 cls[nc].ptcl = model->antiParticle(finalPart[j]);
3638 cls[nc].cltyp = GRCC_AT_Final;
3640 = model->particles[Abs(cls[nc].ptcl)]->cmindeg;
3642 = model->particles[Abs(cls[nc].ptcl)]->cmaxdeg;
3643 if (opt->values[GRCC_OPT_NoExtSelf] > 0 && nExtern > 2) {
3644 cls[nc].cmind = Max(cls[nc].cmind, 3);
3650 for (j = 0; j < nfinal; j++) {
3651 cls[nc].ptcl = model->antiParticle(finalPart[j]);
3655 cls[nc].cltyp = GRCC_AT_Final;
3656 cls[nc].cmind = model->particles[Abs(cls[nc].ptcl)]->cmindeg;
3657 cls[nc].cmaxd = model->particles[Abs(cls[nc].ptcl)]->cmaxdeg;
3658 if (opt->values[GRCC_OPT_NoExtSelf] > 0 && nExtern > 2) {
3659 cls[nc].cmind = Max(cls[nc].cmind, 3);
3664 for (j = 0; j < model->ncplgcp; j++) {
3667 cls[nc].cple = model->cplgcp[j];
3668 cls[nc].cldeg = model->cplglg[j];
3669 cls[nc].clnum = nl[j];
3671 cls[nc].cltyp = (model->cplgcp[j] - model->cplglg[j] + 2)/2;
3678 sproc =
new SProcess(model,
this, opt, nSubproc, clist, nc, cls);
3681 if (nSubproc >= GRCC_MAXSUBPROCS) {
3682 erEnd(
"Subclass: too many sprocesses (GRCC_MAXSUBPROCS)");
3684 sptbl[nSubproc++] = sproc;
3690 nMGraphs += sproc->nMGraphs;
3691 nMOPI += sproc->nMOPI;
3692 wMGraphs.add(sproc->wMGraphs);
3693 wMOPI.add(sproc->wMOPI);
3695 nAGraphs += sproc->nAGraphs;
3696 nAOPI += sproc->nAOPI;
3697 wAGraphs.add(sproc->wAGraphs);
3698 wAOPI.add(sproc->wAOPI);
3711 sec = proct1-proct0;
3721MNodeClass::MNodeClass(
int nnodes,
int nclasses)
3730 nClasses = nclasses;
3732 for (j = 0; j < nNodes; j++) {
3733 for (k = 0; k < nNodes; k++) {
3750MNodeClass::~MNodeClass(
void)
3757void MNodeClass::init(
int *cl,
int mxdeg,
int **adjmat)
3767 for (j = 0; j < nClasses; j++) {
3787 nClasses = mnc->nClasses;
3788 for (k = 0; k < nClasses; k++) {
3789 clist[k] = mnc->clist[k];
3791 maxdeg = mnc->maxdeg;
3792 for (j = 0; j < nNodes; j++) {
3793 ndcl[j] = mnc->ndcl[j];
3794 clord[j] = mnc->clord[j];
3795 for (k = 0; k < nClasses; k++) {
3796 clmat[j][k] = mnc->clmat[j][k];
3799 for (k = 0; k < nClasses+1; k++) {
3800 flist[k] = mnc->flist[k];
3805void MNodeClass::mkFlist(
void)
3813 for (j = 0; j < nClasses; j++) {
3817 flist[nClasses] = f;
3821void MNodeClass::mkNdCl(
void)
3829 for (c = 0; c < nClasses; c++) {
3830 for (k = 0; k < clist[c]; k++) {
3837int MNodeClass::clCmp(
int nd0,
int nd1,
int cn)
3845 cmp = ndcl[nd0] - ndcl[nd1];
3851 cmp = - cmpMNCArray(clmat[nd0], clmat[nd1], cn);
3859void MNodeClass::mkClMat(
int **adjmat)
3866 for (nd = 0; nd < nNodes; nd++) {
3867 for (td = 0; td < nNodes; td++) {
3870 clmat[nd][tc] += CLWIGHTD(adjmat[nd][td]);
3872 clmat[nd][tc] += CLWIGHTO(adjmat[nd][td]);
3879void MNodeClass::incMat(
int nd,
int td,
int val)
3884 clmat[nd][tdc] += val;
3888void MNodeClass::printMat(
void)
3894 grcc_fprintf(GRCC_Stdout,
"\n");
3896 grcc_fprintf(GRCC_Stdout,
"nClasses=%d", nClasses);
3897 grcc_fprintf(GRCC_Stdout,
" clord="); prIntArray(nClasses, clord,
"");
3898 grcc_fprintf(GRCC_Stdout,
" flist="); prIntArray(nClasses, flist,
"\n");
3899 grcc_fprintf(GRCC_Stdout,
"flg = (%d, %d, %d)\n", flg0, flg1, flg2);
3902 grcc_fprintf(GRCC_Stdout,
"nd: cl: ");
3903 for (j2 = 0; j2 < nClasses; j2++) {
3904 grcc_fprintf(GRCC_Stdout,
"%2d ", j2);
3906 grcc_fprintf(GRCC_Stdout,
"\n");
3909 for (j1 = 0; j1 < nNodes; j1++) {
3910 grcc_fprintf(GRCC_Stdout,
"%2d; %2d: [", j1, ndcl[j1]);
3911 for (j2 = 0; j2 < nClasses; j2++) {
3912 grcc_fprintf(GRCC_Stdout,
" %2d", clmat[j1][j2]);
3914 grcc_fprintf(GRCC_Stdout,
"]\n");
3919int MNodeClass::cmpMNCArray(
int *a0,
int *a1,
int ma)
3923 for (k = 0; k < ma; k++) {
3925 if (a0[j] < a1[j]) {
3927 }
else if (a0[j] > a1[j]) {
3935void MNodeClass::reorder(
MGraph *mg)
3937 int flg[GRCC_MAXNODES];
3944 for (co = 0; co < nClasses; co++) {
3946 if (mg->nodes[cn]->freelg < 1) {
3949 }
else if (mg->nodes[cn]->freelg < mg->nodes[cn]->deg) {
3950 flg[co] = mg->nodes[cn]->deg + maxdeg*clist[co];
3953 flg[co] = maxdeg*(mg->nodes[cn]->deg + maxdeg*clist[co]);
3957 if (f0 < nClasses) {
3958 bsort(nClasses, clord, flg);
3962 flg2 = f0 + f1 + f2;
3968MNode::MNode(
int vid,
int vclss,
NCInput *mgi)
3979 freelg = mgi->cldeg;
3980 extloop = mgi->cltyp;
3981 cmindeg = mgi->cmind;
3982 cmaxdeg = mgi->cmaxd;
3986MNode::MNode(
int vid,
int vdeg,
int vextlp,
int vclss,
int cmin,
int cmax)
4006MGraph::MGraph(
int pid,
int ncl,
int *cldeg,
int *clnum,
int *cltyp,
int *cmind,
int *cmaxd,
Options *opts)
4012 clist =
new int[ncl];
4020 for (j = 0; j < nClasses; j++) {
4021 clist[j] = clnum[j];
4023 ne += cldeg[j]*clnum[j];
4027 mindeg = Min(mindeg, cldeg[j]);
4032 maxdeg = Max(maxdeg, cldeg[j]);
4037 grcc_fprintf(GRCC_Stderr,
"Sum of degrees are not even\n");
4038 for (j = 0; j < nClasses; j++) {
4039 grcc_fprintf(GRCC_Stderr,
"class %2d: %2d %2d %2d\n",
4040 j, cldeg[j], clnum[j], cltyp[j]);
4043 erEnd(
"illegal degrees of nodes");
4047 nodes =
new MNode*[nNodes];
4054 nLoops = nEdges - nNodes + 1;
4056 egraph =
new EGraph(nNodes, nEdges, maxdeg);
4059 for (j = 0; j < nClasses; j++) {
4060 for (k = 0; k < clist[j]; k++, nn++) {
4061 nodes[nn] =
new MNode(nn, cldeg[j], cltyp[j], j, cmind[j], cmaxd[j]);
4062 egraph->setExtLoop(nn, cltyp[j]);
4068 egraph->endSetExtLoop();
4080 clist =
new int[ncl];
4088 for (j = 0; j < nClasses; j++) {
4089 clist[j] = mgi[j].clnum;
4091 ne += mgi[j].cldeg*mgi[j].clnum;
4093 mindeg = mgi[j].cldeg;
4095 mindeg = Min(mindeg, mgi[j].cldeg);
4098 maxdeg = mgi[j].cldeg;
4100 maxdeg = Max(maxdeg, mgi[j].cldeg);
4105 grcc_fprintf(GRCC_Stderr,
"Sum of degrees are not even\n");
4106 for (j = 0; j < nClasses; j++) {
4107 grcc_fprintf(GRCC_Stderr,
"class %2d: %2d %2d %2d\n",
4108 j, mgi[j].cldeg, mgi[j].clnum, mgi[j].cltyp);
4111 erEnd(
"illegal degrees of nodes");
4116 grcc_fprintf(GRCC_Stdout,
"*** nNodes = %d\n", nNodes);
4117 erEnd(
"MGraph::MGrap : nNodes < 0");
4119 nodes =
new MNode*[nNodes];
4126 nLoops = nEdges - nNodes + 1;
4128 egraph =
new EGraph(nNodes, nEdges, maxdeg);
4131 for (j = 0; j < nClasses; j++) {
4132 for (k = 0; k < clist[j]; k++, nn++) {
4133 nodes[nn] =
new MNode(nn, j, mgi+j);
4134 egraph->setExtLoop(nn, mgi[j].cltyp);
4135 if (mgi[j].cltyp < 0) {
4140 egraph->endSetExtLoop();
4146void MGraph::init(
void)
4156 adjMat = newMat(nNodes, nNodes, 0);
4166 mconn =
new MConn(nNodes, nEdges);
4180 modmat = newMat(nNodes, nNodes, 0);
4183 bidef = newArray(nNodes, 0);
4184 bilow = newArray(nNodes, 0);
4185 bicol = newArray(nNodes, 0);
4190MGraph::~MGraph(
void)
4196 bicol = deleteArray(bicol);
4197 bilow = deleteArray(bilow);
4198 bidef = deleteArray(bidef);
4201 modmat = deleteMat(modmat, nNodes);
4204 adjMat = deleteMat(adjMat, nNodes);
4210 for (j = 0; j < nNodes; j++) {
4224 grcc_fprintf(GRCC_Stdout,
" ");
4225 for (j2 = 0; j2 < nNodes; j2++) {
4226 grcc_fprintf(GRCC_Stdout,
" %2d", j2);
4228 grcc_fprintf(GRCC_Stdout,
"\n");
4229 for (j1 = 0; j1 < nNodes; j1++) {
4230 grcc_fprintf(GRCC_Stdout,
"%2d : [", j1);
4231 for (j2 = 0; j2 < nNodes; j2++) {
4232 grcc_fprintf(GRCC_Stdout,
" %2d", adjMat[j1][j2]);
4234 grcc_fprintf(GRCC_Stdout,
"] %2d\n", cl->ndcl[j1]);
4239void MGraph::print(
void)
4243 grcc_fprintf(GRCC_Stdout,
"MGraph: pId=%d, cDiag=%ld, c1PI=%ld\n", pId, cDiag, c1PI);
4244 grcc_fprintf(GRCC_Stdout,
" nNodes=%d, nEdges=%d, nExtern=%d, nLoops=%d "
4245 "mindeg=%d, maxdeg=%d, sym=(%ld, %ld)\n",
4246 nNodes, nEdges, nExtern, nLoops, mindeg, maxdeg, nsym, esym);
4247 grcc_fprintf(GRCC_Stdout,
" Nodes=%d\n", nNodes);
4248 for (j = 0; j < nNodes; j++) {
4249 grcc_fprintf(GRCC_Stdout,
" %2d: id=%d, deg=%d, clss=%d, extloop=%d, ",
4250 j, nodes[j]->
id, nodes[j]->deg, nodes[j]->clss,
4252 grcc_fprintf(GRCC_Stdout,
"mind=%d, maxd=%d, freelg=%d\n",
4253 nodes[j]->cmindeg, nodes[j]->cmaxdeg, nodes[j]->freelg);
4257 grcc_fprintf(GRCC_Stdout,
"\n");
4263void MGraph::printPy(FILE *fp,
long mId)
4265 fprintf(fp,
"#TGraph : (gseq, gid, nodes, amat)\n");
4266 fprintf(fp,
"(%ld, (%d, %d),\n", mId, pId, -1);
4267 fprintf(fp,
" [ # nodes: (cid, deg0, loop, part)\n");
4268 for (
int j = 0; j < nNodes; j++) {
4269 fprintf(fp,
" (%2d,%2d,%2d,%2d), # %2d\n",
4270 nodes[j]->clss, nodes[j]->deg, nodes[j]->extloop, 0, j);
4272 fprintf(fp,
" ],\n");
4273 fprintf(fp,
" [\n");
4275 for (
int j2 = 0; j2 < nNodes; j2++) {
4276 fprintf(fp,
" %4d", j2);
4279 for (
int j1 = 0; j1 < nNodes; j1++) {
4281 for (
int j2 = 0; j2 < nNodes; j2++) {
4282 fprintf(fp,
" %2d, ", adjMat[j1][j2]);
4284 fprintf(fp,
"], # %2d\n", j1);
4286 fprintf(fp,
" ], # nodes\n");
4287 fprintf(fp,
" [ # group of 2 elements\n");
4288 fprintf(fp,
" [],\n");
4289 fprintf(fp,
" ], # group\n");
4294Bool MGraph::isConnected(
void)
4302 for (j = 0; j < nNodes; j++) {
4303 nodes[j]->visited = -1;
4309 for (n = 0; n < nNodes; n++) {
4310 if (nodes[n]->visited >= 0) {
4314 return (nv == nNodes);
4318Bool MGraph::visit(
int nd)
4327 if (nodes[nd]->freelg > 0) {
4330 nodes[nd]->visited = 0;
4331 for (td = 0; td < nNodes; td++) {
4332 if ((adjMat[nd][td] > 0) && (nodes[td]->visited < 0)) {
4351 int j1, j2, cmp, nself;
4357 group->newGroup(nNodes, cl->nClasses, cl->clist);
4360 orbits->initPerm(nNodes);
4364 perm = group->genNext();
4369 for (j1 = 0; j1 < nNodes; j1++) {
4370 if (adjMat[j1][j1] > 0) {
4371 nself = adjMat[j1][j1]/2;
4372 esym *= factorial(nself);
4373 esym *= ipow(2, nself);
4375 for (j2 = j1+1; j2 < nNodes; j2++) {
4376 if (adjMat[j1][j2] > 0) {
4377 esym *= factorial(adjMat[j1][j2]);
4384 permMat(nNodes, perm, adjMat, modmat);
4385 cmp = compMat(nNodes, adjMat, modmat);
4388 }
else if (cmp == 0) {
4390 group->addGroup(perm);
4392 nsym = nsym + ToBigInt(1);
4394 orbits->fromPerm(perm);
4401void MGraph::permMat(
int size,
int *perm,
int **mat0,
int **mat1)
4407 for (j1 = 0; j1 < size; j1++) {
4408 for (j2 = 0; j2 < size; j2++) {
4409 mat1[j1][j2] = mat0[perm[j1]][perm[j2]];
4415int MGraph::compMat(
int size,
int **mat0,
int **mat1)
4421 for (j1 = 0; j1 < size; j1++) {
4422 for (j2 = 0; j2 < size; j2++) {
4423 cmp = mat0[j1][j2] - mat1[j1][j2];
4444 int ucl[GRCC_MAXNODES];
4445 int ccn = cl->nClasses;
4454 while (ccl->nClasses != nucl) {
4457 for (td = 1; td < nNodes; td++) {
4462 cmp = ccl->clCmp(td-1, td, ccn);
4471 }
else if (cmp < 0) {
4483 ucl[nucl++] = nce + 1;
4487 if (nucl < ccl->nClasses) {
4488 erEnd(
"refineClasses : smaller number of classes");
4494 xcl->init(ucl, maxdeg, adjMat);
4495 ccn = xcl->nClasses;
4497 nccl = ccl->nClasses;
4503 if (ccl->nClasses == nccl) {
4515void MGraph::biconnME(
void)
4519 int j, j1, root, vr, next, nart;
4527 mconn->initCEdges(
this);
4529 for (j = 0; j < nNodes; j++) {
4541 for (j = 0; j < nNodes; j++) {
4543 if (isExternal(j)) {
4546 }
else if (vr < 0) {
4557 bisearchME(root, -1, 0, 1, &mopi, &mblk, &momset, &next, &nart);
4561 for (j = 0; j < mconn->nopic; j++) {
4562 if (mconn->opics[j].loop > 0) {
4564 }
else if (mconn->opics[j].ctloop > 0) {
4569 mconn->ne0bridges = 0;
4570 mconn->ne1bridges = 0;
4571 for (j = 0; j < mconn->nbridges; j++) {
4572 if (mconn->bridges[j].next == 0) {
4573 mconn->ne0bridges++;
4575 if (mconn->bridges[j].next == 1) {
4576 mconn->ne1bridges++;
4580 mconn->nselfloops = 0;
4581 for (j = 0; j < nNodes; j++) {
4582 mconn->nselfloops += adjMat[j][j]/2;
4585 mconn->nmultiedges = 0;
4586 for (j = 0; j < nNodes; j++) {
4587 for (j1 = j+1; j1 < nNodes; j1++) {
4588 if (adjMat[j][j1] > 1) {
4589 mconn->nmultiedges++;
4594 mconn->na1blocks = 0;
4595 for (j = 0; j < mconn->nblocks; j++) {
4596 if (mconn->blocks[j].nartps == 1) {
4600 mconn->neblocks = mconn->nblocks + mconn->nctopic;
4604 if (isExternal(root)) {
4605 momset |= MASK(root);
4608 for (j = 0; j < nExtern; j++) {
4613 grcc_fprintf(GRCC_Stdout,
"*** momset = %ld != %ld\n", momset, m);
4616 if (mconn->nbacked != nLoops) {
4617 grcc_fprintf(GRCC_Stdout,
"*** nbacked = %d != %c\n", mconn->nbacked, nLoops);
4624 erEnd(
"biconnME: illegal connection");
4634 mconn->ne0bridges = 0;
4635 mconn->ne1bridges = 0;
4640void MGraph::bisearchME(
int nd,
int pd,
int ned,
int col,
4641 MCOpi *mopi,
MCBlock *mblk, ULong *momset,
int *next,
int *nart)
4670 int td, td0, lp, conn, next1, nart1, nart2, nvisit, ndart;
4671 int opit, opit0, blkt, l;
4674 ULong momset1, momset2, momset3;
4686 newv = (bidef[nd] < 0);
4688 grcc_fprintf(GRCC_Stdout,
"*** nd=%d, pd=%d, bidef[nd]=%d\n", nd, pd, bidef[nd]);
4689 erEnd(
"bisearchME: illegal connection");
4692 bidef[nd] = bicount;
4693 bilow[nd] = bicount;
4697 opit0 = mconn->opistkptr;
4700 mconn->pushNode(nd);
4704 if (isExternal(nd)) {
4714 mconn->addCEdge(nd, pd, *momset);
4721 for (td = 0; td < nNodes; td++) {
4722 conn = adjMat[nd][td];
4734 if (isExternal(td) && bidef[td] < 0) {
4738 mconn->addArtic(nd, 1);
4744 mconn->addCEdge(td, nd, momset1);
4747 mconn->pushNode(td);
4751 }
else if (td == nd) {
4757 mconn->addArtic(nd, lp);
4758 mconn->addBlockSelf(nd, lp);
4761 for (l = 0; l < lp; l++) {
4762 int m = nExtern + mconn->nbacked;
4765 mconn->addCEdge(td, nd, momset1);
4769 }
else if (td == pd) {
4771 bilow[nd] = Min(bilow[nd], bidef[pd]);
4772 mopi->loop += ned - 1;
4773 mblk->loop += ned - 1;
4780 }
else if (bidef[td] >= 0) {
4781 bipart &= (bicol[td] == - col);
4782 bilow[nd] = Min(bilow[nd], bidef[td]);
4783 if (bidef[nd] >= bidef[td]) {
4786 mopi->nedges += conn;
4788 mconn->pushEdge(td, nd);
4789 for (l = 0; l < conn; l++) {
4790 int m = nExtern + mconn->nbacked;
4794 mconn->addCEdge(td, nd, momset1);
4800 for (
int ed = 0; ed < mconn->sedges; ed++) {
4801 if (mconn->cedges[ed].nodes[0] == nd &&
4802 mconn->cedges[ed].nodes[1] == td &&
4803 mconn->cedges[ed].momdir > 0) {
4805 *momset &= ~ (mconn->cedges[ed].momset);
4806 }
else if (mconn->cedges[ed].nodes[1] == nd &&
4807 mconn->cedges[ed].nodes[0] == td &&
4808 mconn->cedges[ed].momdir < 0) {
4810 *momset &= ~ (mconn->cedges[ed].momset);
4814 grcc_fprintf(GRCC_Stdout,
"*** revisit back-ed (%d --> %d) cn=%d != conn=%d\n",
4816 grcc_fprintf(GRCC_Stdout,
" sedges = %d\n", mconn->sedges);
4817 for (
int ed = 0; ed < mconn->sedges; ed++) {
4818 grcc_fprintf(GRCC_Stdout,
" %d : (%d --> %d) : [%2d*%2ld]\n", ed,
4819 mconn->cedges[ed].nodes[0],
4820 mconn->cedges[ed].nodes[1],
4821 mconn->cedges[ed].momdir,
4822 mconn->cedges[ed].momset);
4827 erEnd(
"bisearchME: illegal connection");
4832 }
else if (bidef[td] < 0) {
4835 if (pd < 0 && isExternal(nd)) {
4838 opit = mconn->opistkptr;
4840 blkt = mconn->blkstkptr;
4843 mconn->pushEdge(nd, td);
4846 bisearchME(td, nd, adjMat[td][nd], - col,
4847 &mopi1, &mblk1, &momset1, &next1, &nart1);
4851 for (l = 0; l < conn - 1; l++) {
4852 int m = nExtern + mconn->nbacked;
4855 mconn->addCEdge(nd, td, momset3);
4858 mconn->addCEdge(td, nd, momset1 | momset2);
4861 if (bilow[td] > bidef[nd]) {
4864 int mom0lg = (momset1 == 0) ? 1 : 0;
4866 mopi1.mom0lg += mom0lg;
4868 mconn->addOPIc(&mopi1, opit);
4872 mopi1.mom0lg = mom0lg;
4873 if (pd >= 0 || !isExternal(nd)) {
4875 mconn->addBridge(nd, td, next1, nExtern);
4879 mconn->addArtic(nd, 1);
4882 mconn->addBlock(&mblk1, blkt);
4889 }
else if (bilow[td] == bidef[nd]) {
4891 mopi1.nedges += conn;
4892 if (pd >= 0 || !isExternal(nd)) {
4894 mconn->addArtic(nd, 1);
4896 mconn->addBlock(&mblk1, blkt);
4905 mopi1.nedges += conn;
4908 bilow[nd] = Min(bilow[nd], bilow[td]);
4909 mopi->nlegs += mopi1.nlegs;
4910 mopi->next += mopi1.next;
4911 mopi->loop += mopi1.loop;
4912 mopi->nedges += mopi1.nedges;
4913 mopi->ctloop += mopi1.ctloop;
4914 mopi->mom0lg += mopi1.mom0lg;
4916 mblk->nartps += mblk1.nartps;
4917 mblk->loop += mblk1.loop;
4933 mconn->pushNode(nd);
4935 mopi->ctloop += Max(0, nodes[nd]->extloop);
4940 if (isExternal(nd)) {
4942 mconn->addArtic(td0, 1);
4944 if (mconn->nopic > 0) {
4945 (mconn->opics[mconn->nopic-1].next)++;
4948 mconn->addOPIc(mopi, opit0);
4953 mconn->addArtic(nd, -1);
4954 (mconn->blocks[mconn->nblocks-1].nartps)--;
4961BigInt MGraph::generate(
void)
4970 cl->init(clist, maxdeg, adjMat);
4987 xcl = refineClass(cl);
4994 if (xcl->flg0 >= xcl->nClasses) {
4997 sc = xcl->clord[xcl->flg0];
4998 sn = xcl->flist[sc];
4999 connectNode(xcl->flg0, sn, xcl);
5002 if (xcl != cl && xcl != NULL) {
5008void MGraph::connectNode(
int so,
int ss,
MNodeClass *cl)
5013 if (ss >= cl->flist[sc+1]) {
5018 for (sn = ss; sn < cl->flist[sc+1]; sn++) {
5019 connectLeg(so, sn, so, sn, cl);
5025void MGraph::connectLeg(
int so,
int sn,
int to,
int ts,
MNodeClass *cl)
5031 int sc, tc, tn, maxself, nc2, nc, maxcon, ts1, ncm, to1;
5036 if (sn >= cl->flist[sc+1]) {
5037 erEnd(
"connectLeg : illegal control");
5043 if (nodes[sn]->freelg < 1) {
5045 if (!isConnected()) {
5051 connectNode(so, sn+1, cl);
5057 for (to1 = to; to1 < cl->nClasses; to1++) {
5058 tc = cl->clord[to1];
5062 ts1 = cl->flist[tc];
5065 if (ts1 >= nNodes) {
5068 if ((nodes[sn]->cmindeg > 0 && nodes[ts1]->deg < nodes[sn]->cmindeg)
5069 ||(nodes[sn]->cmaxdeg > 0 && nodes[ts1]->deg > nodes[sn]->cmaxdeg)
5070 | (nodes[ts1]->cmindeg > 0 && nodes[sn]->deg < nodes[ts1]->cmindeg)
5071 ||(nodes[ts1]->cmaxdeg > 0 && nodes[sn]->deg > nodes[ts1]->cmaxdeg)) {
5075 for (tn = ts1; tn < cl->flist[tc+1]; tn++) {
5076 if (sc == tc && sn > tn) {
5079 }
else if (nodes[tn]->freelg < 1) {
5083 }
else if (sn == tn) {
5084 if (opt->values[GRCC_OPT_NoSelfLoop] > 0) {
5090 maxself = Min((nodes[sn]->freelg)/2, (nodes[sn]->deg-1)/2);
5093 maxself = nodes[sn]->freelg/2;
5096 if ((nExtern > 1) && (opt->values[GRCC_OPT_1PI] > 0
5097 || opt->values[GRCC_OPT_NoTadpole] > 0)
5101 maxself = Min((nodes[sn]->deg-2)/2, maxself);
5105 for (nc2 = maxself; nc2 > 0; nc2--) {
5108 adjMat[sn][sn] = nc;
5109 nodes[sn]->freelg -= nc;
5110 cl->incMat(sn, tn, ncm);
5113 connectLeg(so, sn, to1, tn+1, cl);
5116 cl->incMat(tn, sn, - ncm);
5118 nodes[sn]->freelg += nc;
5124 maxcon = Min(nodes[sn]->freelg, nodes[tn]->freelg);
5127 if (nNodes > 2 && nodes[sn]->deg == nodes[tn]->deg) {
5128 maxcon = Min(maxcon, nodes[sn]->deg-1);
5131 if (opt->values[GRCC_OPT_NoMultiEdge] > 0) {
5132 maxcon = Min(maxcon, 1);
5136 if ((adjMat[sn][tn] != 0) ||
5137 (adjMat[sn][tn] != adjMat[tn][sn])) {
5138 grcc_fprintf(GRCC_Stdout,
"*** inconsistent connection: sn=%d, tn=%d",
5141 erEnd(
"inconsistent connection ");
5146 for (nc = maxcon; nc > 0; nc--) {
5148 adjMat[sn][tn] = nc;
5149 adjMat[tn][sn] = nc;
5150 nodes[sn]->freelg -= nc;
5151 nodes[tn]->freelg -= nc;
5152 cl->incMat(sn, tn, ncm);
5153 cl->incMat(tn, sn, ncm);
5156 connectLeg(so, sn, to1, tn+1, cl);
5159 cl->incMat(sn, tn, - ncm);
5160 cl->incMat(tn, sn, - ncm);
5163 nodes[sn]->freelg += nc;
5164 nodes[tn]->freelg += nc;
5172Bool MGraph::isOptM(
void)
5176 opi = (mconn->nopic == 1);
5177 opiloop = (mconn->nlpopic <= 1);
5178 tadpole = (mconn->ne0bridges >= ((nExtern == 1) ? 0 : 1));
5179 selfloop = (mconn->nselfloops > 0);
5180 multiedge = (mconn->nmultiedges > 0);
5181 tadblock = (mconn->na1blocks > 0);
5182 block = (mconn->neblocks <= 1);
5183 extself = (mconn->ne1bridges > 0);
5185 if (opt->values[GRCC_OPT_1PI] > 0) {
5187 }
else if (opt->values[GRCC_OPT_1PI] < 0) {
5190 if (opt->values[GRCC_OPT_NoExtSelf] > 0) {
5191 ok = ok && !extself;
5192 }
else if (opt->values[GRCC_OPT_NoExtSelf] < 0) {
5195 if (opt->values[GRCC_OPT_NoTadpole] > 0) {
5197 ok = ok && !tadpole;
5200 ok = ok && !tadpole;
5203 }
else if (opt->values[GRCC_OPT_NoTadpole] < 0) {
5205 ok = ok && !tadpole;
5212 if (opt->values[GRCC_OPT_NoSelfLoop] > 0) {
5213 }
else if (opt->values[GRCC_OPT_NoSelfLoop] < 0) {
5214 ok = ok && selfloop;
5216 if (opt->values[GRCC_OPT_NoMultiEdge] > 0) {
5217 }
else if (opt->values[GRCC_OPT_NoMultiEdge] < 0) {
5218 ok = ok && multiedge;
5220 if (opt->values[GRCC_OPT_No1PtBlock] > 0) {
5222 ok = ok && !tadblock;
5224 }
else if (opt->values[GRCC_OPT_No1PtBlock] < 0) {
5226 ok = ok && tadblock;
5229 if (opt->values[GRCC_OPT_Block] > 0) {
5231 }
else if (opt->values[GRCC_OPT_Block] < 0) {
5257 xcl = refineClass(cl);
5265 connected = isConnected();
5274 if (!isIsomorphic(xcl)) {
5283 egraph->fromMGraph(
this);
5284 if (egraph->isOptE()) {
5290 if (opt->outmg != NULL) {
5291 if (opt->proc != NULL) {
5292 egraph->mId = opt->proc->mgrcount + 1;
5293 egraph->sId = opt->proc->agrcount + 1;
5294 }
else if (opt->sproc != NULL) {
5295 egraph->mId = opt->sproc->mgrcount + 1;
5296 egraph->sId = opt->sproc->agrcount + 1;
5298 egraph->mId = cDiag;
5300 ok = (*(opt->outmg))(egraph, opt->argmg);
5306 wscon.add(1, nsym*esym);
5310 wsopi.add(1, nsym*esym);
5327 grcc_fprintf(GRCC_Stdout,
"\n");
5328 grcc_fprintf(GRCC_Stdout,
"Graph : %ld (%ld)", cDiag, ngen);
5329 grcc_fprintf(GRCC_Stdout,
" sym. factor = (%ld*%ld)\n", nsym, esym);
5337 opt->newMGraph(
this);
5348 if (xcl != cl && xcl != NULL) {
5356MOrbits::MOrbits(
void)
5363MOrbits::~MOrbits(
void)
5368void MOrbits::print(
void)
5372 grcc_fprintf(GRCC_Stdout,
"Orbits : nOrbits=%d: nd2or=", nOrbits);
5373 prIntArray(nNodes, nd2or,
"\n");
5374 for (c = 0; c < nOrbits; c++) {
5375 grcc_fprintf(GRCC_Stdout,
" %2d: (%2d -- %2d) :", c, flist[c], flist[c+1]-1);
5376 prIntArray(flist[c+1]-flist[c], or2nd+flist[c],
"\n");
5381void MOrbits::initPerm(
int nnodes)
5386 for (j = 0; j < nNodes; j++) {
5392void MOrbits::fromPerm(
int *perm)
5398 for (j = 0; j < nNodes; j++) {
5399 if (nd2or[j] != j) {
5405 for (j1 = 0; j1 < nNodes && (k = perm[k]) != j; j1++) {
5406 if (nd2or[k] == j) {
5411 for (l = 0; l < nNodes; l++) {
5412 if (nd2or[l] == k) {
5420 grcc_fprintf(GRCC_Stdout,
"*** fromPerm: illegal control: j=%d, k=%d\n", j, k);
5421 grcc_fprintf(GRCC_Stdout,
"perm=");
5422 prIntArray(nNodes, perm,
" nd2or=");
5423 prIntArray(nNodes, nd2or,
"\n");
5424 erEnd(
"fromPerm: illegal control");
5431void MOrbits::toOrbits(
void)
5434 int nd, nel, k, lor;
5439 for (nd = 0; nd < nNodes; nd++) {
5440 if (nd2or[nd] > lor) {
5443 flist[nOrbits++] = nel;
5444 for (k = nd; k < nNodes; k++) {
5445 if (nd2or[k] == lor) {
5451 flist[nOrbits] = nNodes;
5465MCEdge::~MCEdge(
void)
5485void MCOpi::init(
void)
5500MCBridge::MCBridge(
void)
5505MCBridge::~MCBridge(
void)
5512MCBlock::MCBlock(
void)
5518MCBlock::~MCBlock(
void)
5523void MCBlock::init(
void)
5534MConn::MConn(
int nnod,
int nedg)
5542 cedges =
new MCEdge[sedges];
5543 opics =
new MCOpi[snodes];
5546 articuls =
new int[snodes];
5549 opisp =
new int[snodes];
5550 opistk =
new int[snodes];
5551 blksp =
new Edge2n[sedges];
5552 blkstk =
new Edge2n[sedges];
5572void MConn::init(
void)
5597 for (j = 0; j < snodes; j++) {
5603void MConn::pushNode(
int nd)
5605 if (opistkptr >= snodes) {
5606 erEnd(
"MConn::pushNode: opistkptr >= snodes");
5608 opistk[opistkptr++] = nd;
5612void MConn::pushEdge(
int n0,
int n1)
5614 if (blkstkptr >= sedges) {
5615 erEnd(
"MConn::pushEdge: blkstkptr >= sedges");
5618 blkstk[blkstkptr][0] = n0;
5619 blkstk[blkstkptr][1] = n1;
5621 blkstk[blkstkptr][0] = n1;
5622 blkstk[blkstkptr][1] = n0;
5629void MConn::initCEdges(
MGraph *mg)
5634 for (n0 = 0; n0 < mg->nNodes; n0++) {
5635 for (e = 0; e < mg->adjMat[n0][n0]/2; e++, ed++) {
5636 cedges[ed].nodes[0] = n0;
5637 cedges[ed].nodes[1] = n0;
5638 cedges[ed].momdir = 0;
5640 for (n1 = n0+1; n1 < mg->nNodes; n1++) {
5641 for (e = 0; e < mg->adjMat[n0][n1]; e++, ed++) {
5642 cedges[ed].nodes[0] = n0;
5643 cedges[ed].nodes[1] = n1;
5644 cedges[ed].momdir = 0;
5649 grcc_fprintf(GRCC_Stdout,
"*** ed=%d != sedges=%d\n", ed, sedges);
5650 erEnd(
"MConn::initCEdge: table overflow");
5655void MConn::addCEdge(
int n0,
int n1, ULong momset)
5657 int m0, m1, dir, ed;
5668 for (ed = 0; ed < sedges; ed++) {
5669 if (cedges[ed].nodes[0] == m0 && cedges[ed].nodes[1] == m1 &&
5670 cedges[ed].momdir == 0) {
5671 cedges[ed].momdir = dir;
5672 cedges[ed].momset = momset;
5679void MConn::addOPIc(
MCOpi *mopi,
int stp)
5683 if (nopic >= snodes) {
5684 erEnd(
"MConn::addOPIc: table overflow");
5687 nn = opistkptr - stp;
5688 for (j = 0; j < nn; j++) {
5689 opisp[nopisp+j] = opistk[stp+j];
5691 opics[nopic].nnodes = nn;
5692 opics[nopic].nlegs = mopi->nlegs;
5693 opics[nopic].nedges = mopi->nedges;
5694 opics[nopic].next = mopi->next;
5695 opics[nopic].loop = mopi->loop;
5696 opics[nopic].ctloop = mopi->ctloop;
5697 opics[nopic].mom0lg = mopi->mom0lg;
5698 opics[nopic].nodes = opisp + nopisp;
5706void MConn::addBridge(
int n0,
int n1,
int nex,
int nextot)
5708 if (nbridges >= sedges) {
5709 erEnd(
"MConn::addBridge: table overflow");
5712 bridges[nbridges].nodes[0] = n0;
5713 bridges[nbridges].nodes[1] = n1;
5715 bridges[nbridges].nodes[0] = n1;
5716 bridges[nbridges].nodes[1] = n0;
5718 bridges[nbridges].next = Min(nex, nextot-nex);
5723void MConn::addArtic(
int nd,
int mul)
5725 articuls[nd] += mul;
5726 if (articuls[nd] == mul) {
5728 }
else if (articuls[nd] == 0) {
5734void MConn::addBlock(
MCBlock *mblk,
int stp)
5738 if (nopic >= snodes) {
5739 erEnd(
"MConn::addOPIc: table overflow");
5742 nn = blkstkptr - stp;
5743 for (j = 0; j < nn; j++) {
5744 blksp[nblksp+j][0] = blkstk[stp+j][0];
5745 blksp[nblksp+j][1] = blkstk[stp+j][1];
5747 blocks[nblocks].edges = blksp + nblksp;
5748 blocks[nblocks].nmedges = nn;
5749 blocks[nblocks].nartps = mblk->nartps;
5750 blocks[nblocks].loop = mblk->loop;
5758void MConn::addBlockSelf(
int nd,
int mult)
5768 for (j = 0; j < mult; j++) {
5771 addBlock(&mblk, stp);
5776void MConn::print(
void)
5780 grcc_fprintf(GRCC_Stdout,
"+++ MConn object: snodes=%d, sedges=%d\n", snodes, sedges);
5781 grcc_fprintf(GRCC_Stdout,
" nopic=%d, nlpopic=%d, nbacked=%d, nctopic=%d, "
5782 "nbridges=%d, ne0bridges=%d, ne1bridges=%d\n",
5783 nopic, nlpopic, nbacked, nctopic,
5784 nbridges, ne0bridges, ne1bridges);
5785 grcc_fprintf(GRCC_Stdout,
" nblocks=%d, neblocks=%d, na1blocks=%d, narticuls=%d, nselfloop=%d\n",
5786 nblocks, neblocks, na1blocks, narticuls, nselfloops);
5787 grcc_fprintf(GRCC_Stdout,
" nmultiedges=%d\n", nmultiedges);
5789 grcc_fprintf(GRCC_Stdout,
" cEdges (%d)\n", sedges);
5791 for (j = 0; j < sedges; j++) {
5792 grcc_fprintf(GRCC_Stdout,
" %2d: (%d,%d)[%2d*%2ld] \n", j,
5793 cedges[j].nodes[0], cedges[j].nodes[1],
5794 cedges[j].momdir, cedges[j].momset);
5797 grcc_fprintf(GRCC_Stdout,
"\n");
5799 grcc_fprintf(GRCC_Stdout,
" 1PI components (%d)\n", nopic);
5800 for (j = 0; j < nopic; j++) {
5801 grcc_fprintf(GRCC_Stdout,
" %d: nleg=%d, nnodes=%d, nedge=%d, ",
5802 j, opics[j].nlegs, opics[j].nnodes, opics[j].nedges);
5803 grcc_fprintf(GRCC_Stdout,
"next=%d, loop=%d, ctlp=%d: m0lg=%d [",
5804 opics[j].next, opics[j].loop, opics[j].ctloop,
5806 for (k = 0; k < opics[j].nnodes; k++) {
5807 grcc_fprintf(GRCC_Stdout,
" %d", opics[j].nodes[k]);
5809 grcc_fprintf(GRCC_Stdout,
"]\n");
5812 grcc_fprintf(GRCC_Stdout,
" bridges (%d)\n", nbridges);
5814 grcc_fprintf(GRCC_Stdout,
" ");
5815 for (j = 0; j < nbridges; j++) {
5816 grcc_fprintf(GRCC_Stdout,
"(%d,%d)[mom=%d] ",
5817 bridges[j].nodes[0], bridges[j].nodes[1], bridges[j].next);
5819 grcc_fprintf(GRCC_Stdout,
"\n");
5822 grcc_fprintf(GRCC_Stdout,
" blocks (%d)\n", nblocks);
5823 for (j = 0; j < nblocks; j++) {
5824 grcc_fprintf(GRCC_Stdout,
" %d: nmedges=%d, nartps=%d, loop=%d: [",
5825 j, blocks[j].nmedges, blocks[j].nartps, blocks[j].loop);
5826 for (k = 0; k < blocks[j].nmedges; k++) {
5827 grcc_fprintf(GRCC_Stdout,
" (%d,%d)", blocks[j].edges[k][0], blocks[j].edges[k][1]);
5829 grcc_fprintf(GRCC_Stdout,
"]\n");
5832 grcc_fprintf(GRCC_Stdout,
" articulation points (%d)\n", narticuls);
5833 if (narticuls > 0) {
5834 grcc_fprintf(GRCC_Stdout,
" ");
5835 for (j = 0; j < snodes; j++) {
5836 if (articuls[j] != 0) {
5838 grcc_fprintf(GRCC_Stdout,
"%d(%d) ", j, articuls[j]);
5840 grcc_fprintf(GRCC_Stdout,
"%d ", j);
5844 grcc_fprintf(GRCC_Stdout,
"\n");
5849void MConn::prEdges(
void)
5853 grcc_fprintf(GRCC_Stdout,
" cEdges (%d)", sedges);
5855 for (j = 0; j < sedges; j++) {
5857 grcc_fprintf(GRCC_Stdout,
"\n ");
5859 grcc_fprintf(GRCC_Stdout,
"(%d,%d)[%2d*%3ld] ",
5860 cedges[j].nodes[0], cedges[j].nodes[1],
5861 cedges[j].momdir, cedges[j].momset);
5863 grcc_fprintf(GRCC_Stdout,
"\n");
5883SGroup::~SGroup(
void)
5888 for (j = GRCC_MAXGROUP-1; j >= 0; j--) {
5889 if (elem[j] != NULL) {
5900void SGroup::print(
void)
5904 grcc_fprintf(GRCC_Stdout,
"SGroup : nnodes = %d, nelem=%ld, cgen=%d, csav=%d\n",
5905 nnodes, nelem, cgen, csav);
5906 grcc_fprintf(GRCC_Stdout,
" eclass =");
5907 prIntArray(neclass, eclass,
"\n");
5908 for (j = 0; j < nelem; j++) {
5909 grcc_fprintf(GRCC_Stdout,
" %4d: ", j);
5910 prIntArray(nnodes, elem[j],
"\n");
5915void SGroup::newGroup(
int nnd,
int nclss,
int *clss)
5923 if (neclass > GRCC_MAXNODES) {
5924 erEnd(
"newGroup : too many classes (GRCC_MAXNODES)");
5926 for (j = 0; j < neclass; j++) {
5927 eclass[j] = clss[j];
5932 elem =
new int*[GRCC_MAXGROUP];
5933 for (j = 0; j < GRCC_MAXGROUP; j++) {
5940void SGroup::clearGroup(
void)
5948void SGroup::delGroup(
void)
5952 for (j = 0; j < GRCC_MAXGROUP; j++) {
5953 if (elem[j] != NULL) {
5964int *SGroup::genNext(
void)
5966 cgen = nextPerm(nnodes, neclass, eclass, pgr, pgq, permg, cgen);
5974void SGroup::addGroup(
int *p)
5978 if (nelem >= GRCC_MAXGROUP) {
5979 erEnd(
"addGroup : too many elements (GRCC_MAXGROUP)");
5981 if (elem[nelem] == NULL) {
5982 elem[nelem] =
new int[GRCC_MAXNODES];
5984 for (j = 0; j < nnodes; j++) {
5985 elem[nelem][j] = p[j];
5991BigInt SGroup::nElem(
void)
5997int *SGroup::nextElem(
void)
6000 cgen = nextPerm(nelem, neclass, eclass,
6001 psr, psq, perms, cgen);
6007 if (cgen >= nelem) {
6010 return elem[cgen++];
6018static const char *GRCC_ND_names[] = GRCC_ND_NAMES;
6019static const char *GRCC_ED_names[] = GRCC_ED_NAMES;
6034 ndtype = GRCC_ND_Undef;
6040ENode::ENode(
EGraph *egrph,
int loops,
int sdeg)
6051 ndtype = GRCC_ND_Undef;
6055 edges =
new int[sdeg];
6056 klow =
new int[loops+1];
6071void ENode::copy(
ENode *en)
6076 extloop = en->extloop;
6077 ndtype = en->ndtype;
6078 for (d = 0; d < deg; d++) {
6079 edges[d] = en->edges[d];
6084void ENode::initAss(
EGraph *egrph,
int nid,
int sdg)
6092void ENode::setId(
EGraph *egrph,
const int nid)
6099void ENode::setExtern(
int typ,
int pt)
6102 if (typ == GRCC_AT_Initial || typ == GRCC_AT_Final) {
6105 grcc_fprintf(GRCC_Stderr,
"** ENode::setExternal:: illegal typ=%d\n",
6107 erEnd(
"ENode::setExternal:: illegal typ");
6112void ENode::setType(
int typ)
6115 if (ndtype != GRCC_ND_Undef && ndtype != typ) {
6116 grcc_fprintf(GRCC_Stderr,
"*** ndtype is already defined : old=%d, new = %d\n",
6118 erEnd(
"ndtype is already defined");
6125void ENode::print(
void)
6129 grcc_fprintf(GRCC_Stdout,
"Enode %d deg=%d, extl=%2d, intr=%d ",
6130 id, deg, extloop, intrct);
6131 if (egraph->bicount > 0) {
6132 grcc_fprintf(GRCC_Stdout,
"(%-8s) ", GRCC_ND_names[ndtype]);
6134 grcc_fprintf(GRCC_Stdout,
"edge=[");
6135 for (j = 0; j < deg; j++) {
6136 grcc_fprintf(GRCC_Stdout,
" %2d", edges[j]);
6138 grcc_fprintf(GRCC_Stdout,
"]\n");
6152 ptcl = GRCC_PT_Undef;
6155 edtype = GRCC_ED_Undef;
6166EEdge::EEdge(
EGraph *egrph,
int nedges,
int nloops)
6174 ptcl = GRCC_PT_Undef;
6177 edtype = GRCC_ED_Undef;
6180 emom =
new int[nedges];
6181 lmom =
new int[nloops];
6196void EEdge::copy(
EEdge *ee)
6198 nodes[0] = ee->nodes[0];
6199 nodes[1] = ee->nodes[1];
6203 edtype = ee->edtype;
6208void EEdge::setId(
EGraph *egrph,
const int eid)
6215void EEdge::setType(
int typ)
6218 if (edtype != GRCC_ED_Undef) {
6219 erEnd(
"edtype is already defined");
6226void EEdge::setEMom(
int nedges,
int *em,
int dir)
6230 for (e = 0; e < nedges; e++) {
6238void EEdge::setLMom(
int k,
int dir)
6244void EEdge::print(
void)
6248 grcc_fprintf(GRCC_Stdout,
"Edge %2d ext=%d ptcl=%2d ",
id, ext, ptcl);
6249 if (egraph == NULL || !egraph->assigned) {
6250 grcc_fprintf(GRCC_Stdout,
"[%d, %d]", nodes[0], nodes[1]);
6252 grcc_fprintf(GRCC_Stdout,
"[(%d,%d), (%d,%d))", nodes[0], nlegs[0], nodes[1], nlegs[1]);
6255 grcc_fprintf(GRCC_Stdout,
"[%2d*%2lu]", momdir, momset);
6258 if (egraph != NULL && egraph->bicount > 0) {
6260 grcc_fprintf(GRCC_Stdout,
" %-9s ",
"Cut");
6262 grcc_fprintf(GRCC_Stdout,
" %-9s ", GRCC_ED_names[edtype]);
6264 grcc_fprintf(GRCC_Stdout,
"c%2d: ", opicomp);
6265 grcc_fprintf(GRCC_Stdout,
"%s%d =", (ext)?
"Q":
"p", id);
6267 for (n = 0; n < egraph->nEdges; n++) {
6269 prMomStr(emom[n],
"Q", n);
6273 for (k = 0; k < egraph->nLoops; k++) {
6275 prMomStr(lmom[k],
"k", k);
6280 grcc_fprintf(GRCC_Stdout,
" 0");
6283 if (egraph == NULL) {
6284 grcc_fprintf(GRCC_Stdout,
" egraph=NULL");
6286 grcc_fprintf(GRCC_Stdout,
" egraph->bicount=%d", egraph->bicount);
6289 grcc_fprintf(GRCC_Stdout,
"\n");
6303void EFLine::print(
const char *msg)
6305 if (ftype == FL_Open) {
6306 grcc_fprintf(GRCC_Stdout,
" Open");
6307 }
else if (ftype == FL_Closed) {
6308 grcc_fprintf(GRCC_Stdout,
" Loop");
6310 grcc_fprintf(GRCC_Stdout,
" ?%d", ftype);
6312 if (fkind == GRCC_PT_Dirac) {
6313 grcc_fprintf(GRCC_Stdout,
" Dirac");
6314 }
else if (fkind == GRCC_PT_Majorana) {
6315 grcc_fprintf(GRCC_Stdout,
" Major");
6316 }
else if (fkind == GRCC_PT_Ghost) {
6317 grcc_fprintf(GRCC_Stdout,
" Ghost");
6319 grcc_fprintf(GRCC_Stdout,
" ?%d", fkind);
6321 grcc_fprintf(GRCC_Stdout,
" len=%d ", nlist);
6322 prIntArray(nlist, elist, msg);
6328EGraph::EGraph(
int nnodes,
int nedges,
int mxdeg)
6344 sLoops = sEdges - sNodes + 1;
6365 nodes =
new ENode*[sNodes];
6366 for (j = 0; j < sNodes; j++) {
6367 nodes[j] =
new ENode(
this, sNodes, sMaxdeg);
6369 edges =
new EEdge*[sEdges];
6370 for (j = 0; j < sEdges; j++) {
6371 edges[j] =
new EEdge(
this, sEdges, sLoops);
6374 bidef =
new int[sNodes];
6375 bilow =
new int[sNodes];
6378 extMom =
new int[sEdges+1];
6382 for (j = 0; j < GRCC_MAXFLINES; j++) {
6388EGraph::~EGraph(
void)
6392 for (j = GRCC_MAXFLINES-1; j >= 0; j--) {
6393 if (flines[j] != NULL) {
6404 for (j = 0; j < sEdges; j++) {
6410 for (j = 0; j < sNodes; j++) {
6419void EGraph::copy(
EGraph *eg)
6424 if (sNodes < eg->nNodes || sEdges < eg->nEdges ||
6425 sMaxdeg < eg->sMaxdeg || sLoops < eg->nLoops) {
6426 erEnd(
"EGraph::copy: sizes are too small");
6432 mgraph = eg->mgraph;
6434 econn = mgraph->mconn;
6438 gSubId = eg->gSubId;
6439 nNodes = eg->nNodes;
6440 nEdges = eg->nEdges;
6441 nExtern = eg->nExtern;
6442 nLoops = eg->nLoops;
6444 for (nd = 0; nd < nNodes; nd++) {
6445 nodes[nd]->copy(eg->nodes[nd]);
6447 for (ed = 0; ed < nEdges; ed++) {
6448 edges[ed]->copy(eg->edges[ed]);
6453 extperm = eg->extperm;
6461void EGraph::setExtLoop(
int nd,
int val)
6464 nodes[nd]->extloop = val;
6468void EGraph::endSetExtLoop(
void)
6475 for (n = 0; n < nNodes; n++) {
6476 if (isExternal(n)) {
6483void EGraph::fromDGraph(
DGraph *dg)
6485 int deg[GRCC_MAXNODES];
6486 int maxdeg, nextern, nconn, tc;
6489 if (dg->nnodes > GRCC_MAXNODES) {
6490 grcc_fprintf(GRCC_Stderr,
"*** too many nodes (GRCC_MAXNODES)\n");
6493 if (dg->nedges > GRCC_MAXEDGES) {
6494 grcc_fprintf(GRCC_Stderr,
"*** too many edges (GRCC_MAXEDGES)\n");
6498 for (e = 0; e < dg->nedges; e++) {
6499 if (dg->edges[e][0] >= dg->nnodes || dg->edges[e][0] >= dg->nnodes) {
6500 grcc_fprintf(GRCC_Stderr,
"*** undefined node:");
6501 grcc_fprintf(GRCC_Stderr,
"edge[%d] = {%d, %d}\n",
6502 e, dg->edges[e][0], dg->edges[e][0]);
6507 for (n0 = 0; n0 < dg->nnodes; n0++) {
6510 for (e = 0; e < dg->nedges; e++) {
6511 deg[dg->edges[e][0]]++;
6512 deg[dg->edges[e][1]]++;
6516 for (n0 = 0; n0 < dg->nnodes; n0++) {
6517 maxdeg = Max(maxdeg, deg[n0]);
6520 }
if (deg[n0] < 0) {
6521 grcc_fprintf(GRCC_Stderr,
"+++ node %d is isolated\n", n0);
6531 nNodes = dg->nnodes;
6532 nEdges = dg->nedges;
6547 for (n0 = 0; n0 < dg->nnodes; n0++) {
6549 nodes[n0]->extloop = dg->nodes[n0].extloop;
6551 for (e = 0; e < dg->nedges; e++) {
6552 n0 = dg->edges[e][0];
6553 n1 = dg->edges[e][1];
6554 edges[e]->nodes[0] = n0;
6555 edges[e]->nodes[1] = n1;
6556 edges[e]->ext = (isExternal(n0) || isExternal(n1));
6557 edges[e]->momdir = 0;
6559 nodes[n0]->edges[nodes[n0]->deg++] = I2Vedge(e, -1);
6560 nodes[n1]->edges[nodes[n1]->deg++] = I2Vedge(e, +1);
6562 for (n0 = 0; n0 < dg->nnodes; n0++) {
6563 nodes[n0]->setId(
this, n0);
6564 nodes[n0]->intrct = dg->nodes[n0].intrct;
6567 for (e = 0; e < nEdges; e++) {
6568 edges[e]->setId(
this, e);
6571 for (n0 = 0; n0 < nNodes; n0++) {
6572 if (isExternal(n0)) {
6573 setExtern(n0, 1, GRCC_AT_Initial);
6575 tc += 2*nodes[n0]->extloop + nodes[n0]->deg - 2;
6582 nLoops = nEdges - nNodes + nconn;
6589void EGraph::fromMGraph(
MGraph *mg)
6595 int n0, n1, m0, m1, ed, e;
6601 if (sNodes < mg->nNodes || sEdges < mg->nEdges || sMaxdeg < mg->maxdeg) {
6602 erEnd(
"EGraph::fromMGraph: sizes are too small");
6604 if (sLoops < mg->nLoops) {
6605 grcc_fprintf(GRCC_Stdout,
"too small sLoops : %d < %d\n", sLoops, mg->nLoops);
6606 erEnd(
"too small sLoops");
6611 econn = mgraph->mconn;
6613 if (sproc == NULL) {
6618 model = sproc->model;
6621 nNodes = mg->nNodes;
6622 nEdges = mg->nEdges;
6623 nLoops = mg->nLoops;
6624 nExtern = mg->nExtern;
6628 if (opt->proc != NULL) {
6629 mId = mg->opt->proc->mgrcount;
6630 sId = mg->opt->proc->agrcount;
6631 }
else if (mg->opt->sproc != NULL) {
6632 mId = mg->opt->sproc->mgrcount;
6633 sId = mg->opt->sproc->agrcount;
6644 maxdeg = mg->maxdeg;
6648 for (j = 0; j < proc->model->ncouple; j++) {
6649 totalc += proc->clist[j];
6652 for (ed = 0; ed < nEdges; ed++) {
6653 edges[ed]->edtype = GRCC_ED_Undef;
6657 for (n0 = 0; n0 < nNodes; n0++) {
6660 for (n0 = 0; n0 < nNodes; n0++) {
6661 for (e = 0; e < mg->adjMat[n0][n0]/2; e++, ed++) {
6662 edges[ed]->nodes[0] = n0;
6663 edges[ed]->nodes[1] = n0;
6664 nodes[n0]->edges[nodes[n0]->deg++] = I2Vedge(ed, -1);
6665 nodes[n0]->edges[nodes[n0]->deg++] = I2Vedge(ed, +1);
6666 edges[ed]->ext = isExternal(n0);
6668 for (n1 = n0+1; n1 < nNodes; n1++) {
6669 for (e = 0; e < mg->adjMat[n0][n1]; e++, ed++) {
6670 edges[ed]->nodes[0] = n0;
6671 edges[ed]->nodes[1] = n1;
6672 nodes[n0]->edges[nodes[n0]->deg++] = I2Vedge(ed, -1);
6673 nodes[n1]->edges[nodes[n1]->deg++] = I2Vedge(ed, +1);
6674 edges[ed]->ext = (isExternal(n0) || isExternal(n1));
6680 grcc_fprintf(GRCC_Stdout,
"*** EGraph::init: ed=%d != nEdges=%d\n",
6682 erEnd(
"EGraph::init: illegal connection");
6686 for (j = 0; j < nNodes; j++) {
6687 nodes[j]->setId(
this, j);
6688 nodes[j]->intrct = mg->nodes[j]->extloop;
6691 for (j = 0; j < nEdges; j++) {
6692 edges[j]->setId(
this, j);
6699 for (j = 0; j < ni; j++) {
6700 setExtern(j, proc->initlPart[j], GRCC_AT_Initial);
6704 for (j = 0; j < nf; j++) {
6705 setExtern(j+ni, proc->finalPart[j], GRCC_AT_Final);
6708 }
else if (sproc != NULL) {
6709 pnc = sproc->pnclass;
6710 for (j = 0; j < sproc->pnclass->nclass; j++) {
6711 if (pnc->type[j] == GRCC_AT_Initial
6712 || pnc->type[j] == GRCC_AT_External) {
6713 for (k = pnc->cl2nd[j]; k < pnc->cl2nd[j+1]; k++) {
6714 setExtern(j, k, GRCC_AT_Initial);
6717 }
else if (sproc->pnclass->type[j] == GRCC_AT_Final) {
6718 for (k = pnc->cl2nd[j]; k < pnc->cl2nd[j+1]; k++) {
6719 setExtern(j, k, GRCC_AT_Final);
6728 for (j = 0; j < nNodes; j++) {
6729 if (!isExternal(j)) {
6730 totalc += 2*nodes[j]->extloop + nodes[j]->deg - 2;
6734 for (ed = 0; ed < nEdges; ed++) {
6735 edges[ed]->momdir = 0;
6737 for (e = 0; e < econn->sedges; e++) {
6738 m0 = econn->cedges[e].nodes[0];
6739 m1 = econn->cedges[e].nodes[1];
6742 for (ed = 0; ed < nEdges; ed++) {
6743 if (edges[ed]->momdir != 0) {
6746 n0 = edges[ed]->nodes[0];
6747 n1 = edges[ed]->nodes[1];
6748 if (m0 == n0 && m1 == n1) {
6749 edges[ed]->momdir = econn->cedges[e].momdir;
6750 edges[ed]->momset = econn->cedges[e].momset;
6753 }
else if (m0 == n1 && m1 == n0) {
6754 edges[ed]->momdir = - econn->cedges[e].momdir;
6755 edges[ed]->momset = econn->cedges[e].momset;
6762 grcc_fprintf(GRCC_Stdout,
"*** EGraph::fromMGraph:edge (%d->%d) is not found\n",
6771ENode *EGraph::setExtern(
int n0,
int pt,
int ndtyp)
6779 grcc_fprintf(GRCC_Stderr,
"*** illegal external particle %d, %d, %d\n",
6782 erEnd(
"illegal external particle");
6786 if (ndtyp == GRCC_AT_Initial) {
6788 }
else if (ndtyp == GRCC_AT_Final) {
6793 grcc_fprintf(GRCC_Stderr,
"*** illegal type of an external particle : "
6794 "node = %d, particle = %d, type = %d", n0, pt, ndtyp);
6796 erEnd(
"illegal type of an external particle");
6804 npt = model->normalParticle(npt);
6805 ept = model->normalParticle(npt);
6811 nd->setExtern(ndtyp, npt);
6812 edges[eg]->ptcl = ept;
6818void EGraph::print(
void)
6824 nlp = nEdges - nNodes + 1;
6825 grcc_fprintf(GRCC_Stdout,
"\nEGraph\n");
6826 grcc_fprintf(GRCC_Stdout,
" pId=%d, gSubId=%ld, mId=%ld, aId=%ld, sId=%ld\n",
6827 pId, gSubId, mId, aId, sId);
6828 grcc_fprintf(GRCC_Stdout,
" sNodes=%d, sEdges=%d, sMaxdeg=%d, sLoops=%d\n",
6829 sNodes, sEdges, sMaxdeg, sLoops);
6830 grcc_fprintf(GRCC_Stdout,
" nNodes=%d, nEdges=%d, nExtern=%d, nLoops=%d, totalc=%d\n",
6831 nNodes, nEdges, nExtern, nlp, totalc);
6832 grcc_fprintf(GRCC_Stdout,
" ");
6833 if (model == NULL) {
6834 grcc_fprintf(GRCC_Stdout,
"model=NULL,");
6836 grcc_fprintf(GRCC_Stdout,
"model=\"%s\",", model->name);
6839 grcc_fprintf(GRCC_Stdout,
"proc=NULL,");
6841 grcc_fprintf(GRCC_Stdout,
"proc=%d,", proc->id);
6843 if (sproc == NULL) {
6844 grcc_fprintf(GRCC_Stdout,
"sproc=NULL,");
6846 grcc_fprintf(GRCC_Stdout,
"sproc=%d,", sproc->id);
6848 grcc_fprintf(GRCC_Stdout,
"\n");
6849 grcc_fprintf(GRCC_Stdout,
" assigned=%d, sym = (%ld * %ld) ",
6850 assigned, nsym, esym);
6851 grcc_fprintf(GRCC_Stdout,
"extperm=%ld, nsym1=%ld, multp=%ld\n",
6852 extperm, nsym1, multp);
6854 grcc_fprintf(GRCC_Stdout,
" Nodes\n");
6855 for (nd = 0; nd < nNodes; nd++) {
6856 if (isExternal(nd)) {
6857 grcc_fprintf(GRCC_Stdout,
" %2d Extern ", nd);
6859 grcc_fprintf(GRCC_Stdout,
" %2d Vertex ", nd);
6863 grcc_fprintf(GRCC_Stdout,
" Edges\n");
6864 for (ed = 0; ed < nEdges; ed++) {
6865 if (edges[ed]->ext) {
6866 grcc_fprintf(GRCC_Stdout,
" %2d Extern ", ed);
6868 grcc_fprintf(GRCC_Stdout,
" %2d Intern ", ed);
6873 grcc_fprintf(GRCC_Stdout,
" Biconn: nopicomp=%d, nopi2p=%d, opi2plp=%d, nadj2ptv=%d\n",
6874 nopicomp, nopi2p, opi2plp, nadj2ptv);
6876 grcc_fprintf(GRCC_Stdout,
"\n");
6884void EGraph::printPy(FILE *fp,
long mId)
6887 mgraph->printPy(fp, mId);
6891 fprintf(fp,
"#AGraph : (gseq, {spid, gid, ...}, nodes, node-group)\n");
6892 fprintf(fp,
"(%ld,\n", mId);
6895 fprintf(fp,
" {'proc':%d, 'subproc':%d, ", proc->id, sproc->id);
6896 fprintf(fp,
"'tgraph':%ld, ", mId);
6897 fprintf(fp,
"'OPI':%s, ", BOOLSTR(opt->values[GRCC_OPT_1PI]));
6898 fprintf(fp,
"'agraph':%ld, \n", aId);
6899 fprintf(fp,
" 'nnodes':%d, 'nnsym':%ld, 'nesym':%ld, },\n",
6900 nNodes, nsym, esym);
6902 fprintf(fp,
" [ #nodes\n");
6903 for (
int nd = 0; nd < nNodes; nd++) {
6904 ENode *end = nodes[nd];
6908 fprintf(fp,
" [ #node %d: ", nd);
6912 fprintf(fp,
" {'deg':%d, 'cid':%d, 'extern':%s, ",
6913 end->deg, 0, BOOLSTR(isExternal(nd)));
6914 fprintf(fp,
"'intr':%d, 'loop':%d, },\n",
6915 end->intrct, end->extloop);
6920 for (
int lg = 0; lg < end->deg; lg++) {
6921 int ed = V2Iedge(end->edges[lg]);
6922 int ptcl = edges[ed]->ptcl;
6923 int edlg = (edges[ed]->nodes[0] == nd) ? 1 : 0;
6924 int nnd = edges[ed]->nodes[edlg];
6925 int nlg = edges[ed]->nlegs[edlg];
6926 fprintf(fp,
"(%d, %d, %d), ", nnd, ptcl, nlg);
6928 fprintf(fp,
"],\n");
6929 fprintf(fp,
" ],\n");
6931 fprintf(fp,
" ], #node\n");
6932 fprintf(fp,
" [ # group of 1 elements\n");
6933 fprintf(fp,
" [], #\n");
6934 fprintf(fp,
" ], # group\n");
6939void EGraph::prFLines(
void)
6943 grcc_fprintf(GRCC_Stdout,
" Fermion lines %d, sign=%d (mId=%ld, aId=%ld)\n",
6944 nFlines, fsign, mId, aId);
6945 for (j = 0; j < nFlines; j++) {
6946 grcc_fprintf(GRCC_Stdout,
"%4d ", j);
6947 flines[j]->print(
"\n");
6952int EGraph::dirEdge(
int n,
int e)
6954 if (nodes[n]->edges[e] > 0) {
6962int EGraph::cmpMom(
int *lm0,
int *em0,
int *lm1,
int *em1)
6964 int lp, ed, cmp, sgn = 0;
6966 for (lp = 0; lp < nLoops; lp++) {
6970 }
else if (sgn == 0) {
6971 sgn = lm0[lp]*lm1[lp];
6973 if (Abs(sgn) != 1) {
6974 erEnd(
"cmpMom : illegal element of lmom");
6978 cmp = lm0[lp] - sgn*lm1[lp];
6983 }
else if (lm1[lp] != 0) {
6987 for (ed = 0; ed < nEdges; ed++) {
6991 }
else if (sgn == 0) {
6992 sgn = em0[ed]*em1[ed];
6994 if (Abs(sgn) != 1) {
6995 erEnd(
"cmpMom : illegal element of emom");
6999 cmp = em0[ed] - sgn*em1[ed];
7004 }
else if (em1[ed] != 0) {
7012int EGraph::groupLMom(
int *grp,
int *ed2gr)
7017 int ed0, ed1, cmp, ngrp = -1;
7019 for (ed0 = 0; ed0 < nEdges; ed0++) {
7024 for (ed0 = 0; ed0 < nEdges; ed0++) {
7025 if (ed2gr[ed0] < 0) {
7029 for (ed1 = ed0+1; ed1 < nEdges; ed1++) {
7030 cmp = cmpMom(edges[ed0]->lmom, edges[ed0]->emom,
7031 edges[ed1]->lmom, edges[ed1]->emom);
7045Bool EGraph::optQGrafM(
Options *opt)
7047 int *qgopt = opt->qgopt;
7048 int nopis[GRCC_MAXEDGES];
7050 mext = (~ mext) << nExtern;
7053 grcc_fprintf(GRCC_Stdout,
"optQGrafM:");
7058 grcc_fprintf(GRCC_Stdout,
"optQGrafM: %8ld\n", mId);
7064 for (
int j = 0; j < econn->nopic; j++) {
7065 maxlegs = Max(maxlegs, econn->opics[j].nlegs);
7067 if (maxlegs >= GRCC_MAXEDGES) {
7068 grcc_fprintf(GRCC_Stdout,
"*** table overflow\n");
7071 for (
int k = 0; k < GRCC_MAXEDGES; k++) {
7074 for (
int j = 0; j < econn->nopic; j++) {
7075 nopis[econn->opics[j].nlegs]++;
7079 if (qgopt[GRCC_QGRAF_OPT_ONEPI] > 0) {
7080 if (econn->nopic != 1) {
7083 }
else if (qgopt[GRCC_QGRAF_OPT_ONEPI] < 0) {
7084 if (econn->nopic < 2) {
7089 if (qgopt[GRCC_QGRAF_OPT_ONSHELL] != 0) {
7091 if (qgopt[GRCC_QGRAF_OPT_ONSHELL] > 0) {
7092 if (econn->nopic != 1) {
7095 }
else if (qgopt[GRCC_QGRAF_OPT_ONSHELL] < 0) {
7096 if (econn->nopic == 1) {
7101 if (qgopt[GRCC_QGRAF_OPT_ONSHELL] > 0) {
7102 if (econn->ne1bridges > 0) {
7105 }
else if (qgopt[GRCC_QGRAF_OPT_ONSHELL] < 0) {
7106 if (econn->ne1bridges <= 0) {
7113 if (qgopt[GRCC_QGRAF_OPT_NOSNAIL] != 0) {
7114 if (qgopt[GRCC_QGRAF_OPT_NOSNAIL] > 0) {
7116 if (econn->nblocks != 1) {
7120 if (mgraph->mconn->ne0bridges >= 1 ||
7121 mgraph->mconn->na1blocks >= 1) {
7125 }
else if (qgopt[GRCC_QGRAF_OPT_NOSNAIL] < 0) {
7127 if (econn->nblocks == 1) {
7131 if (mgraph->mconn->ne0bridges < 1 &&
7132 mgraph->mconn->na1blocks < 1) {
7139 if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] != 0) {
7141 if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] > 0) {
7142 if (econn->nopic != 1) {
7145 }
else if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] < 0) {
7146 if (econn->nopic == 1) {
7151 if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] > 0) {
7152 if (mgraph->mconn->ne0bridges != 0) {
7155 }
else if (qgopt[GRCC_QGRAF_OPT_NOTADPOLE] < 0) {
7156 if (mgraph->mconn->ne0bridges == 0) {
7163 if (qgopt[GRCC_QGRAF_OPT_NOSIGMA] != 0) {
7170 for (
int j = 0; j < econn->nopic; j++) {
7171 if (econn->opics[j].nlegs >= 2 &&
7172 econn->opics[j].nlegs == econn->opics[j].mom0lg) {
7178 for (
int j = 0; j < econn->sedges - 1; j++) {
7179 ULong momj = econn->cedges[j].momset;
7184 if (econn->cedges[j].nodes[0] < nExtern) {
7186 }
else if (econn->cedges[j].nodes[1] < nExtern) {
7189 for (
int k = j+1; k < econn->sedges; k++) {
7190 ULong momk = econn->cedges[k].momset;
7195 if (econn->cedges[k].nodes[0] < nExtern) {
7197 }
else if (econn->cedges[k].nodes[1] < nExtern) {
7200 if ( momj == momk) {
7202 if (extj + extk != 2) {
7211 if (qgopt[GRCC_QGRAF_OPT_NOSIGMA] > 0) {
7215 }
else if (qgopt[GRCC_QGRAF_OPT_NOSIGMA] < 0) {
7222 if (qgopt[GRCC_QGRAF_OPT_SIMPLE] > 0) {
7223 if (mgraph->selfloop || mgraph->multiedge) {
7226 }
else if (qgopt[GRCC_QGRAF_OPT_SIMPLE] < 0) {
7227 if (!mgraph->selfloop && !mgraph->multiedge) {
7232 if (qgopt[GRCC_QGRAF_OPT_BIPART] > 0) {
7233 if (! mgraph->bipart) {
7236 }
else if (qgopt[GRCC_QGRAF_OPT_BIPART] < 0) {
7237 if (mgraph->bipart) {
7242 if (qgopt[GRCC_QGRAF_OPT_CYCLI] != 0) {
7244 for (
int k = 0; k < econn->nblocks; k++) {
7245 if (econn->blocks[k].loop > 0) {
7249 if (qgopt[GRCC_QGRAF_OPT_CYCLI] > 0) {
7253 }
else if (qgopt[GRCC_QGRAF_OPT_CYCLI] < 0) {
7264Bool EGraph::optQGrafA(
Options *opt)
7267 grcc_fprintf(GRCC_Stdout,
"optQGrafA: %8ld\n", mId);
7271 if (opt->qgopt[GRCC_QGRAF_OPT_FLOOP] != 0) {
7272 for (
int fl=0; fl < nFlines; fl++) {
7273 if (flines[fl]->ftype == FL_Closed) {
7274 if (flines[fl]->nlist % 2 != 0) {
7281 if (opt->qgopt[GRCC_QGRAF_OPT_FLOOP] == -1) {
7282 retval = (retval == True ? False : True);
7289Bool EGraph::isOptE(
void)
7293 int grp[GRCC_MAXEDGES], ed2gr[GRCC_MAXEDGES];
7304 for (ed = 0; ed < nEdges; ed++) {
7305 edges[ed]->cut = False;
7311 if (! optQGrafM(opt)) {
7315 if (opt->values[GRCC_OPT_NoAdj2PtV] > 0) {
7319 }
else if (opt->values[GRCC_OPT_NoAdj2PtV] < 0) {
7324 if (opt->values[GRCC_OPT_No2PtL1PI] == 0) {
7328 if (nExtern == 2 && nopicomp == 1) {
7334 if (opi2plp > 0 && nopi2p >= minopi2p) {
7335 if (opt->values[GRCC_OPT_No2PtL1PI] > 0) {
7337 }
else if (opt->values[GRCC_OPT_No2PtL1PI] < 0) {
7343 ngrp = groupLMom(grp, ed2gr);
7345 for (g = 0; g < ngrp; g++) {
7350 for (ed = 0; ed < nEdges; ed++) {
7351 if (ed2gr[ed] == g) {
7352 if (edges[ed]->ext) {
7356 edup->edges[ed]->cut = True;
7358 edup->edges[ed]->cut = False;
7365 if (edup->nExtern == 2 && edup->nopicomp == 1) {
7371 if (edup->opi2plp > 0 && edup->nopi2p >= minopi2p) {
7372 if (opt->values[GRCC_OPT_No2PtL1PI] > 0) {
7374 }
else if (opt->values[GRCC_OPT_No2PtL1PI] < 0) {
7381 if (opt->values[GRCC_OPT_No2PtL1PI] > 0) {
7383 }
else if (opt->values[GRCC_OPT_No2PtL1PI] < 0) {
7386 erEnd(
"isOptE: illegal control");
7392int EGraph::findRoot(
void)
7394 int root, vr, er, e, nd0, nd1;
7399 for (e = 0; e < nEdges; e++) {
7400 if (edges[e]->cut || edges[e]->edtype == GRCC_ED_Deleted) {
7403 if (edges[e]->visited) {
7407 if (edges[e]->ext) {
7408 nd0 = edges[e]->nodes[0];
7409 nd1 = edges[e]->nodes[1];
7410 if (isExternal(nd0) && !isExternal(nd1)) {
7413 }
else if (!isExternal(nd0) && isExternal(nd1)) {
7420 er = edges[e]->nodes[0];
7423 vr = edges[e]->nodes[0];
7441int EGraph::connComp(
void)
7449 for (j = 0; j < nNodes; j++) {
7450 nodes[j]->visited = -1;
7454 while (nelem < nNodes) {
7455 for (j = 0; j < nNodes; j++) {
7456 if (nodes[j]->visited < 0) {
7463 nelem += connVisit(j, ncc);
7466 if (nelem != nNodes) {
7467 erEnd(
"EGraph::connComp: illegal control");
7473int EGraph::connVisit(
int nd,
int ncc)
7477 int e, en, nn, j, nelem;
7480 nodes[nd]->visited = ncc;
7481 for (e = 0; e < nodes[nd]->deg; e++) {
7482 en = V2Iedge(nodes[nd]->edges[e]);
7483 for (j = 0; j < 2; j++) {
7484 nn = edges[en]->nodes[j];
7485 if (nodes[nn]->visited < 0) {
7486 nelem += connVisit(nn, ncc);
7494void EGraph::biconnE(
void)
7497 int extlst[GRCC_MAXEDGES], intlst[GRCC_MAXEDGES];
7498 int opiext, opiloop;
7501 for (e = 0; e < nEdges; e++) {
7502 if (nodes[edges[e]->nodes[0]]->deg == 2 &&
7503 nodes[edges[e]->nodes[1]]->deg == 2 &&
7504 edges[e]->nodes[0] != edges[e]->nodes[1]) {
7515 for (e = 0; e < nEdges; e++) {
7524 bisearchE(root, extlst, intlst, &opiext, &opiloop);
7529 if (isExternal(root)) {
7535 opi2plp = Max(opi2plp, opiloop);
7538 for (ie = 0; ie < nEdges; ie++) {
7540 edges[ie]->opicomp = nopicomp;
7558void EGraph::biinitE(
void)
7569 for (n = 0; n < nNodes; n++) {
7572 nodes[n]->ndtype = GRCC_ND_Undef;
7573 for (j = 0; j < nLoops; j++) {
7574 nodes[n]->klow[j] = -1;
7578 for (e = 0; e < nEdges; e++) {
7580 edges[e]->visited = False;
7581 edges[e]->conid = -1;
7582 edges[e]->edtype = GRCC_ED_Undef;
7583 edges[e]->opicomp = -1;
7584 for (j = 0; j < nEdges; j++) {
7585 edges[e]->emom[j] = 0;
7587 for (j = 0; j < nLoops; j++) {
7588 edges[e]->lmom[j] = 0;
7594void EGraph::bisearchE(
int nd,
int *extlst,
int *intlst,
int *opiext,
int *opiloop)
7613 int k, e, ie, ed, td, dir, j;
7614 int opiext1, opiloop1;
7615 int extlst1[GRCC_MAXEDGES];
7616 int intlst1[GRCC_MAXEDGES];
7622 bidef[nd] = bicount;
7623 bilow[nd] = bicount;
7624 for (k = 0; k < nLoops; k++) {
7625 nodes[nd]->klow[k] = bicount;
7628 for (j = 0; j < nEdges; j++) {
7634 for (e = 0; e < nodes[nd]->deg; e++) {
7635 ed = V2Iedge(nodes[nd]->edges[e]);
7638 if (edges[ed]->edtype == GRCC_ED_Deleted) {
7641 }
else if (edges[ed]->visited) {
7644 }
else if (edges[ed]->cut) {
7647 nodes[nd]->setType(GRCC_ND_CPoint);
7649 }
else if (!isExternal(nd) && edges[ed]->ext) {
7651 edges[ed]->setType(GRCC_ED_Extern);
7652 edges[ed]->visited = True;
7653 edges[ed]->emom[ed] = 1;
7656 nodes[nd]->setType(GRCC_ND_CPoint);
7661 edges[ed]->visited = True;
7664 td = edges[ed]->nodes[0];
7667 td = edges[ed]->nodes[1];
7672 if (bidef[td] >= 0) {
7674 bilow[nd] = Min(bilow[nd], bilow[td]);
7675 edges[ed]->setType(GRCC_ED_Back);
7681 nodes[nd]->klow[k] = Min(nodes[nd]->klow[k], nodes[td]->klow[k]);
7682 edges[ed]->setLMom(k, dir);
7686 nodes[nd]->setType(GRCC_ND_CPoint);
7692 bisearchE(td, extlst1, intlst1, &opiext1, &opiloop1);
7695 for (j = 0; j < nEdges; j++) {
7696 extlst[j] = extlst[j] || extlst1[j];
7700 for (k = 0; k < nLoops; k++) {
7701 if (nodes[td]->klow[k] <= nodes[nd]->klow[k]) {
7703 edges[ed]->setLMom(k, dir);
7705 nodes[nd]->klow[k] = Min(nodes[nd]->klow[k], nodes[td]->klow[k]);
7709 if (bilow[td] >= bidef[nd]) {
7711 if (nodes[nd]->ndtype == GRCC_ND_Undef) {
7712 nodes[nd]->setType(GRCC_ND_CPoint);
7714 }
else if (nodes[nd]->ndtype != GRCC_ND_CPoint) {
7716 grcc_fprintf(GRCC_Stderr,
"bisearch: node %d is a cut point ",
7718 grcc_fprintf(GRCC_Stderr,
"(not undef %d)\n",
7726 if (bilow[td] > bidef[nd]) {
7728 if (edges[ed]->edtype == GRCC_ED_Undef) {
7729 edges[ed]->setType(GRCC_ED_Bridge);
7730 nodes[td]->setType(GRCC_ND_CPoint);
7732 }
else if (edges[ed]->edtype != GRCC_ED_Bridge) {
7734 grcc_fprintf(GRCC_Stderr,
"bisearch: edges %d is a bridge ", ed);
7735 grcc_fprintf(GRCC_Stderr,
"(not undef %d)\n", edges[ed]->edtype);
7740 if (!edges[ed]->ext) {
7743 for (ie = 0; ie < nEdges; ie++) {
7746 edges[ie]->opicomp = nopicomp;
7747 intlst1[ie] = False;
7752 opi2plp = Max(opi2plp, opiloop1);
7766 bilow[nd] = Min(bilow[nd], bilow[td]);
7768 if (edges[ed]->edtype == GRCC_ED_Undef) {
7769 edges[ed]->setType(GRCC_ED_Inloop);
7771 }
else if (edges[ed]->edtype != GRCC_ED_Inloop) {
7773 grcc_fprintf(GRCC_Stderr,
"bisearch: ");
7774 grcc_fprintf(GRCC_Stderr,
"edges %d is not undef (%d)\n",
7775 ed, edges[ed]->edtype);
7783 for (ie = 0; ie < nEdges; ie++) {
7784 intlst[ie] = intlst[ie] || intlst1[ie];
7786 (*opiext) += opiext1;
7787 (*opiloop) += opiloop1;
7790 edges[ed]->setEMom(nEdges, extlst1, dir);
7799void EGraph::extMomConsv(
void)
7804 for (e = 0; e < nEdges; e++) {
7805 if (edges[e]->ext) {
7807 extMom[e] = edges[e]->emom[e];
7816 for (e = 0; e < nEdges; e++) {
7817 rs = edges[e]->emom[lex];
7819 for (ex = 0; ex < nEdges; ex++) {
7820 edges[e]->emom[ex] -= rs*extMom[ex];
7827void EGraph::chkMomConsv(
void)
7831 int esum[GRCC_MAXEDGES];
7832 int lsum[GRCC_MAXNODES];
7834 int n, ex, lk, ej, e, dir;
7840 for (n = 0; n < nNodes; n++) {
7841 if (isExternal(n)) {
7844 for (ex = 0; ex < nEdges; ex++) {
7847 for (lk = 0; lk < nLoops; lk++) {
7850 for (ej = 0; ej < nodes[n]->deg; ej++) {
7851 e = V2Iedge(nodes[n]->edges[ej]);
7852 if (edges[e]->nodes[0] == edges[e]->nodes[0]) {
7855 dir = dirEdge(n, ej);
7856 for (ex = 0; ex < nEdges; ex++) {
7857 if (edges[e]->ext) {
7858 esum[ex] += dir*edges[e]->emom[ex];
7861 for (lk = 0; lk < nLoops; lk++) {
7862 lsum[lk] += dir*edges[e]->lmom[lk];
7867 for (ex = 0; ex < nEdges; ex++) {
7868 if (esum[ex] != 0) {
7870 grcc_fprintf(GRCC_Stderr,
"chkMomConsv:n=%d, esum[%d]=%d\n",
7875 for (lk = 0; lk < nLoops; lk++) {
7876 if (lsum[lk] != 0) {
7878 grcc_fprintf(GRCC_Stderr,
"chkMomConsv:n=%d, lsum[%d]=%d\n",
7886 grcc_fprintf(GRCC_Stderr,
"*** Violation of momentum conservation ");
7887 grcc_fprintf(GRCC_Stderr,
"at node =%d\n",n);
7893 erEnd(
"inconsistent momentum");
7898int EGraph::legParticle(
int ed,
int el)
7902 return model->normalParticle((2*el-1)*edges[ed]->ptcl);
7906int EGraph::isFermion(
int ed)
7912 ptcl = edges[ed]->ptcl;
7913 ptype = model->particles[Abs(ptcl)]->ptype;
7915 return (ptype == GRCC_PT_Dirac || ptype == GRCC_PT_Majorana
7916 || ptype == GRCC_PT_Ghost);
7920int EGraph::fltrace(
int fk,
int nd0,
int *fl)
7930 int nfl, k, i, nd, nl, e, ed, el, fgcnt, fkind, lk;
7934 for (k = 0; k < nEdges; k++) {
7936 grcc_fprintf(GRCC_Stdout,
"*** fltrace:illegal contorl: k=%d, nEdges=%d\n", k, nEdges);
7946 grcc_fprintf(GRCC_Stderr,
"*** fltrace: ed=%d > nEdges=%d, fl=",
7948 prIntArray(nNodes, fl,
"\n");
7949 erEnd(
"fltrace: illegal fl");
7952 nd = edges[ed]->nodes[el];
7953 nl = edges[ed]->nlegs[el];
7956 if (isExternal(nd) || nd == nd0) {
7965 for (i = 0; i < nodes[nd]->deg; i++) {
7966 e = nodes[nd]->edges[i];
7968 fkind = model->particles[Abs(edges[ed]->ptcl)]->ptype;
7980 if (!edges[ed]->visited) {
7981 if (!isFermion(ed)) {
7982 edges[ed]->visited = True;
7989 if ((lk == 1 && nl == i+1) || (lk == -1 && nl == i-1)) {
7990 edges[ed]->visited = True;
8003 grcc_fprintf(GRCC_Stdout,
"*** fline: Fermion number is not conserved\n");
8004 grcc_fprintf(GRCC_Stdout,
" nd=%d, e=%d, ed=%d, fgcnt=%d\n",
8007 }
else if (fgcnt > 1) {
8008 grcc_fprintf(GRCC_Stdout,
"+++ fline: more than two fermions: check fsign\n");
8009 grcc_fprintf(GRCC_Stdout,
" nd=%d, e=%d, ed=%d, fgcnt=%d\n",
8013 if (k >= nEdges || k >= nfl) {
8014 grcc_fprintf(GRCC_Stdout,
"*** fline: illegal control\n");
8015 grcc_fprintf(GRCC_Stdout,
" nEdges=%d, nfl=%d, k=%d, ", nEdges, nfl, k);
8016 prIntArray(nfl, fl,
"\n");
8017 erEnd(
"fline: illegal control");
8023void EGraph::getFLines(
void)
8025 int fl[GRCC_MAXNODES];
8026 int nextn, exto[GRCC_MAXNODES];
8027 int e, ed, nd, floop, nswap, nfl, fkind, el, ptcl;
8030 for (ed = 0; ed < nEdges; ed++) {
8031 edges[ed]->visited = False;
8035 for (ed = 0; ed < nEdges; ed++) {
8036 if (edges[ed]->visited) {
8039 if (!edges[ed]->ext) {
8043 nd = edges[ed]->nodes[el];
8044 if (!isExternal(nd)) {
8046 nd = edges[ed]->nodes[el];
8047 if (!isExternal(nd)) {
8048 erEnd(
"*** illegal control");
8051 if (!isFermion(ed)) {
8054 ptcl = legParticle(ed, 1-el);
8064 fl[0] = - I2Vedge(ed, el);
8066 fkind = model->particles[Abs(edges[ed]->ptcl)]->ptype;
8067 edges[ed]->visited = True;
8069 nfl = fltrace(fkind, nd, fl);
8070 addFLine(FL_Open, fkind, nfl, fl);
8072 exto[nextn++] = edges[V2Iedge(e)]->nodes[V2Ileg(e)];
8077 for (ed = 0; ed < nEdges; ed++) {
8078 if (edges[ed]->visited) {
8081 if (!isFermion(ed)) {
8085 ptcl = legParticle(ed, 1-el);
8089 nd = edges[ed]->nodes[el];
8095 fl[0] = - I2Vedge(ed, el);
8097 edges[ed]->visited = True;
8099 fkind = model->particles[Abs(edges[ed]->ptcl)]->ptype;
8100 nfl = fltrace(fkind, nd, fl);
8101 addFLine(FL_Closed, fkind, nfl, fl);
8105 nswap = sortb(nextn, exto);
8109 if ((floop+nswap) % 2 == 0) {
8117void EGraph::addFLine(
const FLType ft,
int fk,
int nfl,
int *fl)
8121 if (nFlines >= GRCC_MAXFLINES) {
8122 erEnd(
"too many Fermion lines (GRCC_MAXEDGES)");
8124 if (flines[nFlines] == NULL) {
8125 flines[nFlines] =
new EFLine();
8127 flines[nFlines]->ftype = ft;
8128 flines[nFlines]->fkind = fk;
8129 flines[nFlines]->nlist = nfl;
8130 for (j = 0; j < nfl; j++) {
8131 flines[nFlines]->elist[j] = fl[j];
8147NCand::NCand(
const NCandSt sta,
const int dega,
const int nilst,
int *ilst)
8156 for (j = 0; j < nilist; j++) {
8161 if (st == AS_Assigned) {
8163 erEnd(
"illegal assignment");
8165 }
else if (st == AS_AssExt) {
8167 erEnd(
"illegal assignment");
8179void NCand::prNCand(
const char* msg)
8181 grcc_fprintf(GRCC_Stdout,
"%d %d ", st, deg);
8182 prIntArray(nilist, ilist, msg);
8187ECand::ECand(
int dt,
int nplst,
int *plst)
8193 for (j = 0; j < nplist; j++) {
8197 if (det && nplist != 1) {
8199 grcc_fprintf(GRCC_Stderr,
"*** ECand : len(plist) != 1 : det=%d ", det);
8200 prIntArrayErr(nplist, plist,
"\n");
8202 erEnd(
"ECand : len(plist) != 1");
8212void ECand::prECand(
const char *msg)
8214 prIntArray(nplist, plist,
"");
8215 grcc_fprintf(GRCC_Stdout,
" (det=%d)%s", det, msg);
8226 anodes =
new int[deg];
8227 aedges =
new int[deg];
8228 aelegs =
new int[deg];
8230 for (j = 0; j < deg; j++) {
8244 if (aelegs != NULL) {
8255int ANode::newleg(
void)
8264 grcc_fprintf(GRCC_Stderr,
"*** ANode::newleg : nlegs = %d > deg = %d\n",
8266 erEnd(
"ANode::newleg : nlegs > deg");
8275AEdge::AEdge(
int n0,
int l0,
int n1,
int l1)
8281 ptcl = GRCC_PT_Undef;
8301 erEnd(
"Assign: sproc == NULL");
8306 model = sproc->model;
8311 egraph = mgraph->egraph;
8312 if (egraph == NULL) {
8313 erEnd(
"Assign: egraph == NULL");
8316 astack = sproc->astack;
8317 if (astack == NULL) {
8318 erEnd(
"Assign: astack == NULL");
8321 nNodes = mgraph->nNodes;
8322 nEdges = mgraph->nEdges;
8323 nExtern = sproc->nExtern;
8332 nodes =
new ANode*[nNodes];
8333 edges =
new AEdge*[nEdges];
8335 for (j = 0; j < nNodes; j++) {
8338 for (j = 0; j < nEdges; j++) {
8341 for (j = 0; j < model->ncouple; j++) {
8342 cplleft[j] = sproc->clist[j];
8346 astack->setAGraph(
this);
8347 astack->checkPoint(checkpoint0);
8352 checkAG(
"Assign::Assign");
8357 assignAllVertices();
8364 astack->restoreMsg(checkpoint0,
"Assign");
8366 astack->restore(checkpoint0);
8371Assign::~Assign(
void)
8375 for (j = 0; j < nEdges; j++) {
8382 for (j = 0; j < nNodes; j++) {
8391void Assign::prCand(
const char *msg)
8397 grcc_fprintf(GRCC_Stdout,
"\n");
8398 grcc_fprintf(GRCC_Stdout,
"+++ Candidate list: %s\n", msg);
8399 grcc_fprintf(GRCC_Stdout,
" Nodes %d\n", nNodes);
8400 for (n = 0; n < nNodes; n++) {
8401 grcc_fprintf(GRCC_Stdout,
"%d: edges=", n);
8402 prIntArray(nodes[n]->deg, nodes[n]->aedges,
": cand=");
8403 if (nodes[n]->cand == NULL) {
8404 grcc_fprintf(GRCC_Stdout,
"NULL\n");
8406 nodes[n]->cand->prNCand(
"\n");
8409 ne = Min(nEdges, nETotal);
8410 grcc_fprintf(GRCC_Stdout,
" Edges %d\n", ne);
8411 for (e = 0; e < ne; e++) {
8414 grcc_fprintf(GRCC_Stdout,
"NULL_Edge\n");
8416 grcc_fprintf(GRCC_Stdout,
"%d: %d->%d: cand=", e, ed->nodes[0], ed->nodes[1]);
8417 if (edges[e]->cand == NULL) {
8418 grcc_fprintf(GRCC_Stdout,
"NULL\n");
8420 edges[e]->cand->prECand(
"\n");
8427void Assign::checkAG(
const char *msg)
8432 for (n = 0; n < nNodes; n++) {
8433 for (lg = 0; lg < nodes[n]->deg; lg++) {
8434 if (nodes[n]->aedges[lg] < 0) {
8435 grcc_fprintf(GRCC_Stdout,
"*** checkAG:%s: n=%d, lg=%d, aedges=%d\n",
8436 msg, n, lg, nodes[n]->aedges[lg]);
8443 erEnd(
"checkAG: failed");
8451Bool Assign::assignAllVertices(
void)
8457 checkCand(
"assignAllVertices:1");
8462 ok = selectVertexSimp(-1);
8464 ok = selectVertex();
8468 checkCand(
"assignAllVertices:2");
8475Bool Assign::selectVertexSimp(
int lastv)
8486 checkCand(
"selectVertexSimp");
8490 v = selUnAssVertexSimp(lastv);
8497 if (!egraph->optQGrafA(opt)) {
8501 opt->newAGraph(egraph);
8508 ok = selectLeg(v, -1);
8514Bool Assign::selectVertex(
void)
8525 checkCand(
"selectVertex");
8529 v = selUnAssVertex();
8536 if (!egraph->optQGrafA(opt)) {
8540 opt->newAGraph(egraph);
8547 ok = selectLeg(v, -1);
8553Bool Assign::selectLeg(
int v,
int lastlg)
8561 int nplist, plist[GRCC_MAXMPARTICLES2];
8564 checkNode(v,
"selectLeg:0");
8565 checkCand(
"selectLeg");
8569 ln = selUnAssLeg(v, lastlg);
8574 ok = assignVertex(v);
8580 astack->checkPoint(sav0);
8581 nplist = candPart(v, ln, plist, GRCC_MAXMPARTICLES2);
8589 astack->checkPoint(sav);
8590 for (j = 0; j < nplist; j++) {
8594 if(assignPLeg(v, ln, pt)) {
8600 astack->restoreMsg(sav,
"selectLeg2");
8602 astack->restore(sav);
8607 astack->restoreMsg(sav0,
"selectLeg2");
8609 astack->restore(sav0);
8615void Assign::saveCouple(
int *sav)
8619 if (model->ncouple > 1) {
8620 for (j = 0; j < model->ncouple; j++) {
8621 sav[j] = cplleft[j];
8627void Assign::restoreCouple(
int *sav)
8631 if (model->ncouple > 1) {
8632 for (j = 0; j < model->ncouple; j++) {
8633 cplleft[j] = sav[j];
8639Bool Assign::subCouple(
int *cpl)
8643 if (model->ncouple > 1) {
8644 for (j = 0; j < model->ncouple; j++) {
8645 cplleft[j] -= cpl[j];
8646 if (cplleft[j] < 0) {
8655Bool Assign::assignVertex(
int v)
8665 int savc0[GRCC_MAXNCPLG], savc1[GRCC_MAXNCPLG];
8668 checkCand(
"assignVertex");
8672 astack->checkPoint(sav);
8677 for (n = 0; n < nNodes; n++) {
8680 cl = pnclass->nd2cl[n];
8681 if (!isATExternal(pnclass->type[cl])) {
8683 nc = nodes[n]->cand;
8684 if (nc->st == AS_UnAssLegs && nc->nilist == 1) {
8688 vst = assignIVertex(n, ia);
8689 if (vst == AS_Impossible) {
8692 astack->restoreMsg(sav,
"assignVertex");
8694 astack->restore(sav);
8696 restoreCouple(savc0);
8699 }
else if (vst == AS_Assigned) {
8700 if (!subCouple(model->interacts[ia]->clist)) {
8702 astack->restoreMsg(sav,
"assignVertex");
8704 astack->restore(sav);
8706 restoreCouple(savc0);
8711 done = (vst == AS_Assigned);
8722 ok = selectVertexSimp(v);
8724 ok = selectVertex();
8728 astack->restoreMsg(sav,
"assignVertex");
8730 astack->restore(sav);
8732 restoreCouple(savc0);
8739 astack->checkPoint(sav1);
8741 for (j = 0; j < nodes[v]->cand->nilist; j++) {
8742 ia = nodes[v]->cand->ilist[j];
8746 vst = assignIVertex(v, ia);
8747 if (vst == AS_Assigned) {
8748 ok1 = subCouple(model->interacts[ia]->clist);
8750 if (ok1 && (vst == AS_Assigned || vst == AS_Assigned0)) {
8753 ok = selectVertexSimp(v);
8755 ok = selectVertex();
8760 astack->restoreMsg(sav1,
"assignVertex");
8762 astack->restore(sav1);
8764 restoreCouple(savc1);
8768 astack->restoreMsg(sav,
"assignVertex");
8770 astack->restore(sav);
8772 restoreCouple(savc0);
8778Bool Assign::allAssigned(
void)
8784 BigInt nsym, esym, nsym1;
8794 checkCand(
"allAssigned");
8799 if (!checkOrderCpl()) {
8800 erEnd(
"allAssigned: checkOrderCpl = False");
8816 ok = isIsomorphic(cl, &nsym, &esym, &nsym1);
8817 if (!ok || nsym < 1 || esym < 1) {
8826 wAGraphs.add(1,nsym*esym);
8829 wAOPI.add(1,nsym*esym);
8833 checkAG(
"allAssigned");
8836 fillEGraph(nAGraphs, nsym, esym, nsym1);
8841 for (e = 0; e < nEdges; e++) {
8842 p = Abs(egraph->edges[e]->ptcl);
8843 ext = egraph->edges[e]->ext;
8845 if (model->particles[p]->extonly) {
8858Bool Assign::fromMGraph(
void)
8865 int j, k, n, nn, n1, deg, nc, ptcl, cl, typ, mcl, cl1, typ1;
8867 int nilist, ilist[GRCC_MAXMINTERACT];
8871 for (j = 0; j < nNodes; j++) {
8872 nodes[j] =
new ANode(mgraph->nodes[j]->deg);
8879 pall = model->allParticles(&npall);
8882 for (n = 0; n < mgraph->nNodes; n++) {
8883 cl = pnclass->nd2cl[n];
8884 typ = pnclass->type[cl];
8887 if (isATExternal(typ)) {
8890 if (mgraph->nodes[n]->deg != 1) {
8891 grcc_fprintf(GRCC_Stdout,
"*** assign:fromMGraph : "
8892 "external but deg[%d] = %d != 1, type=%d\n",
8893 n, mgraph->nodes[n]->deg, typ);
8895 erEnd(
"assign:fromMGraph : illegal external particle");
8898 np[0] = pnclass->particle[cl];
8899 nodes[n]->cand =
new NCand(AS_AssExt, 1, 1, np);
8900 for (n1 = 0; n1 < mgraph->nNodes; n1++) {
8901 nc = mgraph->adjMat[n][n1];
8902 for (k = 0; k < nc; k++) {
8903 addEdge(n, n1, 1, np);
8911 deg = mgraph->nodes[n]->deg;
8913 cl = sproc->pnclass->nd2cl[n];
8914 mcl = sproc->pnclass->cl2mcl[cl];
8916 if (nodes[n]->cand != NULL) {
8917 delete nodes[n]->cand;
8919 lcn = model->ncouple;
8924 for (j = 0; j < model->cplgnvl[mcl]; j++) {
8925 ia = model->cplgvl[mcl][j];
8928 if (leqArray(lcn, model->interacts[ia]->clist, cplleft)) {
8929 ilist[nilist++] = ia;
8932 nodes[n]->cand =
new NCand(st, deg, nilist, ilist);
8937 new NCand(st, deg, model->cplgnvl[mcl], model->cplgvl[mcl]);
8941 for (n1 = n; n1 < mgraph->nNodes; n1++) {
8942 cl1 = pnclass->nd2cl[n1];
8943 typ1 = pnclass->type[cl1];
8944 if (!isATExternal(typ1)) {
8945 nc = mgraph->adjMat[n][n1];
8949 for (k = 0; k < nc; k++) {
8950 addEdge(n, n1, npall, pall);
8957 if (nETotal != nEdges) {
8958 grcc_fprintf(GRCC_Stdout,
"*** Assign::fromMGraph nETotal=%d != nEdges=%d\n",
8960 erEnd(
"Assign::fromMGraph nETotal= != nEdges");
8965 for (n = 0; n < mgraph->nNodes; n++) {
8966 cl = pnclass->nd2cl[n];
8967 typ = pnclass->type[cl];
8968 if (isATExternal(typ)) {
8970 cl = pnclass->nd2cl[n];
8971 ptcl = pnclass->particle[cl];
8978 erEnd(
"fromMGraph: ptcl=0");
8981 ok = assignPLeg(n, 0, -ptcl);
8986 nn = nodes[n]->anodes[0];
8987 ok = updateCandNode(nn);
8996 checkCand(
"fromMGraph");
8997 checkAG(
"fromMGraph");
9004void Assign::addEdge(
int n0,
int n1,
int nplist,
int *plist)
9013 if (n0 >= nNodes || n1 >= nNodes) {
9014 grcc_fprintf(GRCC_Stdout,
"*** Assign::addEdge : undefined nodes %d: [%d, %d]",
9016 erEnd(
"Assign::addEdge : undefined nodes");
9021 lg0 = nodes[n0]->newleg();
9022 lg1 = nodes[n1]->newleg();
9026 if (nETotal >= nEdges) {
9027 erEnd(
"too many edges");
9031 aed =
new AEdge(n0, lg0, n1, lg1);
9032 aed->cand =
new ECand(False, nplist, plist);
9035 if (edges[nETotal] != NULL) {
9036 erEnd(
"edges[nETotal] != NULL");
9042 edges[nETotal++] = aed;
9045 connect(n0, lg0, eid, 0, n1, lg1);
9049void Assign::connect(
int n0,
int l0,
int eg,
int el,
int n1,
int l1)
9062 nd0->anodes[l0] = n1;
9063 nd0->aedges[l0] = eg;
9064 nd0->aelegs[l0] = el;
9066 nd1->anodes[l1] = n0;
9067 nd1->aedges[l1] = eg;
9068 nd1->aelegs[l1] = eo;
9078Bool Assign::fillEGraph(
int aid, BigInt nsym, BigInt esym, BigInt nsym1)
9086 int work[3][GRCC_MAXLEGS];
9094 egraph->assigned = True;
9097 egraph->nsym = nsym;
9098 egraph->esym = esym;
9099 egraph->bicount = -1;
9100 if (sproc != NULL) {
9101 egraph->extperm = sproc->extperm;
9103 egraph->extperm = 1;
9105 egraph->nsym1 = nsym1;
9106 egraph->multp = (egraph->extperm * egraph->nsym1) / egraph->nsym;
9109 checkAG(
"fillEGraph:0");
9112 for (n = 0; n < nNodes; n++) {
9115 cl = pnclass->nd2cl[n];
9116 if (isATExternal(pnclass->type[cl])) {
9118 }
else if (nodes[n]->cand->st != AS_Assigned) {
9119 grcc_fprintf(GRCC_Stdout,
"*** fillEGraph : node %d is not assigned", n);
9120 prCand(
"fillEGraph: node");
9121 erEnd(
"fillEGraph : node is not assigned");
9126 en = egraph->nodes[n];
9128 en->initAss(egraph, n, an->deg);
9130 en->intrct = an->cand->ilist[0];
9131 elist = reordLeg(n, work[0], work[1], work[2]);
9132 for (lr = 0; lr < nodes[n]->deg; lr++) {
9133 if (elist == NULL) {
9139 if (lg < 0 || lg >= nodes[n]->deg) {
9140 grcc_fprintf(GRCC_Stdout,
"*** fillEGraph: n=%d, lr=%d: 0 <= lg=%d < %d\n",
9141 n, lr, lg, nodes[n]->deg);
9142 erEnd(
"fillEGraph: illegal reordering");
9145 ed = an->aedges[lg];
9146 eg = an->aelegs[lg];
9148 en->edges[lr] = I2Vedge(ed, 2*eg-1);
9150 egraph->edges[ed]->nodes[eg] = n;
9151 egraph->edges[ed]->nlegs[eg] = lr;
9155 for (e = 0; e < nEdges; e++) {
9157 if (edges[e]->cand->nplist != 1) {
9158 grcc_fprintf(GRCC_Stdout,
"*** fillEGraph : edge %d is not assigned", e);
9159 prCand(
"fillEGraph: edge");
9160 erEnd(
"fillEGraph : edge is not assigned");
9163 egraph->edges[e]->ptcl = edges[e]->cand->plist[0];
9166 for (e = 0; e < nEdges; e++) {
9167 if (egraph->edges[e]->ptcl < 0) {
9168 egraph->edges[e]->ptcl = - edges[e]->cand->plist[0];
9169 n = egraph->edges[e]->nodes[0];
9170 lr = egraph->edges[e]->nlegs[0];
9171 egraph->edges[e]->nodes[0] = egraph->edges[e]->nodes[1];
9172 egraph->edges[e]->nlegs[0] = egraph->edges[e]->nlegs[1];
9173 egraph->edges[e]->nodes[1] = n;
9174 egraph->edges[e]->nlegs[1] = lr;
9176 egraph->nodes[n]->edges[lr] = - egraph->nodes[n]->edges[lr];
9177 n = egraph->edges[e]->nodes[0];
9178 lr = egraph->edges[e]->nlegs[0];
9179 egraph->nodes[n]->edges[lr] = - egraph->nodes[n]->edges[lr];
9184 for (ed = 0; ed < nEdges; ed++) {
9185 n = egraph->edges[ed]->nodes[0];
9186 lr = egraph->edges[ed]->nlegs[0];
9187 if (egraph->nodes[n]->edges[lr] != -ed-1) {
9188 grcc_fprintf(GRCC_Stderr,
"+++ node[%d][%d]=%d != - (edge[%d][0] + 1) = %d\n",
9189 n, lr, egraph->nodes[n]->edges[lr], e, -ed-1);
9192 n = egraph->edges[ed]->nodes[1];
9193 lr = egraph->edges[ed]->nlegs[1];
9194 if (egraph->nodes[n]->edges[lr] != ed+1) {
9195 grcc_fprintf(GRCC_Stderr,
"+++ node[%d][%d]=%d != + (edge[%d][0] + 1) = %d\n",
9196 n, lr, egraph->nodes[n]->edges[lr], e, ed+1);
9202 erEnd(
"fillEGraph: illegal connection");
9208 egraph->getFLines();
9211 checkAG(
"fillEGraph:0");
9218int *Assign::reordLeg(
int n,
int *reord,
int *plist,
int *used)
9225 int lg, ia, lr, deg, pt, ed;
9232 if (nodes[n]->cand->st == AS_AssExt) {
9238 deg = nodes[n]->deg;
9246 for (lg = 0; lg < deg; lg++) {
9247 ed = nodes[n]->aedges[lg];
9248 pt = edges[ed]->cand->plist[0];
9249 plist[lg] = legEdgeParticle(n, lg, pt);
9254 ia = nodes[n]->cand->ilist[0];
9255 ilegs = model->interacts[ia]->plist;
9261 for (lr = 0; lr < deg; lr++) {
9262 for (lg = 0; lg < deg; lg++) {
9263 if (used[lg] == 0 && plist[lg] == ilegs[lr]) {
9275 grcc_fprintf(GRCC_Stdout,
"*** reordLeg: illegal list of particles:"
9276 "interaction %d ", ia);
9277 prIntArray(deg, ilegs,
"; ");
9278 grcc_fprintf(GRCC_Stdout,
"vertex %d ", n);
9279 prIntArray(deg, plist,
"\n");
9281 erEnd(
"reordLeg: illegal list of particles");
9293int Assign::getLegParticle(
int n,
int ln)
9307 elg = nd->aelegs[ln];
9308 ed = nd->aedges[ln];
9309 pt = edges[ed]->cand->plist[0];
9313 return model->normalParticle(-pt);
9315 return model->normalParticle(pt);
9320int Assign::legEdgeParticle(
int n,
int ln,
int pt)
9332 if (nodes[n]->aelegs[ln] == 0) {
9333 return model->normalParticle(-pt);
9335 return model->normalParticle(pt);
9340int Assign::legPart(
int v,
int lg,
int nplist,
int *plist,
int *rlist,
const int size)
9348 for (j = 0; j < nplist; j++) {
9349 nrlist = intSetAdd(nrlist, rlist,
9350 legEdgeParticle(v, lg, plist[j]), size);
9356int Assign::candPart(
int v,
int ln,
int *plist,
const int size)
9364 en = nodes[v]->aedges[ln];
9367 if (edges[en]->cand->det) {
9368 grcc_fprintf(GRCC_Stdout,
"*** candPart : particle of leg (%d, %d) "
9369 "is assigned to %d\n",
9370 v, ln, edges[en]->cand->plist[0]);
9371 checkCand(
"candPart");
9372 mgraph->printAdjMat(mgraph->curcl);
9374 erEnd(
"candPart : particle of leg is assigned");
9378 ec = edges[en]->cand;
9380 ntplist = legPart(v, ln, ec->nplist, ec->plist, plist, size);
9389int Assign::selUnAssVertexSimp(
int lastv)
9399 v0 = Max(lastv+1, nExtern);
9400 for (v = v0; v < nNodes; v++) {
9401 if (nodes[v]->cand->st == AS_UnAssLegs) {
9409int Assign::selUnAssVertex(
void)
9421 for (v = 0; v < nNodes; v++) {
9422 if (nodes[v]->cand->st == AS_UnAssLegs) {
9423 if (nodes[v]->cand->nilist == 1) {
9425 }
else if (nodes[v]->cand->nilist > 1) {
9427 if (v0 < 0 || nodes[v]->cand->nilist < nl) {
9429 nl = nodes[v]->cand->nilist;
9438int Assign::selUnAssLeg(
int v,
int lastlg)
9451 for (lg = lg0; lg < nodes[v]->deg; lg++) {
9455 n0 = nodes[v]->anodes[lastlg];
9459 n1 = nodes[v]->anodes[lg];
9461 grcc_fprintf(GRCC_Stdout,
"*** selUnAssLeg: n0=%d > n1=%d\n", n0, n1);
9462 grcc_fprintf(GRCC_Stdout,
"*** illegal connection\n");
9463 erEnd(
"selUnAssLeg: n0 > n1");
9467 e = nodes[v]->aedges[lg];
9468 if (!edges[e]->cand->det) {
9477NCandSt Assign::assignIVertex(
int v,
int ia)
9481 int lplist[GRCC_MAXLEGS], iplist[GRCC_MAXLEGS];
9485 if (nodes[v]->cand->st == AS_Assigned || nodes[v]->cand->st == AS_AssExt) {
9486 return AS_Assigned0;
9489 deg = nodes[v]->deg;
9490 if (deg != model->interacts[ia]->nlegs) {
9491 return AS_Impossible;
9495 for (lg = 0; lg < deg; lg++) {
9496 e = nodes[v]->aedges[lg];
9497 if (edges[e]->cand->nplist != 1) {
9498 return AS_UnAssLegs;
9500 pt = edges[e]->cand->plist[0];
9501 lplist[lg] = legEdgeParticle(v, lg, pt);
9503 iplist[lg] = model->interacts[ia]->plist[lg];
9510 if (cmpArray(deg, lplist, deg, iplist) == 0) {
9512 astack->pushNode(v);
9513 nodes[v]->cand->st = AS_Assigned;
9514 nodes[v]->cand->nilist = 1;
9515 nodes[v]->cand->ilist[0] = ia;
9519 return AS_Impossible;
9523Bool Assign::assignPLeg(
int n,
int ln,
int pt)
9532 checkNode(n,
"assignPLeg0");
9533 checkCand(
"assignPLeg");
9541 if (edges[e]->cand->det) {
9545 if (!isOrdPLeg(n, ln, pt)) {
9550 ept = legEdgeParticle(n, ln, pt);
9553 if (!isIn(edges[e]->cand->nplist, edges[e]->cand->plist, ept)) {
9554 grcc_fprintf(GRCC_Stdout,
"*** assignPLeg: particle %d is not in the cand. of e=%d",
9556 edges[e]->cand->prECand(
"\n");
9557 prCand(
"assignPLeg");
9558 erEnd(
"assignPLeg: particle is not in the cand.");
9563 astack->pushEdge(e);
9566 edges[e]->cand->nplist = 1;
9567 edges[e]->cand->plist[0] = ept;
9575Bool Assign::detEdge(
int e)
9579 if (edges[e]->cand->nplist != 1) {
9582 edges[e]->cand->det = True;
9585 ok0 = updateCandNode(edges[e]->nodes[0]);
9587 ok1 = updateCandNode(edges[e]->nodes[1]);
9592 return (ok0 && ok1);
9596Bool Assign::isOrdPLeg(
int n,
int ln,
int pt)
9598 int e0, el, ln0, ep0, pt0, nn, n0, n1, ln1, pt1, e1, ep1;
9600 if (nodes[n]->deg < 2) {
9603 nn = nodes[n]->anodes[ln];
9612 e0 = nodes[n]->aedges[ln];
9613 el = nodes[n]->aelegs[ln];
9614 ln0 = edges[e0]->nlegs[1-el];
9615 ep0 = legEdgeParticle(n, ln, pt);
9616 pt0 = legEdgeParticle(n0, ln0, ep0);
9619 for (ln1 = 0; ln1 < nodes[n0]->deg; ln1++) {
9620 if (ln1 == ln0 || n1 != nodes[n0]->anodes[ln1]) {
9623 e1 = nodes[n0]->aedges[ln1];
9624 if (!edges[e1]->cand->det) {
9627 ep1 = edges[e1]->cand->plist[0];
9628 pt1 = legEdgeParticle(n0, ln1, ep1);
9629 if ((ln0 < ln1 && pt0 > pt1) ||
9630 (ln0 > ln1 && pt0 < pt1)) {
9640Bool Assign::candPartClassify(
int v,
int *npdass,
int *pdass,
int *npuass,
int *puass,
const int size)
9649 int lg, e, npl, ass, jd, ju, j, pt, pts;
9653 for (lg = 0; lg < nodes[v]->deg; lg++) {
9654 e = nodes[v]->aedges[lg];
9655 npl = edges[e]->cand->nplist;
9660 for (j = 0; j < edges[e]->cand->nplist; j++) {
9661 pt = edges[e]->cand->plist[j];
9662 pts = legEdgeParticle(v, lg, pt);
9664 jd = intSListAdd(jd, pdass, pts, size);
9666 ju = intSListAdd(ju, puass, pts, size);
9677Bool Assign::updateCandNode(
int v)
9684 int npdass, pdass[GRCC_MAXPSLIST];
9685 int npuass, puass[GRCC_MAXPSLIST];
9686 int nsub0, sub0[GRCC_MAXPSLIST];
9687 int niplist, iplist[GRCC_MAXPSLIST];
9688 int nd0, d0[GRCC_MAXMPARTICLES2];
9689 int nd1, d1[GRCC_MAXMPARTICLES2];
9690 int ns0, s0[GRCC_MAXMPARTICLES2];
9691 int neplist, eplist[GRCC_MAXMPARTICLES2];
9692 int nitlist, itlist[GRCC_MAXMINTERACT];
9693 int e, it, j, k, i, lg;
9696 if (nodes[v]->cand->st == AS_AssExt) {
9698 e = nodes[v]->aedges[0];
9699 if (e < 0 || edges[e]->cand->nplist != 1) {
9700 grcc_fprintf(GRCC_Stdout,
"*** illegal external node: v=%d e=%d :", v, e);
9701 prCand(
"updateCandNode");
9702 erEnd(
"illegal external node");
9710 if (nodes[v]->cand->nilist < 1) {
9720 if (!candPartClassify(v, &npdass, pdass, &npuass, puass, GRCC_MAXPSLIST)) {
9736 for (j = 0; j < nodes[v]->cand->nilist; j++) {
9737 it = nodes[v]->cand->ilist[j];
9741 if (isSubSList(npdass, pdass,
9742 model->interacts[it]->nslist,
9743 model->interacts[it]->slist) ) {
9744 nsub0 = subtrSet(model->interacts[it]->nslist,
9745 model->interacts[it]->slist,
9746 npdass, pdass, sub0, GRCC_MAXPSLIST);
9749 nitlist = intSetAdd(nitlist, itlist, it, GRCC_MAXMINTERACT);
9752 if (isSubSList(nsub0, sub0, npuass, puass)) {
9753 nitlist = intSetAdd(nitlist, itlist, it, GRCC_MAXMINTERACT);
9754 for (k = 0; k < nsub0; k++) {
9755 niplist = intSetAdd(niplist,iplist,sub0[k],GRCC_MAXPSLIST);
9766 if (nitlist != nodes[v]->cand->nilist) {
9767 astack->pushNode(v);
9768 nodes[v]->cand->nilist = nitlist;
9769 for (i = 0; i < nitlist; i++) {
9770 nodes[v]->cand->ilist[i] = itlist[i];
9773 }
else if (nitlist > nodes[v]->cand->nilist) {
9774 erEnd(
"larger ncand");
9779 for (lg = 0; lg < nodes[v]->deg; lg++) {
9780 e = nodes[v]->aedges[lg];
9781 if (edges[e]->cand->nplist > 1) {
9784 neplist = legPart(v, lg, niplist, iplist, eplist, GRCC_MAXMPARTICLES2);
9785 listDiff(edges[e]->cand->nplist, edges[e]->cand->plist,
9787 &nd0, d0, &ns0, s0, &nd1, d1);
9792 if (ns0 < edges[e]->cand->nplist) {
9793 astack->pushEdge(e);
9794 edges[e]->cand->nplist = ns0;
9795 for (i = 0; i < ns0; i++) {
9796 edges[e]->cand->plist[i] = s0[i];
9801 for (lg = 0; lg < nodes[v]->deg; lg++) {
9802 e = nodes[v]->aedges[lg];
9803 if (edges[e]->cand->nplist == 1 && !edges[e]->cand->det) {
9816Bool Assign::checkOrderCpl(
void)
9820 int ord[GRCC_MAXNCPLG];
9821 int lcn, n, ia, j, cl;
9824 lcn = model->ncouple;
9831 for (j = 0; j < GRCC_MAXNCPLG; j++) {
9834 for (n = 0; n < nNodes; n++) {
9835 cl = pnclass->nd2cl[n];
9836 if (!isATExternal(pnclass->type[cl])) {
9837 ia = nodes[n]->cand->ilist[0];
9838 for (j = 0; j < lcn; j++) {
9839 ord[j] += model->interacts[ia]->clist[j];
9845 for (j = 0; j < lcn; j++) {
9846 if (sproc->clist[j] != ord[j]) {
9854Bool Assign::isOrdLegs(
void)
9860 int v, l0, v0, l1, v1, e0, e1, q0, q1, p0, p1, cl;
9862 for (v = 0; v < nNodes; v++) {
9863 cl = pnclass->nd2cl[v];
9864 if (!isATExternal(pnclass->type[cl])) {
9865 for (l0 = 0; l0 < nodes[v]->deg; l0++) {
9866 v0 = nodes[v]->anodes[l0];
9867 for (l1 = l0+1; l1 < nodes[v]->deg; l1++) {
9868 v1 = nodes[v]->anodes[l1];
9870 if (v0 == v1 && v <= v0) {
9872 e0 = nodes[v]->aedges[l0];
9873 e1 = nodes[v]->aedges[l1];
9874 q0 = edges[e0]->cand->plist[0];
9875 q1 = edges[e1]->cand->plist[0];
9876 p0 = legEdgeParticle(v, l0, q0);
9877 p1 = legEdgeParticle(v, l1, q1);
9892Bool Assign::isIsomorphic(
MNodeClass *cl, BigInt *nsym, BigInt *esym, BigInt *nsym1)
9906 ngelem = mgraph->group->nElem();
9908 if (mgraph->nsym > 1 && ngelem <= 1) {
9909 grcc_fprintf(GRCC_Stdout,
"*** isIsomorphic: illegal group: "
9910 "ngelem=%ld, mgraph->sym=(%ld, %ld)\n",
9911 ngelem, mgraph->nsym, mgraph->esym);
9912 erEnd(
"Assign::isIsomorphic: illegal group");
9921 for (j = 0; j < ngelem; j++) {
9925 p = mgraph->group->elem[j];
9927 cmp = cmpPermGraph(p, cl);
9931 }
else if(cmp == 0) {
9934 if (opt->values[GRCC_OPT_SymmInitial] || opt->values[GRCC_OPT_SymmFinal]) {
9936 for (n = 0; n < nNodes; n++) {
9937 cln = pnclass->nd2cl[n];
9938 if (!isATExternal(pnclass->type[cln])) {
9968int Assign::cmpPermGraph(
int *p,
MNodeClass *cl)
9972 int n, cmp, n1, n2, p1, p2, j;
9975 int njn, jn[GRCC_MAXLEGS];
9976 int njp, jp[GRCC_MAXLEGS];
9980 erEnd(
"Assign::cmpPermGraph: p==NULL");
9983 for (n = 0; n < nNodes; n++) {
9984 if (!isATExternal(pnclass->type[pnclass->nd2cl[n]])) {
9985 cmp = cmpNodes(n, p[n], cl);
9994 for (n1 = 0; n1 < nNodes; n1++) {
9998 for (n2 = n1; n2 < nNodes; n2++) {
9999 if (mgraph->adjMat[n1][n2] == 0) {
10004 if (p1 == n1 && p2 == n2) {
10007 cmp = mgraph->adjMat[n1][n2] - mgraph->adjMat[p1][p2];
10014 for (j = 0; j < nd1->deg; j++) {
10015 if (nd1->anodes[j] == n2) {
10016 jn[njn++] = getLegParticle(n1, j);
10020 for (j = 0; j < np1->deg; j++) {
10021 if (np1->anodes[j] == p2) {
10022 jp[njp++] = getLegParticle(p1, j);
10032 if (mgraph->adjMat[n1][n2] >= 2) {
10037 for (j = 0; j < njn; j++) {
10038 cmp = jn[j] - jp[j];
10050int Assign::cmpNodes(
int nd0,
int nd1,
MNodeClass *cn)
10059 cmp = cn->ndcl[nd0] - cn->ndcl[nd1];
10065 cmp = nodes[nd0]->cand->ilist[0] - nodes[nd1]->cand->ilist[0];
10070BigInt Assign::edgeSym(
void)
10074 int lg[GRCC_MAXLEGS], lt[GRCC_MAXLEGS], lc;
10076 int n1, n2, j, k, mult;
10080 for (n1 = 0; n1 < nNodes; n1++) {
10082 for (n2 = n1; n2 < nNodes; n2++) {
10083 if (mgraph->adjMat[n1][n2] > 1) {
10085 for (j = 0; j < nd1->deg; j++) {
10086 if (nd1->anodes[j] == n2) {
10088 lt[lc] = getLegParticle(n1, j);
10092 for (j = 0; j < lc-1; j++) {
10094 for (k = lc-1; k > j; k--) {
10095 if (lt[j] == lt[k]) {
10104 for (j = 0; j < lc; j++) {
10112 esym *= ipow(2,mult/2)*factorial(mult/2);
10114 esym *= factorial(mult);
10129Bool Assign::checkCand(
const char *msg)
10141 for (n = 0; n < nNodes; n++) {
10146 if (nc->st == AS_UnAssLegs) {
10149 }
else if (nc->st == AS_Assigned) {
10150 if (nc->nilist < 1) {
10151 grcc_fprintf(GRCC_Stdout,
"*** checkCand:7:%s:status (%d) of node %d says"
10152 " interaction is assigned to %d but ilist=",
10153 msg, nc->st, n, nc->st);
10154 prIntArray(nc->nilist, nc->ilist,
"\n");
10158 for (lg = 0; lg < nc->deg; lg++) {
10159 e = na->aedges[lg];
10160 ec = edges[e]->cand;
10161 if (ec->nplist != 1) {
10162 grcc_fprintf(GRCC_Stdout,
"*** checkCand:8:%s:status (%d) of node %d says"
10163 " interaction is assigned "
10164 " but unassigned edge %d is found\n",
10165 msg, nc->st, n, e);
10171 }
else if (nc->st == AS_AssExt) {
10175 grcc_fprintf(GRCC_Stdout,
"*** checkCand:10:%s:illegal status of node %d : %d",
10183 for (e = 0; e < nEdges; e++) {
10184 ec = edges[e]->cand;
10185 if (ec != NULL && ec->nplist < 1) {
10186 grcc_fprintf(GRCC_Stdout,
"*** checkCand:12:%s:illegal edge %d\n", msg, e);
10192 grcc_fprintf(GRCC_Stdout,
"*** checkCand:15:%s:illegal configuration\n", msg);
10193 prCand(
"checkCand");
10194 grcc_fprintf(GRCC_Stdout,
"*** checkCand:16:illegal configuration\n");
10195 erEnd(
"checkCand:16:illegal configuration");
10201void Assign::checkNode(
int n,
const char *msg)
10205 for (j = 0; j < nodes[n]->cand->nilist; j++) {
10206 it = nodes[n]->cand->ilist[j];
10207 if (Abs(it) >= GRCC_MAXMINTERACT) {
10208 grcc_fprintf(GRCC_Stdout,
"*** %s: n=%d, j=%d, it=%d\n", msg, n, j, it);
10209 nodes[n]->cand->prNCand(msg);
10210 grcc_fprintf(GRCC_Stdout,
"\n");
10211 erEnd(
"checkNode:illegal it");
10226void NStack::print(
const char *msg)
10228 grcc_fprintf(GRCC_Stdout,
" node=%d, deg=%d, st=%d, ilist=", noden, deg, st);
10229 prilist(nilist, ilist, msg);
10233void EStack::print(
const char *msg)
10235 grcc_fprintf(GRCC_Stdout,
" edge=%d, det=%d, plist=", edgen, det);
10236 prilist(nplist, plist, msg);
10240AStack::AStack(
int nsize,
int esize)
10249 nStack =
new NStack*[nSize];
10254 eStack =
new EStack*[eSize];
10262 for (j = 0; j < nSize; j++) {
10263 nStack[j] =
new NStack();
10266 for (j = 0; j < eSize; j++) {
10267 eStack[j] =
new EStack();
10272AStack::~AStack(
void)
10276 if (eStack != NULL) {
10277 for (j = eSize-1; j >= 0; j--) {
10278 if (eStack[j] != NULL) {
10287 if (nStack != NULL) {
10288 for (j = nSize-1; j >= 0; j--) {
10289 if (nStack[j] != NULL) {
10300void AStack::setAGraph(
Assign *ag)
10306void AStack::pushNode(
int n)
10313 if (agraph == NULL) {
10314 erEnd(
"pushNode: agraph is not defined");
10317 if (nStackP >= GRCC_MAXNSTACK) {
10318 erEnd(
"N-stack overflow (GRCC_MAXNSTACK)");
10320 nc = agraph->nodes[n]->cand;
10321 if (nc->nilist >= GRCC_MAXMINTERACT) {
10322 erEnd(
"N-stack: too long list (GRCC_MAXMINTERACT)");
10324 ns = nStack[nStackP];
10328 ns->nilist = nc->nilist;
10329 for (j = 0; j < nc->nilist; j++) {
10330 ns->ilist[j] = nc->ilist[j];
10336void AStack::pushEdge(
int e)
10343 if (agraph == NULL) {
10344 erEnd(
"pushEdge: agraph is not defined");
10347 if (eStackP >= GRCC_MAXESTACK) {
10348 erEnd(
"E-stack overflow (GRCC_MAXESTACK)");
10350 ec = agraph->edges[e]->cand;
10351 if (ec->nplist >= GRCC_MAXMPARTICLES2) {
10352 erEnd(
"E-stack: Too long list (GRCC_MAXMPARTICLES)");
10354 es = eStack[eStackP];
10357 es->nplist = ec->nplist;
10358 for (j = 0; j < ec->nplist; j++) {
10359 es->plist[j] = ec->plist[j];
10365void AStack::checkPoint(CheckPt sav)
10372void AStack::restoreNode(
int spr)
10377 for (sp = nStackP-1; sp >= spr; sp--) {
10380 agraph->nodes[n]->cand->deg = stc->deg;
10381 agraph->nodes[n]->cand->st = stc->st;
10382 agraph->nodes[n]->cand->nilist = stc->nilist;
10383 for (j = 0; j < stc->nilist; j++) {
10384 agraph->nodes[stc->noden]->cand->ilist[j] = stc->ilist[j];
10391void AStack::restoreEdge(
int spr)
10396 for (sp = eStackP-1; sp >= spr; sp--) {
10398 e = eStack[sp]->edgen;
10399 agraph->edges[e]->cand->det = stc->det;
10400 agraph->edges[e]->cand->nplist = stc->nplist;
10401 for (j = 0; j < stc->nplist; j++) {
10402 agraph->edges[e]->cand->plist[j] = stc->plist[j];
10409void AStack::restore(CheckPt sav)
10411 restoreNode(sav[0]);
10412 restoreEdge(sav[1]);
10417void AStack::restoreMsg(CheckPt sav,
const char *msg)
10419 if (agraph == NULL) {
10420 erEnd(
"restore: agraph is not defined");
10425 if (!agraph->checkCand(
"restore")) {
10426 grcc_fprintf(GRCC_Stdout,
"restore is called from %s\n", msg);
10432void AStack::prStack(
void)
10436 grcc_fprintf(GRCC_Stdout,
"+++ prStack : (%d, %d)", nStackP, eStackP);
10437 for (j = 0; j < nStackP; j++) {
10438 grcc_fprintf(GRCC_Stdout,
"N:%4d ", j);
10439 nStack[j]->print(
"\n");
10441 for (j = 0; j < eStackP; j++) {
10442 grcc_fprintf(GRCC_Stdout,
"E:%4d ", j);
10443 eStack[j]->print(
"\n");
10450Fraction::Fraction(BigInt n, BigInt d)
10457 ratio = Real(n)/Real(d);
10461void Fraction::print(
const char *msg)
10463 double err = Abs(Real(num)/Real(den) - ratio);
10465 if (err > GRCC_FRACERROR) {
10466 grcc_fprintf(GRCC_Stdout,
"%ld/%ld(%g)(overflow)%s", num, den, ratio, msg);
10468 grcc_fprintf(GRCC_Stdout,
"%ld/%ld(%g)%s", num, den, ratio, msg);
10473void Fraction::setValue(BigInt n, BigInt d)
10477 ratio = Real(n)/Real(d);
10481void Fraction::setValue(
Fraction &f)
10489BigInt Fraction::gcd(BigInt n0, BigInt n1)
10499 r = nn[nc] % nn[1-nc];
10509void Fraction::normal(
void)
10514 erEnd(
"Fraction: den==0");
10526void Fraction::add(BigInt n, BigInt d)
10540 num = num * d1 + n * (den/g);
10545 ratio += Real(n)/Real(d);
10551 double r = ratio + f.ratio;
10560 double r = ratio - f.ratio;
10562 add(-f.num, f.den);
10569 return (num == f.num && den == f.den);
10579static void grcc_fprintf(FILE* out,
const char* fmt, ...)
10582 va_start(args, fmt);
10588 va_copy(args_copy, args);
10590 const int len = vsnprintf(NULL, 0, fmt, args);
10591 char *buffer =
new char[len+1];
10592 vsnprintf(buffer, len+1, fmt, args_copy);
10593 MLOCK(ErrorMessageLock);
10594 MesPrint(
"%s%", buffer);
10595 MUNLOCK(ErrorMessageLock);
10599 vfprintf(out, fmt, args);
10606static void erEnd(
const char *msg)
10608 if (erExit != NULL) {
10609 (*erExit)(msg, erExitArg);
10611 grcc_fprintf(GRCC_Stderr,
"*** Error : %s\n\n", msg);
10617static Bool nextPart(
int nelem,
int nclist,
int *clist,
int *nl,
int *r)
10634 for (j = 0; j < nclist; j++) {
10635 nl[j] = rem/clist[j];
10636 rem -= nl[j]*clist[j];
10644 for (c = 0; c < MAXPART; c++) {
10645 rem += nl[nclist-1]*clist[nclist-1];
10646 for (pn = nclist-2; pn >= 0 && nl[pn] == 0; pn--) {
10654 for (j = pn+1; j < nclist; j++) {
10655 nl[j] = rem/clist[j];
10656 rem -= nl[j]*clist[j];
10662 erEnd(
"*** nextPart : illegal control : too many repetition");
10667static int *intdup(
int n,
int *a)
10672 for (j = 0; j < n; j++) {
10679static int *delintdup(
int *a)
10688static void prilist(
int n,
const int *a,
const char *msg)
10690 grcc_fprintf(GRCC_Stdout,
"[");
10691 for (
int j = 0; j < n; j++) {
10692 if (j!=0) grcc_fprintf(GRCC_Stdout,
", ");
10693 grcc_fprintf(GRCC_Stdout,
"%d", a[j]);
10695 grcc_fprintf(GRCC_Stdout,
"]%s", msg);
10699static int nextPerm(
int nelem,
int nclass,
int *cl,
int *r,
int *q,
int *p,
int count)
10706 for (j = 0; j < nelem; j++) {
10710 for (j = 0; j < nelem; j++) {
10715 for (k = 0; k < nclass; k++) {
10717 for (e = 0; e < n; e++) {
10724 erEnd(
"inconsistent # elements");
10730 for (j = nelem-1; j >= 0; j--) {
10732 for (k = j+1; k < nelem; k++) {
10744 for (j = 0; j < nelem; j++) {
10754static BigInt factorial(
int n)
10761 for (j = 2; j <= n; j++) {
10768static BigInt ipow(
int n,
int p)
10773 for (j = 0; j < p; j++) {
10780static int *newArray(
int size,
int val)
10787 for (j = 0; j < size; j++) {
10794static int *deleteArray(
int *a)
10803static int **newMat(
int n0,
int n1,
int val)
10808 for (j = 0; j < n0; j++) {
10809 m[j] =
new int[n1];
10810 for (k = 0; k < n1; k++) {
10818static int **deleteMat(
int **m,
int n0)
10822 for (j = n0-1; j >= 0; j--) {
10831static void bsort(
int n,
int *ord,
int *a)
10835 for (j = n-1; j > 0; j--) {
10836 for (i = 0; i < j; i++) {
10837 if(a[ord[i+1]] < a[ord[i]]) {
10847static void prMomStr(
int mom,
const char *ms,
int mn)
10851 }
else if (mom == 1) {
10852 grcc_fprintf(GRCC_Stdout,
" + %s%d", ms, mn);
10853 }
else if (mom > 0) {
10854 grcc_fprintf(GRCC_Stdout,
" + %d*%s%d", mom, ms, mn);
10855 }
else if (mom == -1) {
10856 grcc_fprintf(GRCC_Stdout,
" - %s%d", ms, mn);
10858 grcc_fprintf(GRCC_Stdout,
" - %d*%s%d", -mom, ms, mn);
10863static void prIntArray(
int n,
int *p,
const char *msg)
10868 grcc_fprintf(GRCC_Stdout,
"[");
10869 for (j = 0; j < n; j++) {
10870 if (j!=0) grcc_fprintf(GRCC_Stdout,
", ");
10871 grcc_fprintf(GRCC_Stdout,
"%2d", p[j]);
10873 grcc_fprintf(GRCC_Stdout,
"]%s", msg);
10877static void prIntArrayErr(
int n,
int *p,
const char *msg)
10882 grcc_fprintf(GRCC_Stderr,
"[");
10883 for (j = 0; j < n; j++) {
10884 if (j!=0) grcc_fprintf(GRCC_Stderr,
", ");
10885 grcc_fprintf(GRCC_Stderr,
"%2d", p[j]);
10887 grcc_fprintf(GRCC_Stderr,
"]%s", msg);
10891static void sortrec(
int l,
int r,
int *a)
10900 while(a[i] < x) i++;
10901 while(x < a[j]) j--;
10903 w = a[i]; a[i] = a[j]; a[j] = w;
10907 if (l < j) sortrec(l, j, a);
10908 if (i < r) sortrec(i, r, a);
10912static void sorti(
int size,
int *a)
10914 sortrec(0, size-1, a);
10918static int sortb(
int n,
int *a)
10920 int i, j, t, nswap;
10923 for (j = n-1; j > 0; j--) {
10924 for (i = 0; i < j; i++) {
10925 if(a[i+1] < a[i]) {
10940static Bool isIn(
int n,
int *a,
int v)
10944 for (j = 0; j < n; j++) {
10954static int intSetAdd(
int n,
int *a,
int v,
const int size)
10964 grcc_fprintf(GRCC_Stderr,
"*** intSetAdd : array out of range (>%d)\n", size);
10965 erEnd(
"intSetAdd : array out of range (GRCC_MAXPSLIST)");
10967 for (j = 0; j < n; j++) {
10970 }
else if (a[j] == v) {
10976 for (k = n-1; k >= j; k--) {
10984static int intSListAdd(
int n,
int *a,
int v,
const int size)
10994 grcc_fprintf(GRCC_Stderr,
"*** intSListAdd : array out of range (>%d)\n", size);
10995 erEnd(
"intSListAdd : array out of range");
10997 for (j = 0; j < n; j++) {
11004 for (k = n-1; k >= j; k--) {
11012static int cmpArray(
int na,
int *a,
int nb,
int *b)
11020 for (j = 0; j < na; j++) {
11030static Bool leqArray(
int n,
int *a,
int *b)
11034 for (j = 0; j < n; j++) {
11044static int toSList(
int n,
int *a)
11054static void listDiff(
int na,
int *a,
int nb,
int *b,
int *np,
int *p,
int *nq,
int *q,
int *nr,
int *r)
11058 *np = *nq = *nr = 0;
11060 while (ja < na && jb < nb) {
11061 while (ja < na && a[ja] < b[jb]) {
11062 p[(*np)++] = a[ja++];
11064 while (jb < nb && b[jb] < a[ja]) {
11065 r[(*nr)++] = b[jb++];
11067 if (ja < na && jb < nb && a[ja] == b[jb]) {
11068 q[(*nq)++] = a[ja++];
11072 for (; ja < na; ja++) {
11073 p[(*np)++] = a[ja];
11075 for (; jb < nb; jb++) {
11076 r[(*nr)++] = b[jb];
11081static int isSubSList(
int na,
int *a,
int nb,
int *b)
11086 while (ja < na && jb < nb) {
11087 for ( ; jb < nb && b[jb] < a[ja]; jb++) {
11090 if (jb >= nb || b[jb] > a[ja]) {
11093 for ( ; jb < nb && ja < na && b[jb] == a[ja]; jb++, ja++) {
11101static int subtrSet(
int na,
int *a,
int nb,
int *b,
int *c,
int size)
11111 while (ja < na && jb < nb) {
11112 while (ja < na && a[ja] < b[jb]) {
11113 if (jc == 0 || c[jc-1] != a[ja]) {
11115 erEnd(
"subsutSet: too small size of array (GRCC_MAXPSLIST)");
11121 while (jb < nb && b[jb] < a[ja]) {
11124 if (ja < na && jb < nb && a[ja] == b[jb]) {
11129 for (; ja < na; ja++) {
11130 if (jc == 0 || c[jc-1] != a[ja]) {