diff --git a/check/features.frm b/check/features.frm index a1b91032..7ed6a871 100644 --- a/check/features.frm +++ b/check/features.frm @@ -3013,3 +3013,43 @@ TableBase "no212.tbl" open, readonly; .end assert runtime_error?('Trying to open non-existent TableBase in readonly mode: no212.tbl') *--#] tablebase_ro_2 : +*--#[ diagrams_err_1 : +Model PHI3; + Particle phi,1; + Vertex phi,phi,phi:1; +EndModel; +.end +#pend_if mpi? +assert runtime_error?('Invalid coupling constant in vertex statement.') +*--#] diagrams_err_1 : +*--#[ diagrams_err_2 : +Model PHI3; + Particle phi,1; + Vertex phi,phi,phi:g^-1; +EndModel; +.end +#pend_if mpi? +assert runtime_error?('Invalid negative power of coupling constant.') +*--#] diagrams_err_2 : +*--#[ diagrams_err_3 : +Vector q1,q2,p1,p2; +Model PHI3; + Particle phi,1; + Vertex phi,phi,phi:g; +EndModel; +Local test = diagrams_(PHI3,{phi},{phi},{},{p1,p2},1,0); +.end +#pend_if mpi? +assert runtime_error?('Insufficient external momenta in diagrams_') +*--#] diagrams_err_3 : +*--#[ diagrams_err_4 : +Vector q1,q2,p1,p2; +Model PHI3; + Particle phi,1; + Vertex phi,phi,phi:g; +EndModel; +Local test = diagrams_(PHI3,{phi},{phi},{q1,q2},{},1,0); +.end +#pend_if mpi? +assert runtime_error?('Insufficient internal momenta in diagrams_') +*--#] diagrams_err_4 : diff --git a/sources/declare.h b/sources/declare.h index 4d6bfcca..591d1efe 100644 --- a/sources/declare.h +++ b/sources/declare.h @@ -560,11 +560,9 @@ extern int TestPartitions(WORD *, PARTI *); extern int DoPartitions(PHEAD WORD *,WORD); extern int CoCanonicalize(UBYTE *); extern int DoCanonicalize(PHEAD WORD *, WORD *); -extern int GenTopologies(PHEAD WORD *,WORD); extern int GenDiagrams(PHEAD WORD *,WORD); extern int DoTopologyCanonicalize(PHEAD WORD *,WORD,WORD,WORD *); extern int DoShattering(PHEAD WORD *,WORD *,WORD *,WORD); -extern int GenerateTopologies(PHEAD WORD,WORD,WORD,WORD); extern int DoTableExpansion(WORD *,WORD); extern int DoDistrib(PHEAD WORD *,WORD); extern int DoShuffle(WORD *,WORD,WORD,WORD); diff --git a/sources/diagrams.c b/sources/diagrams.c index 68354878..018fb7e2 100644 --- a/sources/diagrams.c +++ b/sources/diagrams.c @@ -328,77 +328,6 @@ int DoCanonicalize(PHEAD WORD *term, WORD *params) /* #] DoCanonicalize : - #[ GenTopologies : - - This function has the syntax - topologies_(nloops, Number of loops - nlegs, Number of legs - setvertexsizes, A set which tells which vertices are allowed like {3,4}. - set_extmomenta, The name of a set with the external momenta - set_intmomenta The name of a set with the internal momenta - [,options]) - The output will be using the built in functions vertex_ and edge_. - - The test for whether this function can be evaluated is in TestSub (inside - file proces.c) (search for the string TOPOLOGIES). - This passes the code -15 in AN.TeInFun to Generator, which then calls - the GenTopologies routine. -*/ - -#ifdef OLDTOPO - -int GenTopologies(PHEAD WORD *term,WORD level) -{ - WORD *t1, *tt1, *tstop, *t, *tt; - WORD *oldworkpointer = AT.WorkPointer; - WORD option1 = 0, option2 = 0, setoption = -1; - int retval; -/* - - We have to go through the testing procedure again, because there could - be more than one topologies_ function and not all have to be expandable. -*/ - tstop = term+*term; tstop -= ABS(tstop[-1]); - tt = term+1; - while ( tt < tstop ) { - t = tt; tt = t+t[1]; - if ( *t != TOPOLOGIES ) continue; - tt = t + t[1]; t1 = t + FUNHEAD; - if ( t1+10 > tt || *t1 != -SNUMBER || t1[1] < 0 || /* loops */ - t1[2] != -SNUMBER || ( t1[3] < 0 && t1[3] != -2 ) ||/* legs */ - t1[4] != -SETSET || Sets[t1[5]].type != CNUMBER || /* set vertices */ - t1[6] != -SETSET || Sets[t1[7]].type != CVECTOR || /* outvectors */ - t1[8] != -SETSET || Sets[t1[9]].type != CVECTOR ) continue; - tt1 = t1 + 10; - if ( tt1+2 <= tt && tt1[0] == -SETSET ) { - if ( Sets[t1[5]].last-Sets[t1[5]].first != - Sets[tt1[1]].last-Sets[tt1[1]].first ) continue; - setoption = tt1[1]; tt1 += 2; - } - if ( tt1+2 <= tt && tt1[0] == -SNUMBER ) { option1 = tt1[1]; tt1 += 2; } - if ( tt1+2 <= tt && tt1[0] == -SNUMBER ) { option2 = tt1[1]; tt1 += 2; } - AT.setinterntopo = t1[9]; - AT.setexterntopo = t1[7]; - AT.TopologiesTerm = term; - AT.TopologiesStart = t; - AT.TopologiesLevel = level; - AT.TopologiesOptions[0] = option1; - AT.TopologiesOptions[1] = option2; - retval = GenerateTopologies(BHEAD t1[1],t1[3],t1[5],setoption); - AT.WorkPointer = oldworkpointer; - return(retval); - } - MLOCK(ErrorMessageLock); - MesPrint("Internal error: topologies_ function not encountered."); - MUNLOCK(ErrorMessageLock); - return(-1); - -} - -#endif - -/* - #] GenTopologies : #[ DoTopologyCanonicalize : term: The term diff --git a/sources/diawrap.cc b/sources/diawrap.cc index e9b0800f..458c6617 100644 --- a/sources/diawrap.cc +++ b/sources/diawrap.cc @@ -237,7 +237,7 @@ void ProcessDiagram(EGraph *eg, void *ti) int i, j, intr; Model *model = (Model *)info->currentModel; MODEL *m = (MODEL *)info->currentMODEL; - int numlegs, vect, edge; + int numlegs, vect, edge, maxmom = 0; newterm = term + *term; for ( i = 1; i < info->diaoffset; i++ ) newterm[i] = term[i]; @@ -301,6 +301,8 @@ void ProcessDiagram(EGraph *eg, void *ti) } else { // Look up in set of internal momenta set *fill++ = SetElements[Sets[info->internalset].first+(vect-eg->nExtern)]; + // determine the number of momenta required from internalset: + maxmom = MaX(maxmom, vect-eg->nExtern); } *fill++ = 1; *fill++ = 1; *fill++ = 3; } @@ -335,6 +337,7 @@ void ProcessDiagram(EGraph *eg, void *ti) } else { // Look up in set of internal momenta set *fill++ = SetElements[Sets[info->internalset].first+(i-eg->nExtern)]; + maxmom = MaX(maxmom, i-eg->nExtern); } *fill++ = 1; *fill++ = 1; *fill++ = 3; // @@ -385,6 +388,7 @@ void ProcessDiagram(EGraph *eg, void *ti) } else { // Look up in set of internal momenta set *fill++ = SetElements[Sets[info->internalset].first+(vect-info->numextern)]; + maxmom = MaX(maxmom, vect-info->numextern); } } funfill[1] = fill-funfill; @@ -445,6 +449,15 @@ void ProcessDiagram(EGraph *eg, void *ti) *fill++ = SNUMBER; *fill++ = 4; *fill++ = (WORD)eg->extperm; *fill++ = 1; } // +// verify internalset has sufficient momenta: +// + if ( maxmom >= Sets[info->internalset].last - Sets[info->internalset].first ) { + MLOCK(ErrorMessageLock); + MesPrint("&Insufficient internal momenta in diagrams_"); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + } +// // finish it off // while ( tail < tend ) *fill++ = *tail++; @@ -518,7 +531,7 @@ Bool ProcessTopology(EGraph *eg, void *ti) Model *model = (Model *)info->currentModel; MODEL *m = (MODEL *)info->currentMODEL; int i, j; - int numlegs, vect, edge; + int numlegs, vect, edge, maxmom = 0; newterm = term + *term; for ( i = 1; i < info->diaoffset; i++ ) newterm[i] = term[i]; @@ -568,6 +581,8 @@ Bool ProcessTopology(EGraph *eg, void *ti) } else { // Look up in set of internal momenta set *fill++ = SetElements[Sets[info->internalset].first+(vect-info->numextern)]; + // determine the number of momenta required from internalset: + maxmom = MaX(maxmom, vect-info->numextern); } } startfill[1] = fill-startfill; @@ -590,6 +605,7 @@ Bool ProcessTopology(EGraph *eg, void *ti) } else { // Look up in set of internal momenta set *fill++ = SetElements[Sets[info->internalset].first+(i-eg->nExtern)]; + maxmom = MaX(maxmom, i-eg->nExtern); } // *fill++ = -SNUMBER; *fill++ = n1+1; // number of the node from @@ -642,6 +658,7 @@ Bool ProcessTopology(EGraph *eg, void *ti) } else { // Look up in set of internal momenta set *fill++ = SetElements[Sets[info->internalset].first+(vect-info->numextern)]; + maxmom = MaX(maxmom, vect-info->numextern); } } funfill[1] = fill-funfill; @@ -704,13 +721,13 @@ Bool ProcessTopology(EGraph *eg, void *ti) *fill++ = 0; *fill++ = 1; *fill++ = 5; } // -// Symmetry factors. We let Normalize do the multiplication. +// verify internalset has sufficient momenta: // - if ( eg->nsym != 1 ) { - *fill++ = SNUMBER; *fill++ = 4; *fill++ = (WORD)eg->nsym; *fill++ = -1; - } - if ( eg->esym != 1 ) { - *fill++ = SNUMBER; *fill++ = 4; *fill++ = (WORD)eg->esym; *fill++ = -1; + if ( maxmom >= Sets[info->internalset].last - Sets[info->internalset].first ) { + MLOCK(ErrorMessageLock); + MesPrint("&Insufficient internal momenta in diagrams_"); + MUNLOCK(ErrorMessageLock); + Terminate(-1); } // // finish it off @@ -847,6 +864,13 @@ int GenDiagrams(PHEAD WORD *term, WORD level) info.legcouple[i+ninitl] = m->vertices[numParticle(m,x)]->couplings; } info.numextern = ninitl + nfinal; + // Check that we have sufficient external momenta in the set: + if ( info.numextern > Sets[info.externalset].last - Sets[info.externalset].first ) { + MLOCK(ErrorMessageLock); + MesPrint("&Insufficient external momenta in diagrams_"); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + } for ( i = 2; i <= MAXLEGS; i++ ) { if ( m->legcouple[i] == 1 ) { for ( j = 0; j < info.numextern; j++ ) { @@ -976,114 +1000,4 @@ int processVertex(TOPOTYPE *TopoInf, int pointsremaining, int level) } // #] processVertex : -// #[ GenTopologies : - -#define TOPO_MAXVERT 10 - -int GenTopologies(PHEAD WORD *term, WORD level) -{ - Options *opt = new Options; - int nlegs, nloops, i, identical; - TERMINFO info; - WORD *t, *t1, *tstop; - TOPOTYPE TopoInf; - SETS s; -// - info.term = term; - info.level = level; - info.diaoffset = AR.funoffset; - info.flags = 0; - - t = term + info.diaoffset; // the function - t1 = t + FUNHEAD; // its arguments - tstop = t + t[1]; - - info.externalset = t1[7]; - info.internalset = t1[9]; - - s = &(Sets[t1[5]]); - TopoInf.nvert = s->last - s->first; - TopoInf.vert = &(SetElements[s->first]); - - nloops = t1[1]; - nlegs = t1[3]; - - info.numextern = nlegs; - - for ( i = 0; i <= MAXLEGS; i++ ) { TopoInf.cmind[i] = TopoInf.cmaxd[i] = 0; } - - t1 += 10; - if ( t1 < tstop && t1[0] == -SETSET ) { - TopoInf.vertmax = &(SetElements[Sets[t1[1]].first]); - t1 += 2; - } - else TopoInf.vertmax = NULL; - - info.flags |= TOPOLOGIESONLY; // this is the topologies_ function after all. - if ( t1 < tstop && t1[0] == -SNUMBER ) { - if ( ( t1[1] & WITHOUTNODES ) == WITHOUTNODES ) info.flags |= WITHOUTNODES; - if ( ( t1[1] & WITHEDGES ) == WITHEDGES ) info.flags |= WITHEDGES; - if ( ( t1[1] & WITHBLOCKS ) == WITHBLOCKS ) info.flags |= WITHBLOCKS; - if ( ( t1[1] & WITHONEPISETS ) == WITHONEPISETS ) info.flags |= WITHONEPISETS; - opt->values[GRCC_OPT_1PI] = ( t1[1] & ONEPARTI ) == ONEPARTI; -// opt->values[GRCC_OPT_NoTadpole] = ( t1[1] & NOTADPOLE ) == NOTADPOLE; - opt->values[GRCC_OPT_NoTadpole] = ( t1[1] & NOSNAIL ) == NOSNAIL; - opt->values[GRCC_OPT_No1PtBlock] = ( t1[1] & NOTADPOLE ) == NOTADPOLE; -// opt->values[GRCC_OPT_NoExtSelf] = ( t1[1] & NOEXTSELF ) == NOEXTSELF; - -// if ( ( t1[1] & WITHINSERTIONS ) == WITHINSERTIONS ) { -// opt->values[GRCC_OPT_No2PtL1PI] = True; -// opt->values[GRCC_OPT_NoAdj2PtV] = True; -// opt->values[GRCC_OPT_No2PtL1PI] = True; -// } - opt->values[GRCC_OPT_SymmInitial] = ( t1[1] & WITHSYMMETRIZEI ) == WITHSYMMETRIZEI; - opt->values[GRCC_OPT_SymmFinal] = ( t1[1] & WITHSYMMETRIZEF ) == WITHSYMMETRIZEF; - } - - info.numdia = 0; - info.numtopo = 1; - - opt->setOutAG(ProcessDiagram, &info); - opt->setOutMG(ProcessTopology, &info); -// -// Now we should sum over all possible vertices and run MGraph for -// each combination. This is done by recursion in the processVertex routine -// First load up the relevant arrays. -// - -// First the external nodes. - - if ( nlegs == -2 ) { - nlegs = 2; - identical = 1; - } - for ( i = 0; i < nlegs; i++ ) { - TopoInf.cldeg[i] = 1; TopoInf.clnum[i] = 1; TopoInf.clext[i] = -1; - } - int points = 2*nloops-2+nlegs; - - if ( identical == 1 ) { /* Only propagator topologies..... */ - nlegs = 1; - TopoInf.clnum[0] = 2; - } - TopoInf.ncl = nlegs; - TopoInf.opt = opt; - - if ( points >= MAXPOINTS ) { - MLOCK(ErrorMessageLock); - MesPrint("GenTopologies: %d loops and %d legs considered excessive",nloops,nlegs); - MUNLOCK(ErrorMessageLock); - Terminate(-1); - } - if ( processVertex(&TopoInf,points,0) != 0 ) { - MLOCK(ErrorMessageLock); - MesPrint("Called from GenTopologies with %d loops and %d legs",nloops,nlegs); - MUNLOCK(ErrorMessageLock); - Terminate(-1); - } - delete opt; - return(0); -} - -// #] GenTopologies : diff --git a/sources/gentopo.cc b/sources/gentopo.cc deleted file mode 100644 index 936e57e8..00000000 --- a/sources/gentopo.cc +++ /dev/null @@ -1,1583 +0,0 @@ -// { ( [ - -#ifdef HAVE_CONFIG_H -#ifndef CONFIG_H_INCLUDED -#define CONFIG_H_INCLUDED -#include -#endif -#endif - -#include -#include -#include -#include -#include - -#include "grcc.h" -#include "gentopo.h" - -#define DUMMYUSE(x) (void)(x) - -using namespace std; - -// The next limitation is imposed by the fact that the latest compilers -// can give warnings on declarations like "int dtcl1[nNodes]" -// With a bit of overkill there should be no real problems. -#define MAXNODES 100 -#define MAXNCLASSES 100 - -// Generate scalar connected Feynman graphs. - -//============================================================== - -typedef int Bool; -const int T_True = 1; -const int T_False = 0; - -// compile options -const int CHECK = T_True; -//const int MONITOR = T_True; -const int MONITOR = T_False; - -const int OPTPRINT = T_False; - -const int DEBUG0 = T_False; -//const int DEBUG1 = T_False; -const int DEBUG = T_False; - -// for debugging memory use -#define DEBUGM T_False - -//============================================================== -class T_MGraph; -class T_EGraph; - -#define Extern -#define Global -#define Local - -// External functions -Extern void toForm(T_EGraph *egraph); -Extern int countPhiCon(int ex, int lp, int v4); - -// Temporal definition: they should be replaced by FROM functions. -Local BigInt factorial(int n); -Local BigInt ipow(int n, int p); - -// Local functions -Local int nextPerm(int nelem, int nclass, int *cl, int *r, int *q, int *p, int count); - -Local void erEnd(const char *msg); -Local int *newArray(int size, int val); -Local void deletArray(int *a); -Local void copyArray(int size, int *a0, int *a1); -Local int **newMat(int n0, int n1, int val); -Local void deleteMat(int **m, int n0, int n1); -Local void printArray(int n, int *p); -Local void printMat(int n0, int n1, int **m); - -#if DEBUGM -static int countNMC = 0; -static int countEG = 0; -static int countMG = 0; -#endif - -//============================================================== -// class T_EGraph - -T_EGraph::T_EGraph(int nnodes, int nedges, int mxdeg) -{ - int j; - - nNodes = nnodes; - nEdges = nedges; - maxdeg = mxdeg; // maximum value of degree of nodes - nExtern = 0; - - nodes = new T_ENode[nNodes]; - edges = new T_EEdge[nEdges+1]; - for (j = 0; j < nNodes; j++) { - nodes[j].deg = 0; - nodes[j].edges = new int[maxdeg]; - } -#if DEBUGM - printf("+++ new T_EGraph %d\n", ++countEG); - if(countEG > 100) { exit(1); } -#endif -} - -T_EGraph::~T_EGraph() -{ - int j; - - for (j = 0; j < nNodes; j++) { - delete[] nodes[j].edges; - } - delete[] nodes; - delete[] edges; -#if DEBUGM - printf("+++ delete T_EGraph %d\n", countEG--); -#endif -} - -// construct T_EGraph from adjacency matrix -void T_EGraph::init(int pid, long gid, int **adjmat, Bool sopi, BigInt nsm, BigInt esm) -{ - int n0, n1, ed, e, eext, eint; -// Bool ok; - - pId = pid; - gId = gid; - opi = sopi; - nsym = nsm; - esym = esm; - - for (n0 = 0; n0 < nNodes; n0++) { - nodes[n0].deg = 0; - } - ed = 1; - for (n0 = 0; n0 < nNodes; n0++) { - for (e = 0; e < adjmat[n0][n0]/2; e++, ed++) { - edges[ed].nodes[0] = n0; - edges[ed].nodes[1] = n0; - nodes[n0].edges[nodes[n0].deg++] = - ed; - nodes[n0].edges[nodes[n0].deg++] = ed; - edges[ed].ext = nodes[n0].ext; - } - for (n1 = n0+1; n1 < nNodes; n1++) { - for (e = 0; e < adjmat[n0][n1]; e++, ed++) { - edges[ed].nodes[0] = n0; - edges[ed].nodes[1] = n1; - nodes[n0].edges[nodes[n0].deg++] = - ed; - nodes[n1].edges[nodes[n1].deg++] = ed; - edges[ed].ext = (nodes[n0].ext || nodes[n1].ext); - } - } - } - if (CHECK) { - if (ed != nEdges+1) { - printf("*** T_EGraph::init: ed=%d != nEdges=%d\n", ed, nEdges); - erEnd("*** T_EGraph::init: illegal connection\n"); - } - } - - // name of momenta - eext = 1; - eint = 1; - for (ed = 1; ed <= nEdges; ed++) { - if (edges[ed].ext) { - edges[ed].momn = eext++; - edges[ed].momc[0] = 'Q'; - edges[ed].momc[1] = ((char)0); - } else { - edges[ed].momn = eint++; - edges[ed].momc[0] = 'P'; - edges[ed].momc[1] = ((char)0); - } - } -} - -// set external particle to node 'nd' -void T_EGraph::setExtern(int nd, Bool val) -{ - nodes[nd].ext = val; -} - -// end of calling 'setExtern' -void T_EGraph::endSetExtern(void) -{ - int n; - - nExtern = 0; - for (n = 0; n < nNodes; n++) { - if (nodes[n].ext) { - nExtern++; - } - } -} - -// print the T_EGraph -void T_EGraph::print() -{ - int nd, lg, ed, nlp; - - nlp = nEdges - nNodes + 1; - printf("\n"); - printf("T_EGraph: pId=%d gId=%ld nExtern=%d nLoops=%d nNodes=%d nEdges=%d\n", - pId, gId, nExtern, nlp, nNodes, nEdges); - printf(" sym = (%ld * %ld) maxdeg=%d\n", nsym, esym, maxdeg); - printf(" Nodes\n"); - for (nd = 0; nd < nNodes; nd++) { - if (nodes[nd].ext) { - printf(" %2d Extern ", nd); - } else { - printf(" %2d Vertex ", nd); - } - printf("deg=%d [", nodes[nd].deg); - for (lg = 0; lg < nodes[nd].deg; lg++) { - printf(" %3d", nodes[nd].edges[lg]); - } - printf("]\n"); - } - printf(" Edges\n"); - for (ed = 1; ed <= nEdges; ed++) { - if (edges[ed].ext) { - printf(" %2d Extern ", ed); - } else { - printf(" %2d Intern ", ed); - } - printf("%s%d", (edges[ed].ext)?"Q":"p", edges[ed].momn); - printf(" [%3d %3d]\n", edges[ed].nodes[1], edges[ed].nodes[0]); - } -} - -//============================================================== -// class T_MNode : nodes in T_MGraph - -T_MNode::T_MNode(int vid, int vdeg, Bool vext, int vclss) -{ - id = vid; // id of the node - deg = vdeg; // degree of the node - freelg = vdeg; // number of free legs - clss = vclss; // class to which the node belongs - ext = vext; // external node or not -} - -//=============================================================== -// class of node-classes for T_MGraph -// -class T_MNodeClass { - public: - int nNodes; // the number of nodes - int nClasses; // the number of classes - int *clist; // the number of nodes in each class - int *ndcl; // node --> class - int **clmat; // matrix used for classification - int *flist; // the first node in each class - int maxdeg; // maximal value of degree(node) - int forallignment; - - T_MNodeClass(int nnodes, int ncl); - ~T_MNodeClass(); - void init(int *cl, int mxdeg, int **adjmat); - void copy(T_MNodeClass* mnc); - int clCmp(int nd0, int nd1, int cn); - void printMat(void); - - void mkFlist(void); - void mkNdCl(void); - void mkClMat(int **adjmat); - void incMat(int nd, int td, int val); - Bool chkOrd(int nd, int ndc, T_MNodeClass *cl, int *dtcl); - int cmpArray(int *a0, int *a1, int ma); -}; - -T_MNodeClass::T_MNodeClass(int nnodes, int ncl) -{ - nNodes = nnodes; - nClasses = ncl; - clist = new int[nClasses]; - ndcl = new int[nNodes]; - clmat = newMat(nNodes, nClasses, 0); - flist = new int[nClasses+1]; - -#if DEBUGM - printf("+++ new T_MNodeClass %d\n", ++countNMC); - if(countNMC > 100) { exit(1); } -#endif -} - -T_MNodeClass::~T_MNodeClass() -{ - delete[] clist; - delete[] ndcl; - deleteMat(clmat, nNodes, nClasses); - delete[] flist; - -#if DEBUGM - printf("+++ delete T_MNodeClass %d\n", countNMC--); -#endif -} - -void T_MNodeClass::init(int *cl, int mxdeg, int **adjmat) -{ - int j; - - for (j = 0; j < nClasses; j++) { - clist[j] = cl[j]; - } - maxdeg = mxdeg; - mkNdCl(); - mkClMat(adjmat); - mkFlist(); -} - -void T_MNodeClass::copy(T_MNodeClass* mnc) -{ - int j, k; - - for (k = 0; k < nClasses; k++) { - clist[k] = mnc->clist[k]; - } - maxdeg = mnc->maxdeg; - for (j = 0; j < nNodes; j++) { - ndcl[j] = mnc->ndcl[j]; - for (k = 0; k < nClasses; k++) { - clmat[j][k] = mnc->clmat[j][k]; - } - } - for (k = 0; k < nClasses+1; k++) { - flist[k] = mnc->flist[k]; - } -} - -// Construct flist -// The set of nodes in class 'cl' is [flist[cl],...,flist[cl+1]-1] -void T_MNodeClass::mkFlist(void) -{ - int j, f; - - f = 0; - for (j = 0; j < nClasses; j++) { - flist[j] = f; - f += clist[j]; - } - flist[nClasses] = f; -} - -// Construct ndcl -// ndcl[nd] = (the class id in which node 'nd' belongs) -void T_MNodeClass::mkNdCl(void) -{ - int c, k; - int nd = 0; - - for (c = 0; c < nClasses; c++) { - for (k = 0; k < clist[c]; k++) { - ndcl[nd++] = c; - } - } -} - -// Comparison of two nodes 'nd0' and 'nd1' -// Ordering is lexicographic (class, connection configuration) -int T_MNodeClass::clCmp(int nd0, int nd1, int cn) -{ - - // Whether two nodes are in a same class or not. - int cmp = ndcl[nd0] - ndcl[nd1]; - if (cmp != 0) { - return cmp; - } - - // Sign '-' signifies the reverse ordering - cmp = - cmpArray(clmat[nd0], clmat[nd1], cn); - if (cmp != 0) { - return cmp; - } - // for particles ??? - return cmp; -} - -// Construct a matrix 'clmat[nd][tc]' which is the number -// of edges connecting 'nd' and all nodes in class 'tc'. -void T_MNodeClass::mkClMat(int **adjmat) -{ - int nd, td, tc; - - for (nd = 0; nd < nNodes; nd++) { - for (td = 0; td < nNodes; td++) { - tc = ndcl[td]; // another node - clmat[nd][tc] += adjmat[nd][td]; // the number of edges 'nd'--'td' - } - - } -} - -// Increase the number of edges by 'val' between 'nd' and 'td'. -void T_MNodeClass::incMat(int nd, int td, int val) -{ - int tdc = ndcl[td]; // class of 'td' - clmat[nd][tdc] += val; // modify matrix 'clmat'. -} - -// Check whether the configuration satisfies the ordering condition or not. -Bool T_MNodeClass::chkOrd(int nd, int ndc, T_MNodeClass *cl, int *dtcl) -{ - Bool tcl[MAXNCLASSES]; -// Bool tcl[cl->nClasses]; - int tn, tc, mxn, cmp, n; - DUMMYUSE(nd); - - for (tc = 0; tc < cl->nClasses; tc++) { - tcl[tc] = T_False; - } - for (tn = 0; tn < nNodes; tn++) { - if (dtcl[tn] == 0) { - tcl[cl->ndcl[tn]] = T_False; - } - } - for (tc = 0; tc < cl->nClasses; tc++) { - if (tcl[tc] && ndc != tc) { - mxn = flist[tc+1]; - for (n = flist[tc]+1; n < mxn; n++) { - cmp = - clmat[n-1][tc] + clmat[n][tc]; - if (cmp > 0) { - return T_False; - } - } - } - } - return T_True; -} - -// Print configuration matrix. -void T_MNodeClass::printMat(void) -{ - int j1, j2; - - cout << endl; - - // the first line - cout << setw(2) << "nd" << ": " << setw(2) << "cl: "; - for (j2 = 0; j2 < nClasses; j2++) { - cout << setw(2) << j2 << " "; - } - cout << endl; - - // print raw - for (j1 = 0; j1 < nNodes; j1++) { - cout << setw(2) << j1 << ": " << setw(2) << ndcl[j1] << ": ["; - for (j2 = 0; j2 < nClasses; j2++) { - cout << " " << setw(2) << clmat[j1][j2]; - } - cout << "] " << endl; - } -} - -int T_MNodeClass::cmpArray(int *a0, int *a1, int ma) -{ - for (int j = 0; j < ma; j++) { - if (a0[j] < a1[j]) { - return -1; - } else if (a0[j] > a1[j]) { - return 1; - } - } - return 0; -} - -//=============================================================== -// class T_MGraph : scalar graph expressed by matrix form - -T_MGraph::T_MGraph(int pid, int ncl, int *cldeg, int *clnum, int *clext, Bool sopi) -{ - int nn, ne, j, k; - - // initial conditions - nClasses = ncl; - clist = new int[ncl]; - - mindeg = -1; - maxdeg = -1; - ne = 0; - nn = 0; - for (j = 0; j < nClasses; j++) { - clist[j] = clnum[j]; - nn += clnum[j]; - ne += cldeg[j]*clnum[j]; - if (mindeg < 0) { - mindeg = cldeg[j]; - } else { - mindeg = min(mindeg, cldeg[j]); - } - if (maxdeg < 0) { - maxdeg = cldeg[j]; - } else { - maxdeg = max(maxdeg, cldeg[j]); - } - } - if (ne % 2 != 0) { - printf("Sum of degrees are not even\n"); - for (j = 0; j < nClasses; j++) { - printf("class %2d: %2d %2d %2d\n", - j, cldeg[j], clnum[j], clext[j]); - } - erEnd("illegal degrees of nodes"); - } - pId = pid; - nNodes = nn; - nodes = new T_MNode*[nNodes]; - - selOPI = sopi; - - nEdges = ne / 2; - nLoops = nEdges - nNodes + 1; - - egraph = new T_EGraph(nNodes, nEdges, maxdeg); - nn = 0; - for (j = 0; j < nClasses; j++) { - for (k = 0; k < clist[j]; k++, nn++) { - nodes[nn] = new T_MNode(nn, cldeg[j], clext[j], j); - egraph->setExtern(nn, clext[j]); - } - } - egraph->endSetExtern(); - - // generated set of graphs - ndiag = 0; - n1PI = 0; - - // the current graph - adjMat = newMat(nNodes, nNodes, 0); - nsym = ToBigInt(0); - esym = ToBigInt(0); - c1PI = 0; - wsum = ToFraction(0, 1); - wsopi = ToFraction(0, 1); - - // current node classification - curcl = new T_MNodeClass(nNodes, nClasses); - - // measures of efficiency - ngen = 0; - ngconn = 0; - - if (MONITOR) { - nCallRefine = 0; - discardOrd = 0; - discardRefine = 0; - discardDisc = 0; - discardIso = 0; - } - - // work space for isIsomorphic - modmat = newMat(nNodes, nNodes, 0); - permp = newArray(nNodes, 0); - permq = newArray(nNodes, 0); - permr = newArray(nNodes, 0); - - // work space for bisearchM - bidef = newArray(nNodes, 0); - bilow = newArray(nNodes, 0); - bicount = 0; - DUMMYUSE(padding); - -#if DEBUGM - printf("+++ new T_MGraph %d\n", ++countMG); - if(countEG > 100) { exit(1); } -#endif -} - -T_MGraph::~T_MGraph() -{ - int j; - - deletArray(bilow); - deletArray(bidef); - deletArray(permr); - deletArray(permq); - deletArray(permp); - - deleteMat(modmat, nNodes, nNodes); - delete curcl; - deleteMat(adjMat, nNodes, nNodes); - delete egraph; - for (j = 0; j < nNodes; j++) { - delete nodes[j]; - } - delete[] nodes; - delete[] clist; - -#if DEBUGM - printf("+++ delete T_MGraph %d\n", countMG++); -#endif -} - - -void T_MGraph::printAdjMat(T_MNodeClass *cl) -{ - int j1, j2; - - cout << " "; - for (j2 = 0; j2 < nNodes; j2++) { - cout << " " << setw(2) << j2; - } - cout << endl; - for (j1 = 0; j1 < nNodes; j1++) { - cout << setw(2) << j1 << ": ["; - for (j2 = 0; j2 < nNodes; j2++) { - cout << " " << setw(2) << adjMat[j1][j2]; - } - cout << "] " << cl->ndcl[j1] << endl; - } -} - -//--------------------------------------------------------------- -// Check graph can be a connected one. -// If a connected component without free leg is not the whole graph then -// return T_False, otherwise return T_True. -Bool T_MGraph::isConnected(void) -{ - int j, n, nv; - - for (j = 0; j < nNodes; j++) { - nodes[j]->visited = -1; - } - if (visit(0)) { - return T_True; - } - nv = 0; - for (n = 0; n < nNodes; n++) { - if (nodes[n]->visited >= 0) { - nv++; - } - } - return (nv == nNodes); -} - -// Visiting connected node used for 'isConnected' -// If child nodes has free legs, then this function returns T_True. -// otherwise it returns T_False. -Bool T_MGraph::visit(int nd) -{ - int td; - - // This node has free legs. - if (nodes[nd]->freelg > 0) { - return T_True; - } - nodes[nd]->visited = 0; - for (td = 0; td < nNodes; td++) { - if ((adjMat[nd][td] > 0) and (nodes[td]->visited < 0)) { - if (visit(td)) { - return T_True; - } - } - } - // all the child nodes has no free legs. - return T_False; -} - -// Check whether the current graph is the Representative of a isomorphic class. -// nsym = symmetry factor by the permutation of nodes. -// esym = symmetry factor by the permutation of edge. -// If this graph is not a representative, then returns T_False. -Bool T_MGraph::isIsomorphic(T_MNodeClass *cl) -{ - int j1, j2, cmp, count, nself; - - nsym = ToBigInt(0); - esym = ToBigInt(1); - - count = 0; - while (T_True) { - count = nextPerm(nNodes, cl->nClasses, cl->clist, - permr, permq, permp, count); - - if (count < 0) { - // calculate permutations of edges - esym = ToBigInt(1); - for (j1 = 0; j1 < nNodes; j1++) { - if (adjMat[j1][j1] > 0) { - nself = adjMat[j1][j1]/2; - esym *= factorial(nself); - esym *= ipow(2, nself); - } - for (j2 = j1+1; j2 < nNodes; j2++) { - if (adjMat[j1][j2] > 0) { - esym *= factorial(adjMat[j1][j2]); - } - } - } - return T_True; - } - - permMat(nNodes, permp, adjMat, modmat); - cmp = compMat(nNodes, adjMat, modmat); - if (cmp < 0) { - return T_False; - } else if (cmp == 0) { - // save permutation ??? - nsym = nsym + ToBigInt(1); - } - } -} - -// apply permutation to matrix 'mat0' and obtain 'mat1' -void T_MGraph::permMat(int size, int *perm, int **mat0, int **mat1) -{ - int j1, j2; - - for (j1 = 0; j1 < size; j1++) { - for (j2 = 0; j2 < size; j2++) { - mat1[j1][j2] = mat0[perm[j1]][perm[j2]]; - } - } -} - -// comparison of matrix -int T_MGraph::compMat(int size, int **mat0, int **mat1) -{ - int j1, j2, cmp; - - for (j1 = 0; j1 < size; j1++) { - for (j2 = 0; j2 < size; j2++) { - cmp = mat0[j1][j2] - mat1[j1][j2]; - if (cmp != 0) { - return cmp; - } - } - } - return 0; -} - - -// Refine the classification -// cl : the current 'T_MNodeClass' object -// cn : the class number -// Returns (the new class number corresponds to 'cn') if OK, or -// 'None' if ordering condition is not satisfied -T_MNodeClass *T_MGraph::refineClass(T_MNodeClass *cl, int cn) -{ - T_MNodeClass *ccl = cl; - int ccn = cn; - T_MNodeClass *ncl = NULL; - int ucl[MAXNODES]; - int nucl, nce; - int td, cmp; - - nucl = 0; - while (ccl->nClasses != nucl) { // repeat refinement. - nce = 0; - nucl = 0; - for (td = 1; td < nNodes; td++) { - // 'td' is the next node and the current node is 'td-1'. - // Count up the number of the elements in the current class - // corresponding to the current node - nce++; - cmp = ccl->clCmp(td-1, td, ccn); - - // the ordering condition is not satisfied. - if (cmp > 0) { - if (DEBUG) { - cout << "refine: cls = " << cl->clist << endl; - cl->printMat(); - cout << "clmat" << endl; - ccl->printMat(); - cout << "refine: discard: cls = " << ccl->clist - << "ucl = " << ucl << endl; - } - if (ccl != cl) { - delete ccl; - } - return NULL; - - } else if (cmp < 0) { - // 'td' is in the next class to the current node. - ucl[nucl++] = nce; // close the current class - - // start new class - nce = 0; - } - // nothing to do for the case of 'cmp == 0'. - - } - - // close array 'ucl'. - ucl[nucl++] = nce + 1; - - // class is not modified - if (nucl == ccl->nClasses) { - ncl = ccl; - break; - - // inconsistent - } else if (nucl < ccl->nClasses) { - erEnd("refineClasses : smaller number of classes"); - - // preparation of the next repetition. - } else { - if (cn == ccl->nClasses) { - td = nNodes; - } else { - td = ccl->flist[cn+1]-1; - } - if (ccl != cl) { - delete ccl; - } - ccl = new T_MNodeClass(nNodes, nucl); - ccl->init(ucl, maxdeg, adjMat); - if (td == nNodes) { - ccn = ccl->nClasses; - } else { - ccn = ccl->ndcl[td]; - } - nucl = 0; - } - } - - if (DEBUG) { - cout << "refine: ucl = " << ucl << endl; - cout << "refine: ncl = " << ncl->clist << endl; - ncl->printMat(); - } - - return ncl; -} - -// Search biconnected component -// visit : pd --> nd --> td -// ne : the number of edges between pd and nd. -void T_MGraph::bisearchM(int nd, int pd, int ne) -{ - int td; - - bidef[nd] = bicount; - bilow[nd] = bicount; - bicount++; - - for (td = 0; td < nNodes; td++) { - if (nodes[td]->ext) { // ignore external node - continue; - - } else if (td == nd) { // pd --> nd --> nd - continue; - - } else if (td == pd) { // pd --> nd --> pd - if (ne > 1) { - bilow[nd] = min(bilow[nd], bidef[pd]); - } - - } else if (adjMat[td][nd] < 1) { // td is not adjacent to nd - continue; - - } else if (bidef[td] >= 0) { // back edge - bilow[nd] = min(bilow[nd], bidef[td]); - - // new node - } else { - - bisearchM(td, nd, adjMat[td][nd]); - - // ordinary case - if (bilow[td] > bidef[nd]) { - nBridges++; - } - bilow[nd] = min(bilow[nd], bilow[td]); - } - } - - // nd is the starting point and not an external line - // if (pd < 0 && (!nodes[nd]->ext || nodes[nd]->deg > 1)) { - // if (pd < 0 && !(nodes[nd]->ext && nodes[nd]->deg ==1)) { - if (pd < 0 && nodes[nd]->deg != 1) { - // nBridges += nodes[nd]->deg; - nBridges++; - } -} - -// Count the number of 1PI components. -// The algorithm is described in -// A.V. Aho, J.E. Hopcroft and J.D. Ullman -// 'The Design and Analysis of Computer Algorithms', Chap. 5 -// 1974, Addison-Wesley. - -int T_MGraph::count1PI(void) -{ - int j; - - if (nLoops < 0) { - return 1; - } - - // initialization - bicount = 0; - nBridges = 0; - for (j = 0; j < nNodes; j++) { - bidef[j] = -1; - bilow[j] = -1; - } - - bisearchM(0, -1, 0); - - return nBridges; -} - -//--------------------------------------------------------------- -// Generate graphs -// the generation process starts from 'connectClass'. - -long T_MGraph::generate(void) -{ - T_MNodeClass *cl; - int dscl[MAXNODES]; - int n; - - for (n = 0; n < nNodes; n++) { - dscl[n] = T_False; - } - - // Initial classification of nodes. - cl = new T_MNodeClass(nNodes, nClasses); - cl->init(clist, maxdeg, adjMat); - connectClass(cl, dscl); - - // Print the result. - - if (MONITOR) { - cout << endl; - cout << endl; - cout << "* Total " << ndiag << " Graphs."; - cout << "(" << n1PI << " 1PI)"; - // cout << " wsum = " << wsum << "(" << wsopi << "1PI)" << endl; - cout << endl; - } - - if (MONITOR) { - cout << "* refine: " << nCallRefine << endl; - cout << "* discard for ordering: " << discardOrd << endl; - cout << "* discard for refinement: " << discardRefine << endl; - cout << "* discard for disconnected: " << discardDisc << endl; - cout << "* discard for duplication: " << discardIso << endl; - } - delete cl; - - return ndiag; -} - -int T_MGraph::findNextCl(T_MNodeClass *cl, int *dscl) -{ - int mine, cr, c, n, me; - - mine = -1; - cr = -1; - for (c = 0; c < cl->nClasses; c++) { - n = cl->flist[c]; - if (!dscl[n]) { - me = nodes[n]->freelg; - if (me > 0) { - if (nodes[n]->freelg < nodes[n]->deg) { - return c; - } - if (mine < 0 || mine > me) { - mine = me; - cr = c; - } - } - } - } - return cr; -} - -int T_MGraph::findNextTCl(T_MNodeClass *cl, int *dtcl) -{ - int c, n, mine, cr, me; - - mine = -1; - cr = -1; - for (c = 0; c < cl->nClasses; c++) { - for (n = cl->flist[c]; n < cl->flist[c+1]; n++) { - if (!dtcl[n]) { - me = nodes[n]->freelg; - if (me > 0) { - if (nodes[n]->freelg < nodes[n]->deg) { - return c; - } - if (mine < 0 || mine > me) { - mine = me; - cr = c; - } - } - } - } - } - return cr; -} - -// Connect nodes in a class to others -void T_MGraph::connectClass(T_MNodeClass *cl, int *dscl) -{ - int sc, sn; - T_MNodeClass *xcl; - - if (DEBUG0) { - printf("connectClass:begin:"); - printf(" dscl="); printArray(nNodes, dscl); - printf("\n"); - } - xcl = refineClass(cl, cl->nClasses); - - if (xcl == NULL) { - if (MONITOR) { - discardRefine++; - } - } else { - sc = findNextCl(xcl, dscl); - if (sc < 0) { - newGraph(cl); - } else { - sn = xcl->flist[sc]; - connectNode(sc, sn, xcl, dscl); - } - } - if (xcl != cl && xcl != NULL) { - delete xcl; - } - if (DEBUG0) { - printf("connectClass:end\n"); - } -} - -//------------------------------------------------------------ -void T_MGraph::connectNode(int sc, int ss, T_MNodeClass *cl, int *dscl) -{ - int sn; - int dtcl[MAXNODES]; - - if (DEBUG0) { - printf("connectNode:begin:(%d,%d)", sc, ss); - printf(" dscl="); printArray(nNodes, dscl); - printf("\n"); - } - if (ss >= cl->flist[sc+1]) { - connectClass(cl, dscl); - if (DEBUG0) { - printf("connectNode:end1\n"); - } - return; - } - - copyArray(nNodes, dscl, dtcl); - for (sn = ss; sn < cl->flist[sc+1]; sn++) { - if (!dscl[sn]) { - connectLeg(sc, sn, sc, sn, cl, dscl, dtcl); - if (DEBUG0) { - printf("connectNode:end2\n"); - } - return; - } - } - if (DEBUG0) { - printf("connectNode:end3\n"); - } -} - -// Add one connection between two legs. -// 1. select another node to connect -// 2. determine multiplicity of the connection -// -// Arguments -// cn : the current class -// nd : the current node to be connected. -// nextnd : {nextnd, ...} is the possible target node of the connection. -// cl : the current node class - -void T_MGraph::connectLeg(int sc, int sn, int tc, int ts, T_MNodeClass *cl, int *dscl, int* dtcl) -{ - int tn, maxself, nc2, nc, maxcon, ts1, wc, ncm; - int dtcl1[MAXNODES]; - - if (DEBUG0) { - printf("connectLeg:begin:(%d,%d,%d,%d)", sc, sn, tc, ts); - printf(" dscl="); printArray(nNodes, dscl); - printf(" dtcl="); printArray(nNodes, dtcl); - printf("\n"); - printAdjMat(cl); - } - - if (sn >= cl->flist[sc+1]) { - erEnd("*** connectLeg : illegal control"); - return; - - // There remains no free legs in the node 'sn' : move to next node. - } else if (nodes[sn]->freelg < 1) { - - if (!isConnected()) { - if (DEBUG0) { - printf("connectLeg:disconnected\n"); - } - if (MONITOR) { - discardDisc++; - } - } else { - if (DEBUG0) { - printf("connectLeg: call conNode\n"); - } - dscl[sn] = T_True; - - // next node in the current class. - connectNode(sc, sn+1, cl, dscl); - - dscl[sn] = T_False; - } - - // connect a free leg of the current node 'sn'. - } else { - copyArray(nNodes, dtcl, dtcl1); - - ts1 = ts; - for (wc = 0; wc < nNodes; wc++) { - if (ts1 >= cl->flist[tc+1]) { - tc = findNextTCl(cl, dtcl1); - if (DEBUG0) { - printf("connectLeg:1:tc=%d\n", tc); - } - if (tc < 0) { - if (DEBUG0) { - printf("connectLeg:end1:(%d,%d,%d,%d)\n", sc, sn, tc, ts); - } - return; - } - if (tc != sc && !cl->chkOrd(sn, sc, cl, dtcl)) { - if (MONITOR) { - discardOrd++; - } - if (DEBUG0) { - printf("connectLeg:end2:(%d,%d,%d,%d)\n", sc, sn, tc, ts); - } - return; - } - } - - ts1 = -1; - - // repeat for all possible target nodes - // for all tn in tc - if (DEBUG0) { - printf("connectLeg:2:tc=%d, fl:%d-->%d\n", tc, cl->flist[tc], cl->flist[tc+1]); - } - for (tn = cl->flist[tc]; tn < cl->flist[tc+1]; tn++) { - if (DEBUG0) { - printf("connectLeg:2:%d=>try %d\n", sn, tn); - } - - if (dtcl1[tn]) { - if (DEBUG0) { - printf("connectLeg:3\n"); - } - continue; - } - dtcl1[tn] = T_True; - - if (sc == tc && sn > tn) { - if (DEBUG0) { - printf("connectLeg:4\n"); - } - continue; - - // self-loop - } else if (sn == tn) { - if (nNodes > 1) { - // there are two or more nodes in the graph : - // avoid disconnected graph - maxself = min((nodes[sn]->freelg)/2, (nodes[sn]->deg-1)/2); - } else { - // there is only one node the graph. - maxself = nodes[sn]->freelg/2; - } - - // If we can assume no tadpole, the following line can be used. - if (selOPI and nNodes > 2) { - maxself = min((nodes[sn]->deg-2)/2, maxself); - } - - // vary the number of connection. - for (nc2 = maxself; nc2 > 0; nc2--) { - // if nNodes > 1 and ndg == nc2: - // // disconnected - // continue - nc = 2*nc2; - if (DEBUG0) { - printf("connectLeg: call conLeg: (same) %d=>%d(%d)\n", sn, tn,nc); - } - ncm = 5*nc; - adjMat[sn][sn] = nc; - nodes[sn]->freelg -= nc; - cl->incMat(sn, tn, ncm); - ts1 = tn + 1; - - // next connection - connectLeg(sc, sn, tc, ts1, cl, dscl, dtcl1); - - // restore the configuration - cl->incMat(tn, sn, - ncm); - adjMat[sn][sn] = 0; - nodes[sn]->freelg += nc; - if (DEBUG0) { - printf("connectLeg: ret conLeg: (same) %d=>%d(%d)\n", sn, tn,nc); - } - } - - // connections between different nodes. - } else { - // maximum possible connection number - maxcon = min(nodes[sn]->freelg, nodes[tn]->freelg); - - // avoid disconnected graphs. - if (nNodes > 2 && nodes[sn]->deg == nodes[tn]->deg) { - maxcon = min(maxcon, nodes[sn]->deg-1); - } - - if (CHECK) { - if ((adjMat[sn][tn] != 0) || (adjMat[sn][tn] != adjMat[tn][sn])) { - printf("*** inconsistent connection: sn=%d, tn=%d", sn, tn); - printf(" dscl="); printArray(nNodes, dscl); - printf(" dtcl="); printArray(nNodes, dtcl); - printf("\n"); - printAdjMat(cl); - erEnd("*** inconsistent connection "); - } - } - - // vary number of connections - for (nc = maxcon; nc > 0; nc--) { - if (DEBUG0) { - printf("connectLeg: call conLeg: (diff) %d=>%d(%d)", sn, tn,nc); - printf(" dtcl="); printArray(nNodes, dtcl); - printf("\n"); - } - adjMat[sn][tn] = nc; - adjMat[tn][sn] = nc; - nodes[sn]->freelg -= nc; - nodes[tn]->freelg -= nc; - cl->incMat(sn, tn, nc); - cl->incMat(tn, sn, nc); - ts1 = tn + 1; - - // next connection - connectLeg(sc, sn, tc, ts1, cl, dscl, dtcl1); - - // restore configuration - cl->incMat(sn, tn, - nc); - cl->incMat(tn, sn, - nc); - adjMat[sn][tn] = 0; - adjMat[tn][sn] = 0; - nodes[sn]->freelg += nc; - nodes[tn]->freelg += nc; - if (DEBUG0) { - printf("connectLeg: ret conLeg: (diff) %d=>%d(%d)\n", sn, tn,nc); - printf(" dtcl="); printArray(nNodes, dtcl); - printAdjMat(cl); - printf("\n"); - } - } - } - } - if (ts1 < 0) { - ts1 = cl->flist[tc+1]; - } - } // for wc - } - if (DEBUG0) { - printf("connectLeg:end3:(%d,%d,%d,%d)\n", sc, sn, tc, ts); - } -} - -// A new candidate diagram is obtained. -// It may be a disconnected graph - -void T_MGraph::newGraph(T_MNodeClass *cl) -{ - int connected; - T_MNodeClass *xcl; - Bool sopi; - - ngen++; - // printf("newGraph: %d\n", ngen); - - // refine class and check ordering condition - xcl = refineClass(cl, cl->nClasses); - if (xcl == NULL) { - // printf("newGraph: fail refine\n"); - if (MONITOR) { - discardRefine++; - } - - // check whether this is a connected graph or not. - } else { - connected = isConnected(); - if (!connected) { - // printf("newGraph: not connected\n"); - if (DEBUG0) { - cout << "+++ disconnected graph" << ngen << endl; - xcl->printMat(); - } - - if (MONITOR) { - discardDisc++; - } - - // connected graph : check isomorphism of the graph - } else { - ngconn++; - if (!isIsomorphic(xcl)) { - // printf("newGraph: not isomorphic\n"); - if (DEBUG) { - cout << "+++ duplicated graph" << ngen << ngconn << endl; - xcl->printMat(); - } - - if (MONITOR) { - discardIso++; - } - - // We got a new connected and a unique representation of a class. - } else { - - c1PI = count1PI(); - if (!selOPI || c1PI == 1) { - wsum = wsum + ToFraction(1, nsym*esym); - if (c1PI == 1) { - wsopi = wsopi + ToFraction(1, nsym*esym); - } - curcl->copy(xcl); - ndiag++; - sopi = (c1PI == 1); - if (sopi) { - n1PI++; - } - - if (OPTPRINT) { - cout << endl; - cout << "Graph :" << ndiag - << "(" << ngen << ")" - << " 1PI comp. = " << c1PI - << " sym. factor = (" << nsym << "*" << esym << ")" - << endl; - printAdjMat(cl); - // cl->printMat(); - if (MONITOR) { - cout << "refine: " << nCallRefine << endl; - cout << "discard for ordering: " << discardOrd << endl; - cout << "discard for refinement: " << discardRefine << endl; - cout << "discard for disconnected: " << discardDisc << endl; - cout << "discard for duplication: " << discardIso << endl; - } - } - - // go to next step - egraph->init(pId, ndiag, adjMat, sopi, nsym, esym); - toForm(egraph); - } - } - } - } -// printf("in newGraph cl=%p, xcl=%p\n", cl, xcl); - if (xcl != cl && xcl != NULL) { - delete xcl; - } -} - - -// go to next step -// 1. convert data format which treat the edges as objects. -// 2. definition of loop momenta to the edges. -// 3. analysis of loop structure in a graph. -// 4. assignment of particles to edges and vertices to nodes -// 5. recalculate symmetry factor - -// Permutation -Local int nextPerm(int nelem, int nclass, int *cl, int *r, int *q, int *p, int count) -{ - int j, k, n, e, t; - Bool b; - - for (j = 0; j < nelem; j++) { - p[j] = j; - } - if (count < 1) { - for (j = 0; j < nelem; j++) { - q[j] = 0; - r[j] = 0; - } - j = 0; - for (k = 0; k < nclass; k++) { - n = cl[k]; - for (e = 0; e < n; e++) { - r[j] = n - e - 1; - j++; - } - } - if (j != nelem) { - erEnd("*** inconsistent # elements"); - } - return 1; - } - b = T_False; - for (j = nelem-1; j >= 0; j--) { - if (q[j] < r[j]) { - for (k = j+1; k < nelem; k++) { - q[k] = 0; - } - q[j]++; - b = T_True; - break; - } - } - if (!b) { - return (-count); - } - - for (j = 0; j < nelem; j++) { - k = j + q[j]; - t = p[j]; - p[j] = p[k]; - p[k] = t; - } - return count + 1; -} - -Local BigInt factorial(int n) -/* return (n < 1) ? 1 : n! - */ -{ - int r, j; - - r = 1; - for (j = 2; j <= n; j++) { - r *= j; - } - return r; -} - -Local BigInt ipow(int n, int p) -{ - int r, j; - DUMMYUSE(p); - - r = 1; - for (j = 2; j <= n; j++) { - r *= n; - } - return r; -} - - -Local void erEnd(const char *msg) -{ - printf("*** Error : %s\n", msg); - exit(1); -} - -/* memory allocation */ -Local int *newArray(int size, int val) -{ - int *a, j; - - a = new int[size]; - for (j = 0; j < size; j++) { - a[j] = val; - } - return a; -} - -Local void deletArray(int *a) -{ - delete[] a; -} - -Local void copyArray(int size, int *a0, int *a1) -{ - int j; - for (j = 0; j < size; j++) { - a1[j] = a0[j]; - } -} - -Local int **newMat(int n0, int n1, int val) -{ - int **m, j; - - m = new int*[n0]; - for (j = 0; j < n0; j++) { - m[j] = newArray(n1, val); - } - return m; -} - -Local void deleteMat(int **m, int n0, int n1) -{ - int j; - DUMMYUSE(n1); - - for (j = 0; j < n0; j++) { - deletArray(m[j]); - } - delete[] m; -} - -//============================================================== -// Functions for testing the program. -// -//-------------------------------------------------------------- -// Testing function nextPerm -// -Local void printArray(int n, int *p) -{ - int j; - - printf("["); - for (j = 0; j < n; j++) { - printf(" %2d", p[j]); - } - printf("]"); -} - -Local void printMat(int n0, int n1, int **m) -{ - int j; - - for (j = 0; j < n0; j++) { - printArray(n1, m[j]); - printf("\n"); - } -} - -Global void testPerm() -{ - int nelem, nperm, nclist, n, count; - int *p, *q, *r; - int clist[] = {1, 2, 2, 3}; - - nclist = sizeof(clist)/sizeof(int); - nelem = 0; - nperm = 1; - for (n = 0; n < nclist; n++) { - nelem += clist[n]; - nperm *= factorial(clist[n]); - } - - printf("+++ clist = (%d) ", nclist); - printArray(nclist, clist); - printf("\n"); - printf("+++ nelem = %d, nperm = %d\n", nelem, nperm); - - p = new int[nelem]; - q = new int[nelem]; - r = new int[nelem]; - count = 0; - while (T_True) { - count = nextPerm(nelem, nclist, clist, r, q, p, count); - if (count < 0) { - count = - count; - break; - } - printf("%4d:", count); - printArray(nelem, p); - printf("\n"); - - if (count > nperm) { - break; - } - } - if (count != nperm) { - printf("*** %d != %d\n", count, nperm); - } - delete[] p; - delete[] q; - delete[] r; -} - -// } ) ] - diff --git a/sources/gentopo.h b/sources/gentopo.h deleted file mode 100644 index ea8c32c9..00000000 --- a/sources/gentopo.h +++ /dev/null @@ -1,162 +0,0 @@ -#pragma once - -// Temporal definition of big integer -#define BigInt long -#define ToBigInt(x) ((BigInt) (x)) -// -// Temporal definition of ratio of two big integers -#define Fraction double -#define ToFraction(x, y) (((double) (x))/((double) (y))) - -//============================================================== -// Output to FROM - -//============================================================== -// class of nodes for T_EGraph - -class T_ENode { - public: - int deg; - int ext; - int *edges; -}; - -class T_EEdge { - public: - int nodes[2]; - int ext; - int momn; - char momc[2]; // no more used - // momentum is printed like: ("%s%d", (enode.ext)?"Q":"p", enode.momn) - char padding[6]; -}; - -class T_EGraph { - public: - long gId; - int pId; - int nNodes; - int nEdges; - int maxdeg; - int nExtern; - - int opi; - BigInt nsym, esym; - - T_ENode *nodes; - T_EEdge *edges; - - T_EGraph(int nnodes, int nedges, int mxdeg); - ~T_EGraph(); - void print(void); - void init(int pid, long gid, int **adjmat, int sopi, BigInt nsym, BigInt esym); - - void setExtern(int nd, int val); - void endSetExtern(void); -}; - -//============================================================== -// class of nodes for T_MGraph - -class T_MNode { - public: - int id; // node id - int deg; // degree(node) = the number of legs - int clss; // initial class number in which the node belongs - int ext; - - int freelg; // the number of free legs - int visited; - - T_MNode(int id, int deg, int ext, int clss); -}; - -//=============================================================== -// class of node-classes for T_MGraph -class T_MNodeClass; - -//=============================================================== -// class of scalar graph expressed by matrix form -// -// Input : the classified set of nodes. -// Output : control passed to 'T_EGraph(self)' - -class T_MGraph { - - public: - - // initial conditions - int pId; // process/subprocess ID - int nNodes; // the number of nodes - int nEdges; // the number of edges - int nLoops; // the number of loops - T_MNode **nodes; // table of T_MNode object - int *clist; // list of initial classes - int nClasses; // the number of initial classes - // int ndcl; // node --> initial class number - int mindeg; // minimum value of degree of nodes - int maxdeg; // maximum value of degree of nodes - - int selOPI; // flag to select 1PI graphs - - // generated set of graphs - long ndiag; // the total number of generated graphs - long n1PI; // the total number of 1PI graphs - int nBridges; // the number of bridges - - // the current graph - int c1PI; // the number of 1PI components - int **adjMat; // adjacency matrix - BigInt nsym; // symmetry factor from nodes - BigInt esym; // symmetry factor from edges - Fraction wsum; // weighted sum of graphs - Fraction wsopi; // weighted sum of 1PI graphs - T_MNodeClass *curcl; // the current 'T_MNodeClass' object - T_EGraph *egraph; - - // measures of efficiency - long ngen; // generated graph before check - long ngconn; // generated connected graph before check - - long nCallRefine; - long discardOrd; - long discardRefine; - long discardDisc; - long discardIso; - - // functions */ - T_MGraph(int pid, int ncl, int *cldeg, int *clnum, int *clext, int sopi); - ~T_MGraph(); - long generate(void); - - private: - // work space for isomorphism - int **modmat; // permutated adjacency matrix - int *permp; - int *permq; - int *permr; - - // work space for biconnected component - int *bidef; - int *bilow; - int bicount; - int padding; - - void printAdjMat(T_MNodeClass *cl); - int isConnected(void); - int visit(int nd); - int isIsomorphic(T_MNodeClass *cl); - void permMat(int size, int *perm, int **mat0, int **mat1); - int compMat(int size, int **mat0, int **mat1); - T_MNodeClass *refineClass(T_MNodeClass *cl, int cn); - void bisearchM(int nd, int pd, int ne); - int count1PI(void); - - int findNextCl(T_MNodeClass *cl, int *dscl); - int findNextTCl(T_MNodeClass *cl, int *dcl); - void connectClass(T_MNodeClass *cl, int *dscl); - void connectNode(int sc, int ss, T_MNodeClass *cl, int *dscl); - void connectLeg(int sc, int sn, int tc, int ts, T_MNodeClass *cl, int *dscl, int* dtcl); - void newGraph(T_MNodeClass *cl); - -}; diff --git a/sources/grcc.cc b/sources/grcc.cc index e2d2829a..a2c86a9a 100644 --- a/sources/grcc.cc +++ b/sources/grcc.cc @@ -138,6 +138,7 @@ static void *erExitArg = NULL; // Utility functions //============================================================== +static void grcc_fprintf(FILE* out, const char* fmt, ...); static void erEnd(const char *msg); static Bool nextPart(int nelem, int nclist, int *clist, int *nl, int *r); static void prilist(int n, const int *a, const char *msg); @@ -176,10 +177,10 @@ static Bool isIn(int n, int *a, int v); Options::Options(void) { if (nOptDef != GRCC_OPT_Size) { - fprintf(GRCC_Stderr, "*** Options: inconsistent default values\n"); - fprintf(GRCC_Stderr, "nOptDef=%d, GRCC_OPT_Size=%d\n", + grcc_fprintf(GRCC_Stderr, "*** Options: inconsistent default values\n"); + grcc_fprintf(GRCC_Stderr, "nOptDef=%d, GRCC_OPT_Size=%d\n", nOptDef, GRCC_OPT_Size); - exit(1); + GRCC_ABORT(); } model = NULL; @@ -198,10 +199,10 @@ Options::Options(void) // QGraf options if (nOptQGDef != GRCC_QGRAF_OPT_Size) { - fprintf(GRCC_Stderr, "*** Options: inconsistent default values\n"); - fprintf(GRCC_Stderr, "nOptQGDef=%d, GRCC_QGRAF_OPT_Size=%d\n", + grcc_fprintf(GRCC_Stderr, "*** Options: inconsistent default values\n"); + grcc_fprintf(GRCC_Stderr, "nOptQGDef=%d, GRCC_QGRAF_OPT_Size=%d\n", nOptQGDef, GRCC_QGRAF_OPT_Size); - exit(1); + GRCC_ABORT(); } nqgopt = 0; for (int j = 0; j < nOptQGDef; j++) { @@ -318,9 +319,9 @@ void Options::setValue(int ind, int val) values[ind] = val; } else { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** Options::setValue : invalid index=%d ", + grcc_fprintf(GRCC_Stderr, "*** Options::setValue : invalid index=%d ", ind); - fprintf(GRCC_Stderr, "(val=%d)\n", val); + grcc_fprintf(GRCC_Stderr, "(val=%d)\n", val); } erEnd("Options::setValue : invalid index"); } @@ -333,7 +334,7 @@ int Options::getValue(int ind) return values[ind]; } else { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** Options::getValue : invalid index=%d\n", ind); + grcc_fprintf(GRCC_Stderr, "*** Options::getValue : invalid index=%d\n", ind); } erEnd("Options::getValue : invalid index"); } @@ -420,18 +421,18 @@ void Options::print(void) { int j; - printf("Options\n"); - printf("+++ GRCC_OPT_Size=%d, print level=%d: ", + grcc_fprintf(GRCC_Stdout, "Options\n"); + grcc_fprintf(GRCC_Stdout, "+++ GRCC_OPT_Size=%d, print level=%d: ", GRCC_OPT_Size, prlevel); - printf("symbol = value (default)\n"); + grcc_fprintf(GRCC_Stdout, "symbol = value (default)\n"); for (j=0; j < GRCC_OPT_Size; j++) { - printf(" %4d GRCC_OPT_%-15s = %2d (%2d)\n", + grcc_fprintf(GRCC_Stdout, " %4d GRCC_OPT_%-15s = %2d (%2d)\n", j, optDef[j].name, values[j], optDef[j].defaultv); } - printf(" outgrf=%s, outgrp=%s\n", out->outgrf, out->outgrp); - printf(" GRCC_QGRAF_OPT_Size=%d:\n", GRCC_QGRAF_OPT_Size); + grcc_fprintf(GRCC_Stdout, " outgrf=%s, outgrp=%s\n", out->outgrf, out->outgrp); + grcc_fprintf(GRCC_Stdout, " GRCC_QGRAF_OPT_Size=%d:\n", GRCC_QGRAF_OPT_Size); for (j=0; j < GRCC_QGRAF_OPT_Size; j++) { - printf(" %4d %-10s = %2d\n", j, optQGDef[j].name, qgopt[j]); + grcc_fprintf(GRCC_Stdout, " %4d %-10s = %2d\n", j, optQGDef[j].name, qgopt[j]); } } @@ -474,7 +475,7 @@ void Options::printModel(void) if (model != NULL) { model->prModel(); } else { - printf("*** model is not defined\n"); + grcc_fprintf(GRCC_Stdout, "*** model is not defined\n"); } } } @@ -512,23 +513,23 @@ void Options::begin(Model *mdl) void Options::end(void) { if (prlevel > 1) { - printf("Optimization: "); + grcc_fprintf(GRCC_Stdout, "Optimization: "); #ifdef SIMPSEL - printf("SIMPSEL=1 "); + grcc_fprintf(GRCC_Stdout, "SIMPSEL=1 "); #else - printf("SIMPSEL=0 "); + grcc_fprintf(GRCC_Stdout, "SIMPSEL=0 "); #endif #ifdef MINMAXLEG - printf("MINMAXLEG=1 "); + grcc_fprintf(GRCC_Stdout, "MINMAXLEG=1 "); #else - printf("MINMAXLEG=0 "); + grcc_fprintf(GRCC_Stdout, "MINMAXLEG=0 "); #endif #ifdef OPTEXTONLY - printf("OPTEXTONLY=1 "); + grcc_fprintf(GRCC_Stdout, "OPTEXTONLY=1 "); #else - printf("OPTEXTONLY=0 "); + grcc_fprintf(GRCC_Stdout, "OPTEXTONLY=0 "); #endif - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); } if (out != NULL && values[GRCC_OPT_Outgrf]) { @@ -571,68 +572,68 @@ void Options::endProc(void) return; } if (prlevel > 0) { - printf("\n"); - printf("+++ Proc %d: ext=%d, loop=%d, ", + grcc_fprintf(GRCC_Stdout, "\n"); + grcc_fprintf(GRCC_Stdout, "+++ Proc %d: ext=%d, loop=%d, ", proc->id, proc->nExtern, proc->loop); if (model != NULL) { - printf("order="); + grcc_fprintf(GRCC_Stdout, "order="); prIntArray(model->ncouple, proc->clist, ": "); model->prParticleArray(proc->ninitl, proc->initlPart, "-->"); model->prParticleArray(proc->nfinal, proc->finalPart, ""); } - printf(" (%8.2f sec)\n", proc->sec); + grcc_fprintf(GRCC_Stdout, " (%8.2f sec)\n", proc->sec); - printf(" Proc %d: Total M-Graphs=%ld, M-Graphs=", + grcc_fprintf(GRCC_Stdout, " Proc %d: Total M-Graphs=%ld, M-Graphs=", proc->id, proc->nMGraphs); proc->wMGraphs.print(" (Conn)\n"); - printf(" Proc %d: Total M-Graphs=%ld, M-Graphs=", + grcc_fprintf(GRCC_Stdout, " Proc %d: Total M-Graphs=%ld, M-Graphs=", proc->id, proc->nMOPI); proc->wMOPI.print(" (1PI)\n"); - printf(" Proc %d: Total A-Graphs=%ld, A-Graphs=", + grcc_fprintf(GRCC_Stdout, " Proc %d: Total A-Graphs=%ld, A-Graphs=", proc->id, proc->nAGraphs); proc->wAGraphs.print(" (Conn)\n"); - printf(" Proc %d: Total A-Graphs=%ld, A-Graphs=", + grcc_fprintf(GRCC_Stdout, " Proc %d: Total A-Graphs=%ld, A-Graphs=", proc->id, proc->nAOPI); proc->wAOPI.print(" (1PI)\n"); - printf("# { %d,{", proc->ninitl); + grcc_fprintf(GRCC_Stdout, "# { %d,{", proc->ninitl); for (k = 0; k < proc->ninitl; k++) { if (k != 0) { - printf(", "); + grcc_fprintf(GRCC_Stdout, ", "); } if (proc->model != NULL) { - printf("\"%s\"", proc->model->particleName(proc->initlPart[k])); + grcc_fprintf(GRCC_Stdout, "\"%s\"", proc->model->particleName(proc->initlPart[k])); } else { - printf("%d", proc->initlPart[k]); + grcc_fprintf(GRCC_Stdout, "%d", proc->initlPart[k]); } } - printf("}, %d,{", proc->nfinal); + grcc_fprintf(GRCC_Stdout, "}, %d,{", proc->nfinal); for (k = 0; k < proc->nfinal; k++) { if (k != 0) { - printf(", "); + grcc_fprintf(GRCC_Stdout, ", "); } if (proc->model != NULL) { - printf("\"%s\"", proc->model->particleName(proc->finalPart[k])); + grcc_fprintf(GRCC_Stdout, "\"%s\"", proc->model->particleName(proc->finalPart[k])); } else { - printf("%d", proc->initlPart[k]); + grcc_fprintf(GRCC_Stdout, "%d", proc->initlPart[k]); } } - printf("}, "); + grcc_fprintf(GRCC_Stdout, "}, "); if (model != NULL) { - printf("{"); + grcc_fprintf(GRCC_Stdout, "{"); for (k = 0; k < model->ncouple; k++) { if (k != 0) { - printf(", "); + grcc_fprintf(GRCC_Stdout, ", "); } - printf("%d", proc->clist[k]); + grcc_fprintf(GRCC_Stdout, "%d", proc->clist[k]); } - printf("},"); + grcc_fprintf(GRCC_Stdout, "},"); } - printf("},%6ldL,%6ldL,%3ldL, -1.0, %4.2f},\n", + grcc_fprintf(GRCC_Stdout, "},%6ldL,%6ldL,%3ldL, -1.0, %4.2f},\n", proc->nAOPI, proc->wAOPI.num, proc->wAOPI.den, proc->sec); } @@ -687,39 +688,39 @@ void Options::endSubProc(void) } if (prlevel > 1) { - printf("\n"); - printf("+++ Subproc %d: ext=%d, loop=%d, nodes=%d, edges=%d\n", + grcc_fprintf(GRCC_Stdout, "\n"); + grcc_fprintf(GRCC_Stdout, "+++ Subproc %d: ext=%d, loop=%d, nodes=%d, edges=%d\n", sproc->id, sproc->nExtern, sproc->loop, sproc->nNodes, sproc->nEdges); - printf(" Subproc %d: Total M-Graphs=%ld, M-Wsum=", + grcc_fprintf(GRCC_Stdout, " Subproc %d: Total M-Graphs=%ld, M-Wsum=", sproc->id, sproc->nMGraphs); sproc->wMGraphs.print(" (Conn)\n"); - printf(" Subproc %d: Total M-Graphs=%ld, M-Wsum=", + grcc_fprintf(GRCC_Stdout, " Subproc %d: Total M-Graphs=%ld, M-Wsum=", sproc->id, sproc->nMOPI); sproc->wMOPI.print(" (1PI)\n"); - printf(" Subproc %d: Total A-Graphs=%ld, A-Wsum=", + grcc_fprintf(GRCC_Stdout, " Subproc %d: Total A-Graphs=%ld, A-Wsum=", sproc->id, sproc->nAGraphs); sproc->wAGraphs.print(" (Conn)\n"); - printf(" Subproc %d: Total A-Graphs=%ld, A-Wsum=", + grcc_fprintf(GRCC_Stdout, " Subproc %d: Total A-Graphs=%ld, A-Wsum=", sproc->id, sproc->nAOPI); sproc->wAOPI.print(" (1PI)\n"); } if (prlevel > 0) { - printf("\n"); - printf("* Total %ld MGraphs; %ld 1PI", mgraph->cDiag, mgraph->c1PI); - printf(" wscon = "); + grcc_fprintf(GRCC_Stdout, "\n"); + grcc_fprintf(GRCC_Stdout, "* Total %ld MGraphs; %ld 1PI", mgraph->cDiag, mgraph->c1PI); + grcc_fprintf(GRCC_Stdout, " wscon = "); mgraph->wscon.print(" ( "); mgraph->wsopi.print(" 1PI)\n"); #ifdef MONITOR - printf("* refine: %ld\n", mgraph->nCallRefine); - printf("* discarded for refinement: %ld\n", mgraph->discardRefine); - printf("* discarded for disconnected: %ld\n", mgraph->discardDisc); - printf("* discarded for duplication: %ld\n", mgraph->discardIso); + grcc_fprintf(GRCC_Stdout, "* refine: %ld\n", mgraph->nCallRefine); + grcc_fprintf(GRCC_Stdout, "* discarded for refinement: %ld\n", mgraph->discardRefine); + grcc_fprintf(GRCC_Stdout, "* discarded for disconnected: %ld\n", mgraph->discardDisc); + grcc_fprintf(GRCC_Stdout, "* discarded for duplication: %ld\n", mgraph->discardIso); #endif } @@ -825,7 +826,7 @@ Output::~Output(void) { if (outgrfp != NULL) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** file has not been closed : \"%s\"\n", + grcc_fprintf(GRCC_Stderr, "*** file has not been closed : \"%s\"\n", outgrf); } fclose(outgrfp); @@ -838,7 +839,7 @@ Output::~Output(void) } if (outgrpp != NULL) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** file has not been closed : \"%s\"\n", + grcc_fprintf(GRCC_Stderr, "*** file has not been closed : \"%s\"\n", outgrp); } fclose(outgrpp); @@ -900,7 +901,7 @@ Bool Output::outBeginF(Model *mdl, Bool pr) return True; } else if (outgrfp != NULL) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** outBegin: \"%s\" is already opened\n", + grcc_fprintf(GRCC_Stderr, "*** outBegin: \"%s\" is already opened\n", outgrf); } erEnd("outBegin: file is already opened\n"); @@ -911,7 +912,7 @@ Bool Output::outBeginF(Model *mdl, Bool pr) } if ((outgrfp = fopen(outgrf, "w")) == NULL) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** outBegin: cannot open %s\n", outgrf); + grcc_fprintf(GRCC_Stderr, "*** outBegin: cannot open %s\n", outgrf); } return False; } @@ -946,7 +947,7 @@ Bool Output::outBeginP(Model *mdl, Bool pr) return True; } else if (outgrpp != NULL) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** outBegin: \"%s\" is already opened\n", + grcc_fprintf(GRCC_Stderr, "*** outBegin: \"%s\" is already opened\n", outgrp); } erEnd("outBegin: file is already opened\n"); @@ -957,7 +958,7 @@ Bool Output::outBeginP(Model *mdl, Bool pr) } if ((outgrpp = fopen(outgrp, "w")) == NULL) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** outBegin: cannot open %s\n", outgrp); + grcc_fprintf(GRCC_Stderr, "*** outBegin: cannot open %s\n", outgrp); } return False; } @@ -1341,11 +1342,11 @@ void Output::outEGraphP(EGraph *eg) } if (eg->assigned) { - printf("outEGraphP:egraph:egraph[%ld]\n", eg->mId-1); + grcc_fprintf(GRCC_Stdout, "outEGraphP:egraph:egraph[%ld]\n", eg->mId-1); eg->printPy(outgrpp, eg->mId); return; } else { - printf("outEGraphP:mgraph:egraph[%ld]\n", eg->mId-1); + grcc_fprintf(GRCC_Stdout, "outEGraphP:mgraph:egraph[%ld]\n", eg->mId-1); eg->mgraph->printPy(outgrpp, eg->mId); return; } @@ -1475,7 +1476,7 @@ void Output::outModelF(void) if (model == NULL) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** Output::outModel : model is not defined\n"); + grcc_fprintf(GRCC_Stderr, "*** Output::outModel : model is not defined\n"); } return; } @@ -1485,7 +1486,7 @@ void Output::outModelF(void) if ((mdlfp = fopen(mdlfn, "w")) == NULL) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** Output::outModel : cannot open \"%s\"\n", + grcc_fprintf(GRCC_Stderr, "*** Output::outModel : cannot open \"%s\"\n", mdlfn); } return; @@ -1556,8 +1557,8 @@ void Output::outModelP(void) if (model == NULL) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** Output::outModel : "); - fprintf(GRCC_Stderr, "model is not defined\n"); + grcc_fprintf(GRCC_Stderr, "*** Output::outModel : "); + grcc_fprintf(GRCC_Stderr, "model is not defined\n"); } return; } @@ -1567,7 +1568,7 @@ void Output::outModelP(void) if ((mdlfp = fopen(mdlfn, "w")) == NULL) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** Output::outModel : cannot open \"%s\"\n", + grcc_fprintf(GRCC_Stderr, "*** Output::outModel : cannot open \"%s\"\n", mdlfn); } return; @@ -1687,7 +1688,7 @@ Particle::Particle(Model *modl, int pid, PInput *pinp) if (Abs(mdl->defpart) == GRCC_DEFBYCODE) { if (pinp->ptypec < GRCC_PT_Undef || pinp->ptypec > GRCC_PT_Size) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** particle type is not defined: %d\n", + grcc_fprintf(GRCC_Stderr, "*** particle type is not defined: %d\n", pinp->ptypec); } erEnd("particle type is not defined (illegal code)"); @@ -1706,7 +1707,7 @@ Particle::Particle(Model *modl, int pid, PInput *pinp) } if (ptype < 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** particle type \"%s\" is not defined\n", + grcc_fprintf(GRCC_Stderr, "*** particle type \"%s\" is not defined\n", pinp->ptypen); } erEnd("particle type is not defined (name)"); @@ -1782,11 +1783,11 @@ char *Particle::aparticle(void) void Particle::prParticle(void) { if (isNeutral()) { - printf("pid=%d, name=\"%s\", real_field, ", id, name); + grcc_fprintf(GRCC_Stdout, "pid=%d, name=\"%s\", real_field, ", id, name); } else { - printf("pid=%d, name=\"%s\", aname=\"%s\", ", id, name, aname); + grcc_fprintf(GRCC_Stdout, "pid=%d, name=\"%s\", aname=\"%s\", ", id, name, aname); } - printf("ptype=%s, pcode=%d, acode=%d, extonly=%d, cdeg=(%d,%d)\n", + grcc_fprintf(GRCC_Stdout, "ptype=%s, pcode=%d, acode=%d, extonly=%d, cdeg=(%d,%d)\n", ptypenames[ptype], pcode, acode, extonly, cmindeg, cmaxdeg); } @@ -1858,7 +1859,7 @@ Interaction::Interaction(Model *modl, int iid, const char *nam, int icd, int *cp if (Abs(sdir) == 1) { jdir = j; } else if ((Abs(sdir) == 0 && j != jdir + 1) || Abs(sdir) > 1) { - fprintf(GRCC_Stderr, "*** Interaction: Dirac and anti-Dirac should " + grcc_fprintf(GRCC_Stderr, "*** Interaction: Dirac and anti-Dirac should " "arranged in pairs in an interaction.\n"); ok = False; } @@ -1872,7 +1873,7 @@ Interaction::Interaction(Model *modl, int iid, const char *nam, int icd, int *cp if (Abs(sgho) == 1) { jgho = j; } else if ((Abs(sgho) == 0 && j != jgho + 1) || Abs(sgho) > 1) { - fprintf(GRCC_Stderr, "*** Interaction: Ghost and anti-Ghost should " + grcc_fprintf(GRCC_Stderr, "*** Interaction: Ghost and anti-Ghost should " "arranged in pairs in an interaction.\n"); ok = False; } @@ -1881,32 +1882,32 @@ Interaction::Interaction(Model *modl, int iid, const char *nam, int icd, int *cp if (nmaj % 2 == 1) { jmaj = j; } else if (j != jmaj + 1) { - fprintf(GRCC_Stderr, "*** Interaction: two Majoranas should " + grcc_fprintf(GRCC_Stderr, "*** Interaction: two Majoranas should " "arranged in pairs in an interaction.\n"); ok = False; } } } if (sdir != 0) { - fprintf(GRCC_Stderr, "*** Interaction: Dirac number is not conserved\n"); + grcc_fprintf(GRCC_Stderr, "*** Interaction: Dirac number is not conserved\n"); ok = False; } if (sgho != 0) { - fprintf(GRCC_Stderr, "*** Interaction: Ghost number is not conserved\n"); + grcc_fprintf(GRCC_Stderr, "*** Interaction: Ghost number is not conserved\n"); ok = False; } if (nmaj % 2 != 0) { - fprintf(GRCC_Stderr, "*** Interaction: odd number of Majorana particles\n"); + grcc_fprintf(GRCC_Stderr, "*** Interaction: odd number of Majorana particles\n"); ok = False; } if (ndir + nmaj + ngho > 2 && prmsg) { - fprintf(GRCC_Stderr, "+++ Interaction: more than 2 " + grcc_fprintf(GRCC_Stderr, "+++ Interaction: more than 2 " "Dirac/Majorana/Ghost " "particles are interacting.\n"); - fprintf(GRCC_Stderr, " Interaction has %d Dirac, %d Ghost " + grcc_fprintf(GRCC_Stderr, " Interaction has %d Dirac, %d Ghost " "and %d Majorana particles.\n", ndir, nmaj, ngho); - fprintf(GRCC_Stderr, " Sign factors related to fermions may " + grcc_fprintf(GRCC_Stderr, " Sign factors related to fermions may " "be inconsistent with your calculation method.\n"); prmsg = False; } @@ -1936,16 +1937,16 @@ void Interaction::prInteraction(void) { int j; - printf("vid=%d, icode=%d, name=\"%s\", loop=%d, csum=%d, cpl=", + grcc_fprintf(GRCC_Stdout, "vid=%d, icode=%d, name=\"%s\", loop=%d, csum=%d, cpl=", id, icode, name, loop, csum); prIntArray(mdl->ncouple, clist, ", legs=["); for (j = 0; j < nplist; j++) { if (j != 0) { - printf(", "); + grcc_fprintf(GRCC_Stdout, ", "); } - printf("%s", mdl->particleName(plist[j])); + grcc_fprintf(GRCC_Stdout, "%s", mdl->particleName(plist[j])); } - printf("];\n"); + grcc_fprintf(GRCC_Stdout, "];\n"); } @@ -2061,46 +2062,46 @@ void Model::prModel(void) static char hd1[] = "#-------------------------------------------------\n"; int j; - printf("%s", hdr); - printf("Model=\"%s\", ", name); - printf("coupling=["); + grcc_fprintf(GRCC_Stdout, "%s", hdr); + grcc_fprintf(GRCC_Stdout, "Model=\"%s\", ", name); + grcc_fprintf(GRCC_Stdout, "coupling=["); for (j = 0; j < ncouple; j++) { if (j != 0) { - printf(", "); + grcc_fprintf(GRCC_Stdout, ", "); } - printf("\"%s\"", cnlist[j]); + grcc_fprintf(GRCC_Stdout, "\"%s\"", cnlist[j]); } - printf("];\n"); - printf("%s", hd1); - printf("Particles=%d;\n", nParticles); + grcc_fprintf(GRCC_Stdout, "];\n"); + grcc_fprintf(GRCC_Stdout, "%s", hd1); + grcc_fprintf(GRCC_Stdout, "Particles=%d;\n", nParticles); for (j = 1; j < nParticles; j++) { particles[j]->prParticle(); } - printf("EndParticle;\n"); - printf("allPart (%d) = [", nallPart); + grcc_fprintf(GRCC_Stdout, "EndParticle;\n"); + grcc_fprintf(GRCC_Stdout, "allPart (%d) = [", nallPart); for (j = 0; j < nallPart; j++) { if (j != 0) { - printf(", "); + grcc_fprintf(GRCC_Stdout, ", "); } - printf("%d", allPart[j]); + grcc_fprintf(GRCC_Stdout, "%d", allPart[j]); } - printf("]\n"); + grcc_fprintf(GRCC_Stdout, "]\n"); - printf("%s", hd1); - printf("Interactions=%d;\n", nInteracts); + grcc_fprintf(GRCC_Stdout, "%s", hd1); + grcc_fprintf(GRCC_Stdout, "Interactions=%d;\n", nInteracts); for (j = 0; j < nInteracts; j++) { interacts[j]->prInteraction(); } - printf("EndInteraction;\n"); - printf("%s", hd1); - printf("InteractionTable;\n"); - printf("# %d class in (total coupling, degree)\n", ncplgcp); + grcc_fprintf(GRCC_Stdout, "EndInteraction;\n"); + grcc_fprintf(GRCC_Stdout, "%s", hd1); + grcc_fprintf(GRCC_Stdout, "InteractionTable;\n"); + grcc_fprintf(GRCC_Stdout, "# %d class in (total coupling, degree)\n", ncplgcp); for (j = 0; j < ncplgcp; j++) { - printf(" class=%d: cp=%d, lg=%d, vl=", j, cplgcp[j], cplglg[j]); + grcc_fprintf(GRCC_Stdout, " class=%d: cp=%d, lg=%d, vl=", j, cplgcp[j], cplglg[j]); prIntArray(cplgnvl[j], cplgvl[j], "\n"); } - printf("EndInteractionTable;\n"); - printf("%s", hdr); + grcc_fprintf(GRCC_Stdout, "EndInteractionTable;\n"); + grcc_fprintf(GRCC_Stdout, "%s", hdr); } //-------------------------------------------------------------- @@ -2118,9 +2119,9 @@ void Model::addParticle(PInput *pinp) acd = findParticleCode(pinp->acode); if (pcd > 0 || acd > 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** particle code [%d, %d] ", + grcc_fprintf(GRCC_Stderr, "*** particle code [%d, %d] ", pinp->pcode, pinp->acode); - fprintf(GRCC_Stderr, "is already used.\n"); + grcc_fprintf(GRCC_Stderr, "is already used.\n"); } erEnd("particle code is already defined"); } @@ -2133,9 +2134,9 @@ void Model::addParticle(PInput *pinp) aid = findParticleName(pinp->aname); if (nid > 0 || aid > 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** particle name [%s, %s] ", + grcc_fprintf(GRCC_Stderr, "*** particle name [%s, %s] ", pinp->name, pinp->aname); - fprintf(GRCC_Stderr, "is already used.\n"); + grcc_fprintf(GRCC_Stderr, "is already used.\n"); } erEnd("particle name is already defined."); } @@ -2205,7 +2206,7 @@ void Model::addInteraction(IInput *iinp) if (defpart == GRCC_DEFBYCODE) { if (iinp->icode < 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** vertex code %d should be positive\n", + grcc_fprintf(GRCC_Stderr, "*** vertex code %d should be positive\n", iinp->icode); } erEnd("vertex code should be positive"); @@ -2213,7 +2214,7 @@ void Model::addInteraction(IInput *iinp) vid = findInteractionCode(iinp->icode); if (vid >= 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** vertex code %d is already used\n", + grcc_fprintf(GRCC_Stderr, "*** vertex code %d is already used\n", iinp->icode); } erEnd("vertex code is already used"); @@ -2226,7 +2227,7 @@ void Model::addInteraction(IInput *iinp) vid = findInteractionName(iinp->name); if (vid >= 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** vertex name %s is already used\n", + grcc_fprintf(GRCC_Stderr, "*** vertex name %s is already used\n", iinp->name); erEnd("vertex name is already used"); } @@ -2247,9 +2248,9 @@ void Model::addInteraction(IInput *iinp) plist[j] = findParticleCode(iinp->plistc[j]); if (plist[j] == 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** particle code %d ", + grcc_fprintf(GRCC_Stderr, "*** particle code %d ", iinp->plistc[j]); - fprintf(GRCC_Stderr, "is not defined\n"); + grcc_fprintf(GRCC_Stderr, "is not defined\n"); } erEnd("particle is not defined (code)"); } @@ -2257,9 +2258,9 @@ void Model::addInteraction(IInput *iinp) plist[j] = findParticleName(iinp->plistn[j]); if (plist[j] == 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** particle name %s ", + grcc_fprintf(GRCC_Stderr, "*** particle name %s ", iinp->plistn[j]); - fprintf(GRCC_Stderr, "is not defined\n"); + grcc_fprintf(GRCC_Stderr, "is not defined\n"); } erEnd("particle is not defined (name)"); } @@ -2271,9 +2272,9 @@ void Model::addInteraction(IInput *iinp) c = iinp->cvallist[j]; if (c < 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** illegal value of c-constants \n"); + grcc_fprintf(GRCC_Stderr, "*** illegal value of c-constants \n"); for (k = 0; k < ncouple; k++) { - fprintf(GRCC_Stderr, " %s = %d\n", + grcc_fprintf(GRCC_Stderr, " %s = %d\n", cnlist[k], iinp->cvallist[k]); } } @@ -2283,10 +2284,10 @@ void Model::addInteraction(IInput *iinp) } if (csum < 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** illegal total coupling-constants %d\n", + grcc_fprintf(GRCC_Stderr, "*** illegal total coupling-constants %d\n", csum); for (k = 0; k < ncouple; k++) { - fprintf(GRCC_Stderr, " %s = %d\n", + grcc_fprintf(GRCC_Stderr, " %s = %d\n", cnlist[k], iinp->cvallist[k]); } } @@ -2297,8 +2298,8 @@ void Model::addInteraction(IInput *iinp) lp2 = csum - nlegs + 2; if (lp2 % 2 != 0 || lp2 < 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** illegal coupling const : "); - fprintf(GRCC_Stderr, "nlegs - 2 + 2*loop: 2*loop = %d\n", lp2); + grcc_fprintf(GRCC_Stderr, "*** illegal coupling const : "); + grcc_fprintf(GRCC_Stderr, "nlegs - 2 + 2*loop: 2*loop = %d\n", lp2); } } lp = lp2/2; @@ -2421,8 +2422,8 @@ char *Model::particleName(int p) if (Abs(p) >= nParticles) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "\n*** Model::particleName: "); - fprintf(GRCC_Stderr, "illegal particle id=%d\n", p); + grcc_fprintf(GRCC_Stderr, "\n*** Model::particleName: "); + grcc_fprintf(GRCC_Stderr, "illegal particle id=%d\n", p); } erEnd("Model::particleName: illegal particle"); } @@ -2440,7 +2441,7 @@ int Model::particleCode(int p) int q; if (Abs(p) >= nParticles) { - fprintf(GRCC_Stderr, "\n*** Model::particleCode: illegal particle id=%d\n", + grcc_fprintf(GRCC_Stderr, "\n*** Model::particleCode: illegal particle id=%d\n", p); erEnd("Model::particleCode: illegal particle"); } @@ -2457,14 +2458,14 @@ void Model::prParticleArray(int n, int *a, const char *msg) { int j; - printf("["); + grcc_fprintf(GRCC_Stdout, "["); for (j = 0; j < n; j++) { if (j != 0) { - printf(", "); + grcc_fprintf(GRCC_Stdout, ", "); } - printf("%s", particleName(a[j])); + grcc_fprintf(GRCC_Stdout, "%s", particleName(a[j])); } - printf("]%s", msg); + grcc_fprintf(GRCC_Stdout, "]%s", msg); } //-------------------------------------------------------------- @@ -2547,27 +2548,27 @@ int Model::findInteractionCode(int icd) //-------------------------------------------------------------- void Model::printMInput(MInput *min) { - printf(" \"Model\": [\"%s\", %d, [", + grcc_fprintf(GRCC_Stdout, " \"Model\": [\"%s\", %d, [", min->name, min->ncouple); for (int j = 0; j < min->ncouple; j++) { if (j != 0) { - printf(", "); + grcc_fprintf(GRCC_Stdout, ", "); } - printf("\"%s\"", min->cnamlist[j]); + grcc_fprintf(GRCC_Stdout, "\"%s\"", min->cnamlist[j]); } - printf("], "); + grcc_fprintf(GRCC_Stdout, "], "); if (min->defpart == GRCC_DEFBYCODE) { - printf("\"ByCode\""); + grcc_fprintf(GRCC_Stdout, "\"ByCode\""); } else { - printf("\"ByName\""); + grcc_fprintf(GRCC_Stdout, "\"ByName\""); } - printf("],\n"); + grcc_fprintf(GRCC_Stdout, "],\n"); } //-------------------------------------------------------------- void Model::printPInput(PInput *pin) { - printf(" [\"%s\"(%d), \"%s\"(%d), \"%s\"(%d), %d],\n", + grcc_fprintf(GRCC_Stdout, " [\"%s\"(%d), \"%s\"(%d), \"%s\"(%d), %d],\n", pin->name, pin->pcode, pin->aname, pin->acode, pin->ptypen, pin->ptypec, pin->extonly); } @@ -2577,21 +2578,21 @@ void Model::printIInput(IInput *iin) { int j; - printf(" [\"%s\"(%d), %d, [", iin->name, iin->icode, iin->nplistn); + grcc_fprintf(GRCC_Stdout, " [\"%s\"(%d), %d, [", iin->name, iin->icode, iin->nplistn); for (j = 0; j < iin->nplistn; j++) { if (j != 0) { - printf(", "); + grcc_fprintf(GRCC_Stdout, ", "); } - printf("\"%s\"(%d)", iin->plistn[j], iin->plistc[j]); + grcc_fprintf(GRCC_Stdout, "\"%s\"(%d)", iin->plistn[j], iin->plistc[j]); } - printf("], ["); + grcc_fprintf(GRCC_Stdout, "], ["); for (j = 0; j < GRCC_MAXNCPLG; j++) { if (j != 0) { - printf(", "); + grcc_fprintf(GRCC_Stdout, ", "); } - printf("%d", iin->cvallist[j]); + grcc_fprintf(GRCC_Stdout, "%d", iin->cvallist[j]); } - printf("],\n"); + grcc_fprintf(GRCC_Stdout, "],\n"); } //************************************************************** @@ -2631,13 +2632,13 @@ PNodeClass::PNodeClass(SProcess *spc, int nnods, int nclss, NCInput *cls) lp2 = (couple[j]-deg[j]+2); if (!isATExternal(type[j]) && (lp2 % 2 != 0 || lp2 < 0)) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** PNodeClass: illegal loop: " + grcc_fprintf(GRCC_Stderr, "*** PNodeClass: illegal loop: " ": 2*loop=cpl[%d](%d)-deg[%d](%d)+2 =2*loop = %d\n", j, couple[j], j, deg[j], lp2); for (k = 0; k < nclss; k++) { - fprintf(GRCC_Stderr, "k=%d, deg=%d, typ=%d, ptcl=%d, ", + grcc_fprintf(GRCC_Stderr, "k=%d, deg=%d, typ=%d, ptcl=%d, ", k, deg[k], type[k], particle[k]); - fprintf(GRCC_Stderr, "cpl=%d, cnt=%d\n", + grcc_fprintf(GRCC_Stderr, "cpl=%d, cnt=%d\n", couple[k], count[k]); } } @@ -2663,8 +2664,8 @@ PNodeClass::PNodeClass(SProcess *spc, int nnods, int nclss, NCInput *cls) cl2mcl[j] = spc->model->findMClass(couple[j], deg[j]); if (cl2mcl[j] < 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** PNodeClass : no vertex : "); - fprintf(GRCC_Stderr, "coupling=%d, degree=%d\n", + grcc_fprintf(GRCC_Stderr, "*** PNodeClass : no vertex : "); + grcc_fprintf(GRCC_Stderr, "coupling=%d, degree=%d\n", couple[j], deg[j]); } erEnd("PNodeClass : no vertex"); @@ -2709,13 +2710,13 @@ PNodeClass::PNodeClass(SProcess *spc, int nnods, int nclss, int *dgs, int *typ, lp2 = (cpl[j]-dgs[j]+2); if (!isATExternal(type[j]) && (lp2 % 2 != 0 || lp2 < 0)) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** PNodeClass: illegal loop: " + grcc_fprintf(GRCC_Stderr, "*** PNodeClass: illegal loop: " ": 2*loop=cpl[%d](%d)-deg[%d](%d)+2 =2*loop = %d\n", j, cpl[j], j, dgs[j], lp2); for (k = 0; k < nclss; k++) { - fprintf(GRCC_Stderr, "k=%d, dgs=%d, typ=%d, ptcl=%d, ", + grcc_fprintf(GRCC_Stderr, "k=%d, dgs=%d, typ=%d, ptcl=%d, ", k, dgs[k], typ[k], ptcl[k]); - fprintf(GRCC_Stderr, "cpl=%d, cnt=%d\n", cpl[k], cnt[k]); + grcc_fprintf(GRCC_Stderr, "cpl=%d, cnt=%d\n", cpl[k], cnt[k]); } } ok = False; @@ -2740,8 +2741,8 @@ PNodeClass::PNodeClass(SProcess *spc, int nnods, int nclss, int *dgs, int *typ, cl2mcl[j] = spc->model->findMClass(couple[j], deg[j]); if (cl2mcl[j] < 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** PNodeClass : no vertex : "); - fprintf(GRCC_Stderr, "coupling=%d, degree=%d\n", + grcc_fprintf(GRCC_Stderr, "*** PNodeClass : no vertex : "); + grcc_fprintf(GRCC_Stderr, "coupling=%d, degree=%d\n", couple[j], deg[j]); } erEnd("PNodeClass : no vertex"); @@ -2771,11 +2772,11 @@ void PNodeClass::prPNodeClass(void) { int j; - printf("+++ PNodeClass: nclass=%d, nnodes=%d\n", nclass, nnodes); + grcc_fprintf(GRCC_Stdout, "+++ PNodeClass: nclass=%d, nnodes=%d\n", nclass, nnodes); for (j = 0; j < nclass; j++) { prElem(j); } - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); } //-------------------------------------------------------------- @@ -2783,18 +2784,18 @@ void PNodeClass::prElem(int j) { int tp; - printf("%3d: %-7s(%2d), ", j, GRCC_AT_NdStr(type[j]), type[j]); - printf("deg=%d, count=%d, nodes[%d--%d], couple=%d, cmindeg=%d, cmaxdeg=%d", + grcc_fprintf(GRCC_Stdout, "%3d: %-7s(%2d), ", j, GRCC_AT_NdStr(type[j]), type[j]); + grcc_fprintf(GRCC_Stdout, "deg=%d, count=%d, nodes[%d--%d], couple=%d, cmindeg=%d, cmaxdeg=%d", deg[j], count[j], cl2nd[j], cl2nd[j+1], couple[j], cmindeg[j], cmaxdeg[j]); tp = type[j]; if (isATExternal(tp)) { if (sproc->model != NULL) { - printf(", ptcl=%s ", sproc->model->particleName(particle[j])); + grcc_fprintf(GRCC_Stdout, ", ptcl=%s ", sproc->model->particleName(particle[j])); } else { - printf(", ptcl=%d ", particle[j]); + grcc_fprintf(GRCC_Stdout, ", ptcl=%d ", particle[j]); } } - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); } @@ -2821,7 +2822,7 @@ SProcess::SProcess(Model *mdl, Process *prc, Options *opts, int sid, int *clst, bool ok; if (ncls < 1) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: no class %d\n", ncls); + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: no class %d\n", ncls); erEnd("SProcess::SProcess: no Class"); } id = sid; @@ -2873,7 +2874,7 @@ SProcess::SProcess(Model *mdl, Process *prc, Options *opts, int sid, int *clst, lp2 = cpl[j] - cdeg[j] + 2; if (lp2 % 2 != 0 || lp2 < 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: illegal loop:" + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: illegal loop:" "class %d, cp=%d, deg=%d, lp2=%d\n", j, cp, cdeg[j], lp2); } @@ -2888,9 +2889,9 @@ SProcess::SProcess(Model *mdl, Process *prc, Options *opts, int sid, int *clst, } else if (isATExternal(ctyp[j])) { if (model->normalParticle(ptcl[j]) == 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); - fprintf(GRCC_Stderr, "illegal particle code: "); - fprintf(GRCC_Stderr, "code: class %d, ptcl=%d\n", j, ptcl[j]); + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); + grcc_fprintf(GRCC_Stderr, "illegal particle code: "); + grcc_fprintf(GRCC_Stderr, "code: class %d, ptcl=%d\n", j, ptcl[j]); } ok = False; } else { @@ -2900,25 +2901,25 @@ SProcess::SProcess(Model *mdl, Process *prc, Options *opts, int sid, int *clst, } else { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: illegal type: "); - fprintf(GRCC_Stderr, "class %d, type=%d\n", j, ctyp[j]); + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: illegal type: "); + grcc_fprintf(GRCC_Stderr, "class %d, type=%d\n", j, ctyp[j]); } ok = False; } } if (tcpl0 != tCouple) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); - fprintf(GRCC_Stderr, "illegal coupling constants:"); - fprintf(GRCC_Stderr, " %d != %d\n", tcpl0, tCouple); + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); + grcc_fprintf(GRCC_Stderr, "illegal coupling constants:"); + grcc_fprintf(GRCC_Stderr, " %d != %d\n", tcpl0, tCouple); } ok = False; } if (ndeg == nExtern) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); - fprintf(GRCC_Stderr, "no vertices: "); - fprintf(GRCC_Stderr, "nExtern=%d, ndeg=%d\n", nExtern, ndeg); + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); + grcc_fprintf(GRCC_Stderr, "no vertices: "); + grcc_fprintf(GRCC_Stderr, "nExtern=%d, ndeg=%d\n", nExtern, ndeg); } ok = False; } @@ -2926,17 +2927,17 @@ SProcess::SProcess(Model *mdl, Process *prc, Options *opts, int sid, int *clst, nEdges = ndeg/2; } else { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); - fprintf(GRCC_Stderr, "illegal total deg = %d (not even)\n", ndeg); + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); + grcc_fprintf(GRCC_Stderr, "illegal total deg = %d (not even)\n", ndeg); } ok = False; } nNodes = nvrt + nExtern; if (nNodes >= GRCC_MAXNODES) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); - fprintf(GRCC_Stderr, "too many nodes = %d\n", nNodes); - fprintf(GRCC_Stderr, " nExtern=%d, nvert=%d (GRCC_MAXNODES)\n", + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); + grcc_fprintf(GRCC_Stderr, "too many nodes = %d\n", nNodes); + grcc_fprintf(GRCC_Stderr, " nExtern=%d, nvert=%d (GRCC_MAXNODES)\n", nExtern, nvert); } ok = False; @@ -2952,7 +2953,7 @@ SProcess::SProcess(Model *mdl, Process *prc, Options *opts, int sid, int *clst, lp2 = tCouple - nExtern + 2; if (lp2 % 2 != 0 || lp2 < 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: illegal loop: " + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: illegal loop: " "tCouple=%d, nExtern=%d, 2*loop=%d\n", tCouple, nExtern, lp2); } @@ -2984,7 +2985,7 @@ SProcess::SProcess(Model *mdl, Process *prc, Options *opts, int sid, int *clst, bool ok; if (ncls < 1) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: no class %d\n", ncls); + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: no class %d\n", ncls); erEnd("SProcess::SProcess: no Class"); } id = sid; @@ -3036,7 +3037,7 @@ SProcess::SProcess(Model *mdl, Process *prc, Options *opts, int sid, int *clst, lp2 = cls[j].cple - cls[j].cldeg + 2; if (lp2 % 2 != 0 || lp2 < 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: illegal loop: " + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: illegal loop: " "class %d, cp=%d, deg=%d, lp2=%d\n", j, cp, cls[j].cldeg, lp2); } @@ -3051,9 +3052,9 @@ SProcess::SProcess(Model *mdl, Process *prc, Options *opts, int sid, int *clst, } else if (isATExternal(cls[j].cltyp)) { if (model->normalParticle(cls[j].ptcl) == 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); - fprintf(GRCC_Stderr, "illegal particle code: "); - fprintf(GRCC_Stderr, "code: class %d, ptcl=%d\n", j, cls[j].ptcl); + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); + grcc_fprintf(GRCC_Stderr, "illegal particle code: "); + grcc_fprintf(GRCC_Stderr, "code: class %d, ptcl=%d\n", j, cls[j].ptcl); } ok = False; } else { @@ -3063,25 +3064,25 @@ SProcess::SProcess(Model *mdl, Process *prc, Options *opts, int sid, int *clst, } else { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: illegal type: "); - fprintf(GRCC_Stderr, "class %d, type=%d\n", j, cls[j].cltyp); + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: illegal type: "); + grcc_fprintf(GRCC_Stderr, "class %d, type=%d\n", j, cls[j].cltyp); } ok = False; } } if (tcpl0 != tCouple) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); - fprintf(GRCC_Stderr, "illegal coupling constants:"); - fprintf(GRCC_Stderr, " %d != %d\n", tcpl0, tCouple); + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); + grcc_fprintf(GRCC_Stderr, "illegal coupling constants:"); + grcc_fprintf(GRCC_Stderr, " %d != %d\n", tcpl0, tCouple); } ok = False; } if (ndeg == nExtern) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); - fprintf(GRCC_Stderr, "no vertices: "); - fprintf(GRCC_Stderr, "nExtern=%d, ndeg=%d\n", nExtern, ndeg); + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); + grcc_fprintf(GRCC_Stderr, "no vertices: "); + grcc_fprintf(GRCC_Stderr, "nExtern=%d, ndeg=%d\n", nExtern, ndeg); } ok = False; } @@ -3089,17 +3090,17 @@ SProcess::SProcess(Model *mdl, Process *prc, Options *opts, int sid, int *clst, nEdges = ndeg/2; } else { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); - fprintf(GRCC_Stderr, "illegal total deg = %d (not even)\n", ndeg); + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); + grcc_fprintf(GRCC_Stderr, "illegal total deg = %d (not even)\n", ndeg); } ok = False; } nNodes = nvrt + nExtern; if (nNodes >= GRCC_MAXNODES) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); - fprintf(GRCC_Stderr, "too many nodes = %d\n", nNodes); - fprintf(GRCC_Stderr, " nExtern=%d, nvert=%d (GRCC_MAXNODES)\n", + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: "); + grcc_fprintf(GRCC_Stderr, "too many nodes = %d\n", nNodes); + grcc_fprintf(GRCC_Stderr, " nExtern=%d, nvert=%d (GRCC_MAXNODES)\n", nExtern, nvert); } ok = False; @@ -3114,7 +3115,7 @@ SProcess::SProcess(Model *mdl, Process *prc, Options *opts, int sid, int *clst, lp2 = tCouple - nExtern + 2; if (lp2 % 2 != 0 || lp2 < 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::SProcess: illegal loop: " + grcc_fprintf(GRCC_Stderr, "*** SProcess::SProcess: illegal loop: " "tCouple=%d, nExtern=%d, 2*loop=%d\n", tCouple, nExtern, lp2); } @@ -3143,12 +3144,12 @@ SProcess::~SProcess(void) //-------------------------------------------------------------- void SProcess::prSProcess(void) { - printf("\n"); - printf("+++ Subprocess %d, class=%d\n", id, nclass); + grcc_fprintf(GRCC_Stdout, "\n"); + grcc_fprintf(GRCC_Stdout, "+++ Subprocess %d, class=%d\n", id, nclass); if (pnclass != NULL) { pnclass->prPNodeClass(); } else { - printf(" pnclass = NULL\n"); + grcc_fprintf(GRCC_Stdout, " pnclass = NULL\n"); } } @@ -3227,7 +3228,7 @@ PNodeClass *SProcess::match(MGraph *mgr) if (mgr->nNodes != nNodes) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::match: different nNodes=%d != %d\n", + grcc_fprintf(GRCC_Stderr, "*** SProcess::match: different nNodes=%d != %d\n", nNodes, mgr->nNodes); } erEnd("SProcess::match: different nNodes"); @@ -3249,10 +3250,10 @@ PNodeClass *SProcess::match(MGraph *mgr) mc = mgr->nodes[nd]->clss; if (mc != n2m[nd]) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** SProcess::match: "); - fprintf(GRCC_Stderr, "inconsistent class\n"); + grcc_fprintf(GRCC_Stderr, "*** SProcess::match: "); + grcc_fprintf(GRCC_Stderr, "inconsistent class\n"); for (l = 0; l < nNodes; l++) { - fprintf(GRCC_Stderr, " %d: %d, %d\n", + grcc_fprintf(GRCC_Stderr, " %d: %d, %d\n", j, mgr->nodes[l]->clss, n2m[l]); } } @@ -3378,14 +3379,14 @@ Process::Process(int pid, Model *modl, Options *optn, int nini, int *intlPrt, in lp2 = ctotal - nExtern + 2; if (lp2 % 2 != 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** cannot generate : 2*loop is odd : " + grcc_fprintf(GRCC_Stderr, "*** cannot generate : 2*loop is odd : " "2*loop = %d, ctotal=%d, nExtern=%d\n", lp2, ctotal, nExtern); } ok = False; } else if (lp2 < 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** cannot make a connected graph : " + grcc_fprintf(GRCC_Stderr, "*** cannot make a connected graph : " "2*loop=%d, ctotal=%d, nExtern=%d\n", lp2, ctotal, nExtern); } @@ -3396,7 +3397,7 @@ Process::Process(int pid, Model *modl, Options *optn, int nini, int *intlPrt, in if (!ok) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** Process: illegal input: lp2 = %d\n", lp2); + grcc_fprintf(GRCC_Stderr, "*** Process: illegal input: lp2 = %d\n", lp2); } erEnd("Process: illegal input"); } @@ -3424,7 +3425,7 @@ Process::~Process(void) //-------------------------------------------------------------- void Process::prProcess(void) { - printf("+++ process options : OPI = %d, (Step) = %d, coupling=%d: ", + grcc_fprintf(GRCC_Stdout, "+++ process options : OPI = %d, (Step) = %d, coupling=%d: ", opt->values[GRCC_OPT_1PI], opt->values[GRCC_OPT_Step], ctotal); prIntArray(model->ncouple, clist, "\n"); } @@ -3435,7 +3436,7 @@ void Process::prProcessP(const char *fname) FILE *fp; if ((fp = fopen(fname, "w")) == NULL) { - fprintf(GRCC_Stderr, "*** cannot open \"%s\"\n", fname); + grcc_fprintf(GRCC_Stderr, "*** cannot open \"%s\"\n", fname); return; } fprintf(fp, "Process = {\n"); @@ -3866,27 +3867,27 @@ void MNodeClass::printMat(void) int j1, j2; - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); - printf("nClasses=%d", nClasses); - printf(" clord="); prIntArray(nClasses, clord, ""); - printf(" flist="); prIntArray(nClasses, flist, "\n"); - printf("flg = (%d, %d, %d)\n", flg0, flg1, flg2); + grcc_fprintf(GRCC_Stdout, "nClasses=%d", nClasses); + grcc_fprintf(GRCC_Stdout, " clord="); prIntArray(nClasses, clord, ""); + grcc_fprintf(GRCC_Stdout, " flist="); prIntArray(nClasses, flist, "\n"); + grcc_fprintf(GRCC_Stdout, "flg = (%d, %d, %d)\n", flg0, flg1, flg2); // the first line - printf("nd: cl: "); + grcc_fprintf(GRCC_Stdout, "nd: cl: "); for (j2 = 0; j2 < nClasses; j2++) { - printf("%2d ", j2); + grcc_fprintf(GRCC_Stdout, "%2d ", j2); } - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); // print raw for (j1 = 0; j1 < nNodes; j1++) { - printf("%2d; %2d: [", j1, ndcl[j1]); + grcc_fprintf(GRCC_Stdout, "%2d; %2d: [", j1, ndcl[j1]); for (j2 = 0; j2 < nClasses; j2++) { - printf(" %2d", clmat[j1][j2]); + grcc_fprintf(GRCC_Stdout, " %2d", clmat[j1][j2]); } - printf("]\n"); + grcc_fprintf(GRCC_Stdout, "]\n"); } } @@ -4009,9 +4010,9 @@ MGraph::MGraph(int pid, int ncl, int *cldeg, int *clnum, int *cltyp, int *cmind, } if (ne % 2 != 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "Sum of degrees are not even\n"); + grcc_fprintf(GRCC_Stderr, "Sum of degrees are not even\n"); for (j = 0; j < nClasses; j++) { - fprintf(GRCC_Stderr, "class %2d: %2d %2d %2d\n", + grcc_fprintf(GRCC_Stderr, "class %2d: %2d %2d %2d\n", j, cldeg[j], clnum[j], cltyp[j]); } } @@ -4077,9 +4078,9 @@ MGraph::MGraph(int pid, int ncl, NCInput *mgi, Options *opts) } if (ne % 2 != 0) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "Sum of degrees are not even\n"); + grcc_fprintf(GRCC_Stderr, "Sum of degrees are not even\n"); for (j = 0; j < nClasses; j++) { - fprintf(GRCC_Stderr, "class %2d: %2d %2d %2d\n", + grcc_fprintf(GRCC_Stderr, "class %2d: %2d %2d %2d\n", j, mgi[j].cldeg, mgi[j].clnum, mgi[j].cltyp); } } @@ -4088,7 +4089,7 @@ MGraph::MGraph(int pid, int ncl, NCInput *mgi, Options *opts) pId = pid; nNodes = nn; if (nNodes < 0) { - printf("*** nNodes = %d\n", nNodes); + grcc_fprintf(GRCC_Stdout, "*** nNodes = %d\n", nNodes); erEnd("MGraph::MGrap : nNodes < 0"); } nodes = new MNode*[nNodes]; @@ -4196,17 +4197,17 @@ void MGraph::printAdjMat(MNodeClass *cl) { int j1, j2; - printf(" "); + grcc_fprintf(GRCC_Stdout, " "); for (j2 = 0; j2 < nNodes; j2++) { - printf(" %2d", j2); + grcc_fprintf(GRCC_Stdout, " %2d", j2); } - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); for (j1 = 0; j1 < nNodes; j1++) { - printf("%2d : [", j1); + grcc_fprintf(GRCC_Stdout, "%2d : [", j1); for (j2 = 0; j2 < nNodes; j2++) { - printf(" %2d", adjMat[j1][j2]); + grcc_fprintf(GRCC_Stdout, " %2d", adjMat[j1][j2]); } - printf("] %2d\n", cl->ndcl[j1]); + grcc_fprintf(GRCC_Stdout, "] %2d\n", cl->ndcl[j1]); } } @@ -4215,21 +4216,21 @@ void MGraph::print(void) { int j; - printf("MGraph: pId=%d, cDiag=%ld, c1PI=%ld\n", pId, cDiag, c1PI); - printf(" nNodes=%d, nEdges=%d, nExtern=%d, nLoops=%d " + grcc_fprintf(GRCC_Stdout, "MGraph: pId=%d, cDiag=%ld, c1PI=%ld\n", pId, cDiag, c1PI); + grcc_fprintf(GRCC_Stdout, " nNodes=%d, nEdges=%d, nExtern=%d, nLoops=%d " "mindeg=%d, maxdeg=%d, sym=(%ld, %ld)\n", nNodes, nEdges, nExtern, nLoops, mindeg, maxdeg, nsym, esym); - printf(" Nodes=%d\n", nNodes); + grcc_fprintf(GRCC_Stdout, " Nodes=%d\n", nNodes); for (j = 0; j < nNodes; j++) { - printf(" %2d: id=%d, deg=%d, clss=%d, extloop=%d, ", + grcc_fprintf(GRCC_Stdout, " %2d: id=%d, deg=%d, clss=%d, extloop=%d, ", j, nodes[j]->id, nodes[j]->deg, nodes[j]->clss, nodes[j]->extloop); - printf("mind=%d, maxd=%d, freelg=%d\n", + grcc_fprintf(GRCC_Stdout, "mind=%d, maxd=%d, freelg=%d\n", nodes[j]->cmindeg, nodes[j]->cmaxdeg, nodes[j]->freelg); } curcl->printMat(); - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); printAdjMat(curcl); } @@ -4585,11 +4586,11 @@ void MGraph::biconnME(void) } bool ok = True; if (momset != m) { - printf("*** momset = %ld != %ld\n", momset, m); + grcc_fprintf(GRCC_Stdout, "*** momset = %ld != %ld\n", momset, m); ok = False; } if (mconn->nbacked != nLoops) { - printf("*** nbacked = %d != %c\n", mconn->nbacked, nLoops); + grcc_fprintf(GRCC_Stdout, "*** nbacked = %d != %c\n", mconn->nbacked, nLoops); ok = False; } if (! ok) { @@ -4660,7 +4661,7 @@ void MGraph::bisearchME(int nd, int pd, int ned, int col, newv = (bidef[nd] < 0); if (! newv) { - printf("*** nd=%d, pd=%d, bidef[nd]=%d\n", nd, pd, bidef[nd]); + grcc_fprintf(GRCC_Stdout, "*** nd=%d, pd=%d, bidef[nd]=%d\n", nd, pd, bidef[nd]); erEnd("bisearchME: illegal connection"); } @@ -4786,11 +4787,11 @@ void MGraph::bisearchME(int nd, int pd, int ned, int col, } } if (cn != conn) { - printf("*** revisit back-ed (%d --> %d) cn=%d != conn=%d\n", + grcc_fprintf(GRCC_Stdout, "*** revisit back-ed (%d --> %d) cn=%d != conn=%d\n", nd, td, cn, conn); - printf(" sedges = %d\n", mconn->sedges); + grcc_fprintf(GRCC_Stdout, " sedges = %d\n", mconn->sedges); for (int ed = 0; ed < mconn->sedges; ed++) { - printf(" %d : (%d --> %d) : [%2d*%2ld]\n", ed, + grcc_fprintf(GRCC_Stdout, " %d : (%d --> %d) : [%2d*%2ld]\n", ed, mconn->cedges[ed].nodes[0], mconn->cedges[ed].nodes[1], mconn->cedges[ed].momdir, @@ -5110,7 +5111,7 @@ void MGraph::connectLeg(int so, int sn, int to, int ts, MNodeClass *cl) #ifdef CHECK if ((adjMat[sn][tn] != 0) || (adjMat[sn][tn] != adjMat[tn][sn])) { - printf("*** inconsistent connection: sn=%d, tn=%d", + grcc_fprintf(GRCC_Stdout, "*** inconsistent connection: sn=%d, tn=%d", sn, tn); printAdjMat(cl); erEnd("inconsistent connection "); @@ -5299,9 +5300,9 @@ void MGraph::newGraph(MNodeClass *cl) } #ifdef MONITOR - printf("\n"); - printf("Graph : %ld (%ld)", cDiag, ngen); - printf(" sym. factor = (%ld*%ld)\n", nsym, esym); + grcc_fprintf(GRCC_Stdout, "\n"); + grcc_fprintf(GRCC_Stdout, "Graph : %ld (%ld)", cDiag, ngen); + grcc_fprintf(GRCC_Stdout, " sym. factor = (%ld*%ld)\n", nsym, esym); printAdjMat(cl); // cl->printMat(); # ifdef ORBITS @@ -5344,10 +5345,10 @@ void MOrbits::print(void) { int c; - printf("Orbits : nOrbits=%d: nd2or=", nOrbits); + grcc_fprintf(GRCC_Stdout, "Orbits : nOrbits=%d: nd2or=", nOrbits); prIntArray(nNodes, nd2or, "\n"); for (c = 0; c < nOrbits; c++) { - printf(" %2d: (%2d -- %2d) :", c, flist[c], flist[c+1]-1); + grcc_fprintf(GRCC_Stdout, " %2d: (%2d -- %2d) :", c, flist[c], flist[c+1]-1); prIntArray(flist[c+1]-flist[c], or2nd+flist[c], "\n"); } } @@ -5392,8 +5393,8 @@ void MOrbits::fromPerm(int *perm) } #ifdef CHECK if (j1 >= nNodes) { - printf("*** fromPerm: illegal control: j=%d, k=%d\n", j, k); - printf("perm="); + grcc_fprintf(GRCC_Stdout, "*** fromPerm: illegal control: j=%d, k=%d\n", j, k); + grcc_fprintf(GRCC_Stdout, "perm="); prIntArray(nNodes, perm, " nd2or="); prIntArray(nNodes, nd2or, "\n"); erEnd("fromPerm: illegal control"); @@ -5621,7 +5622,7 @@ void MConn::initCEdges(MGraph *mg) } } if (ed != sedges) { - printf("*** ed=%d != sedges=%d\n", ed, sedges); + grcc_fprintf(GRCC_Stdout, "*** ed=%d != sedges=%d\n", ed, sedges); erEnd("MConn::initCEdge: table overflow"); } } @@ -5752,71 +5753,71 @@ void MConn::print(void) { int j, k; - printf("+++ MConn object: snodes=%d, sedges=%d\n", snodes, sedges); - printf(" nopic=%d, nlpopic=%d, nbacked=%d, nctopic=%d, " + grcc_fprintf(GRCC_Stdout, "+++ MConn object: snodes=%d, sedges=%d\n", snodes, sedges); + grcc_fprintf(GRCC_Stdout, " nopic=%d, nlpopic=%d, nbacked=%d, nctopic=%d, " "nbridges=%d, ne0bridges=%d, ne1bridges=%d\n", nopic, nlpopic, nbacked, nctopic, nbridges, ne0bridges, ne1bridges); - printf(" nblocks=%d, neblocks=%d, na1blocks=%d, narticuls=%d, nselfloop=%d\n", + grcc_fprintf(GRCC_Stdout, " nblocks=%d, neblocks=%d, na1blocks=%d, narticuls=%d, nselfloop=%d\n", nblocks, neblocks, na1blocks, narticuls, nselfloops); - printf(" nmultiedges=%d\n", nmultiedges); + grcc_fprintf(GRCC_Stdout, " nmultiedges=%d\n", nmultiedges); - printf(" cEdges (%d)\n", sedges); + grcc_fprintf(GRCC_Stdout, " cEdges (%d)\n", sedges); if (sedges > 0) { for (j = 0; j < sedges; j++) { - printf(" %2d: (%d,%d)[%2d*%2ld] \n", j, + grcc_fprintf(GRCC_Stdout, " %2d: (%d,%d)[%2d*%2ld] \n", j, cedges[j].nodes[0], cedges[j].nodes[1], cedges[j].momdir, cedges[j].momset); } } - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); - printf(" 1PI components (%d)\n", nopic); + grcc_fprintf(GRCC_Stdout, " 1PI components (%d)\n", nopic); for (j = 0; j < nopic; j++) { - printf(" %d: nleg=%d, nnodes=%d, nedge=%d, ", + grcc_fprintf(GRCC_Stdout, " %d: nleg=%d, nnodes=%d, nedge=%d, ", j, opics[j].nlegs, opics[j].nnodes, opics[j].nedges); - printf("next=%d, loop=%d, ctlp=%d: m0lg=%d [", + grcc_fprintf(GRCC_Stdout, "next=%d, loop=%d, ctlp=%d: m0lg=%d [", opics[j].next, opics[j].loop, opics[j].ctloop, opics[j].mom0lg); for (k = 0; k < opics[j].nnodes; k++) { - printf(" %d", opics[j].nodes[k]); + grcc_fprintf(GRCC_Stdout, " %d", opics[j].nodes[k]); } - printf("]\n"); + grcc_fprintf(GRCC_Stdout, "]\n"); } - printf(" bridges (%d)\n", nbridges); + grcc_fprintf(GRCC_Stdout, " bridges (%d)\n", nbridges); if (nbridges > 0) { - printf(" "); + grcc_fprintf(GRCC_Stdout, " "); for (j = 0; j < nbridges; j++) { - printf("(%d,%d)[mom=%d] ", + grcc_fprintf(GRCC_Stdout, "(%d,%d)[mom=%d] ", bridges[j].nodes[0], bridges[j].nodes[1], bridges[j].next); } - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); } - printf(" blocks (%d)\n", nblocks); + grcc_fprintf(GRCC_Stdout, " blocks (%d)\n", nblocks); for (j = 0; j < nblocks; j++) { - printf(" %d: nmedges=%d, nartps=%d, loop=%d: [", + grcc_fprintf(GRCC_Stdout, " %d: nmedges=%d, nartps=%d, loop=%d: [", j, blocks[j].nmedges, blocks[j].nartps, blocks[j].loop); for (k = 0; k < blocks[j].nmedges; k++) { - printf(" (%d,%d)", blocks[j].edges[k][0], blocks[j].edges[k][1]); + grcc_fprintf(GRCC_Stdout, " (%d,%d)", blocks[j].edges[k][0], blocks[j].edges[k][1]); } - printf("]\n"); + grcc_fprintf(GRCC_Stdout, "]\n"); } - printf(" articulation points (%d)\n", narticuls); + grcc_fprintf(GRCC_Stdout, " articulation points (%d)\n", narticuls); if (narticuls > 0) { - printf(" "); + grcc_fprintf(GRCC_Stdout, " "); for (j = 0; j < snodes; j++) { if (articuls[j] != 0) { #if 0 - printf("%d(%d) ", j, articuls[j]); + grcc_fprintf(GRCC_Stdout, "%d(%d) ", j, articuls[j]); #else - printf("%d ", j); + grcc_fprintf(GRCC_Stdout, "%d ", j); #endif } } - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); } } @@ -5825,17 +5826,17 @@ void MConn::prEdges(void) { int j; - printf(" cEdges (%d)", sedges); + grcc_fprintf(GRCC_Stdout, " cEdges (%d)", sedges); if (sedges > 0) { for (j = 0; j < sedges; j++) { if (j % 5 == 0) { - printf("\n "); + grcc_fprintf(GRCC_Stdout, "\n "); } - printf("(%d,%d)[%2d*%3ld] ", + grcc_fprintf(GRCC_Stdout, "(%d,%d)[%2d*%3ld] ", cedges[j].nodes[0], cedges[j].nodes[1], cedges[j].momdir, cedges[j].momset); } - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); } } @@ -5876,12 +5877,12 @@ void SGroup::print(void) { int j; - printf("SGroup : nnodes = %d, nelem=%ld, cgen=%d, csav=%d\n", + grcc_fprintf(GRCC_Stdout, "SGroup : nnodes = %d, nelem=%ld, cgen=%d, csav=%d\n", nnodes, nelem, cgen, csav); - printf(" eclass ="); + grcc_fprintf(GRCC_Stdout, " eclass ="); prIntArray(neclass, eclass, "\n"); for (j = 0; j < nelem; j++) { - printf(" %4d: ", j); + grcc_fprintf(GRCC_Stdout, " %4d: ", j); prIntArray(nnodes, elem[j], "\n"); } } @@ -6077,7 +6078,7 @@ void ENode::setExtern(int typ, int pt) if (typ == GRCC_AT_Initial || typ == GRCC_AT_Final) { extloop = typ; } else { - fprintf(GRCC_Stderr, "** ENode::setExternal:: illegal typ=%d\n", + grcc_fprintf(GRCC_Stderr, "** ENode::setExternal:: illegal typ=%d\n", typ); erEnd("ENode::setExternal:: illegal typ"); } @@ -6088,7 +6089,7 @@ void ENode::setType(int typ) { #ifdef CHECK if (ndtype != GRCC_ND_Undef && ndtype != typ) { - fprintf(GRCC_Stderr, "*** ndtype is already defined : old=%d, new = %d\n", + grcc_fprintf(GRCC_Stderr, "*** ndtype is already defined : old=%d, new = %d\n", ndtype, typ); erEnd("ndtype is already defined"); } @@ -6101,16 +6102,16 @@ void ENode::print(void) { int j; - printf("Enode %d deg=%d, extl=%2d, intr=%d ", + grcc_fprintf(GRCC_Stdout, "Enode %d deg=%d, extl=%2d, intr=%d ", id, deg, extloop, intrct); if (egraph->bicount > 0) { - printf("(%-8s) ", GRCC_ND_names[ndtype]); + grcc_fprintf(GRCC_Stdout, "(%-8s) ", GRCC_ND_names[ndtype]); } - printf("edge=["); + grcc_fprintf(GRCC_Stdout, "edge=["); for (j = 0; j < deg; j++) { - printf(" %2d", edges[j]); + grcc_fprintf(GRCC_Stdout, " %2d", edges[j]); } - printf("]\n"); + grcc_fprintf(GRCC_Stdout, "]\n"); } //============================================================== @@ -6220,24 +6221,24 @@ void EEdge::print(void) { int zero, n, k; - printf("Edge %2d ext=%d ptcl=%2d ", id, ext, ptcl); + grcc_fprintf(GRCC_Stdout, "Edge %2d ext=%d ptcl=%2d ", id, ext, ptcl); if (egraph == NULL || !egraph->assigned) { - printf("[%d, %d]", nodes[0], nodes[1]); + grcc_fprintf(GRCC_Stdout, "[%d, %d]", nodes[0], nodes[1]); } else { - printf("[(%d,%d), (%d,%d))", nodes[0], nlegs[0], nodes[1], nlegs[1]); + grcc_fprintf(GRCC_Stdout, "[(%d,%d), (%d,%d))", nodes[0], nlegs[0], nodes[1], nlegs[1]); } if (momdir != 0) { - printf("[%2d*%2lu]", momdir, momset); + grcc_fprintf(GRCC_Stdout, "[%2d*%2lu]", momdir, momset); } if (egraph != NULL && egraph->bicount > 0) { if (cut) { - printf(" %-9s ", "Cut"); + grcc_fprintf(GRCC_Stdout, " %-9s ", "Cut"); } else { - printf(" %-9s ", GRCC_ED_names[edtype]); + grcc_fprintf(GRCC_Stdout, " %-9s ", GRCC_ED_names[edtype]); } - printf("c%2d: ", opicomp); - printf("%s%d =", (ext)?"Q":"p", id); + grcc_fprintf(GRCC_Stdout, "c%2d: ", opicomp); + grcc_fprintf(GRCC_Stdout, "%s%d =", (ext)?"Q":"p", id); zero = True; for (n = 0; n < egraph->nEdges; n++) { if (emom[n] != 0) { @@ -6252,16 +6253,16 @@ void EEdge::print(void) } } if (zero) { - printf(" 0"); + grcc_fprintf(GRCC_Stdout, " 0"); } } else { if (egraph == NULL) { - printf(" egraph=NULL"); + grcc_fprintf(GRCC_Stdout, " egraph=NULL"); } else { - printf(" egraph->bicount=%d", egraph->bicount); + grcc_fprintf(GRCC_Stdout, " egraph->bicount=%d", egraph->bicount); } } - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); } //============================================================== @@ -6278,22 +6279,22 @@ EFLine::EFLine(void) void EFLine::print(const char *msg) { if (ftype == FL_Open) { - printf(" Open"); + grcc_fprintf(GRCC_Stdout, " Open"); } else if (ftype == FL_Closed) { - printf(" Loop"); + grcc_fprintf(GRCC_Stdout, " Loop"); } else { - printf(" ?%d", ftype); + grcc_fprintf(GRCC_Stdout, " ?%d", ftype); } if (fkind == GRCC_PT_Dirac) { - printf(" Dirac"); + grcc_fprintf(GRCC_Stdout, " Dirac"); } else if (fkind == GRCC_PT_Majorana) { - printf(" Major"); + grcc_fprintf(GRCC_Stdout, " Major"); } else if (fkind == GRCC_PT_Ghost) { - printf(" Ghost"); + grcc_fprintf(GRCC_Stdout, " Ghost"); } else { - printf(" ?%d", fkind); + grcc_fprintf(GRCC_Stdout, " ?%d", fkind); } - printf(" len=%d ", nlist); + grcc_fprintf(GRCC_Stdout, " len=%d ", nlist); prIntArray(nlist, elist, msg); } @@ -6462,20 +6463,20 @@ void EGraph::fromDGraph(DGraph *dg) int n0, n1, e; if (dg->nnodes > GRCC_MAXNODES) { - fprintf(GRCC_Stderr, "*** too many nodes (GRCC_MAXNODES)\n"); - exit(1); + grcc_fprintf(GRCC_Stderr, "*** too many nodes (GRCC_MAXNODES)\n"); + GRCC_ABORT(); } if (dg->nedges > GRCC_MAXEDGES) { - fprintf(GRCC_Stderr, "*** too many edges (GRCC_MAXEDGES)\n"); - exit(1); + grcc_fprintf(GRCC_Stderr, "*** too many edges (GRCC_MAXEDGES)\n"); + GRCC_ABORT(); } for (e = 0; e < dg->nedges; e++) { if (dg->edges[e][0] >= dg->nnodes || dg->edges[e][0] >= dg->nnodes) { - fprintf(GRCC_Stderr, "*** undefined node:"); - fprintf(GRCC_Stderr, "edge[%d] = {%d, %d}\n", + grcc_fprintf(GRCC_Stderr, "*** undefined node:"); + grcc_fprintf(GRCC_Stderr, "edge[%d] = {%d, %d}\n", e, dg->edges[e][0], dg->edges[e][0]); - exit(1); + GRCC_ABORT(); } } @@ -6493,7 +6494,7 @@ void EGraph::fromDGraph(DGraph *dg) if (deg[n0] == 1) { nextern++; } if (deg[n0] < 0) { - fprintf(GRCC_Stderr, "+++ node %d is isolated\n", n0); + grcc_fprintf(GRCC_Stderr, "+++ node %d is isolated\n", n0); } } @@ -6577,7 +6578,7 @@ void EGraph::fromMGraph(MGraph *mg) erEnd("EGraph::fromMGraph: sizes are too small"); } if (sLoops < mg->nLoops) { - printf("too small sLoops : %d < %d\n", sLoops, mg->nLoops); + grcc_fprintf(GRCC_Stdout, "too small sLoops : %d < %d\n", sLoops, mg->nLoops); erEnd("too small sLoops"); } #endif @@ -6652,7 +6653,7 @@ void EGraph::fromMGraph(MGraph *mg) } #ifdef CHECK if (ed != nEdges) { - printf("*** EGraph::init: ed=%d != nEdges=%d\n", + grcc_fprintf(GRCC_Stdout, "*** EGraph::init: ed=%d != nEdges=%d\n", ed, nEdges+1); erEnd("EGraph::init: illegal connection"); } @@ -6734,7 +6735,7 @@ void EGraph::fromMGraph(MGraph *mg) } if (! found) { - printf("*** EGraph::fromMGraph:edge (%d->%d) is not found\n", + grcc_fprintf(GRCC_Stdout, "*** EGraph::fromMGraph:edge (%d->%d) is not found\n", m0, m1); } } @@ -6751,7 +6752,7 @@ ENode *EGraph::setExtern(int n0, int pt, int ndtyp) nd = nodes[n0]; if (nd->deg != 1) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** illegal external particle %d, %d, %d\n", + grcc_fprintf(GRCC_Stderr, "*** illegal external particle %d, %d, %d\n", n0,pt,ndtyp); } erEnd("illegal external particle"); @@ -6765,7 +6766,7 @@ ENode *EGraph::setExtern(int n0, int pt, int ndtyp) } else { npt = 0; if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** illegal type of an external particle : " + grcc_fprintf(GRCC_Stderr, "*** illegal type of an external particle : " "node = %d, particle = %d, type = %d", n0, pt, ndtyp); } erEnd("illegal type of an external particle"); @@ -6797,58 +6798,58 @@ void EGraph::print(void) int nd, ed, nlp; nlp = nEdges - nNodes + 1; - printf("\nEGraph\n"); - printf(" pId=%d, gSubId=%ld, mId=%ld, aId=%ld, sId=%ld\n", + grcc_fprintf(GRCC_Stdout, "\nEGraph\n"); + grcc_fprintf(GRCC_Stdout, " pId=%d, gSubId=%ld, mId=%ld, aId=%ld, sId=%ld\n", pId, gSubId, mId, aId, sId); - printf(" sNodes=%d, sEdges=%d, sMaxdeg=%d, sLoops=%d\n", + grcc_fprintf(GRCC_Stdout, " sNodes=%d, sEdges=%d, sMaxdeg=%d, sLoops=%d\n", sNodes, sEdges, sMaxdeg, sLoops); - printf(" nNodes=%d, nEdges=%d, nExtern=%d, nLoops=%d, totalc=%d\n", + grcc_fprintf(GRCC_Stdout, " nNodes=%d, nEdges=%d, nExtern=%d, nLoops=%d, totalc=%d\n", nNodes, nEdges, nExtern, nlp, totalc); - printf(" "); + grcc_fprintf(GRCC_Stdout, " "); if (model == NULL) { - printf("model=NULL,"); + grcc_fprintf(GRCC_Stdout, "model=NULL,"); } else { - printf("model=\"%s\",", model->name); + grcc_fprintf(GRCC_Stdout, "model=\"%s\",", model->name); } if (proc == NULL) { - printf("proc=NULL,"); + grcc_fprintf(GRCC_Stdout, "proc=NULL,"); } else { - printf("proc=%d,", proc->id); + grcc_fprintf(GRCC_Stdout, "proc=%d,", proc->id); } if (sproc == NULL) { - printf("sproc=NULL,"); + grcc_fprintf(GRCC_Stdout, "sproc=NULL,"); } else { - printf("sproc=%d,", sproc->id); + grcc_fprintf(GRCC_Stdout, "sproc=%d,", sproc->id); } - printf("\n"); - printf(" assigned=%d, sym = (%ld * %ld) ", + grcc_fprintf(GRCC_Stdout, "\n"); + grcc_fprintf(GRCC_Stdout, " assigned=%d, sym = (%ld * %ld) ", assigned, nsym, esym); - printf("extperm=%ld, nsym1=%ld, multp=%ld\n", + grcc_fprintf(GRCC_Stdout, "extperm=%ld, nsym1=%ld, multp=%ld\n", extperm, nsym1, multp); - printf(" Nodes\n"); + grcc_fprintf(GRCC_Stdout, " Nodes\n"); for (nd = 0; nd < nNodes; nd++) { if (isExternal(nd)) { - printf(" %2d Extern ", nd); + grcc_fprintf(GRCC_Stdout, " %2d Extern ", nd); } else { - printf(" %2d Vertex ", nd); + grcc_fprintf(GRCC_Stdout, " %2d Vertex ", nd); } nodes[nd]->print(); } - printf(" Edges\n"); + grcc_fprintf(GRCC_Stdout, " Edges\n"); for (ed = 0; ed < nEdges; ed++) { if (edges[ed]->ext) { - printf(" %2d Extern ", ed); + grcc_fprintf(GRCC_Stdout, " %2d Extern ", ed); } else { - printf(" %2d Intern ", ed); + grcc_fprintf(GRCC_Stdout, " %2d Intern ", ed); } edges[ed]->print(); } if (bicount > 0) { - printf(" Biconn: nopicomp=%d, nopi2p=%d, opi2plp=%d, nadj2ptv=%d\n", + grcc_fprintf(GRCC_Stdout, " Biconn: nopicomp=%d, nopi2p=%d, opi2plp=%d, nadj2ptv=%d\n", nopicomp, nopi2p, opi2plp, nadj2ptv); } - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); if (nFlines > 0) { prFLines(); @@ -6915,10 +6916,10 @@ void EGraph::prFLines(void) { int j; - printf(" Fermion lines %d, sign=%d (mId=%ld, aId=%ld)\n", + grcc_fprintf(GRCC_Stdout, " Fermion lines %d, sign=%d (mId=%ld, aId=%ld)\n", nFlines, fsign, mId, aId); for (j = 0; j < nFlines; j++) { - printf("%4d ", j); + grcc_fprintf(GRCC_Stdout, "%4d ", j); flines[j]->print("\n"); } } @@ -7025,12 +7026,12 @@ Bool EGraph::optQGrafM(Options *opt) mext = (~ mext) << nExtern; #ifdef PRINT - printf("optQGrafM:"); + grcc_fprintf(GRCC_Stdout, "optQGrafM:"); print(); #endif #ifdef PRINT - printf("optQGrafM: %8ld\n", mId); + grcc_fprintf(GRCC_Stdout, "optQGrafM: %8ld\n", mId); econn->print(); #endif @@ -7040,8 +7041,8 @@ Bool EGraph::optQGrafM(Options *opt) maxlegs = Max(maxlegs, econn->opics[j].nlegs); } if (maxlegs >= GRCC_MAXEDGES) { - printf("*** table overflow\n"); - abort(); + grcc_fprintf(GRCC_Stdout, "*** table overflow\n"); + GRCC_ABORT(); } for (int k = 0; k < GRCC_MAXEDGES; k++) { nopis[k] = 0; @@ -7239,7 +7240,7 @@ Bool EGraph::optQGrafM(Options *opt) Bool EGraph::optQGrafA(Options *opt) { #ifdef PRINT - printf("optQGrafA: %8ld\n", mId); + grcc_fprintf(GRCC_Stdout, "optQGrafA: %8ld\n", mId); econn->print(); #endif Bool retval = True; @@ -7688,9 +7689,9 @@ void EGraph::bisearchE(int nd, int *extlst, int *intlst, int *opiext, int *opilo #ifdef CHECK } else if (nodes[nd]->ndtype != GRCC_ND_CPoint) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "bisearch: node %d is a cut point ", + grcc_fprintf(GRCC_Stderr, "bisearch: node %d is a cut point ", nd); - fprintf(GRCC_Stderr, "(not undef %d)\n", + grcc_fprintf(GRCC_Stderr, "(not undef %d)\n", nodes[nd]->ndtype); } #endif @@ -7706,8 +7707,8 @@ void EGraph::bisearchE(int nd, int *extlst, int *intlst, int *opiext, int *opilo #ifdef CHECK } else if (edges[ed]->edtype != GRCC_ED_Bridge) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "bisearch: edges %d is a bridge ", ed); - fprintf(GRCC_Stderr, "(not undef %d)\n", edges[ed]->edtype); + grcc_fprintf(GRCC_Stderr, "bisearch: edges %d is a bridge ", ed); + grcc_fprintf(GRCC_Stderr, "(not undef %d)\n", edges[ed]->edtype); } #endif } @@ -7745,8 +7746,8 @@ void EGraph::bisearchE(int nd, int *extlst, int *intlst, int *opiext, int *opilo #ifdef CHECK } else if (edges[ed]->edtype != GRCC_ED_Inloop) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "bisearch: "); - fprintf(GRCC_Stderr, "edges %d is not undef (%d)\n", + grcc_fprintf(GRCC_Stderr, "bisearch: "); + grcc_fprintf(GRCC_Stderr, "edges %d is not undef (%d)\n", ed, edges[ed]->edtype); } #endif @@ -7842,7 +7843,7 @@ void EGraph::chkMomConsv(void) for (ex = 0; ex < nEdges; ex++) { if (esum[ex] != 0) { okn = False; - fprintf(GRCC_Stderr, "chkMomConsv:n=%d, esum[%d]=%d\n", + grcc_fprintf(GRCC_Stderr, "chkMomConsv:n=%d, esum[%d]=%d\n", n, ex, esum[ex]); } } @@ -7850,7 +7851,7 @@ void EGraph::chkMomConsv(void) for (lk = 0; lk < nLoops; lk++) { if (lsum[lk] != 0) { okn = False; - fprintf(GRCC_Stderr, "chkMomConsv:n=%d, lsum[%d]=%d\n", + grcc_fprintf(GRCC_Stderr, "chkMomConsv:n=%d, lsum[%d]=%d\n", n, lk, lsum[lk]); } @@ -7858,8 +7859,8 @@ void EGraph::chkMomConsv(void) if (!okn) { ok = False; - fprintf(GRCC_Stderr, "*** Violation of momentum conservation "); - fprintf(GRCC_Stderr, "at node =%d\n",n); + grcc_fprintf(GRCC_Stderr, "*** Violation of momentum conservation "); + grcc_fprintf(GRCC_Stderr, "at node =%d\n",n); } } @@ -7908,7 +7909,7 @@ int EGraph::fltrace(int fk, int nd0, int *fl) // maximal possible length of a fline is nEdges for (k = 0; k < nEdges; k++) { if (k >= nfl) { - printf("*** fltrace:illegal contorl: k=%d, nEdges=%d\n", k, nEdges); + grcc_fprintf(GRCC_Stdout, "*** fltrace:illegal contorl: k=%d, nEdges=%d\n", k, nEdges); break; } @@ -7918,7 +7919,7 @@ int EGraph::fltrace(int fk, int nd0, int *fl) el = V2Ileg(e); #ifdef CHECK if (ed > nEdges) { - fprintf(GRCC_Stderr, "*** fltrace: ed=%d > nEdges=%d, fl=", + grcc_fprintf(GRCC_Stderr, "*** fltrace: ed=%d > nEdges=%d, fl=", ed, nEdges); prIntArray(nNodes, fl, "\n"); erEnd("fltrace: illegal fl"); @@ -7975,19 +7976,19 @@ int EGraph::fltrace(int fk, int nd0, int *fl) } // not visited } // end of for k if (fgcnt == 0) { - printf("*** fline: Fermion number is not conserved\n"); - printf(" nd=%d, e=%d, ed=%d, fgcnt=%d\n", + grcc_fprintf(GRCC_Stdout, "*** fline: Fermion number is not conserved\n"); + grcc_fprintf(GRCC_Stdout, " nd=%d, e=%d, ed=%d, fgcnt=%d\n", nd, e, ed, fgcnt); // erEnd("fline: Fermion number is not conserved"); } else if (fgcnt > 1) { - printf("+++ fline: more than two fermions: check fsign\n"); - printf(" nd=%d, e=%d, ed=%d, fgcnt=%d\n", + grcc_fprintf(GRCC_Stdout, "+++ fline: more than two fermions: check fsign\n"); + grcc_fprintf(GRCC_Stdout, " nd=%d, e=%d, ed=%d, fgcnt=%d\n", nd, e, ed, fgcnt); } } if (k >= nEdges || k >= nfl) { - printf("*** fline: illegal control\n"); - printf(" nEdges=%d, nfl=%d, k=%d, ", nEdges, nfl, k); + grcc_fprintf(GRCC_Stdout, "*** fline: illegal control\n"); + grcc_fprintf(GRCC_Stdout, " nEdges=%d, nfl=%d, k=%d, ", nEdges, nfl, k); prIntArray(nfl, fl, "\n"); erEnd("fline: illegal control"); } @@ -8153,7 +8154,7 @@ NCand::~NCand(void) //--------------------------------------------------------------- void NCand::prNCand(const char* msg) { - printf("%d %d ", st, deg); + grcc_fprintf(GRCC_Stdout, "%d %d ", st, deg); prIntArray(nilist, ilist, msg); } @@ -8171,7 +8172,7 @@ ECand::ECand(int dt, int nplst, int *plst) if (det && nplist != 1) { if (prlevel > 0) { - fprintf(GRCC_Stderr, "*** ECand : len(plist) != 1 : det=%d ", det); + grcc_fprintf(GRCC_Stderr, "*** ECand : len(plist) != 1 : det=%d ", det); prIntArrayErr(nplist, plist, "\n"); } erEnd("ECand : len(plist) != 1"); @@ -8187,7 +8188,7 @@ ECand::~ECand(void) void ECand::prECand(const char *msg) { prIntArray(nplist, plist, ""); - printf(" (det=%d)%s", det, msg); + grcc_fprintf(GRCC_Stdout, " (det=%d)%s", det, msg); } //=============================================================== @@ -8236,7 +8237,7 @@ int ANode::newleg(void) nlegs++; #ifdef CHECK if (nlegs > deg) { - fprintf(GRCC_Stderr, "*** ANode::newleg : nlegs = %d > deg = %d\n", + grcc_fprintf(GRCC_Stderr, "*** ANode::newleg : nlegs = %d > deg = %d\n", nlegs, deg); erEnd("ANode::newleg : nlegs > deg"); } @@ -8369,28 +8370,28 @@ void Assign::prCand(const char *msg) int n, e, ne; AEdge *ed; - printf("\n"); - printf("+++ Candidate list: %s\n", msg); - printf(" Nodes %d\n", nNodes); + grcc_fprintf(GRCC_Stdout, "\n"); + grcc_fprintf(GRCC_Stdout, "+++ Candidate list: %s\n", msg); + grcc_fprintf(GRCC_Stdout, " Nodes %d\n", nNodes); for (n = 0; n < nNodes; n++) { - printf("%d: edges=", n); + grcc_fprintf(GRCC_Stdout, "%d: edges=", n); prIntArray(nodes[n]->deg, nodes[n]->aedges, ": cand="); if (nodes[n]->cand == NULL) { - printf("NULL\n"); + grcc_fprintf(GRCC_Stdout, "NULL\n"); } else { nodes[n]->cand->prNCand("\n"); } } ne = Min(nEdges, nETotal); - printf(" Edges %d\n", ne); + grcc_fprintf(GRCC_Stdout, " Edges %d\n", ne); for (e = 0; e < ne; e++) { ed = edges[e]; if (ed == NULL) { - printf("NULL_Edge\n"); + grcc_fprintf(GRCC_Stdout, "NULL_Edge\n"); } else { - printf("%d: %d->%d: cand=", e, ed->nodes[0], ed->nodes[1]); + grcc_fprintf(GRCC_Stdout, "%d: %d->%d: cand=", e, ed->nodes[0], ed->nodes[1]); if (edges[e]->cand == NULL) { - printf("NULL\n"); + grcc_fprintf(GRCC_Stdout, "NULL\n"); } else { edges[e]->cand->prECand("\n"); } @@ -8407,7 +8408,7 @@ void Assign::checkAG(const char *msg) for (n = 0; n < nNodes; n++) { for (lg = 0; lg < nodes[n]->deg; lg++) { if (nodes[n]->aedges[lg] < 0) { - printf("*** checkAG:%s: n=%d, lg=%d, aedges=%d\n", + grcc_fprintf(GRCC_Stdout, "*** checkAG:%s: n=%d, lg=%d, aedges=%d\n", msg, n, lg, nodes[n]->aedges[lg]); ok = False; } @@ -8863,7 +8864,7 @@ Bool Assign::fromMGraph(void) #ifdef CHECK if (mgraph->nodes[n]->deg != 1) { - printf("*** assign:fromMGraph : " + grcc_fprintf(GRCC_Stdout, "*** assign:fromMGraph : " "external but deg[%d] = %d != 1, type=%d\n", n, mgraph->nodes[n]->deg, typ); mgraph->print(); @@ -8930,7 +8931,7 @@ Bool Assign::fromMGraph(void) } #ifdef CHECK if (nETotal != nEdges) { - printf("*** Assign::fromMGraph nETotal=%d != nEdges=%d\n", + grcc_fprintf(GRCC_Stdout, "*** Assign::fromMGraph nETotal=%d != nEdges=%d\n", nETotal, nEdges); erEnd("Assign::fromMGraph nETotal= != nEdges"); } @@ -8986,7 +8987,7 @@ void Assign::addEdge(int n0, int n1, int nplist, int *plist) #ifdef CHECK if (n0 >= nNodes || n1 >= nNodes) { - printf("*** Assign::addEdge : undefined nodes %d: [%d, %d]", + grcc_fprintf(GRCC_Stdout, "*** Assign::addEdge : undefined nodes %d: [%d, %d]", nETotal, n0, n1); erEnd("Assign::addEdge : undefined nodes"); } @@ -9091,7 +9092,7 @@ Bool Assign::fillEGraph(int aid, BigInt nsym, BigInt esym, BigInt nsym1) if (isATExternal(pnclass->type[cl])) { ; } else if (nodes[n]->cand->st != AS_Assigned) { - printf("*** fillEGraph : node %d is not assigned", n); + grcc_fprintf(GRCC_Stdout, "*** fillEGraph : node %d is not assigned", n); prCand("fillEGraph: node"); erEnd("fillEGraph : node is not assigned"); } @@ -9112,7 +9113,7 @@ Bool Assign::fillEGraph(int aid, BigInt nsym, BigInt esym, BigInt nsym1) } #ifdef CHECK if (lg < 0 || lg >= nodes[n]->deg) { - printf("*** fillEGraph: n=%d, lr=%d: 0 <= lg=%d < %d\n", + grcc_fprintf(GRCC_Stdout, "*** fillEGraph: n=%d, lr=%d: 0 <= lg=%d < %d\n", n, lr, lg, nodes[n]->deg); erEnd("fillEGraph: illegal reordering"); } @@ -9130,7 +9131,7 @@ Bool Assign::fillEGraph(int aid, BigInt nsym, BigInt esym, BigInt nsym1) for (e = 0; e < nEdges; e++) { #ifdef CHECK if (edges[e]->cand->nplist != 1) { - printf("*** fillEGraph : edge %d is not assigned", e); + grcc_fprintf(GRCC_Stdout, "*** fillEGraph : edge %d is not assigned", e); prCand("fillEGraph: edge"); erEnd("fillEGraph : edge is not assigned"); } @@ -9160,14 +9161,14 @@ Bool Assign::fillEGraph(int aid, BigInt nsym, BigInt esym, BigInt nsym1) n = egraph->edges[ed]->nodes[0]; lr = egraph->edges[ed]->nlegs[0]; if (egraph->nodes[n]->edges[lr] != -ed-1) { - fprintf(GRCC_Stderr, "+++ node[%d][%d]=%d != - (edge[%d][0] + 1) = %d\n", + grcc_fprintf(GRCC_Stderr, "+++ node[%d][%d]=%d != - (edge[%d][0] + 1) = %d\n", n, lr, egraph->nodes[n]->edges[lr], e, -ed-1); ok = False; } n = egraph->edges[ed]->nodes[1]; lr = egraph->edges[ed]->nlegs[1]; if (egraph->nodes[n]->edges[lr] != ed+1) { - fprintf(GRCC_Stderr, "+++ node[%d][%d]=%d != + (edge[%d][0] + 1) = %d\n", + grcc_fprintf(GRCC_Stderr, "+++ node[%d][%d]=%d != + (edge[%d][0] + 1) = %d\n", n, lr, egraph->nodes[n]->edges[lr], e, ed+1); ok = False; } @@ -9247,10 +9248,10 @@ int *Assign::reordLeg(int n, int *reord, int *plist, int *used) #ifdef CHECK if (!found) { - printf("*** reordLeg: illegal list of particles:" + grcc_fprintf(GRCC_Stdout, "*** reordLeg: illegal list of particles:" "interaction %d ", ia); prIntArray(deg, ilegs, "; "); - printf("vertex %d ", n); + grcc_fprintf(GRCC_Stdout, "vertex %d ", n); prIntArray(deg, plist, "\n"); prCand("reordLeg"); erEnd("reordLeg: illegal list of particles"); @@ -9340,7 +9341,7 @@ int Assign::candPart(int v, int ln, int *plist, const int size) #ifdef CHECK if (edges[en]->cand->det) { - printf("*** candPart : particle of leg (%d, %d) " + grcc_fprintf(GRCC_Stdout, "*** candPart : particle of leg (%d, %d) " "is assigned to %d\n", v, ln, edges[en]->cand->plist[0]); checkCand("candPart"); @@ -9433,8 +9434,8 @@ int Assign::selUnAssLeg(int v, int lastlg) } n1 = nodes[v]->anodes[lg]; if (n0 > n1) { - printf("*** selUnAssLeg: n0=%d > n1=%d\n", n0, n1); - printf("*** illegal connection\n"); + grcc_fprintf(GRCC_Stdout, "*** selUnAssLeg: n0=%d > n1=%d\n", n0, n1); + grcc_fprintf(GRCC_Stdout, "*** illegal connection\n"); erEnd("selUnAssLeg: n0 > n1"); } #endif @@ -9526,7 +9527,7 @@ Bool Assign::assignPLeg(int n, int ln, int pt) #ifdef CHECK if (!isIn(edges[e]->cand->nplist, edges[e]->cand->plist, ept)) { - printf("*** assignPLeg: particle %d is not in the cand. of e=%d", + grcc_fprintf(GRCC_Stdout, "*** assignPLeg: particle %d is not in the cand. of e=%d", ept, e); edges[e]->cand->prECand("\n"); prCand("assignPLeg"); @@ -9672,7 +9673,7 @@ Bool Assign::updateCandNode(int v) #ifdef CHECK e = nodes[v]->aedges[0]; if (e < 0 || edges[e]->cand->nplist != 1) { - printf("*** illegal external node: v=%d e=%d :", v, e); + grcc_fprintf(GRCC_Stdout, "*** illegal external node: v=%d e=%d :", v, e); prCand("updateCandNode"); erEnd("illegal external node"); } @@ -9881,7 +9882,7 @@ Bool Assign::isIsomorphic(MNodeClass *cl, BigInt *nsym, BigInt *esym, BigInt *ns ngelem = mgraph->group->nElem(); #ifdef CHECK if (mgraph->nsym > 1 && ngelem <= 1) { - printf("*** isIsomorphic: illegal group: " + grcc_fprintf(GRCC_Stdout, "*** isIsomorphic: illegal group: " "ngelem=%ld, mgraph->sym=(%ld, %ld)\n", ngelem, mgraph->nsym, mgraph->esym); erEnd("Assign::isIsomorphic: illegal group"); @@ -10123,7 +10124,7 @@ Bool Assign::checkCand(const char *msg) // check assigned vertex } else if (nc->st == AS_Assigned) { if (nc->nilist < 1) { - printf("*** checkCand:7:%s:status (%d) of node %d says" + grcc_fprintf(GRCC_Stdout, "*** checkCand:7:%s:status (%d) of node %d says" " interaction is assigned to %d but ilist=", msg, nc->st, n, nc->st); prIntArray(nc->nilist, nc->ilist, "\n"); @@ -10134,7 +10135,7 @@ Bool Assign::checkCand(const char *msg) e = na->aedges[lg]; ec = edges[e]->cand; if (ec->nplist != 1) { - printf("*** checkCand:8:%s:status (%d) of node %d says" + grcc_fprintf(GRCC_Stdout, "*** checkCand:8:%s:status (%d) of node %d says" " interaction is assigned " " but unassigned edge %d is found\n", msg, nc->st, n, e); @@ -10147,7 +10148,7 @@ Bool Assign::checkCand(const char *msg) // pt = nc->ilist[0]; // *** pte = legEdgeParticle(n, 0, - pt); } else { - printf("*** checkCand:10:%s:illegal status of node %d : %d", + grcc_fprintf(GRCC_Stdout, "*** checkCand:10:%s:illegal status of node %d : %d", msg, n, nc->st); ok = False; } @@ -10158,15 +10159,15 @@ Bool Assign::checkCand(const char *msg) for (e = 0; e < nEdges; e++) { ec = edges[e]->cand; if (ec != NULL && ec->nplist < 1) { - printf("*** checkCand:12:%s:illegal edge %d\n", msg, e); + grcc_fprintf(GRCC_Stdout, "*** checkCand:12:%s:illegal edge %d\n", msg, e); ok = False; } } if (!ok) { - printf("*** checkCand:15:%s:illegal configuration\n", msg); + grcc_fprintf(GRCC_Stdout, "*** checkCand:15:%s:illegal configuration\n", msg); prCand("checkCand"); - printf("*** checkCand:16:illegal configuration\n"); + grcc_fprintf(GRCC_Stdout, "*** checkCand:16:illegal configuration\n"); erEnd("checkCand:16:illegal configuration"); } return ok; @@ -10180,9 +10181,9 @@ void Assign::checkNode(int n, const char *msg) for (j = 0; j < nodes[n]->cand->nilist; j++) { it = nodes[n]->cand->ilist[j]; if (Abs(it) >= GRCC_MAXMINTERACT) { - printf("*** %s: n=%d, j=%d, it=%d\n", msg, n, j, it); + grcc_fprintf(GRCC_Stdout, "*** %s: n=%d, j=%d, it=%d\n", msg, n, j, it); nodes[n]->cand->prNCand(msg); - printf("\n"); + grcc_fprintf(GRCC_Stdout, "\n"); erEnd("checkNode:illegal it"); } } @@ -10200,14 +10201,14 @@ void Assign::checkNode(int n, const char *msg) //-------------------------------------------------------------- void NStack::print(const char *msg) { - printf(" node=%d, deg=%d, st=%d, ilist=", noden, deg, st); + grcc_fprintf(GRCC_Stdout, " node=%d, deg=%d, st=%d, ilist=", noden, deg, st); prilist(nilist, ilist, msg); } //-------------------------------------------------------------- void EStack::print(const char *msg) { - printf(" edge=%d, det=%d, plist=", edgen, det); + grcc_fprintf(GRCC_Stdout, " edge=%d, det=%d, plist=", edgen, det); prilist(nplist, plist, msg); } @@ -10398,7 +10399,7 @@ void AStack::restoreMsg(CheckPt sav, const char *msg) restore(sav); if (!agraph->checkCand("restore")) { - printf("restore is called from %s\n", msg); + grcc_fprintf(GRCC_Stdout, "restore is called from %s\n", msg); } } #endif @@ -10408,13 +10409,13 @@ void AStack::prStack(void) { int j; - printf("+++ prStack : (%d, %d)", nStackP, eStackP); + grcc_fprintf(GRCC_Stdout, "+++ prStack : (%d, %d)", nStackP, eStackP); for (j = 0; j < nStackP; j++) { - printf("N:%4d ", j); + grcc_fprintf(GRCC_Stdout, "N:%4d ", j); nStack[j]->print("\n"); } for (j = 0; j < eStackP; j++) { - printf("E:%4d ", j); + grcc_fprintf(GRCC_Stdout, "E:%4d ", j); eStack[j]->print("\n"); } } @@ -10438,9 +10439,9 @@ void Fraction::print(const char *msg) double err = Abs(Real(num)/Real(den) - ratio); if (err > GRCC_FRACERROR) { - printf("%ld/%ld(%g)(overflow)%s", num, den, ratio, msg); + grcc_fprintf(GRCC_Stdout, "%ld/%ld(%g)(overflow)%s", num, den, ratio, msg); } else { - printf("%ld/%ld(%g)%s", num, den, ratio, msg); + grcc_fprintf(GRCC_Stdout, "%ld/%ld(%g)%s", num, den, ratio, msg); } } @@ -10547,12 +10548,43 @@ Bool Fraction::isEq(Fraction f) //************************************************************** // common.cc //============================================================== + +// Wrapper function for printing messages. This allows the use +// of FORM MesPrint when compiled as part of FORM, and the +// usual fprintf to GRCC_Stdout or GRCC_Stderr otherwise. +static void grcc_fprintf(FILE* out, const char* fmt, ...) +{ + va_list args; + va_start(args, fmt); + +#ifndef NOFORM + DUMMYUSE(out); + // the second call of vsnprintf requires a copy of args + va_list args_copy; + va_copy(args_copy, args); + // determine the required buffer size for the formatted string: + const int len = vsnprintf(NULL, 0, fmt, args); + char *buffer = new char[len+1]; + vsnprintf(buffer, len+1, fmt, args_copy); + MLOCK(ErrorMessageLock); + MesPrint("%s", buffer); + MUNLOCK(ErrorMessageLock); + delete[] buffer; + va_end(args_copy); +#else + vfprintf(out, fmt, args); +#endif + + va_end(args); + return; +} + static void erEnd(const char *msg) { if (erExit != NULL) { (*erExit)(msg, erExitArg); } - fprintf(GRCC_Stderr, "*** Error : %s\n\n", msg); + grcc_fprintf(GRCC_Stderr, "*** Error : %s\n\n", msg); GRCC_ABORT(); } @@ -10631,12 +10663,12 @@ static int *delintdup(int *a) //------------------------------------------------------------ static void prilist(int n, const int *a, const char *msg) { - printf("["); + grcc_fprintf(GRCC_Stdout, "["); for (int j = 0; j < n; j++) { - if (j!=0) printf(", "); - printf("%d", a[j]); + if (j!=0) grcc_fprintf(GRCC_Stdout, ", "); + grcc_fprintf(GRCC_Stdout, "%d", a[j]); } - printf("]%s", msg); + grcc_fprintf(GRCC_Stdout, "]%s", msg); } //------------------------------------------------------------ @@ -10793,13 +10825,13 @@ static void prMomStr(int mom, const char *ms, int mn) if (mom == 0) { return; } else if (mom == 1) { - printf(" + %s%d", ms, mn); + grcc_fprintf(GRCC_Stdout, " + %s%d", ms, mn); } else if (mom > 0) { - printf(" + %d*%s%d", mom, ms, mn); + grcc_fprintf(GRCC_Stdout, " + %d*%s%d", mom, ms, mn); } else if (mom == -1) { - printf(" - %s%d", ms, mn); + grcc_fprintf(GRCC_Stdout, " - %s%d", ms, mn); } else { - printf(" - %d*%s%d", -mom, ms, mn); + grcc_fprintf(GRCC_Stdout, " - %d*%s%d", -mom, ms, mn); } } @@ -10809,12 +10841,12 @@ static void prIntArray(int n, int *p, const char *msg) int j; - printf("["); + grcc_fprintf(GRCC_Stdout, "["); for (j = 0; j < n; j++) { - if (j!=0) printf(", "); - printf("%2d", p[j]); + if (j!=0) grcc_fprintf(GRCC_Stdout, ", "); + grcc_fprintf(GRCC_Stdout, "%2d", p[j]); } - printf("]%s", msg); + grcc_fprintf(GRCC_Stdout, "]%s", msg); } //-------------------------------------------------------------- @@ -10823,12 +10855,12 @@ static void prIntArrayErr(int n, int *p, const char *msg) int j; - fprintf(GRCC_Stderr, "["); + grcc_fprintf(GRCC_Stderr, "["); for (j = 0; j < n; j++) { - if (j!=0) fprintf(GRCC_Stderr, ", "); - fprintf(GRCC_Stderr, "%2d", p[j]); + if (j!=0) grcc_fprintf(GRCC_Stderr, ", "); + grcc_fprintf(GRCC_Stderr, "%2d", p[j]); } - fprintf(GRCC_Stderr, "]%s", msg); + grcc_fprintf(GRCC_Stderr, "]%s", msg); } //-------------------------------------------------------------- @@ -10905,7 +10937,7 @@ static int intSetAdd(int n, int *a, int v, const int size) int j, k; if (n >= size) { - fprintf(GRCC_Stderr, "*** intSetAdd : array out of range (>%d)\n", size); + grcc_fprintf(GRCC_Stderr, "*** intSetAdd : array out of range (>%d)\n", size); erEnd("intSetAdd : array out of range (GRCC_MAXPSLIST)"); } for (j = 0; j < n; j++) { @@ -10935,7 +10967,7 @@ static int intSListAdd(int n, int *a, int v, const int size) int j, k; if (n >= size) { - fprintf(GRCC_Stderr, "*** intSListAdd : array out of range (>%d)\n", size); + grcc_fprintf(GRCC_Stderr, "*** intSListAdd : array out of range (>%d)\n", size); erEnd("intSListAdd : array out of range"); } for (j = 0; j < n; j++) { diff --git a/sources/grccparam.h b/sources/grccparam.h index 88d03e5d..7ab17eb2 100644 --- a/sources/grccparam.h +++ b/sources/grccparam.h @@ -26,15 +26,15 @@ #define GRCC_FRACERROR (1.0e-10) -/* error message */ -/* -#define GRCC_Stderr stderr -*/ +/* Standard and error message destination */ +#define GRCC_Stdout stdout #define GRCC_Stderr stdout -/* -#define GRCC_ABORT() exit(1) -*/ -#define GRCC_ABORT() abort() + +#ifndef NOFORM + #define GRCC_ABORT() Terminate(-1); +#else + #define GRCC_ABORT() exit(1) +#endif /*--------------------------------------------------------------- * Constants diff --git a/sources/model.c b/sources/model.c index 81ce8606..18026e0c 100644 --- a/sources/model.c +++ b/sources/model.c @@ -318,6 +318,11 @@ int CoVertex(UBYTE *s) return(1); } ss = s; s = SkipName(s); c = *s; *s = 0; name = ConstructName(ss,0); + // Forbid couplings which have no symbol or an overal numerical factor. + if ( *name == 0 ) { + v->error = 1; + MesPrint("&Invalid coupling constant in vertex statement."); + } if ( GetVar(name,&type,&v->couplings[v->ncouplings],CSYMBOL,WITHAUTO) == NAMENOTFOUND ) { WORD minpow = -MAXPOWER; @@ -337,6 +342,10 @@ int CoVertex(UBYTE *s) if ( *s == '-' ) sign = -sign; s++; } + if ( sign == -1 ) { + v->error = 1; + MesPrint("&Invalid negative power of coupling constant."); + } while ( FG.cTable[*s] == 1 ) x = 10*x + *s++ - '0'; v->couplings[v->ncouplings++] = x; } diff --git a/sources/proces.c b/sources/proces.c index 40fb5cf5..832f5726 100644 --- a/sources/proces.c +++ b/sources/proces.c @@ -1183,28 +1183,9 @@ Important: we may not have enough spots here } } else if ( *t == TOPOLOGIES ) { -/* - Syntax: - topologies_(nloops,nlegs,setvertexsizes,setext,setint[,options]) -*/ - t1 = t+FUNHEAD; t2 = t+t[1]; - if ( *t1 == -SNUMBER && t1[1] >= 0 && - t1[2] == -SNUMBER && ( t1[3] >= 0 || t1[3] == -2 ) && - t1[4] == -SETSET && Sets[t1[5]].type == CNUMBER && - t1[6] == -SETSET && Sets[t1[7]].type == CVECTOR && - t1[8] == -SETSET && Sets[t1[9]].type == CVECTOR && - t1+10 <= t2 ) { - if ( t1+10 == t2 || ( t1+12 <= t2 && ( t1[10] == -SNUMBER || - ( t1[10] == -SETSET && - Sets[t1[5]].last-Sets[t1[5]].first == - Sets[t1[11]].last-Sets[t1[11]].first ) ) ) ) { - AN.TeInFun = -15; - AN.TeSuOut = 0; - AR.TePos = -1; - AR.funoffset = t - term; - DONE(1) - } - } + MesPrint("&The topologies_ function was removed in FORM 5.0."); + MesPrint("&See the TopologiesOnly_ option of diagrams_."); + Terminate(-1); } else if ( *t == DIAGRAMS ) { /* @@ -1220,8 +1201,10 @@ Important: we may not have enough spots here || ( Sets[t1[3]].type == ANYTYPE && ( Sets[t1[3]].first == Sets[t1[3]].last ) ) ) && t1[4] == -SETSET && ( Sets[t1[5]].type == CFUNCTION || ( Sets[t1[5]].type == ANYTYPE && ( Sets[t1[5]].first == Sets[t1[5]].last ) ) ) && - t1[6] == -SETSET && Sets[t1[7]].type == CVECTOR && - t1[8] == -SETSET && Sets[t1[9]].type == CVECTOR && + t1[6] == -SETSET && ( Sets[t1[7]].type == CVECTOR + || ( Sets[t1[7]].type == ANYTYPE && ( Sets[t1[7]].first == Sets[t1[7]].last ) ) ) && + t1[8] == -SETSET && ( Sets[t1[9]].type == CVECTOR + || ( Sets[t1[9]].type == ANYTYPE && ( Sets[t1[9]].first == Sets[t1[9]].last ) ) ) && t1+12 <= t2 ) { /* Test that the sets of particles correspond to particles @@ -1260,7 +1243,7 @@ Important: we may not have enough spots here tt1 += 2; } } - AN.TeInFun = -16; + AN.TeInFun = -15; AN.TeSuOut = 0; AR.TePos = -1; AR.funoffset = t - term; @@ -1286,7 +1269,7 @@ Important: we may not have enough spots here tt1 += 2; } } - AN.TeInFun = -16; + AN.TeInFun = -15; AN.TeSuOut = 0; AR.TePos = -1; AR.funoffset = t - term; @@ -4078,9 +4061,6 @@ AutoGen: i = *AT.TMout; if ( DIVfunction(BHEAD term,level,3) < 0 ) goto GenCall; break; case -15: - if ( GenTopologies(BHEAD term,level) < 0 ) goto GenCall; - break; - case -16: if ( GenDiagrams(BHEAD term,level) < 0 ) goto GenCall; break; } diff --git a/sources/topowrap.cc b/sources/topowrap.cc deleted file mode 100644 index 804021c4..00000000 --- a/sources/topowrap.cc +++ /dev/null @@ -1,283 +0,0 @@ -/** @file topowrap.cc - * - * routines for conversion of topology and diagram output to FORM notation - */ - -/* #[ License : */ -/* - * Copyright (C) 1984-2023 J.A.M. Vermaseren - * When using this file you are requested to refer to the publication - * J.A.M.Vermaseren "New features of FORM" math-ph/0010025 - * This is considered a matter of courtesy as the development was paid - * for by FOM the Dutch physics granting agency and we would like to - * be able to track its scientific use to convince FOM of its value - * for the community. - * - * This file is part of FORM. - * - * FORM is free software: you can redistribute it and/or modify it under the - * terms of the GNU General Public License as published by the Free Software - * Foundation, either version 3 of the License, or (at your option) any later - * version. - * - * FORM is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - * details. - * - * You should have received a copy of the GNU General Public License along - * with FORM. If not, see . - */ -/* - #] License : - #[ includes : -*/ - -extern "C" { -#include "form3.h" -} - -#include "gentopo.h" - -//using namespace std; - -#define MAXPOINTS 120 - -typedef struct ToPoTyPe { - int cldeg[MAXPOINTS], clnum[MAXPOINTS], clext[MAXPOINTS]; - int ncl, nloops, nlegs, npadding; - WORD *vert; - WORD *vertmax; - WORD nvert; - WORD sopi; -} TOPOTYPE; - -/* - #] includes : - #[ GenerateVertices : - - Routine is to be used recursively to work its way through a list - of possible vertices. The array of vertices is in TopoInf->vert - with TopoInf->nvert the number of possible vertices. - Currently we allow in TopoInf->vert only vertices with 3 or more edges. - - We work with a point system. Each n-point vertex contributes n-2 points. - When all points are assigned, we can call mgraph->generate(). - - The number of vertices of a given number of edges is stored in - TopoInf->clnum[..] but the loop that determines how many there are - may be limited by the corresponding element in TopoInf->vertmax[level] -*/ - -int GenerateVertices(TOPOTYPE *TopoInf, int pointsremaining, int level) -{ - int i, j; - - for ( i = pointsremaining, j = 0; i >= 0; i -= TopoInf->vert[level]-2, j++ ) { - if ( TopoInf->vertmax && TopoInf->vertmax[level] >= 0 - && j > TopoInf->vertmax[level] ) break; - if ( i == 0 ) { // We got one! - T_MGraph *mgraph; - TopoInf->cldeg[TopoInf->ncl] = TopoInf->vert[level]; - TopoInf->clnum[TopoInf->ncl] = j; - TopoInf->clext[TopoInf->ncl] = 0; - TopoInf->ncl++; - - mgraph = new T_MGraph(0, TopoInf->ncl, TopoInf->cldeg, - TopoInf->clnum, TopoInf->clext, TopoInf->sopi); - - mgraph->generate(); - - delete mgraph; - - TopoInf->ncl--; - - break; - } - if ( level < TopoInf->nvert-1 ) { - if ( j > 0 ) { - TopoInf->cldeg[TopoInf->ncl] = TopoInf->vert[level]; - TopoInf->clnum[TopoInf->ncl] = j; - TopoInf->clext[TopoInf->ncl] = 0; - TopoInf->ncl++; - } - if ( GenerateVertices(TopoInf,i,level+1) < 0 ) return(-1); - if ( j > 0 ) { TopoInf->ncl--; } - } - } - return(0); -} - -/* - #] GenerateVertices : - #[ GenerateTopologies : - - Note that setmax, option1 and option2 are optional. - Default values are -1,0,0 - - vert is a pointer to a set of numbers indicating the types of vertices - that are allowed. - nvert is the number of elements in vert. -*/ - -int GenerateTopologies(PHEAD WORD nloops, WORD nlegs, WORD setvert, WORD setmax) -{ - TOPOTYPE TopoInf; - int i, points, identical = 0; - DUMMYUSE(AT.nfac); - - if ( nlegs == -2 ) { - nlegs = 2; - identical = 1; - } - TopoInf.vert = &(SetElements[Sets[setvert].first]); - TopoInf.nvert = Sets[setvert].last-Sets[setvert].first; - - if ( setmax >= 0 ) TopoInf.vertmax = &(SetElements[Sets[setmax].first]); - else TopoInf.vertmax = 0; - -// point counting: an n-point vertex counts for n-2 points. - - points = 2*nloops-2+nlegs; - if ( points >= MAXPOINTS ) { - MLOCK(ErrorMessageLock); - MesPrint("GenerateTopologies: %d loops and %d legs considered excessive",nloops,nlegs); - MUNLOCK(ErrorMessageLock); - Terminate(-1); - } - -// First the external nodes. - - for ( i = 0; i < nlegs; i++ ) { - TopoInf.cldeg[i] = 1; TopoInf.clnum[i] = 1; TopoInf.clext[i] = 1; - } - if ( identical == 1 ) { /* Only propagator topologies..... */ - nlegs = 1; - TopoInf.clnum[0] = 2; - } - TopoInf.ncl = nlegs; - TopoInf.sopi = 1; // For now - if ( GenerateVertices(&TopoInf,points,0) != 0 ) { - MLOCK(ErrorMessageLock); - MesPrint("Called from GenerateTopologies with %d loops and %d legs considered excessive",nloops,nlegs); - MUNLOCK(ErrorMessageLock); - Terminate(-1); - } - return(0); -} - -/* - #] GenerateTopologies : - #[ toForm : -*/ - -//============================================================== -// Output for FROM called by genTopo each time a new graph is -// generated -// -void toForm(T_EGraph *egraph) -{ - GETIDENTITY; -// -// skipext : boolean; whether to skip external node/line in the output. -// -// const int skipext = 1; - - int n, lg, ed, i, fromset; - - WORD *termout = AT.WorkPointer; - WORD *t, *tt, *ttstop, *ttend, *ttail; - - t = termout+1; - -// First pick up part of the original term - - tt = AT.TopologiesTerm + 1; - i = AT.TopologiesStart - tt; - NCOPY(t,tt,i) - ttail = tt+tt[1]; - -// Now the vertices/nodes -// Options are in AT.TopologiesOptions[] - - for ( n = 0; n < egraph->nNodes; n++ ) { - if ( ( AT.TopologiesOptions[1] &1 ) == 1 && egraph->nodes[n].ext ) continue; - tt = t; - *t++ = NODEFUNCTION; t++; FILLFUN(t); - *t++ = -SNUMBER; *t++ = n; - for ( lg = 0; lg < egraph->nodes[n].deg; lg++ ) { - ed = egraph->nodes[n].edges[lg]; - if ( ed >= 0 ) { *t++ = -VECTOR; } - else { ed = -ed; *t++ = -MINVECTOR; } - -// Now we need to pick up the proper set element. - - fromset = egraph->edges[ed].ext ? AT.setexterntopo : AT.setinterntopo; - *t++ = SetElements[Sets[fromset].first+egraph->edges[ed].momn-1]; - } - tt[1] = t-tt; - } - - if ( ( AT.TopologiesOptions[0] & 1 ) == 1 ) { -// Note that the edges count from 1. - for ( n = 1; n <= egraph->nEdges; n++ ) { - if ( ( AT.TopologiesOptions[1] & 1 ) == 1 && egraph->edges[n].ext ) continue; - tt = t; - *t++ = EDGE; t++; FILLFUN(t); - *t++ = -SNUMBER; *t++ = egraph->edges[n].nodes[0]; - *t++ = -SNUMBER; *t++ = egraph->edges[n].nodes[1]; - *t++ = -VECTOR; - fromset = egraph->edges[n].ext ? AT.setexterntopo : AT.setinterntopo; - *t++ = SetElements[Sets[fromset].first+egraph->edges[n].momn-1]; - tt[1] = t-tt; - } - } - -// Now the tail end - - ttend = AT.TopologiesTerm; ttend = ttend+ttend[0]; - ttstop = ttend - ABS(ttend[-1]); - i = ttstop - ttail; - NCOPY(t,ttail,i) - -// Finally the coefficient -// The topological coefficient should be in egraph->wsum, egraph->nwsum -// as a FORM Long number - -#ifdef WITHFACTOR - if ( egraph->nwsum == 1 && egraph->wsum[0] == 1 ) { - while (ttstop < ttend ) *t++ = *ttstop++; - } - else { - WORD newsize; - newsize = ttend[-1]; - newsize = REDLENG(newsize); - if ( Divvy(BHEAD (UWORD *)ttstop,&newsize,egraph->wsum,egraph->nwsum) ) - goto OnError; - newsize = INCLENG(newsize); - i = ABS(newsize)-1; - NCOPY(t,ttstop,i) - *t++ = newsize; - } -#else - while (ttstop < ttend ) *t++ = *ttstop++; -#endif - *termout = t - termout; - - AT.WorkPointer = t; - - if ( Generator(BHEAD termout,AT.TopologiesLevel) < 0 ) { -// OnError: - MLOCK(ErrorMessageLock); - MesPrint("Called from the topologies routine toForm"); - MUNLOCK(ErrorMessageLock); - Terminate(-1); - } - - AT.WorkPointer = termout; -} - -/* - #] toForm : -*/ -